雷電產(chǎn)生的海洋電磁場研究中Levinson-Durbin算法實現(xiàn)

崔培，李琳琳，李沛劍，岳瑞永

（1.大連測控技術(shù)研究所，遼寧 大連 116013；2.海軍駐航天一院軍代表室，北京 100076）

海洋中存在多種電磁效應(yīng)，根據(jù)場源的形式可分為自然電磁場和人為因素形成的電磁場。

自然電磁輻射是非人為因素產(chǎn)生的電磁波輻射。在自然電磁環(huán)境中，靜電、雷電和地磁場等自然輻射是幾種最主要的電磁輻射[1]。地球大氣中雷電現(xiàn)象是海洋自然電磁場的主要場源之一。雷電產(chǎn)生的地磁變異，呈短暫的脈動形式，具有較強的地域性和季節(jié)性。發(fā)生雷擊時，雷電通道周圍會產(chǎn)生強烈的電磁輻射場，當雷擊發(fā)生在海面上的時候，會對水下工作的電磁場傳感器以及水中兵器產(chǎn)生影響[2]。由于海水對高頻電磁波具有強烈的吸收作用，雷電產(chǎn)生的海洋電磁場能量主要集中在極低頻。

使用經(jīng)典算法進行信號譜估計時，為獲得必需的頻率分辨率需要記錄較長的數(shù)據(jù)，同時，由于加窗必然存在功率泄露和頻率混疊使接受的弱信號被掩蓋。近代譜估計就是為了克服以上缺點而提出的，它的基本思想是對所觀察的有限數(shù)據(jù)以外的數(shù)據(jù)不作任何確定性假設(shè)。在一定先驗知識的基礎(chǔ)上，采用外推或預(yù)測的方法，從觀測的數(shù)據(jù)中推導(dǎo)以后的數(shù)據(jù)[3]。由于雷電現(xiàn)象持續(xù)的時間較短，能量又主要集中在極低頻，屬于短數(shù)據(jù)記錄。利用傳統(tǒng)的譜估計方法難以克服頻率分辨率低和譜泄露的缺點，因此，為了解決這些問題，在處理雷電短記錄數(shù)據(jù)時引入了現(xiàn)代譜估計技術(shù)。

1 Levinson-Durbin算法

以AR模型為基礎(chǔ)的現(xiàn)代譜估計，需要知道模型的階和AR 系數(shù)，以及模型激勵源的方差。為此，必須把這些參數(shù)和已知的自相關(guān)函數(shù)聯(lián)系起來，這就是著名的Yule-Walker 方程。用線性方程組的常用解法求解Yule-Walker 方程，需要的運算量較大，但若利用系數(shù)矩陣的對稱性和Toeplitz性質(zhì)，則可構(gòu)成一些高效算法。Levinson-Durbin（簡稱L-D）算法是其中最著名、應(yīng)用最廣泛的一種[4]。

L-D算法的運算步驟如下：

2 仿真信號的比較

由于雷電為短暫的脈動形式，并且海水對高頻電磁波的衰減較快。雷電產(chǎn)生的海洋電磁場能量主要集中在10 Hz以下的極低頻，信號的時間長度甚至不足1個周期。由于傳統(tǒng)的譜估計方法對短記錄數(shù)據(jù)的頻率分辨率較低，會出現(xiàn)嚴重的譜泄漏現(xiàn)象，從而導(dǎo)致估計的線譜頻率不準，而L-D 算法甚至可以從部分周期中提取出線譜信息，但FFT 如果沒有幾個周期就無法做到這一點[6]。

為了比較研究FFT 與L-D 算法的性能，首先將其應(yīng)用在仿真信號中。仿真信號的頻率為5 Hz，利用 MATLAB 軟件進行繪圖，得到 0.5T，0.9T，2.5T 及4T 的頻譜圖（T 為信號1 個周期的長度）。為了方便比較，對頻譜分析結(jié)果進行了歸一化處理，得到的頻譜如圖1所示。

從圖1 中可以看出，當信號不足1 個周期時，經(jīng)FFT 得到的分析結(jié)果線譜較寬，無法準確判斷信號頻率；L-D 算法的頻譜結(jié)果表明，當信號記錄不足1個周期時（如圖1b 所示）就可以提取出線譜頻率。隨著信號周期的增加，F(xiàn)FT的線譜特性越來越明顯，但與L-D算法分析結(jié)果相比，其線譜仍較寬。

3 實測數(shù)據(jù)的比較

利用大連測控技術(shù)研究所自主研制的海洋環(huán)境電磁場觀測系統(tǒng)，多次采集到海面雷電現(xiàn)象。下面以一次雷電產(chǎn)生的垂直分量電磁場為例，比較分析FFT與L-D算法的應(yīng)用結(jié)果。雷電產(chǎn)生的電磁場時域信號如圖2所示，數(shù)據(jù)截取的時間長度為1 s，雷電信號的持續(xù)時間約為0.5 s。

以參數(shù)模型為基礎(chǔ)的譜估計方法一般按下列3個步驟進行。

1）為被估計的隨機過程確定或選擇一個合理的模型，文中選擇的是AR模型。

2）根據(jù)已知觀測數(shù)據(jù)估計模型參數(shù)。

3）用估計得到的模型參數(shù)計算功率譜[4]。

圖1 不同時間長度的頻譜Fig.1 Spectrum of different time length

圖2 雷電產(chǎn)生的電磁場時域波形Fig. 2 Time domain waveform of electromagnetic field in sea generated by lightning

采用AR 模型譜估計方法，既要估計AR 模型參數(shù)，又要估計模型的階。文中以最終預(yù)測誤差（FPE）準則和Akaike 信息準則（AIC）結(jié)合實測數(shù)據(jù)進行討論。

FPE準則函數(shù)：

AIC準則函數(shù)：

圖3 FPE與模型階的關(guān)系曲線Fig.3 The relationship between FPE and order

圖4 AIC與模型階的關(guān)系曲線Fig.4 The relationship between AIC and order

根據(jù)模型階數(shù)和L-D 遞推算法得到的AR 模型參數(shù)可以估計信號的功率譜。為了方便研究，對功率譜進行了幅度歸一化，得到各自的頻譜如圖5、圖6所示。

圖5 FFT直接輸出結(jié)果和L-D算法比較Fig.5 Comparison between FFT（primary）and L-D

圖6 FFT插值后和L-D算法比較Fig.6 Comparison between FFT（interpolation）and L-D

由分析處理結(jié)果可知：L-D 算法具有明顯的線譜特性，線譜較窄，檢測出的信號頻率為6.2 Hz；FFT的線譜較寬，線譜分辨率低，線譜特性不明顯，檢測出的信號頻率在5.9 Hz附近。

4 結(jié)語

針對雷電產(chǎn)生的海洋電磁場記錄時間短的特點，研究了仿真信號FFT 與L-D 算法的頻譜分析結(jié)果。結(jié)果表明，L-D算法利用不足1個周期的數(shù)據(jù)便可以準確提取出信號的線譜頻率；FFT 在數(shù)據(jù)量不足幾個周期時，其頻譜分辨率仍然較低，無法得到確切的線譜頻率信息。基于仿真信號的研究，選擇AR模型作為雷電隨機過程的模型。結(jié)合FPE和AIC準則確定AR 模型階數(shù)，利用L-D 算法求取AR 模型參數(shù)，從而計算出其功率譜。對雷電實測數(shù)據(jù)的處理結(jié)果達到了預(yù)期效果，體現(xiàn)了L-D算法的應(yīng)用價值。

利用現(xiàn)代譜估計技術(shù)準確識別雷電環(huán)境，可以為海軍裝備設(shè)計、研制提供參考，而且也是合理使用、充分發(fā)揮武器裝備性能的基礎(chǔ)信息。

[1]李楠，張雪飛.戰(zhàn)場復(fù)雜電磁環(huán)境構(gòu)成分析[J].裝備環(huán)境工程，2008，5（1）：16—19.

[2]李松，龔沈光.空氣-海水兩層模型中雷電在海水中產(chǎn)生的電磁場計算[J].海軍工程大學(xué)學(xué)報，2008，20（6）：40—44.

[3]劉茵，李誠人. 自回歸濾波器的研究[J]. 計算機仿真，2006，23（6）：107—108.

[4]姚天任，孫洪.現(xiàn)代數(shù)字信號處理[M].武漢：華中科技大學(xué)出版社，1996：123—131.

[5]張永春，王明江.線性預(yù)測編碼中Levinson-Durbin 算法的ASIC實現(xiàn)[J].計算機與數(shù)字工程，2005，33（9）：118—119.

[6]IAN B Murray，ALASTAIR D McAulay. Magnetic Detection and Localization Using Multichannel Levinson-Durbin Algorithm[J]. Proceedings of SPIE，2004：561—566.（余不詳）

[7]陳國強，趙俊偉，黃俊杰，等.基于Matlab的AR模型參數(shù)估計[J].工具技術(shù)，2005，39（4）：39—40.