雷電產生的海洋電磁場研究中Levinson-Durbin算法實現

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

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

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

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

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

1 Levinson-Durbin算法

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

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

2 仿真信號的比較

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

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

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

3 實測數據的比較

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

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

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

2）根據已知觀測數據估計模型參數。

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

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

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

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

FPE準則函數：

AIC準則函數：

圖3 FPE與模型階的關系曲線Fig.3 The relationship between FPE and order

圖4 AIC與模型階的關系曲線Fig.4 The relationship between AIC and order

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

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

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

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

4 結語

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

利用現代譜估計技術準確識別雷電環境，可以為海軍裝備設計、研制提供參考，而且也是合理使用、充分發揮武器裝備性能的基礎信息。

[1]李楠，張雪飛.戰場復雜電磁環境構成分析[J].裝備環境工程，2008，5（1）：16—19.

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

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

[4]姚天任，孫洪.現代數字信號處理[M].武漢：華中科技大學出版社，1996：123—131.

[5]張永春，王明江.線性預測編碼中Levinson-Durbin 算法的ASIC實現[J].計算機與數字工程，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模型參數估計[J].工具技術，2005，39（4）：39—40.