一種基于Savitzky-Golay濾波器的航磁校準方法

竇振家, 葛 林, 韓 琦

(1 中國船舶集團有限公司 第七一五研究所，杭州 310023; 2 哈爾濱工業大學 計算機科學與技術學院， 哈爾濱 150001)

0 引 言

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

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

本文提出了一種基于Savitzky-Golay濾波器的航磁校準方法，通過級聯兩個Savitzky-Golay濾波器實現帶通濾波。實驗結果表明，與目前廣泛采用的FIR帶通濾波器相比，本文所提出的方法能夠取得更好的補償結果。

1 Savitzky-Golay濾波器基本原理及性質

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

(1)

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

(2)

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

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

(3)

式中，y[0]為x[0]的濾波結果。在濾波時，寬度為2M+1的窗口沿逐個樣本點向前滑動，并計算對應窗口中心位置處信號樣本點的濾波結果。實際上，在滑動窗口的過程中不必實時計算多項式系數ak。這是因為窗口中心位置樣本點的濾波結果為a0可表示為式(4)：

(4)

即窗口中心位置樣本點的濾波結果是窗口中所有樣本點的一個線性組合，相應的系數為h0,m，這組系數可根據多項式的階數和窗口寬度來預先計算得到。

已知多項式的階數N及其滑動窗口的寬度2M+1，構造矩陣AT：

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

不難看出，多項式的階數N以及滑動窗口的寬度2M+1對濾波結果具有重要影響。可將Savitzky-Golay濾波器近似為一個低通濾波器[2]。當滑動窗口的寬度一定時，階數N越高，濾波器的截止頻率越大。當階數一定時，滑動窗口的寬度越大，濾波器的截止頻率越小。歸一化截止頻率與多項式及滑動窗口寬度之間的近似關系可以表示為式(5)：

(5)

當M<10時，歸一化截止頻率與多項式及滑動窗口寬度之間的關系較為復雜，式(5)不再成立。可通過組合多組N和M來得到某一通帶內的信號，類似于經典濾波方法中級聯多個濾波器。由于fc是歸一化頻率，在已知采樣率fs時不難得到實際頻率的估計值f，式(6)：

(6)

本文將通過“級聯”兩個Savitzky-Golay濾波器實現帶通濾波器。

2 算法流程及分析

2.1 算法流程

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

表1 算法流程

Step 3令N=2，Ml=134，計算相應Savitzky-Golay濾波器的濾波系數向量ωl；同理，令N=2，Mr=4，計算相應的濾波系數向量ωr；

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

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

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

Step 7返回β。

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

(7)

(8)

(9)

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

(10)

HI=aβT,

(11)

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

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

在Step 5中，由于hf是hr和hl的差，Af是Ar和Al的差，根據前文分析可知hf和Af中主要包含與航空平臺磁干擾有關的信號分量。

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

2.2 算法分析

不難看出，算法的核心在于設置Ml和Mr的最佳取值。本文的出發點是有效去除信號中的地磁場分量，并對信號中的噪聲進行抑制，以提取出更準確的航空平臺磁干擾分量。

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

3 測試結果及分析

為了驗證所提方法的有效性，選取兩組校準飛行期間采集的數據用于對所提方法進行測試，第1組用于估計航空平臺磁干擾補償系數，第2組用于驗證航磁干擾補償效果。需要說明的是，這里兩種方法均采用了考慮地磁梯度因素影響的擴展Tolles-Lawson模型[3]。

本文將0.04-1.0 Hz通帶內剩余磁干擾的標準差作為評價指標。

第2組測試數據補償前總場信號波形如圖1所示，對其通帶濾波后所得結果如圖2所示，通帶內磁干擾的標準差為0.192 nT，峰峰值約1.68 nT。

利用基于FIR帶通濾波器的航磁校準方法所得補償系數和利用本文方法所得補償系數對第2組測試數據總場信號分別進行補償，所得補償結果經帶通濾波后如圖3所示。為了便于分析，將基于FIR帶通濾波器的航磁補償結果沿縱軸向下平移了0.5 nT。

圖2 補償前通帶內磁干擾波形

(a) 通帶內剩余磁干擾

(b) 磁場三分量信號

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

4 結束語

本文提出了一種基于Savitzky-Golay濾波器的航磁校準方法，該方法通過對兩個Savitzky-Golay濾波器級聯來實現“帶通”濾波，從而提取與航空平臺磁干擾有關的信號分量。實驗結果表明，作為航磁校準算法的預處理環節，該方法有助于提高航空平臺磁干擾補償系數的估計精度。