基于功率譜估計的航磁補償優(yōu)化處理方法*

吳佩霖，張群英，陳路昭，費春嬌，朱萬華，方廣有

(1 中國科學(xué)院電子學(xué)研究所電磁輻射與探測技術(shù)重點實驗室， 北京 100190； 2 中國科學(xué)院大學(xué)， 北京 100049) (2017年4月10日收稿； 2017年7月13日收修改稿)

航空磁法測量是磁異常檢測的重要方法之一，廣泛應(yīng)用于地球物理和軍事領(lǐng)域，尤其在礦產(chǎn)資源勘查和航空反潛中獲得成功應(yīng)用[1-3]。航空磁法測量具有以下優(yōu)點：

1) 高效：航空磁法測量采用飛機搭載磁力儀實現(xiàn)空間磁場的測量，可以快速實現(xiàn)大面積的勘探作業(yè)。

2) 安全：通過機載平臺進行空間磁場測量，可以避免人員進入高危區(qū)域，降低作業(yè)風(fēng)險。

3) 可靠：相比于地面磁場測量方式，航空磁法測量可以避免地面局部磁異常體及地表噪聲的影響，提高數(shù)據(jù)的可靠性。

在航空磁法測量過程中，設(shè)備搭載平臺通常是固定翼、直升機、無人機等飛行器。由于搭載平臺含有鐵磁性材料，在航磁測量時，必然會對飛機搭載的磁力儀產(chǎn)生干擾磁場，影響數(shù)據(jù)的準(zhǔn)確性，因此在進行航空磁法測量時，有效地補償飛機產(chǎn)生的磁干擾具有重要的意義。

航磁補償技術(shù)起源于二戰(zhàn)時期磁異常反潛技術(shù)，為提高對潛艇的檢測能力，Tolles[4-5]提出航磁補償模型，將飛機的干擾磁場分成3部分：恒定干擾磁場、感應(yīng)干擾磁場和渦流干擾磁場，并提出硬補償方案。其后，Leliak[6]進一步完善航磁補償模型，針對模型中的復(fù)共線性問題，給出正弦機動飛行的方式實現(xiàn)標(biāo)定飛行，顯著地提升了補償器的補償效果。隨著電子技術(shù)的發(fā)展，傳統(tǒng)的硬補償方法逐漸發(fā)展為軟補償方法，Leach[7]將航磁補償模型作為最小二乘問題，給出航磁補償問題的軟補償方法。21世紀(jì)，Nelson[8-11]在該領(lǐng)域進行了大量研究，包括補償后剩余磁場的功率譜分析、地面補償方法以及飛機噪聲源分析等。其后Noriega[12-14]對航磁補償結(jié)果的穩(wěn)定性進行研究，提出有效評估標(biāo)定飛行結(jié)果的數(shù)值方法，引入評估補償系數(shù)泛化能力的交叉標(biāo)定系數(shù)。在國內(nèi)相關(guān)領(lǐng)域，也有眾多學(xué)者開展相關(guān)研究并且有豐富的成果涌現(xiàn)，相關(guān)單位主要有中船重工第715研究所、海軍工程大學(xué)、航遙中心、國防科學(xué)技術(shù)大學(xué)、哈爾濱工業(yè)大學(xué)、吉林大學(xué)、中國科學(xué)院遙感應(yīng)用研究所、中國科學(xué)院電子學(xué)研究所等[15-25]。

在航磁補償過程中，飛機首先需要進行標(biāo)定飛行，在標(biāo)定飛行時，飛機需要以一定的規(guī)范進行機動飛行，由于飛機航速不同，導(dǎo)致飛機機動飛行頻率不同，最終每次用于求解補償系數(shù)的頻帶會相應(yīng)地發(fā)生變化。因此在航磁補償過程中有必要對補償系數(shù)求解的最優(yōu)頻帶進行估計。

本文針對機動飛行導(dǎo)致磁補償系數(shù)求解頻帶發(fā)生變化的問題，提出采用功率譜估計的方法實現(xiàn)最優(yōu)補償系數(shù)求解頻帶的估計，結(jié)合航磁補償模型，達到提高補償效果的目的。該方法成功應(yīng)用于搭載固定伸桿的直升機航磁勘探系統(tǒng)，并通過野外標(biāo)定飛行實驗驗證該方法的有效性。

1 航磁補償

1.1 航磁補償模型

航磁補償模型的基本原理是利用三軸磁通門磁力儀測量飛機的姿態(tài)，從而補償光泵總場測量磁力儀受到的干擾磁場。在航磁補償模型中，通過建立三軸磁通門磁力儀輸出和地磁場矢量在飛機載體坐標(biāo)系中方向余弦的關(guān)系實現(xiàn)飛機姿態(tài)的測量，二者關(guān)系可表示為

式中：u1、u2和u3是地磁場矢量和飛機載體坐標(biāo)系3個方向軸正方向夾角的余弦值；T(t)、L(t)和V(t)是三軸磁通門在3個方向上的輸出。

航磁補償模型中光泵探頭處受到的飛機磁干擾可以分成以下3類：

1) 恒定干擾磁場：飛機的制作材料中存在永磁性材料，這些材料在光泵磁力儀探頭處產(chǎn)生的干擾磁場即是恒定干擾磁場。其表達式為

HPERM(t) =c1u1+c2u2+c3u3

=A1c1+A2c2+A3c3，

(2)

式中：ci,i=1,2，3為恒定干擾磁場的補償系數(shù)；Ai,i=1,2,3為方向余弦組成的特征項。

2) 感應(yīng)干擾磁場：飛機中的鐵磁性材料會受到地磁場的磁化，產(chǎn)生的磁化磁場在光泵磁力儀探頭處產(chǎn)生的干擾磁場即是感應(yīng)干擾磁場，其表達式為

=A4c4+A5c5+A6c6+A7c7+A8c8+A9c9，

(3)

式中：He(t)為地磁場的大??；ci,i=4,…,9為感應(yīng)干擾磁場的補償系數(shù)；Ai,i=4,…,9為對應(yīng)的方向余弦和He(t)組成的特征項。

3) 渦流干擾磁場：飛機在地磁場中做切割地磁場磁力線運動時，飛機的金屬蒙皮等金屬部件會感生渦旋電流，該渦旋電流會進一步產(chǎn)生感生磁場，進而對光泵探頭處產(chǎn)生干擾磁場，其表達式為

=A10c10+A11c11+A12c12+A13c13+

A14c14+A15c15+A16c16+A17c17+A18c18，

(4)

飛機在光泵探頭處的干擾磁場可以表示為

Hd(t) =HPERM(t)+HIND(t)+HEDDY(t)

對應(yīng)的矩陣表達式可以表示為

Hd=AC+z,

(6)

式(6)中補償系數(shù)的最小二乘解可以表示為

C=(ATA)-1ATHd.

(7)

由于航磁補償模型中存在一定程度的復(fù)共線性，因此采用嶺回歸算法來實現(xiàn)補償系數(shù)C的求解，將提升補償參數(shù)求解的準(zhǔn)確性，其表達式為

Ck=(ATA+kI)-1ATHd.

(8)

通過標(biāo)定飛行獲取的數(shù)據(jù)求取補償系數(shù)Ck，進而對飛機的干擾磁場實現(xiàn)補償，其中補償系數(shù)Ck作用于飛機機動飛行存在的頻帶。在標(biāo)準(zhǔn)標(biāo)定飛行過程中，由于在不同頻帶上機動飛行引起的干擾磁場強度不同，而與機動飛行不相關(guān)的噪聲卻大致相同，因此計算獲得的補償系數(shù)Ck跟飛機的機動飛行相關(guān)，不同的頻帶會影響補償系數(shù)Ck計算的準(zhǔn)確性。在求解補償參數(shù)之前，需要先預(yù)估機動飛行的頻率范圍，從而獲取最優(yōu)的補償系數(shù)求解頻帶。

1.2 航磁補償模型頻帶估計

航磁補償模型中地磁場矢量在飛機載體坐標(biāo)系的方向角的余弦cosX(t)，cosY(t)和cosZ(t)會隨著飛機機動飛行的動作發(fā)生周期性的變化。為了方便書寫，將3個方向角的余弦簡記為cosX，cosY和cosZ，因此可以將3個方向角的余弦假定為下式：

式中：d1，d2和d3分別為常數(shù)項；k1，k2和k3為方向角的變化幅度；θ1，θ2和θ3為姿態(tài)角間的相差；f為飛機的機動頻率。

式(3)感應(yīng)干擾磁場可以變化為

(10)

式(4)渦流干擾磁場可以變化為

(11)

將式(2)，式(10)和式(11)代入式(5)中，可知飛機干擾磁場的頻率分量包含方向余弦cosX，cosY和cosZ的基頻f和其二次諧頻2f。因此為了獲得最優(yōu)的補償頻帶，從而有效地求取磁干擾補償系數(shù)，需要對地磁場矢量在飛機載體坐標(biāo)系中的方向余弦進行頻譜分析。本文采用功率譜估計的方法來實現(xiàn)最優(yōu)補償系數(shù)求解頻帶的估計。

2 功率譜估計

功率譜估計可以獲得信號的功率譜密度，用來實現(xiàn)對信號各頻率成分的估計[26]。常見的功率譜估計方法可以分為非參數(shù)功率譜估計和參數(shù)功率譜估計，針對不同類型的信號，不同的方法效果有所不同。

2.1 Welch功率譜估計

Welch法屬于非參數(shù)譜估計方法，其特點在于兩個方面：1)允許子序列xi(n)相互疊加；2)對各子序列可以進行加窗處理，產(chǎn)生的是平均的修正周期圖。

假設(shè)相鄰的子序列偏移D個點，各子序列長度為L，則第i個序列是

xi(n)=x(n+iD),n=0,1,…,L-1,

(12)

其中xi(n)和xi+1(n)的重疊(L-D)個點，若整個樣本點由K個子序列覆蓋，則樣本點長度N和K的關(guān)系可表示為

N=L+D(K-1).

(13)

在使用Welch功率譜估計時，為改善旁瓣較大導(dǎo)致的譜失真，對子序列信號進行加窗處理，窗函數(shù)為w(n)，信號x(n)的Welch譜估計可表達為

Welch譜估計是功率譜的漸進無偏估計。

2.2 AR模型功率譜估計

AR模型功率譜估計屬于參數(shù)化的功率譜估計方法，該方法假定信號是一個均方誤差為σ2的高斯白噪聲序列u(n)激勵一個線性時不變系統(tǒng)H(z)所得到的。該系統(tǒng)H(z)即是對應(yīng)的AR模型。輸入信號u(n)和分析信號x(n)間關(guān)系滿足

所分析信號x(n)的功率譜估計為

其中H(z)的表達式為

航磁補償中通過功率譜估計能有效地獲得飛機干擾場頻帶，進而獲得最優(yōu)的補償系數(shù)Ck。下面采用兩種不同的功率譜估計方法實現(xiàn)最優(yōu)補償系數(shù)求解頻帶的估計，并實現(xiàn)飛機干擾磁場的補償和補償質(zhì)量評估。

3 野外飛行實驗

3.1 飛行實驗設(shè)計

航磁補償實驗采用直升機搭載固定伸桿來實現(xiàn)，搭載固定伸桿的直升機在3 000 m高空進行標(biāo)定飛行，完成標(biāo)定飛行后，直升機進一步完成驗證飛行，以便對數(shù)據(jù)的質(zhì)量進行評估。搭載設(shè)備的直升機見圖1(a)所示,GPS記錄的飛行航跡如圖1(b)所示。

圖1 野外直升機飛行實驗Fig.1 Helicopter experiment

3.2 最優(yōu)補償頻帶估計

通過與飛機固聯(lián)的三軸磁通門實現(xiàn)飛機姿態(tài)的測量，進而獲得地磁場矢量在飛機載體坐標(biāo)系中的方向余弦，在飛機機動飛行過程中，該方向余弦與飛機機動保持一致的變化，因此采用該方向余弦作為譜估計的對象能獲得飛機機動飛行的頻帶范圍。圖2為磁通門的輸出計算得到的方向余弦信號。

圖2 地磁場在載體坐標(biāo)系中的方向余弦Fig.2 Direction cosine of the geomagnetic field in aircraft coordinate system

直接傅里葉變換后信號頻域的結(jié)果見圖3所示。

圖3 方向余弦的FFTFig.3 FFT of direction cosine

從圖3可見，直接對信號做傅里葉變換，雖然飛機的機動頻率也以尖峰的形式出現(xiàn)在頻譜上，但是由于傅里葉變換結(jié)果中存在較多的尖峰，在沒有任何先驗知識的情況下，難以分辨出哪個尖峰對應(yīng)飛機的機動頻率，因此直接FFT變換難以正確區(qū)分最優(yōu)補償系數(shù)求解頻帶。

Welch功率譜估計的結(jié)果見圖4(a)所示；AR模型功率譜估計的結(jié)果見圖4(b)所示。

從圖4(a)和圖4(b)可見，方向余弦的基頻的最大值不超過0.1 Hz(峰值點位于0.078 Hz)，因此二次諧頻不超過0.2 Hz，可見Welch功率譜估計和AR模型功率譜估計均可以很好地將最優(yōu)補償系數(shù)求解頻帶提取出來，兩種不同的方法具有高度的一致性，結(jié)果可互相驗證。

3.3 航磁補償質(zhì)量評估

圖4 方向余弦功率譜估計Fig.4 Power spectrum estimation of direction cosine

圖5 不同截止頻率濾波器預(yù)處理后補償結(jié)果Fig.5 Compensation results of filters with different cut-off frequencies

在飛機干擾磁場補償過程中，采用包含最優(yōu)補償系數(shù)求解頻帶和不包含最優(yōu)補償系數(shù)求解頻帶的兩種不同的高通濾波器對數(shù)據(jù)進行預(yù)處理，并比較最終補償效果。其中包含最優(yōu)補償系數(shù)求解頻帶的濾波器的通帶應(yīng)該低于方向余弦的基頻，因此選用截止頻率為0.065 Hz的高通濾波器；不包含最優(yōu)補償系數(shù)求解頻帶的濾波器的通帶應(yīng)該高于二次諧頻，因此選用截止頻率為0.22 Hz的高通濾波器。采用式(8)從標(biāo)定飛行中獲得飛機磁干擾補償系數(shù)，并實現(xiàn)對標(biāo)定飛行的補償。補償結(jié)果如圖5所示。其中圖5(a)和圖5(b)為截止頻率為0.065 Hz的高通濾波器預(yù)處理后再進行補償?shù)慕Y(jié)果；圖5(c)和圖5(d)為截止頻率為0.22 Hz的高通濾波器預(yù)處理后再進行補償?shù)慕Y(jié)果。

圖5采用雙縱軸顯示，虛線為補償前的飛機干擾磁場，對應(yīng)左縱軸；實線為補償后的剩余磁場對應(yīng)右縱軸?？梢妰煞N不同截止帶寬的濾波器均可以對飛機干擾磁場實現(xiàn)較好的預(yù)處理，直觀上難以直接分辨哪種濾波器處理后的補償結(jié)果最優(yōu)，因此采用定量的方式對補償結(jié)果進行評估。

航磁補償中的評估準(zhǔn)則采用補償提升比IR，其表達式為

式中：σu為補償前信號的標(biāo)準(zhǔn)差，σc為補償后信號的標(biāo)準(zhǔn)差。兩種不同截止頻率的濾波器處理后的補償信號的提升比見表1所示。其中濾波器1為截止頻率為0.065Hz的濾波器，濾波器2為截止頻率為0.22Hz的濾波器。

表1 標(biāo)準(zhǔn)差與補償提升比Table 1 Standard deviation and IR

表中前3行為標(biāo)定飛行數(shù)據(jù)結(jié)果，后3行為驗證飛行數(shù)據(jù)結(jié)果。可見采用截止頻率為0.065 Hz的高通濾波器處理后獲得最優(yōu)補償系數(shù)再進行補償，對于標(biāo)定飛行和驗證飛行均可以獲得更好的補償效果。

4 結(jié)論

本文提出采用功率譜估計的方法對航磁補償最優(yōu)補償系數(shù)求解頻帶進行估計，并設(shè)計標(biāo)定飛行實驗，對航磁補償結(jié)果進行驗證。實驗結(jié)果表明，參數(shù)化譜估計和非參數(shù)化譜估計對飛機機動飛行頻率的估計均優(yōu)于傅里葉變換。在航磁補償數(shù)據(jù)預(yù)處理階段，采用功率譜估計獲得最優(yōu)補償系數(shù)求解頻帶后計算得到的最優(yōu)補償系數(shù)，能有效地提高航磁補償?shù)难a償提升比。