張 堅(jiān),林春生,羅 青,龔沈光

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

0 引言

航空磁性探測(cè)廣泛應(yīng)用于考古、礦物勘探和水中磁性目標(biāo)探測(cè)等領(lǐng)域。但是在飛機(jī)上對(duì)地磁異常場(chǎng)進(jìn)行測(cè)量時(shí),必須對(duì)飛機(jī)干擾磁場(chǎng)進(jìn)行補(bǔ)償。這就需要利用測(cè)量數(shù)據(jù)求解飛機(jī)干擾磁場(chǎng)模型系數(shù),使用所求系數(shù)計(jì)算干擾磁場(chǎng),最后從測(cè)量的總磁場(chǎng)中減去干擾磁場(chǎng),完成飛機(jī)磁補(bǔ)償。因此,飛機(jī)干擾磁場(chǎng)模型的求解是飛機(jī)磁補(bǔ)償中的關(guān)鍵技術(shù)。

美國(guó)學(xué)者Tolles和Lawson在20世紀(jì)40年代最早對(duì)飛機(jī)干擾磁場(chǎng)的數(shù)學(xué)模型進(jìn)行了研究,建立了經(jīng)典的Tolles和Lawson模型[1],但沒(méi)有給出模型的求解方法。其后,Leliak和Bickel均對(duì)該模型進(jìn)行了改進(jìn)[2-3],但其求解精度仍然不高。文獻(xiàn)[4]和文獻(xiàn)[5]分別研究了基于神經(jīng)網(wǎng)絡(luò)和截?cái)嗥娈愔捣纸獾娘w機(jī)干擾磁場(chǎng)模型求解算法。神經(jīng)網(wǎng)絡(luò)的計(jì)算量較大,收斂速度慢,不適合在實(shí)際工程應(yīng)用中進(jìn)行實(shí)時(shí)求解。截?cái)嗥娈愔捣纸夥▽?duì)奇異值截?cái)辔恢帽容^敏感,截?cái)辔恢貌缓侠頃?huì)使均方誤差迅速增大。同時(shí),這兩種方法都沒(méi)有利用求解結(jié)果進(jìn)行飛機(jī)磁補(bǔ)償實(shí)驗(yàn),其求解效果也就無(wú)法從根本上進(jìn)行檢驗(yàn)。

飛機(jī)干擾磁場(chǎng)模型屬線性回歸模型,并且存在較強(qiáng)的復(fù)共線性。線性模型求解常用方法有最小二乘估計(jì)(LS估計(jì))、截?cái)嗥娈愔捣纸夂蛶X估計(jì)等[6]。截?cái)嗥娈愔捣纸獾牟蛔阋言谇懊嫣岬?LS估計(jì)運(yùn)算量小,但在模型存在復(fù)共線性時(shí)均方誤差會(huì)很大。針對(duì)上述問(wèn)題,本文提出了將嶺估計(jì)的預(yù)測(cè)殘差平方和(PRESS)法引入飛機(jī)干擾磁場(chǎng)模型求解。

1 飛機(jī)干擾磁場(chǎng)模型與嶺估計(jì)

1.1 飛機(jī)干擾磁場(chǎng)模型

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

圖1 飛機(jī)干擾磁場(chǎng)模型坐標(biāo)系Fig.1 Coordinates of aircraft magnetic interference field model

建模時(shí)首先分別建立恒定磁場(chǎng)、感應(yīng)磁場(chǎng)和渦流磁場(chǎng)的數(shù)學(xué)模型,將三者疊加后再進(jìn)行合并與化簡(jiǎn),最終可以得到總的飛機(jī)干擾磁場(chǎng)ΔH為[1]:

在式(1)兩邊同時(shí)加上地磁場(chǎng)值He后,方程左邊就變成了磁場(chǎng)總強(qiáng)度的測(cè)量值。利用磁探設(shè)備的測(cè)量數(shù)據(jù)和計(jì)算出來(lái)的方向余弦及其導(dǎo)數(shù),就可以求解式(1)中的16項(xiàng)未知系數(shù),完成飛機(jī)干擾磁場(chǎng)模型的參數(shù)估計(jì)。

1.2 模型復(fù)共線性分析

考慮線性回歸模型

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

對(duì)于飛機(jī)干擾磁場(chǎng)模型,X由方向余弦、方向余弦的導(dǎo)數(shù)及其乘積項(xiàng)構(gòu)成,Y由磁場(chǎng)總強(qiáng)度的測(cè)量值構(gòu)成,e為測(cè)量誤差向量。用LS估計(jì)求解模型時(shí)均方誤差(MSE)很大,原因在于設(shè)計(jì)矩陣XTX的奇異性。該矩陣產(chǎn)生奇異性的原因是:在保證飛機(jī)穩(wěn)定飛行的情況下,飛機(jī)的機(jī)動(dòng)角度不可能很大,這樣X(jué)TX每一行的對(duì)應(yīng)元素值非常接近,每行之間有比較強(qiáng)的相關(guān)性,這就導(dǎo)致飛機(jī)干擾磁場(chǎng)模型存在嚴(yán)重的復(fù)共線性。

1.3 模型參數(shù)的嶺估計(jì)

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

文獻(xiàn)[6]提出了基于預(yù)測(cè)殘差平方和(Prediction residual error square sum,PRESS)的嶺參數(shù)確定方法,其過(guò)程如下:

模型(2)中,第 i次觀測(cè)的數(shù)據(jù)點(diǎn)為(X1i,X2i,…,),記為 XTi,此時(shí)

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

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

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

本文將PRESS法引入飛機(jī)干擾磁場(chǎng)模型求解,并從航空磁性探測(cè)的實(shí)際應(yīng)用出發(fā),對(duì)PRESS的求解過(guò)程作了一定改進(jìn)。

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

選取最小預(yù)測(cè)殘差平方和意義下的最優(yōu)嶺參數(shù)后,采用該嶺參數(shù)計(jì)算飛機(jī)干擾磁場(chǎng)模型的嶺估計(jì),就完成了模型參數(shù)的求解?；赑RESS法的模型求解原理框圖如圖2所示。

在航空磁性探測(cè)的實(shí)際應(yīng)用中,求解模型系數(shù)的目的是在后續(xù)的飛行過(guò)程中,利用所求系數(shù)預(yù)測(cè)飛機(jī)干擾磁場(chǎng)的估計(jì)值,再?gòu)目偞艌?chǎng)測(cè)量值中減去干擾磁場(chǎng)的估計(jì)值,完成飛機(jī)磁補(bǔ)償。其中,為第n+1個(gè)觀測(cè)點(diǎn)處磁場(chǎng)總強(qiáng)度的實(shí)際測(cè)量值,為利用前n個(gè)測(cè)量值對(duì)其預(yù)測(cè)的結(jié)果。預(yù)測(cè)殘差平方和越小,說(shuō)明利用所求模型系數(shù)計(jì)算的干擾磁場(chǎng)與實(shí)際干擾磁場(chǎng)越接近,模型求解效果越好,即選用的嶺參數(shù)k越好。因此,PRESS法符合航空磁性探測(cè)的實(shí)際物理意義。

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

3 模型求解仿真驗(yàn)證

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

圖3 飛機(jī)干擾磁場(chǎng)信號(hào)仿真Fig.3 Simulation of aircraft magnetic interference field signal

衡量矩陣病態(tài)程度的方法之一是求它的條件數(shù)。計(jì)算得到法矩陣XTX的條件數(shù)為2.37×1010,呈嚴(yán)重病態(tài)性。采用LS估計(jì)和嶺估計(jì)的計(jì)算結(jié)果如表1所示,嶺估計(jì)分別采用嶺跡法、GCV法和PRESS法進(jìn)行k值估計(jì),結(jié)果分別為:0.087,不收斂,0.316。采用均方誤差(MSE)作為求解效果評(píng)判標(biāo)準(zhǔn)。

表1 LS估計(jì)和嶺估計(jì)的計(jì)算結(jié)果Tab.1 Calculation result of LS estimation and ridge estimation

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

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

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

4 飛機(jī)磁補(bǔ)償實(shí)驗(yàn)驗(yàn)證

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

圖4 磁偶極子目標(biāo)信號(hào)Fig.4 Target signal of magnetic dipole

圖5 目標(biāo)與飛機(jī)干擾磁場(chǎng)的混合信號(hào)及補(bǔ)償后的結(jié)果Fig.5 Mixed signal of target and aircraft magnetic interference field and the result after compensation

5 結(jié)論

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

[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]龐學(xué)亮,林春生,張寧.飛機(jī)磁場(chǎng)模型系數(shù)的截?cái)嗥娈愔捣纸夥ü烙?jì)[J].探測(cè)與控制學(xué)報(bào),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]張金槐.線性模型參數(shù)估計(jì)及其改進(jìn)[M].第2版.長(zhǎng)沙:國(guó)防科技大學(xué)出版社,1999.

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

[8]葉仁玉,曾建軍.特殊數(shù)據(jù)的廣義嶺估計(jì)的迭代算法及其實(shí)現(xiàn)[J].合肥工業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版),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ì)的嶺參數(shù)確定方法[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào),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.