潘惠坤, 李勝, 錢(qián)晨, 季蔡娟, 徐振

(南京理工大學(xué)自動(dòng)化學(xué)院, 210094, 南京)

地球磁場(chǎng)是地球系統(tǒng)的基本物理場(chǎng)之一,因?yàn)槠錈o(wú)源、穩(wěn)定且與地理位置有關(guān)的優(yōu)點(diǎn),利用其進(jìn)行導(dǎo)航的地磁匹配導(dǎo)航技術(shù)受到了學(xué)者們的關(guān)注[1-4]。與其他輔助導(dǎo)航技術(shù)相比,地磁匹配導(dǎo)航技術(shù)具有無(wú)源自主、誤差不隨時(shí)間累積、全天候、隱蔽性強(qiáng)、連續(xù)導(dǎo)航、適用范圍廣的優(yōu)點(diǎn)[5]。獲取高精度的地磁數(shù)據(jù)是地磁匹配導(dǎo)航的關(guān)鍵之一,決定著導(dǎo)航的精度。

三軸磁傳感器進(jìn)行地磁場(chǎng)強(qiáng)度測(cè)量時(shí),磁場(chǎng)數(shù)據(jù)不可避免的會(huì)包含誤差和異常值。誤差可以分為磁測(cè)系統(tǒng)誤差和外界磁干擾兩類(lèi)。磁測(cè)系統(tǒng)誤差包括三軸磁傳感器自身的誤差(三軸磁傳感器在制造過(guò)程中的零偏誤差、非正交誤差、刻度因子誤差等)以及磁測(cè)系統(tǒng)的安裝誤差[6],外界磁干擾包括傳感器周?chē)妮d體磁場(chǎng)及磁日變磁場(chǎng)等干擾信息。這些誤差使得磁傳感器測(cè)得的數(shù)據(jù)無(wú)法直接用來(lái)進(jìn)行地磁匹配,需要對(duì)測(cè)量值進(jìn)行校準(zhǔn)。因此,需要研究高精度的磁傳感器誤差補(bǔ)償算法。

常用的磁傳感器誤差校準(zhǔn)方法分為輔助矢量校正方法和獨(dú)立標(biāo)量校正方法兩類(lèi)。輔助矢量校正方法通過(guò)與已知的磁矢量場(chǎng)比較從而對(duì)磁傳感器數(shù)據(jù)進(jìn)行輔助校正參考[7-8]。然而,在實(shí)際應(yīng)用中,難以獲取高精度的已知磁場(chǎng)矢量。獨(dú)立標(biāo)量校正法則可以避免該缺點(diǎn),通過(guò)三軸磁傳感器在恒定磁場(chǎng)內(nèi)旋轉(zhuǎn),以合成總場(chǎng)是固定值作為約束條件對(duì)其進(jìn)行校正。獨(dú)立標(biāo)量校正方法因其在實(shí)際環(huán)境中易操作的優(yōu)點(diǎn)吸引了大量學(xué)者的研究,典型代表為基于橢球假設(shè)的磁補(bǔ)償方法[9-10]。文獻(xiàn)[11]中提出了基于最小二乘法的誤差標(biāo)定方法并應(yīng)用于磁航向誤差補(bǔ)償中,測(cè)試結(jié)果滿(mǎn)足精度要求;文獻(xiàn)[12]利用自適應(yīng)最小二乘估計(jì)法解決了橢球體擬合問(wèn)題,取得了良好的校正效果;文獻(xiàn)[13]設(shè)計(jì)了一種基于自適應(yīng)遺傳算法的空間橢球磁強(qiáng)計(jì)校準(zhǔn)方法,能夠同時(shí)兼顧三軸磁傳感器誤差及其安裝誤差;文獻(xiàn)[14]提出了基于遞推最小二乘法的三軸磁傳感器誤差在線(xiàn)自校正方法。基于橢球假設(shè)的磁補(bǔ)償方法是最常見(jiàn)的磁傳感器標(biāo)定和補(bǔ)償方法,并廣泛應(yīng)用于實(shí)際測(cè)量中。

這些算法均假設(shè)三軸磁傳感器測(cè)量時(shí)不包含誤差較大的采樣點(diǎn)或異常值。然而,在實(shí)際測(cè)量中,該假設(shè)不成立,并會(huì)對(duì)最終精度造成較大影響。M估計(jì)法是一種有效提高系統(tǒng)魯棒性的技術(shù)[15]。通過(guò)構(gòu)造權(quán)重代價(jià)函數(shù),獲得與采樣點(diǎn)殘差相對(duì)應(yīng)的權(quán)重值,降低甚至消除異常點(diǎn)造成的影響,從而增強(qiáng)系統(tǒng)魯棒性提高精度[16]。

基于這一思路,本文提出了M估計(jì)補(bǔ)充策略下的三軸磁傳感器誤差補(bǔ)償算法。在最小二乘法的橢球擬合算法的基礎(chǔ)上,引入M估計(jì)的思想,根據(jù)其擬合的殘差數(shù)據(jù),構(gòu)造Huber權(quán)重目標(biāo)函數(shù),確定各個(gè)采樣點(diǎn)的權(quán)重,為偏離較大的采樣點(diǎn)設(shè)置低權(quán)重,為偏離較小的點(diǎn)設(shè)置高權(quán)重,對(duì)數(shù)據(jù)進(jìn)行靈活處理,從而增強(qiáng)系統(tǒng)魯棒性,降低異常點(diǎn)對(duì)擬合的影響,最終提高地磁測(cè)量的精度。

1 地磁傳感器誤差分析與建模

1.1 地磁傳感器誤差分析

地磁傳感器誤差主要由載體磁場(chǎng)誤差和自身誤差組成。

載體磁場(chǎng)誤差是由磁傳感器周?chē)母鞣N磁性材料造成的。載體磁場(chǎng)誤差是地磁場(chǎng)測(cè)量中所特有的一種誤差,可以分為硬磁誤差、軟磁誤差以及隨機(jī)磁場(chǎng)誤差共3類(lèi)主要誤差。

三軸磁傳感器在實(shí)際使用過(guò)程中,存在安裝差異,即自身誤差。傳感器誤差來(lái)源主要有3類(lèi):三軸靈敏度不一致造成的靈敏度誤差、傳感器三軸未完全正交造成的非正交誤差以及傳感器的零位偏移誤差。自身誤差屬于機(jī)械誤差,比較固定,出廠(chǎng)后不易改變。

1.2 地磁傳感器誤差模型建模

綜合考慮載體磁場(chǎng)誤差和自身誤差,參考文獻(xiàn)[17]建立誤差模型

Hm=CsiCn(CsHe+bh)+b0+ε

(1)

將式(1)簡(jiǎn)化可得

Hm=CHe+b+ε

(2)

式中:C=CsiCnCs;b=CsiCnbh+b0。顯然,C是可逆矩陣,因此可以得到地磁場(chǎng)誤差校正模型

He=C-1(Hm-b-ε)

(3)

在求出校準(zhǔn)矩陣C-1和偏移量b的基礎(chǔ)上,當(dāng)獲得磁傳感器量測(cè)數(shù)據(jù)后,通過(guò)式(3)可以計(jì)算出當(dāng)?shù)氐恼鎸?shí)地磁強(qiáng)度,式中ε可以通過(guò)合理的硬件技術(shù)手段減小,因此計(jì)算過(guò)程中可以忽略。

2 橢球擬合算法

在磁場(chǎng)強(qiáng)度恒定的區(qū)域內(nèi)對(duì)磁傳感器進(jìn)行任意旋轉(zhuǎn),測(cè)量出的磁場(chǎng)應(yīng)該是一個(gè)定值,即磁場(chǎng)矢量的軌跡應(yīng)該是一個(gè)正球體,可表達(dá)為

‖He‖2=Cconst

(4)

式中Cconst代表磁場(chǎng)強(qiáng)度的定值。結(jié)合式(2)(4)得

‖He‖2=HeTHe=

(Hm-b)T(C-1)TC-1(Hm-b)

(5)

進(jìn)一步轉(zhuǎn)化,得到

(6)

對(duì)式(6)進(jìn)行化簡(jiǎn),得到

(7)

橢球方程的一般形式為

a1X2+b1Y2+c1Z2+2f1XY+2g1XZ+

2h1YZ+2p1X+2q1Y+2r1Z+d1=0

(8)

為方便后續(xù)計(jì)算,將式(8)歸一化為

aX2+bY2+cZ2+2fXY+2gXZ+
2hYZ+2pX+2qY+2rZ=1

(9)

將(9)式寫(xiě)成向量形式,則有

Hiξ=1

(10)

至此,磁場(chǎng)量測(cè)值擬合的橢球面參數(shù)已經(jīng)獲得。下一步求解校準(zhǔn)矩陣C-1和偏移量b。將式(9)轉(zhuǎn)化為矩陣形式

HmTEHm+(2F)THm+G=0

(11)

建立式(11)與式(7)之間的聯(lián)系,對(duì)式(11)進(jìn)行轉(zhuǎn)換可得

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)

結(jié)合式(6)(7)(12),得到

b=x0=-E-1F

(13)

(14)

對(duì)式(14)進(jìn)一步轉(zhuǎn)化可得

(15)

通過(guò)式(13)(15)得到校準(zhǔn)矩陣C-1和偏移矢量b,結(jié)合式(3)可求出實(shí)際地磁場(chǎng)數(shù)據(jù),完成誤差補(bǔ)償。

3 M估計(jì)補(bǔ)充策略下的誤差補(bǔ)償算法

當(dāng)數(shù)據(jù)中包含離群點(diǎn)的時(shí)候,式(10)中ξ的估計(jì)將受到嚴(yán)重的影響,利用最小二乘法對(duì)磁場(chǎng)數(shù)據(jù)進(jìn)行橢球擬合時(shí),無(wú)法有效消除異常點(diǎn)帶來(lái)的影響,其精度無(wú)法得到保證。為此,本文引入M估計(jì)算法解決上述問(wèn)題。魯棒M估計(jì)通過(guò)自適應(yīng)地為樣本分配不同的權(quán)值(為離群點(diǎn)分配接近于0的權(quán)值),消除離群點(diǎn)對(duì)模型參數(shù)估計(jì)結(jié)果的影響。因此ξ的M估計(jì)相當(dāng)于橢球面參數(shù)的優(yōu)化問(wèn)題。

3.1 M估計(jì)算法

在數(shù)據(jù)采集過(guò)程中,會(huì)不可避免地出現(xiàn)一些異常值。最小二乘法即

(16)

M估計(jì)是1964年由Huber提出的,是最常用的穩(wěn)健估計(jì)算法[18]?；舅枷胧遣捎玫訖?quán)最小二乘估計(jì)回歸系數(shù),根據(jù)回歸殘差的大小確定各點(diǎn)的權(quán)值,從而達(dá)到穩(wěn)健的目的。M估計(jì)一般定義為

(17)

(18)

式中med為取中位數(shù)函數(shù)。

對(duì)于式(17),不同的目標(biāo)函數(shù)最后的效果都差不多。本文使用Huber法,其目標(biāo)函數(shù)為

(19)

(20)

式中ψ0(u)是ρ(u)的導(dǎo)數(shù),公式為

(21)

權(quán)重代價(jià)函數(shù)定義為ω(u)=ψ0(u)/u,則式(20)變?yōu)?/p>

(22)

(23)

其中k為可調(diào)參數(shù)。不同的k對(duì)各采樣點(diǎn)的權(quán)重會(huì)有一定的影響。通常k取為1.345,此時(shí)估計(jì)算法既是穩(wěn)健的,又有較高的估計(jì)效率。后文4.1小節(jié)實(shí)驗(yàn)也證明了該數(shù)值下本文算法的有效性。

M估計(jì)的估計(jì)方程寫(xiě)成矩陣形式為

XTWXβ=XTWY

(24)

迭代公式為

(25)

式中:W是以Wi為對(duì)角元的權(quán)矩陣,i=1,2,…,n;X是解釋變量矩陣,X=[x1,x2,…,xn]T;Y是因變量向量,Y=[y1,y2,…,yn]T。

3.2 M估計(jì)補(bǔ)充策略下的磁傳感器誤差補(bǔ)償算法

在磁場(chǎng)內(nèi)部旋轉(zhuǎn)磁傳感器,得到采樣數(shù)據(jù),將采樣數(shù)據(jù)代入地磁傳感器誤差模型中,然后設(shè)計(jì)M估計(jì)補(bǔ)充策略下的磁傳感器誤差補(bǔ)償算法,步驟如下。

步驟1 初始化。通過(guò)最小二乘法得到的橢球面參數(shù)的估計(jì)值ξ(0)作為回歸系數(shù)初始值。

步驟3 權(quán)重獲取。將標(biāo)準(zhǔn)化殘差代入權(quán)函數(shù)中,得到各采樣點(diǎn)權(quán)重,并構(gòu)成權(quán)重函數(shù)矩陣W(k)。

步驟4 橢球參數(shù)更新。由式(25)進(jìn)行拓展,得到橢球面參數(shù)更新公式為

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

4 數(shù)值分析與實(shí)驗(yàn)驗(yàn)證

4.1 仿真驗(yàn)證

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

圖1 模擬的測(cè)量數(shù)據(jù)Fig.1 Simulated measurement data

對(duì)模擬出的測(cè)量數(shù)據(jù)分別用最小二乘法、遞推最小二乘法和M估計(jì)補(bǔ)充策略下的誤差補(bǔ)償算法進(jìn)行橢球擬合。3種算法擬合出的橢球面如圖2所示?？梢钥闯?遞推最小二乘法擬合的橢球面受到異常值的影響最大,這是因?yàn)樵谶f推最小二乘法中,越晚采集到的數(shù)據(jù)所占權(quán)重越大,如果這些數(shù)據(jù)中有偏離值,則會(huì)嚴(yán)重影響擬合效果。

圖2 橢球擬合結(jié)果對(duì)比Fig.2 Comparison of results of ellipsoid fitting algorithm

將測(cè)量、最小二乘法校正后、遞推最小二乘法校正后、M估計(jì)法校正后的三軸磁場(chǎng)數(shù)據(jù)轉(zhuǎn)化為總磁場(chǎng)強(qiáng)度,結(jié)果如圖3所示。可以看出:未校正數(shù)據(jù)波動(dòng)較大;3種算法校正后的數(shù)據(jù)總體波動(dòng)幅度明顯降低;M估計(jì)法受到異常值的影響明顯小于最小二乘法和遞推最小二乘法。

圖3 補(bǔ)償前后總磁場(chǎng)強(qiáng)度對(duì)比Fig.3 Comparison of total magnetic field intensity before and after compensation

在仿真中,3種算法去除異常測(cè)量值后的數(shù)據(jù)統(tǒng)計(jì)特性如表1所示?？梢钥闯?M估計(jì)法的校正效果明顯好于最小二乘法和遞推最小二乘法的,M估計(jì)法有效降低了異常測(cè)量值對(duì)誤差補(bǔ)償?shù)挠绊?提高了系統(tǒng)的魯棒性。

表1 仿真中3種算法的評(píng)估指標(biāo)

4.2 實(shí)驗(yàn)驗(yàn)證

在南京某開(kāi)闊地帶進(jìn)行地磁數(shù)據(jù)采集實(shí)驗(yàn),使用美國(guó)Honeywell公司的HMR2300磁力計(jì)作為三軸磁傳感器,如圖4所示。將磁力計(jì)分別繞x、y、z軸旋轉(zhuǎn)一周采集數(shù)據(jù)。對(duì)數(shù)據(jù)分別用最小二乘法和M估計(jì)法和遞推最小二乘法進(jìn)行橢球擬合,結(jié)果如圖5所示?？梢钥闯?采集的數(shù)據(jù)中包含了離群點(diǎn),明顯脫離橢球面,因此造成3種算法擬合的橢球面有明顯差異。

圖4 實(shí)驗(yàn)設(shè)備Fig.4 Experimental magnetometer

圖5 擬合模型Fig.5 Fitting model

圖6 3種補(bǔ)償算法的總場(chǎng)強(qiáng) Fig.6 Total field intensity of three algorithms for compensation

將校正前和3種算法校正后的三軸磁場(chǎng)數(shù)據(jù)轉(zhuǎn)化為總磁場(chǎng)強(qiáng)度,結(jié)果如圖6所示。可以看出,在第1 000個(gè)采樣點(diǎn)附近,磁場(chǎng)強(qiáng)度發(fā)生明顯波動(dòng),這可能是因?yàn)榇嬖谕饨绱艌?chǎng)干擾,因此造成采樣出現(xiàn)異常值。本文認(rèn)為單位化殘差(Huber目標(biāo)函數(shù)中的k)大于1.345的為異常值。

測(cè)量數(shù)據(jù)中去掉異常值后,各采樣點(diǎn)與3種算法擬合橢球面的殘差如圖7所示。可以看出,M估計(jì)法殘差波動(dòng)明顯小于最小二乘法和遞推最小二乘法。

圖7 實(shí)驗(yàn)中3種算法的擬合殘差Fig.7 Residual error of the fitting results of three algorithms

在實(shí)驗(yàn)中,3種算法去除異常測(cè)量值后的數(shù)據(jù)統(tǒng)計(jì)特性如表2所示。M估計(jì)法的指標(biāo)均優(yōu)于傳統(tǒng)最小二乘法和遞推最小二乘法的,這表明M估計(jì)法有效克服了最小二乘法和遞推最小二乘法對(duì)異常值敏感的缺點(diǎn),抑制了磁測(cè)噪聲,提高了系統(tǒng)魯棒性。

表2 實(shí)驗(yàn)中3種算法的評(píng)估指標(biāo)

5 結(jié) 論

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