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

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

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

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

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

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

1 Levinson-Durbin算法

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

L-D算法的運(yùn)算步驟如下：

2 仿真信號的比較

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

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

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

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

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

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

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

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

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

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

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

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

FPE準(zhǔn)則函數(shù)：

AIC準(zhǔn)則函數(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ù)可以估計(jì)信號的功率譜。為了方便研究，對功率譜進(jìn)行了幅度歸一化，得到各自的頻譜如圖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)生的海洋電磁場記錄時(shí)間短的特點(diǎn)，研究了仿真信號FFT 與L-D 算法的頻譜分析結(jié)果。結(jié)果表明，L-D算法利用不足1個(gè)周期的數(shù)據(jù)便可以準(zhǔn)確提取出信號的線譜頻率；FFT 在數(shù)據(jù)量不足幾個(gè)周期時(shí)，其頻譜分辨率仍然較低，無法得到確切的線譜頻率信息。基于仿真信號的研究，選擇AR模型作為雷電隨機(jī)過程的模型。結(jié)合FPE和AIC準(zhǔn)則確定AR 模型階數(shù)，利用L-D 算法求取AR 模型參數(shù)，從而計(jì)算出其功率譜。對雷電實(shí)測數(shù)據(jù)的處理結(jié)果達(dá)到了預(yù)期效果，體現(xiàn)了L-D算法的應(yīng)用價(jià)值。

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

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

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

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

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

[5]張永春，王明江.線性預(yù)測編碼中Levinson-Durbin 算法的ASIC實(shí)現(xiàn)[J].計(jì)算機(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]陳國強(qiáng)，趙俊偉，黃俊杰，等.基于Matlab的AR模型參數(shù)估計(jì)[J].工具技術(shù)，2005，39（4）：39—40.