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

張堅，林春生，周建軍

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

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

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

1 飛機干擾磁場建模

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

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

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

則剩余磁場系數有3項.

其中，K為感應磁場系數矩陣，有9項系數.

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

上述模型是傳統的飛機干擾磁場模型，共有21項模型系數.利用感應磁場系數矩陣K的對稱性[3]，以及方向余弦cos2X0+cos2Y0+cos2Z0=1的性質，可將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項系數.飛機干擾磁場的補償過程是，利用測量數據估計線性方程組(1)中的16項未知系數，再根據估計出的系數計算飛機干擾磁場，從磁場測量值中減去飛機干擾磁場，得到補償后的磁場值.

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

(2)

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

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

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

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

3 算法驗證

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

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

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

表1 LS估計、嶺估計和WTRE的計算結果

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

由表1的計算結果可以得出以下結論：

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

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

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

圖3 磁偶極子目標信號

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

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

(b)LS估計的補償結果

(c)嶺估計的補償結果

(d)WTRE的補償結果

4 結論

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

[1] 熊盛青. 我國航空重磁勘探技術現狀與發展趨勢[J].地球物理學進展， 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] 陳希孺， 王松桂. 近代回歸分析-原理方法及應用[M].合肥:安徽教育出版社， 1987:294-320.

[6] 劉莉. NQD樣本下部分線性模型中估計的強相合性[J].湖北大學學報：自然科學版， 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-曲線法確定嶺估計中的嶺參數[J].武漢大學學報：信息科學版， 2004， 29(3):235-238.

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