竇振家, 葛 林, 韓 琦

(1 中國船舶集團(tuán)有限公司 第七一五研究所，杭州 310023; 2 哈爾濱工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院， 哈爾濱 150001)

0 引 言

航空磁探是指利用機(jī)載高靈敏度磁力儀對(duì)地磁場進(jìn)行測(cè)量，并通過地磁異常現(xiàn)象發(fā)現(xiàn)感興趣的磁異常目標(biāo)，目前已廣泛用于礦產(chǎn)勘探、地質(zhì)研究、磁力反潛等領(lǐng)域，具有較高的應(yīng)用價(jià)值。然而，機(jī)載磁力儀不可避免地會(huì)受到航空平臺(tái)自身磁干擾的影響，并且航空平臺(tái)磁干擾通常和感興趣的磁異常信號(hào)所處頻帶重疊，故難以通過傳統(tǒng)濾波方法對(duì)其予以去除。工程中通常采用所謂的航磁干擾補(bǔ)償技術(shù)，即：預(yù)先建立航空平臺(tái)磁干擾的數(shù)學(xué)模型并估計(jì)其待定系數(shù)，在實(shí)際探測(cè)過程中利用航磁干擾模型及其系數(shù)實(shí)時(shí)計(jì)算航空平臺(tái)產(chǎn)生的磁干擾并將其從磁力儀輸出信號(hào)中減除，從而得到不含航空平臺(tái)磁干擾的磁場信號(hào)。其中，航空平臺(tái)磁干擾模型系數(shù)估計(jì)通常也被稱為航磁校準(zhǔn)，是高精度航磁干擾補(bǔ)償技術(shù)的關(guān)鍵。近年來，隨著高靈敏度磁力儀技術(shù)水平的不斷進(jìn)步，磁力儀精度越來越高，航磁干擾補(bǔ)償，特別是航磁校準(zhǔn)已經(jīng)成為制約航空磁探性能提升的瓶頸，并日漸引起研究人員重視。

Tolles-Lawson模型是目前應(yīng)用最廣泛的航空平臺(tái)磁干擾模型，該模型將航空平臺(tái)磁干擾建模為與飛機(jī)機(jī)體三軸與地磁場方向所成夾角的余弦值及其導(dǎo)數(shù)有關(guān)的函數(shù)[1]。為了求解Tolles-Lawson模型中的待定系數(shù)，首先需要令飛機(jī)在高空校準(zhǔn)飛行，即沿4個(gè)標(biāo)準(zhǔn)航向分別執(zhí)行俯仰、搖擺和偏航3種機(jī)動(dòng)動(dòng)作，然后根據(jù)Tolles-Lawson模型及校準(zhǔn)飛行期間采集的磁總場數(shù)據(jù)和磁場三分量數(shù)據(jù)建立線性方程組，最后通過求解線性方程組來估計(jì)Tolles-Lawson模型中的待定系數(shù)。由于校準(zhǔn)飛行期間機(jī)載磁力儀測(cè)量的磁場是地球磁場與飛機(jī)干擾磁場的疊加，在求解方程組前需要對(duì)采集到的總場信號(hào)及線性方程組方向余弦矩陣的各列進(jìn)行帶通濾波，以提取與飛機(jī)磁干擾有關(guān)的信號(hào)成分。鑒于航空平臺(tái)磁干擾與飛機(jī)機(jī)動(dòng)動(dòng)作有關(guān)，工程中通常根據(jù)校準(zhǔn)飛行機(jī)動(dòng)動(dòng)作的頻帶范圍來設(shè)置帶通濾波器。不難看出，所用帶通濾波器將影響線性方程組的性質(zhì)，從而間接影響Tolles-Lawson模型參數(shù)的估計(jì)精度，并最終影響航空平臺(tái)磁干擾補(bǔ)償效果。因此，尋找更好的帶通濾波器是提升最終航磁干擾補(bǔ)償性能的一種有效途徑。

本文提出了一種基于Savitzky-Golay濾波器的航磁校準(zhǔn)方法，通過級(jí)聯(lián)兩個(gè)Savitzky-Golay濾波器實(shí)現(xiàn)帶通濾波。實(shí)驗(yàn)結(jié)果表明，與目前廣泛采用的FIR帶通濾波器相比，本文所提出的方法能夠取得更好的補(bǔ)償結(jié)果。

1 Savitzky-Golay濾波器基本原理及性質(zhì)

Savitzky-Golay濾波器又被稱為最小二乘多項(xiàng)式平滑濾波器，由Savitzky和Golay于1964年提出。設(shè)信號(hào)x[n]，其中n=…-2,-1,0,1,2,…。若以n=0為中心，兩側(cè)分別有M個(gè)信號(hào)樣本點(diǎn)，則可以將這段信號(hào)用多項(xiàng)式來擬合，即式(1)：

(1)

并且要求該多項(xiàng)式與這段信號(hào)之間的擬合誤差取最小值，這里擬合誤差取均方差，即式(2)：

(2)

取得最小值。Savitzky-Golay濾波器取p(0)作為x[0]的濾波結(jié)果，即式(3)：

y[0]=p(0)=a0.

(3)

式中，y[0]為x[0]的濾波結(jié)果。在濾波時(shí)，寬度為2M+1的窗口沿逐個(gè)樣本點(diǎn)向前滑動(dòng)，并計(jì)算對(duì)應(yīng)窗口中心位置處信號(hào)樣本點(diǎn)的濾波結(jié)果。實(shí)際上，在滑動(dòng)窗口的過程中不必實(shí)時(shí)計(jì)算多項(xiàng)式系數(shù)ak。這是因?yàn)榇翱谥行奈恢脴颖军c(diǎn)的濾波結(jié)果為a0可表示為式(4)：

(4)

即窗口中心位置樣本點(diǎn)的濾波結(jié)果是窗口中所有樣本點(diǎn)的一個(gè)線性組合，相應(yīng)的系數(shù)為h0,m，這組系數(shù)可根據(jù)多項(xiàng)式的階數(shù)和窗口寬度來預(yù)先計(jì)算得到。

已知多項(xiàng)式的階數(shù)N及其滑動(dòng)窗口的寬度2M+1，構(gòu)造矩陣AT：

進(jìn)而計(jì)算矩陣H=(ATA)-1AT，其第0行中的2M+1個(gè)元素的取值便是h0,m，這里m的取值為同(4)。

不難看出，多項(xiàng)式的階數(shù)N以及滑動(dòng)窗口的寬度2M+1對(duì)濾波結(jié)果具有重要影響?？蓪avitzky-Golay濾波器近似為一個(gè)低通濾波器[2]。當(dāng)滑動(dòng)窗口的寬度一定時(shí)，階數(shù)N越高，濾波器的截止頻率越大。當(dāng)階數(shù)一定時(shí)，滑動(dòng)窗口的寬度越大，濾波器的截止頻率越小。歸一化截止頻率與多項(xiàng)式及滑動(dòng)窗口寬度之間的近似關(guān)系可以表示為式(5)：

(5)

當(dāng)M<10時(shí)，歸一化截止頻率與多項(xiàng)式及滑動(dòng)窗口寬度之間的關(guān)系較為復(fù)雜，式(5)不再成立?？赏ㄟ^組合多組N和M來得到某一通帶內(nèi)的信號(hào)，類似于經(jīng)典濾波方法中級(jí)聯(lián)多個(gè)濾波器。由于fc是歸一化頻率，在已知采樣率fs時(shí)不難得到實(shí)際頻率的估計(jì)值f，式(6)：

(6)

本文將通過“級(jí)聯(lián)”兩個(gè)Savitzky-Golay濾波器實(shí)現(xiàn)帶通濾波器。

2 算法流程及分析

2.1 算法流程

本文所提方法的算法流程見表1。算法輸入是校準(zhǔn)飛行期間利用標(biāo)量磁力儀測(cè)量到的磁總場信號(hào)和利用三分量磁力儀測(cè)量到的磁場三分量信號(hào)，采樣率為10 Hz。算法輸出為Tolles-Lawson模型的系數(shù)向量，即航磁干擾補(bǔ)償系數(shù)。

表1 算法流程

Step 3令N=2，Ml=134，計(jì)算相應(yīng)Savitzky-Golay濾波器的濾波系數(shù)向量ωl；同理，令N=2，Mr=4，計(jì)算相應(yīng)的濾波系數(shù)向量ωr；

Step 4根據(jù)式(4) ，利用濾波器系數(shù)向量ωl對(duì)總場信號(hào)h和方向余弦矩陣A的各列進(jìn)行濾波，分別得到hl和Al，利用濾波器系數(shù)向量ωr對(duì)總場信號(hào)h和方向余弦矩陣A的各列進(jìn)行濾波，分別得到hr和Ar；

Step 5令hf=hr-hl，Af=Ar-Al；

Step 6利用遞推最小二乘法求解線性方程組hf=Afβ；

Step 7返回β。

在Step 1中，方向余弦的計(jì)算公式為(7)～(9)：

(7)

(8)

(9)

在Step 2中，Tolles-Lawson模型如式(10)所示：

(10)

HI=aβT,

(11)

在Step 3中，考慮到在校準(zhǔn)飛行期間各個(gè)標(biāo)準(zhǔn)航向上地磁場總體變化趨勢(shì)是線性的，故令N=2即用一次函數(shù)來擬合信號(hào)。當(dāng)Ml=134時(shí)根據(jù)式(6)可知濾波器系數(shù)向量ωl的等效截止頻率約為0.001 Hz，可以濾除除地磁場以外的信號(hào)成分；當(dāng)Mr=4時(shí)雖然不能根據(jù)式(6)估計(jì)等效截止頻率，但由于Savitzky-Golay本質(zhì)上是低通濾波器，故濾波系數(shù)向量ωr的主要目的是對(duì)信號(hào)進(jìn)行降噪。

在Step 4中，利用ωl對(duì)總場信號(hào)h和方向余弦矩陣A的各列進(jìn)行濾波后所得結(jié)果hl和Al主要包含地磁場成分，利用ωr總場信號(hào)h和方向余弦矩陣A的各列進(jìn)行濾波后所得結(jié)果hr和Ar中的噪聲得到了抑制，且主要包含地磁場和航空平臺(tái)磁干擾成分。

在Step 5中，由于hf是hr和hl的差，Af是Ar和Al的差，根據(jù)前文分析可知hf和Af中主要包含與航空平臺(tái)磁干擾有關(guān)的信號(hào)分量。

在Step 6中，建立線性方程組并利用遞推最小二乘法求解Tolles-Lawson模型系數(shù)向量β。

2.2 算法分析

不難看出，算法的核心在于設(shè)置Ml和Mr的最佳取值。本文的出發(fā)點(diǎn)是有效去除信號(hào)中的地磁場分量，并對(duì)信號(hào)中的噪聲進(jìn)行抑制，以提取出更準(zhǔn)確的航空平臺(tái)磁干擾分量。

同利用FIR帶通濾波器提取航空平臺(tái)磁干擾的方式相比，由于采用了2階多項(xiàng)式對(duì)信號(hào)進(jìn)行擬合，并以“相減”的方式提取航空平臺(tái)磁干擾分量，本文所提出的方法能夠更好的保持航空平臺(tái)磁干擾信號(hào)的原有波形特征。這將有助于提高Tolles-Lawson模型系數(shù)的估計(jì)精度。

3 測(cè)試結(jié)果及分析

為了驗(yàn)證所提方法的有效性，選取兩組校準(zhǔn)飛行期間采集的數(shù)據(jù)用于對(duì)所提方法進(jìn)行測(cè)試，第1組用于估計(jì)航空平臺(tái)磁干擾補(bǔ)償系數(shù)，第2組用于驗(yàn)證航磁干擾補(bǔ)償效果。需要說明的是，這里兩種方法均采用了考慮地磁梯度因素影響的擴(kuò)展Tolles-Lawson模型[3]。

本文將0.04-1.0 Hz通帶內(nèi)剩余磁干擾的標(biāo)準(zhǔn)差作為評(píng)價(jià)指標(biāo)。

第2組測(cè)試數(shù)據(jù)補(bǔ)償前總場信號(hào)波形如圖1所示，對(duì)其通帶濾波后所得結(jié)果如圖2所示，通帶內(nèi)磁干擾的標(biāo)準(zhǔn)差為0.192 nT，峰峰值約1.68 nT。

利用基于FIR帶通濾波器的航磁校準(zhǔn)方法所得補(bǔ)償系數(shù)和利用本文方法所得補(bǔ)償系數(shù)對(duì)第2組測(cè)試數(shù)據(jù)總場信號(hào)分別進(jìn)行補(bǔ)償，所得補(bǔ)償結(jié)果經(jīng)帶通濾波后如圖3所示。為了便于分析，將基于FIR帶通濾波器的航磁補(bǔ)償結(jié)果沿縱軸向下平移了0.5 nT。

圖2 補(bǔ)償前通帶內(nèi)磁干擾波形

(a) 通帶內(nèi)剩余磁干擾

(b) 磁場三分量信號(hào)

從圖3中可以看出，基于Savitzky-Golay濾波器的航磁補(bǔ)償結(jié)果通帶內(nèi)剩余磁干擾波形的波動(dòng)幅度明顯小于基于FIR濾波器的航磁補(bǔ)償結(jié)果，前者通帶內(nèi)剩余磁干擾標(biāo)準(zhǔn)差為0.03 nT，改善比(即補(bǔ)償前和補(bǔ)償后通帶內(nèi)剩余磁干擾標(biāo)準(zhǔn)差之比，越大表示補(bǔ)償效果越好[4])為6.4，后者通帶內(nèi)剩余磁干擾標(biāo)準(zhǔn)差約為0.06 nT，改善比為3.2。由此可見本文所提方法的有效性。

4 結(jié)束語

本文提出了一種基于Savitzky-Golay濾波器的航磁校準(zhǔn)方法，該方法通過對(duì)兩個(gè)Savitzky-Golay濾波器級(jí)聯(lián)來實(shí)現(xiàn)“帶通”濾波，從而提取與航空平臺(tái)磁干擾有關(guān)的信號(hào)分量。實(shí)驗(yàn)結(jié)果表明，作為航磁校準(zhǔn)算法的預(yù)處理環(huán)節(jié)，該方法有助于提高航空平臺(tái)磁干擾補(bǔ)償系數(shù)的估計(jì)精度。