小波閾值去噪與嶺估計在飛機磁補償中的應(yīng)用

張堅，林春生，周建軍

(1.海軍青島雷達聲納修理廠，山東 青島 266100；2.海軍工程大學(xué)兵器工程系，湖北 武漢 430033)

航空磁性探測廣泛應(yīng)用于考古、礦物勘探和水中磁性目標探測等領(lǐng)域[1].但是在飛機上對地磁異常場進行測量時，飛機干擾磁場疊加在地磁異常場上，掩蓋了地磁異常場原有的特性，所以必須對飛機干擾磁場進行補償，利用測量數(shù)據(jù)求解飛機干擾磁場模型系數(shù)，使用所求系數(shù)計算干擾磁場，最后從測量的總磁場中減去干擾磁場，完成飛機磁補償.因此，飛機干擾磁場模型的求解是飛機磁補償中的關(guān)鍵技術(shù).

美國學(xué)者Tolles和Lawson在上世紀40年代最早對飛機干擾磁場的數(shù)學(xué)模型進行了研究，建立了經(jīng)典的Tolles和Lawson模型，但沒有給出模型的求解方法[2].其后，Leliak和Bickel均對該模型進行了改進[3-4]，但其求解精度仍然不高.飛機干擾磁場模型屬于線性模型，并且存在較強的復(fù)共線性.線性模型求解的常用方法有最小二乘估計LSE和嶺估計RE等[5-6].最小二乘估計運算量小，但在模型存在復(fù)共線性時均方誤差很大；嶺估計可以在一定程度上克服模型復(fù)共線性，但在實際探測中，飛機上電氣設(shè)備(如繼電器、電機等)產(chǎn)生較強的與飛機姿態(tài)變化無關(guān)的干擾磁場，這種測量噪聲疊加在觀測數(shù)據(jù)上無法進行補償，并會對嶺估計的求解效果產(chǎn)生較大影響.針對飛機干擾磁場模型求解時均方誤差較大的問題，本文中提出小波閾值去噪結(jié)合嶺估計(WTRE)的模型參數(shù)求解方法，并通過參數(shù)估計和飛機磁補償?shù)姆抡鎸嶒瀸Ψ椒ㄟM行了驗證.

1 飛機干擾磁場建模

圖1 飛機干擾磁場模型坐標系

飛機干擾磁場包括剩余磁場、感應(yīng)磁場和渦流磁場，本文中采用經(jīng)典的Tolles和Lawson模型進行飛機干擾磁場的建模，并對模型進行簡化，模型坐標系如圖1所示.圖中，機載磁傳感器為坐標原點，地磁場矢量He與X，Y，Z軸的夾角為(X0，Y0，Z0)，其方向余弦定義為(cosX0，cosY0，cosZ0).

1)剩余磁場Hp是由飛機上鐵磁性物質(zhì)的剩余磁化產(chǎn)生的，在短時間內(nèi)可以視為穩(wěn)定不變，建模時將其分解為與X，Y，Z軸平行的三分量Hpx，Hpy，Hpz，即Hp=(HpxHpyHpz)T,

則剩余磁場系數(shù)有3項.

其中，K為感應(yīng)磁場系數(shù)矩陣，有9項系數(shù).

其中，L為渦流磁場系數(shù)矩陣，有9項系數(shù).

上述模型是傳統(tǒng)的飛機干擾磁場模型，共有21項模型系數(shù).利用感應(yīng)磁場系數(shù)矩陣K的對稱性[3]，以及方向余弦cos2X0+cos2Y0+cos2Z0=1的性質(zhì)，可將3種磁場相加后再進行合并與化簡，最終得到總的飛機干擾磁場ΔH為ΔH=ABTHe+kzzHe

(1)

A=(Hpx，Hpy，Hpz，kxx-kzz，kxy+kyx，kzx+kxz，kyz+kzy，kyy-kzz，lxx-lzz，lyx，lzx，lxy，lyy-lzz，lzy，lxz，lyz)

B= (cosX0/He，cosY0/He，cosZ0/He，cos2X0，cosX0cosY0，cosX0cosZ0，cosY0cosZ0，cos2Y0，

(1)式是簡化后的飛機干擾磁場模型，共有16項系數(shù).飛機干擾磁場的補償過程是，利用測量數(shù)據(jù)估計線性方程組(1)中的16項未知系數(shù)，再根據(jù)估計出的系數(shù)計算飛機干擾磁場，從磁場測量值中減去飛機干擾磁場，得到補償后的磁場值.

2 小波閾值去噪結(jié)合嶺估計的模型求解方法

(2)

式中，Y為n×1觀測向量；X為n×m列滿秩設(shè)計矩陣，為已知量；β為m×1待求的模型參數(shù)；e為n×1隨機誤差向量；In為n階單位陣；σ2為方差參數(shù).

對于飛機干擾磁場模型，自變量X由方向余弦、方向余弦的導(dǎo)數(shù)及其乘積項構(gòu)成，Y由磁場總強度的測量值構(gòu)成，e為測量誤差向量.因為在一定范圍內(nèi)地磁偏角和傾角可看做不變，顯然Y是飛機姿態(tài)角變化的函數(shù).用LS估計求解模型時均方誤差很大，原因在于設(shè)計矩陣XTX的奇異性.該矩陣產(chǎn)生奇異性的原因是，在保證飛機穩(wěn)定飛行的情況下，飛機的機動角度不可能很大，以致XTX每一行的對應(yīng)元素值非常接近，每行之間有比較強的相關(guān)性，導(dǎo)致飛機干擾磁場模型存在嚴重的復(fù)共線性.

2.3小波閾值去噪與嶺估計的結(jié)合小波閾值去噪是Donoho提出的一種簡潔有效的去噪方法，通過對小波系數(shù)進行非線性閾值處理，恢復(fù)噪聲中的信號.算法的基本思想如下：屬于Besov空間的信號在小波域的能量主要集中在有限的幾個小波系數(shù)中，而白噪聲的能量卻分布在整個小波域內(nèi).因此對wj，k，可以找到一個合適的閾值λ，當(dāng)wj，k<λ時，認為此時的wj，k主要由噪聲引起，將其置為零；當(dāng)wj，k>λ時，認為此時的wj，k主要由信號引起，將這部分wj，k直接保留下來(硬閾值法)或者按某一固定量向零收縮(軟閾值法).最后將處理過的小波系數(shù)進行小波重構(gòu)而得到去噪后的信號.

在航空磁性探測過程中，飛機上電氣設(shè)備(如繼電器、電機等)會產(chǎn)生較強的干擾磁場，這種測量噪聲與飛行姿態(tài)無關(guān)，疊加在觀測數(shù)據(jù)上將對模型求解產(chǎn)生較大的影響.由于測量噪聲在有限頻帶內(nèi)可以近似看作白噪聲，因此本文中采用小波閾值去噪作為嶺估計前的預(yù)處理手段，對測量噪聲進行濾除.小波閾值去噪結(jié)合嶺估計(WTRE)的實施過程如下：求解模型參數(shù)之前，采用小波閾值去噪對觀測數(shù)據(jù)進行預(yù)處理，降低電氣磁場等測量噪聲對模型求解的影響，閾值函數(shù)選取軟閾值；然后利用嶺估計求解干擾磁場模型系數(shù)，嶺參數(shù)的選取采用L曲線法.

3 算法驗證

3.1參數(shù)估計仿真實驗采用Matlab7.0軟件仿真驗證WTRE在飛機干擾磁場模型求解中的有效性.設(shè)地磁場大小為5×104nT，地磁場方向余弦可以用地磁場方位角θ和俯仰角φ的函數(shù)來表示，且假設(shè)dφ/dt=1，dθ/dt=cos(t/2)，φ∈[0，2π]，θ∈[0，π].觀測數(shù)據(jù)采樣點數(shù)為400，即N=400.設(shè)定模型系數(shù)的真值如下：剩余磁場系數(shù)[10，35，90]，5個感應(yīng)磁場系數(shù)均為0.5，8個渦流磁場系數(shù)均為0.2.利用設(shè)定的系數(shù)真值，仿真產(chǎn)生信噪比分別為SNR1=21.47 dB和SNR2=9.13 dB時的飛機干擾磁場信號，如圖2所示.

圖2 不同信噪比下飛機干擾磁場信號仿真

衡量矩陣病態(tài)程度的方法之一是求它的條件數(shù).計算得到法矩陣XTX的條件數(shù)為1.35×1015，呈嚴重病態(tài)性.分別采用LS估計(LSE)、嶺估計(RE)和WTRE 3種方法進行求解，計算結(jié)果如表1所示，

表1 LS估計、嶺估計和WTRE的計算結(jié)果

嶺估計采用L曲線法確定k值.采用均方誤差(MSE)作為求解效果的評價標準：

由表1的計算結(jié)果可以得出以下結(jié)論：

1)LS估計的均方誤差遠大于嶺估計，原因是法矩陣XTX的條件數(shù)很大，模型復(fù)共線性嚴重.

2)利用嶺估計進行飛機干擾磁場模型參數(shù)的求解可以較好地克服復(fù)共線性的影響，且基于L曲線法的嶺參數(shù)選取方法適用性好、穩(wěn)定性高，能夠滿足航空磁性探測的實際需求，估計精度較高.

3)在不同信噪比下，WTRE的均方誤差均小于RE，特別是在低信噪比條件下，WTRE對RE的改進更加明顯，這是由于WTRE降低了測量噪聲對模型求解的影響.

圖3 磁偶極子目標信號

式中，μ0=4π×10-7H/m為真空磁導(dǎo)率，M為磁性目標的磁矩矢量，r為目標到磁傳感器的距離.通過設(shè)定一組參數(shù)，可以由上述磁偶極子模型仿真生成目標磁場信號，如圖3所示.目標信號與飛機干擾磁場疊加后的混合信號如圖4(a)所示.由圖4(a)可以看出，飛機干擾磁場的變化幅度超過10 nT，目標信號完全被干擾磁場淹沒，無法辨識.分別采用LS估計、嶺估計和WTRE所求系數(shù)進行磁補償后的信號如圖4(b)～(d)所示.可以看出，LS估計補償后的信噪比最低，WTRE補償后的信噪比最高.3種方法中只有WTRE的補償結(jié)果可以準確檢測到目標信號，其他2種方法都在不同程度上受到噪聲的干擾.這進一步驗證了WTRE在飛機磁補償中的有效性.

(a)目標與飛機干擾磁場的混合信號

(b)LS估計的補償結(jié)果

(c)嶺估計的補償結(jié)果

(d)WTRE的補償結(jié)果

4 結(jié)論

針對飛機干擾磁場模型求解時均方誤差較大的問題，本文中提出小波閾值去噪結(jié)合嶺估計(WTRE)的模型參數(shù)求解方法.首先采用小波閾值去噪對觀測數(shù)據(jù)進行預(yù)處理，降低測量噪聲對模型求解的影響；然后利用嶺估計求解干擾磁場模型系數(shù)，嶺參數(shù)的選取采用L曲線法.仿真結(jié)果表明，WTRE的均方誤差明顯小于LS估計和嶺估計，特別是在低信噪比情況下對嶺估計的改進尤為明顯，采用WTRE所求系數(shù)補償后的信噪比最高，在飛機磁補償中與其它方法相比具有明顯的優(yōu)越性.

[1] 熊盛青. 我國航空重磁勘探技術(shù)現(xiàn)狀與發(fā)展趨勢[J].地球物理學(xué)進展， 2009，24(1):113-117.

[2] Tolles W E， Lawson J D. Magnetic compensation of MAD equipped aircraft[R].New York:Airborne Instruments Lab Inc，1950.

[3] Leach B W. Aeromagnetic compensation as a linear regression problem[C]//Information linkage between applied mathematics and industry II. New York:Academic Press Inc， 1980：367-371．

[4] Bickel S H. Small signal compensation of magnetic fields resulting from aircraft maneuvers[J].IEEE Trans on AES，1979， 15(4):515-525.

[2] 陳希孺， 王松桂. 近代回歸分析-原理方法及應(yīng)用[M].合肥:安徽教育出版社， 1987:294-320.

[6] 劉莉. NQD樣本下部分線性模型中估計的強相合性[J].湖北大學(xué)學(xué)報：自然科學(xué)版， 2004， 26(4):290-293.

[7] Hansen P C. Analysis of discrete Ill-posed problems by means of the L-curve[J].SIAM Review，1992，34(4):561-580.

[8] 王振杰， 歐吉坤. 用L-曲線法確定嶺估計中的嶺參數(shù)[J].武漢大學(xué)學(xué)報：信息科學(xué)版， 2004， 29(3):235-238.

[9] 林春生， 龔沈光. 艦船物理場[M].北京:兵器工業(yè)出版社， 2007.