潘惠坤, 李勝, 錢晨, 季蔡娟, 徐振
(南京理工大學自動化學院, 210094, 南京)
地球磁場是地球系統的基本物理場之一,因為其無源、穩定且與地理位置有關的優點,利用其進行導航的地磁匹配導航技術受到了學者們的關注[1-4]。與其他輔助導航技術相比,地磁匹配導航技術具有無源自主、誤差不隨時間累積、全天候、隱蔽性強、連續導航、適用范圍廣的優點[5]。獲取高精度的地磁數據是地磁匹配導航的關鍵之一,決定著導航的精度。
三軸磁傳感器進行地磁場強度測量時,磁場數據不可避免的會包含誤差和異常值。誤差可以分為磁測系統誤差和外界磁干擾兩類。磁測系統誤差包括三軸磁傳感器自身的誤差(三軸磁傳感器在制造過程中的零偏誤差、非正交誤差、刻度因子誤差等)以及磁測系統的安裝誤差[6],外界磁干擾包括傳感器周圍的載體磁場及磁日變磁場等干擾信息。這些誤差使得磁傳感器測得的數據無法直接用來進行地磁匹配,需要對測量值進行校準。因此,需要研究高精度的磁傳感器誤差補償算法。
常用的磁傳感器誤差校準方法分為輔助矢量校正方法和獨立標量校正方法兩類。輔助矢量校正方法通過與已知的磁矢量場比較從而對磁傳感器數據進行輔助校正參考[7-8]。然而,在實際應用中,難以獲取高精度的已知磁場矢量。獨立標量校正法則可以避免該缺點,通過三軸磁傳感器在恒定磁場內旋轉,以合成總場是固定值作為約束條件對其進行校正。獨立標量校正方法因其在實際環境中易操作的優點吸引了大量學者的研究,典型代表為基于橢球假設的磁補償方法[9-10]。文獻[11]中提出了基于最小二乘法的誤差標定方法并應用于磁航向誤差補償中,測試結果滿足精度要求;文獻[12]利用自適應最小二乘估計法解決了橢球體擬合問題,取得了良好的校正效果;文獻[13]設計了一種基于自適應遺傳算法的空間橢球磁強計校準方法,能夠同時兼顧三軸磁傳感器誤差及其安裝誤差;文獻[14]提出了基于遞推最小二乘法的三軸磁傳感器誤差在線自校正方法。基于橢球假設的磁補償方法是最常見的磁傳感器標定和補償方法,并廣泛應用于實際測量中。
這些算法均假設三軸磁傳感器測量時不包含誤差較大的采樣點或異常值。然而,在實際測量中,該假設不成立,并會對最終精度造成較大影響。M估計法是一種有效提高系統魯棒性的技術[15]。通過構造權重代價函數,獲得與采樣點殘差相對應的權重值,降低甚至消除異常點造成的影響,從而增強系統魯棒性提高精度[16]。
基于這一思路,本文提出了M估計補充策略下的三軸磁傳感器誤差補償算法。在最小二乘法的橢球擬合算法的基礎上,引入M估計的思想,根據其擬合的殘差數據,構造Huber權重目標函數,確定各個采樣點的權重,為偏離較大的采樣點設置低權重,為偏離較小的點設置高權重,對數據進行靈活處理,從而增強系統魯棒性,降低異常點對擬合的影響,最終提高地磁測量的精度。
地磁傳感器誤差主要由載體磁場誤差和自身誤差組成。
載體磁場誤差是由磁傳感器周圍的各種磁性材料造成的。載體磁場誤差是地磁場測量中所特有的一種誤差,可以分為硬磁誤差、軟磁誤差以及隨機磁場誤差共3類主要誤差。
三軸磁傳感器在實際使用過程中,存在安裝差異,即自身誤差。傳感器誤差來源主要有3類:三軸靈敏度不一致造成的靈敏度誤差、傳感器三軸未完全正交造成的非正交誤差以及傳感器的零位偏移誤差。自身誤差屬于機械誤差,比較固定,出廠后不易改變。
綜合考慮載體磁場誤差和自身誤差,參考文獻[17]建立誤差模型
Hm=CsiCn(CsHe+bh)+b0+ε
(1)
將式(1)簡化可得
Hm=CHe+b+ε
(2)
式中:C=CsiCnCs;b=CsiCnbh+b0。顯然,C是可逆矩陣,因此可以得到地磁場誤差校正模型
He=C-1(Hm-b-ε)
(3)
在求出校準矩陣C-1和偏移量b的基礎上,當獲得磁傳感器量測數據后,通過式(3)可以計算出當地的真實地磁強度,式中ε可以通過合理的硬件技術手段減小,因此計算過程中可以忽略。
在磁場強度恒定的區域內對磁傳感器進行任意旋轉,測量出的磁場應該是一個定值,即磁場矢量的軌跡應該是一個正球體,可表達為
‖He‖2=Cconst
(4)
式中Cconst代表磁場強度的定值。結合式(2)(4)得
‖He‖2=HeTHe=
(Hm-b)T(C-1)TC-1(Hm-b)
(5)
進一步轉化,得到
(6)
對式(6)進行化簡,得到
(7)

橢球方程的一般形式為
a1X2+b1Y2+c1Z2+2f1XY+2g1XZ+
2h1YZ+2p1X+2q1Y+2r1Z+d1=0
(8)
為方便后續計算,將式(8)歸一化為
aX2+bY2+cZ2+2fXY+2gXZ+
2hYZ+2pX+2qY+2rZ=1
(9)
將(9)式寫成向量形式,則有
Hiξ=1
(10)

至此,磁場量測值擬合的橢球面參數已經獲得。下一步求解校準矩陣C-1和偏移量b。將式(9)轉化為矩陣形式
HmTEHm+(2F)THm+G=0
(11)
建立式(11)與式(7)之間的聯系,對式(11)進行轉換可得
HmTEHm+2FTHm+G=0
?HmTEHm+FTHm+HmTF+G=0
?(HmT+FTE-1)EHm+HmTF=-G
?(HmT+FTE-1)(EHm+F)=FTE-1F-G
?(HmT+FTE-1)E(Hm+E-1F)=FTE-1F-G
?(Hm+E-1F)TE(Hm+E-1F)=FTE-1F-G
(12)
結合式(6)(7)(12),得到
b=x0=-E-1F
(13)
(14)
對式(14)進一步轉化可得
(15)
通過式(13)(15)得到校準矩陣C-1和偏移矢量b,結合式(3)可求出實際地磁場數據,完成誤差補償。
當數據中包含離群點的時候,式(10)中ξ的估計將受到嚴重的影響,利用最小二乘法對磁場數據進行橢球擬合時,無法有效消除異常點帶來的影響,其精度無法得到保證。為此,本文引入M估計算法解決上述問題。魯棒M估計通過自適應地為樣本分配不同的權值(為離群點分配接近于0的權值),消除離群點對模型參數估計結果的影響。因此ξ的M估計相當于橢球面參數的優化問題。
在數據采集過程中,會不可避免地出現一些異常值。最小二乘法即
(16)

M估計是1964年由Huber提出的,是最常用的穩健估計算法[18]。基本思想是采用迭代加權最小二乘估計回歸系數,根據回歸殘差的大小確定各點的權值,從而達到穩健的目的。M估計一般定義為
(17)

(18)
式中med為取中位數函數。
對于式(17),不同的目標函數最后的效果都差不多。本文使用Huber法,其目標函數為
(19)

(20)
式中ψ0(u)是ρ(u)的導數,公式為
(21)
權重代價函數定義為ω(u)=ψ0(u)/u,則式(20)變為
(22)

(23)
其中k為可調參數。不同的k對各采樣點的權重會有一定的影響。通常k取為1.345,此時估計算法既是穩健的,又有較高的估計效率。后文4.1小節實驗也證明了該數值下本文算法的有效性。
M估計的估計方程寫成矩陣形式為
XTWXβ=XTWY
(24)
迭代公式為
(25)
式中:W是以Wi為對角元的權矩陣,i=1,2,…,n;X是解釋變量矩陣,X=[x1,x2,…,xn]T;Y是因變量向量,Y=[y1,y2,…,yn]T。
在磁場內部旋轉磁傳感器,得到采樣數據,將采樣數據代入地磁傳感器誤差模型中,然后設計M估計補充策略下的磁傳感器誤差補償算法,步驟如下。
步驟1 初始化。通過最小二乘法得到的橢球面參數的估計值ξ(0)作為回歸系數初始值。

步驟3 權重獲取。將標準化殘差代入權函數中,得到各采樣點權重,并構成權重函數矩陣W(k)。
步驟4 橢球參數更新。由式(25)進行拓展,得到橢球面參數更新公式為

步驟6 校準矩陣和偏移矢量求取。根據式(12)~(15)得到校準矩陣C-1和偏移矢量b,完成誤差補償。


在該橢球面上選取1 200個點并疊加高斯白噪聲ε作為磁傳感器量測值,其中ε=[εx,εy,εz]T;εx,εy,εz~N(0,150 nT2)。隨機插入100個異常值ηi,ηi=[ηix,ηiy,ηiz]T;ηix,ηiy,ηiz~U(-50 000 nT,50 000 nT)。模擬出1 300個測量數據,分布見圖1。

圖1 模擬的測量數據Fig.1 Simulated measurement data
對模擬出的測量數據分別用最小二乘法、遞推最小二乘法和M估計補充策略下的誤差補償算法進行橢球擬合。3種算法擬合出的橢球面如圖2所示。可以看出,遞推最小二乘法擬合的橢球面受到異常值的影響最大,這是因為在遞推最小二乘法中,越晚采集到的數據所占權重越大,如果這些數據中有偏離值,則會嚴重影響擬合效果。

圖2 橢球擬合結果對比Fig.2 Comparison of results of ellipsoid fitting algorithm
將測量、最小二乘法校正后、遞推最小二乘法校正后、M估計法校正后的三軸磁場數據轉化為總磁場強度,結果如圖3所示。可以看出:未校正數據波動較大;3種算法校正后的數據總體波動幅度明顯降低;M估計法受到異常值的影響明顯小于最小二乘法和遞推最小二乘法。

圖3 補償前后總磁場強度對比Fig.3 Comparison of total magnetic field intensity before and after compensation
在仿真中,3種算法去除異常測量值后的數據統計特性如表1所示。可以看出,M估計法的校正效果明顯好于最小二乘法和遞推最小二乘法的,M估計法有效降低了異常測量值對誤差補償的影響,提高了系統的魯棒性。

表1 仿真中3種算法的評估指標
在南京某開闊地帶進行地磁數據采集實驗,使用美國Honeywell公司的HMR2300磁力計作為三軸磁傳感器,如圖4所示。將磁力計分別繞x、y、z軸旋轉一周采集數據。對數據分別用最小二乘法和M估計法和遞推最小二乘法進行橢球擬合,結果如圖5所示。可以看出,采集的數據中包含了離群點,明顯脫離橢球面,因此造成3種算法擬合的橢球面有明顯差異。

圖4 實驗設備Fig.4 Experimental magnetometer

圖5 擬合模型Fig.5 Fitting model

圖6 3種補償算法的總場強 Fig.6 Total field intensity of three algorithms for compensation
將校正前和3種算法校正后的三軸磁場數據轉化為總磁場強度,結果如圖6所示。可以看出,在第1 000個采樣點附近,磁場強度發生明顯波動,這可能是因為存在外界磁場干擾,因此造成采樣出現異常值。本文認為單位化殘差(Huber目標函數中的k)大于1.345的為異常值。
測量數據中去掉異常值后,各采樣點與3種算法擬合橢球面的殘差如圖7所示。可以看出,M估計法殘差波動明顯小于最小二乘法和遞推最小二乘法。

圖7 實驗中3種算法的擬合殘差Fig.7 Residual error of the fitting results of three algorithms
在實驗中,3種算法去除異常測量值后的數據統計特性如表2所示。M估計法的指標均優于傳統最小二乘法和遞推最小二乘法的,這表明M估計法有效克服了最小二乘法和遞推最小二乘法對異常值敏感的缺點,抑制了磁測噪聲,提高了系統魯棒性。

表2 實驗中3種算法的評估指標
地磁測量誤差校正是獲取高精度真實地磁數據的關鍵。但是,之前的誤差校正算法多假設數據采集過程中無異常點的出現。因此,本文提出了一種M估計補充策略下的三軸磁傳感器誤差的補償算法,能夠降低異常點對擬合的影響,提高地磁測量的精度。首先,分析了地磁測量中各種誤差源,建立了磁傳感器誤差模型;其次,利用最小二乘法擬合橢球的參數;然后,在最小二乘法的基礎上引入了M估計法對橢球參數進行迭代更新;最終,求取校準矩陣和偏移矢量,完成誤差補償過程。仿真結果表明,M估計法優于最小二乘法和遞推最小二乘法。實驗結果表明,用M估計補償獲得的數據均方差比使用遞推最小二乘法的低78.12%,比使用最小二乘法的低47.74%。M估計在評估指標方面均優于傳統的最小二乘法和遞推最小二乘法,體現出較好的魯棒性,具有很好的工程應用前景。