基于預測殘差平方和的飛機干擾磁場模型求解

張 堅,林春生,羅 青,龔沈光

(1.海軍工程大學兵器工程系,湖北武漢 430033;2.海軍92038部隊裝備處,山東青島 266041)

0 引言

航空磁性探測廣泛應用于考古、礦物勘探和水中磁性目標探測等領域。但是在飛機上對地磁異常場進行測量時,必須對飛機干擾磁場進行補償。這就需要利用測量數據求解飛機干擾磁場模型系數,使用所求系數計算干擾磁場,最后從測量的總磁場中減去干擾磁場,完成飛機磁補償。因此,飛機干擾磁場模型的求解是飛機磁補償中的關鍵技術。

美國學者Tolles和Lawson在20世紀40年代最早對飛機干擾磁場的數學模型進行了研究,建立了經典的Tolles和Lawson模型[1],但沒有給出模型的求解方法。其后,Leliak和Bickel均對該模型進行了改進[2-3],但其求解精度仍然不高。文獻[4]和文獻[5]分別研究了基于神經網絡和截斷奇異值分解的飛機干擾磁場模型求解算法。神經網絡的計算量較大,收斂速度慢,不適合在實際工程應用中進行實時求解。截斷奇異值分解法對奇異值截斷位置比較敏感,截斷位置不合理會使均方誤差迅速增大。同時,這兩種方法都沒有利用求解結果進行飛機磁補償實驗,其求解效果也就無法從根本上進行檢驗。

飛機干擾磁場模型屬線性回歸模型,并且存在較強的復共線性。線性模型求解常用方法有最小二乘估計(LS估計)、截斷奇異值分解和嶺估計等[6]。截斷奇異值分解的不足已在前面提到,LS估計運算量小,但在模型存在復共線性時均方誤差會很大。針對上述問題,本文提出了將嶺估計的預測殘差平方和(PRESS)法引入飛機干擾磁場模型求解。

1 飛機干擾磁場模型與嶺估計

1.1 飛機干擾磁場模型

飛機干擾磁場包括恒定磁場、感應磁場和渦流磁場。本文用經典Tolles和Lawson模型來描述飛機干擾磁場,模型坐標系如圖1所示。圖中,機載磁探設備為坐標原點,地磁場矢量He與X,Y,Z軸的夾角為(X0,Y0,Z0),其方向余弦定義為(cosX0,cosY0,cosZ0)。

圖1 飛機干擾磁場模型坐標系Fig.1 Coordinates of aircraft magnetic interference field model

建模時首先分別建立恒定磁場、感應磁場和渦流磁場的數學模型,將三者疊加后再進行合并與化簡,最終可以得到總的飛機干擾磁場ΔH為[1]:

在式(1)兩邊同時加上地磁場值He后,方程左邊就變成了磁場總強度的測量值。利用磁探設備的測量數據和計算出來的方向余弦及其導數,就可以求解式(1)中的16項未知系數,完成飛機干擾磁場模型的參數估計。

1.2 模型復共線性分析

考慮線性回歸模型

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

對于飛機干擾磁場模型,X由方向余弦、方向余弦的導數及其乘積項構成,Y由磁場總強度的測量值構成,e為測量誤差向量。用LS估計求解模型時均方誤差(MSE)很大,原因在于設計矩陣XTX的奇異性。該矩陣產生奇異性的原因是:在保證飛機穩定飛行的情況下,飛機的機動角度不可能很大,這樣XTX每一行的對應元素值非常接近,每行之間有比較強的相關性,這就導致飛機干擾磁場模型存在嚴重的復共線性。

1.3 模型參數的嶺估計

為了克服LS估計在模型存在復共線性時的不足,統計學家們提出了許多改進方法,嶺估計就是其中一種。嶺估計由 Hoerl等人在 1970年首次提出[7],是克服模型復共線性的有效方法。對于線性模型(2),β的嶺估計定義為[7]:

文獻[6]提出了基于預測殘差平方和(Prediction residual error square sum,PRESS)的嶺參數確定方法,其過程如下:

模型(2)中,第 i次觀測的數據點為(X1i,X2i,…,),記為 XTi,此時

對于這n組數據,如果將第i組數據(Xi,yi)放在一旁,利用其余n-1組數據作線性回歸,計算在Xi處的預測值,記為

對特定的嶺參數k,可以計算出嶺估計的PRESS(k),選取使PRESS(k)最小的k值作為嶺參數,就得到了最小預測殘差平方和意義下的最優嶺參數。

2 基于預測殘差平方和(PRESS)法的模型求解

本文將PRESS法引入飛機干擾磁場模型求解,并從航空磁性探測的實際應用出發,對PRESS的求解過程作了一定改進。

利用實時采集的n組數據(Xi,yi)(i=1,…,n)作線性回歸,預測第 n+1組數據(),將處的預測值記為=。其中,為利用第1到第n組數據求出的回歸系數估計值。記=-為第n+1個觀測點處的預測殘差。同樣,利用第2到第n+1組數據來預測第n+2組數據,得到預測殘差。對時間序列中連續的L個觀測點進行預測,得到預測殘差平方和為:

選取最小預測殘差平方和意義下的最優嶺參數后,采用該嶺參數計算飛機干擾磁場模型的嶺估計,就完成了模型參數的求解。基于PRESS法的模型求解原理框圖如圖2所示。

在航空磁性探測的實際應用中,求解模型系數的目的是在后續的飛行過程中,利用所求系數預測飛機干擾磁場的估計值,再從總磁場測量值中減去干擾磁場的估計值,完成飛機磁補償。其中,為第n+1個觀測點處磁場總強度的實際測量值,為利用前n個測量值對其預測的結果。預測殘差平方和越小,說明利用所求模型系數計算的干擾磁場與實際干擾磁場越接近,模型求解效果越好,即選用的嶺參數k越好。因此,PRESS法符合航空磁性探測的實際物理意義。

圖2 基于PRESS法的模型求解原理框圖Fig.2 Principle of model solving based on PRESS method

3 模型求解仿真驗證

采用Matlab7.0軟件仿真驗證本文方法在飛機干擾磁場模型求解中的有效性。設地磁場大小為5×104nT。地磁場方向余弦可用地磁場方位角θ和俯仰角φ的函數來表示,且假設dφ/dt=1,dθ/dt=cos(t/2),φ∈ [0,2π],θ∈ [0,π]。觀測數據采樣點數為576。設定模型系數的真值如下:恒定磁場系數[30,40,80],5個感應磁場系數均為0.005,8個渦流磁場系數均為0.002。利用設定的系數真值,仿真產生飛機干擾磁場信號,如圖3所示。

圖3 飛機干擾磁場信號仿真Fig.3 Simulation of aircraft magnetic interference field signal

衡量矩陣病態程度的方法之一是求它的條件數。計算得到法矩陣XTX的條件數為2.37×1010,呈嚴重病態性。采用LS估計和嶺估計的計算結果如表1所示,嶺估計分別采用嶺跡法、GCV法和PRESS法進行k值估計,結果分別為:0.087,不收斂,0.316。采用均方誤差(MSE)作為求解效果評判標準。

表1 LS估計和嶺估計的計算結果Tab.1 Calculation result of LS estimation and ridge estimation

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

1)LS估計的均方誤差遠大于嶺估計(除了GCV法不收斂以外),這是由于法矩陣XTX的條件數很大,模型復共線性嚴重。

2)在本例中,利用GCV法估計k值時,曲線是逐漸衰減的,所以無法確定最小的k值;利用嶺跡法估計k值時,需要根據經驗來估計,帶有主觀性。幾種方法中PRESS法的均方誤差最小,這也驗證了該方法符合航空磁性探測應用的實際物理意義。

4 飛機磁補償實驗驗證

利用上節LS估計、嶺跡法和PRESS法求解出的模型系數,對仿真產生的飛機干擾磁場信號進行磁補償實驗。利用磁偶極子模型生成磁性目標信號如圖4所示。將目標信號與飛機干擾磁場疊加后的混合信號如圖5(a)所示,可以看出,飛機干擾磁場的變化幅度超過10 nT,目標信號完全被干擾磁場所淹沒,無法檢測到目標。分別采用LS估計、嶺跡法和PRESS法所求系數進行磁補償后的信號如圖5(b)—圖5(d)所示。可以看出:LS估計補償后的信噪比最低,PRESS法的信噪比最高;三種方法中只有PRESS法的補償結果可以準確檢測到目標,其他兩種方法都不同程度上會受到噪聲的干擾。這也進一步驗證了本文方法在飛機磁補償中的有效性。

圖4 磁偶極子目標信號Fig.4 Target signal of magnetic dipole

圖5 目標與飛機干擾磁場的混合信號及補償后的結果Fig.5 Mixed signal of target and aircraft magnetic interference field and the result after compensation

5 結論

本文提出了將嶺估計的預測殘差平方和(PRESS)方法引入飛機干擾磁場模型求解。PRESS法通過計算時間序列的預測殘差平方和,選取最小預測殘差平方和意義下的最優嶺參數。仿真表明:嶺估計的均方誤差明顯小于LS估計,并且PRESS法的均方誤差小于嶺跡法和GCV法。實驗表明:采用PRESS法的嶺估計補償后的信噪比最高,在飛機磁補償中與其他方法相比具有明顯的優越性。

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

[2]Leach B W.Aeromagnetic compensation as a linear regression problem[J].Information Linkage Between Applied Mathematics and Industry,1980(2):139-161.

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

[4]Peter M Williams.Aeromagnetic compensation using neural networks[J].Neural Computing &Applications,1993,(1):207-214.

[5]龐學亮,林春生,張寧.飛機磁場模型系數的截斷奇異值分解法估計[J].探測與控制學報,2009,31(5):48-51.PANG Xueliang,LIN Chunsheng,ZHANG Ning.Parameter estimation of airplane magnetic model based on TSVD[J].Journal of Detection&Control,2009,31(5):48-51.

[6]張金槐.線性模型參數估計及其改進[M].第2版.長沙:國防科技大學出版社,1999.

[7]Hoerl A E,Kennard R W.Ridge regression:biased estimation for non-orthogonal problems[J].Technometrics,1970(12):55-88.

[8]葉仁玉,曾建軍.特殊數據的廣義嶺估計的迭代算法及其實現[J].合肥工業大學學報(自然科學版),2007,30(8):1 065-1 068.YE Renyu,ZENG Jianjun.Iterative algorithm of generalized ridge estimates for some special data and its realization[J].Journal of Hefei University of Technology,2007,30(8):1 065-1 068.

[9]劉雁雨,李建偉,劉曉剛.部分嶺估計的嶺參數確定方法[J].測繪科學技術學報,2007,24(6):413-414.LIU Yanyu,LI Jianwei,LIU Xiaogang.Determination of the partial ridge parameters[J].Journal of Zhengzhou Institute of Surveying and Mapping,2007,24(6):413-414.