潘惠坤, 李勝, 錢晨, 季蔡娟, 徐振
(南京理工大學(xué)自動(dòng)化學(xué)院, 210094, 南京)
地球磁場是地球系統(tǒng)的基本物理場之一,因?yàn)槠錈o源、穩(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ù)具有無源自主、誤差不隨時(shí)間累積、全天候、隱蔽性強(qiáng)、連續(xù)導(dǎo)航、適用范圍廣的優(yōu)點(diǎn)[5]。獲取高精度的地磁數(shù)據(jù)是地磁匹配導(dǎo)航的關(guān)鍵之一,決定著導(dǎo)航的精度。
三軸磁傳感器進(jìn)行地磁場強(qiáng)度測量時(shí),磁場數(shù)據(jù)不可避免的會(huì)包含誤差和異常值。誤差可以分為磁測系統(tǒng)誤差和外界磁干擾兩類。磁測系統(tǒng)誤差包括三軸磁傳感器自身的誤差(三軸磁傳感器在制造過程中的零偏誤差、非正交誤差、刻度因子誤差等)以及磁測系統(tǒng)的安裝誤差[6],外界磁干擾包括傳感器周圍的載體磁場及磁日變磁場等干擾信息。這些誤差使得磁傳感器測得的數(shù)據(jù)無法直接用來進(jìn)行地磁匹配,需要對測量值進(jìn)行校準(zhǔn)。因此,需要研究高精度的磁傳感器誤差補(bǔ)償算法。
常用的磁傳感器誤差校準(zhǔn)方法分為輔助矢量校正方法和獨(dú)立標(biāo)量校正方法兩類。輔助矢量校正方法通過與已知的磁矢量場比較從而對磁傳感器數(shù)據(jù)進(jìn)行輔助校正參考[7-8]。然而,在實(shí)際應(yīng)用中,難以獲取高精度的已知磁場矢量。獨(dú)立標(biāo)量校正法則可以避免該缺點(diǎn),通過三軸磁傳感器在恒定磁場內(nèi)旋轉(zhuǎn),以合成總場是固定值作為約束條件對其進(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ǔ)償中,測試結(jié)果滿足精度要求;文獻(xiàn)[12]利用自適應(yīng)最小二乘估計(jì)法解決了橢球體擬合問題,取得了良好的校正效果;文獻(xiàn)[13]設(shè)計(jì)了一種基于自適應(yīng)遺傳算法的空間橢球磁強(qiáng)計(jì)校準(zhǔn)方法,能夠同時(shí)兼顧三軸磁傳感器誤差及其安裝誤差;文獻(xiàn)[14]提出了基于遞推最小二乘法的三軸磁傳感器誤差在線自校正方法?;跈E球假設(shè)的磁補(bǔ)償方法是最常見的磁傳感器標(biāo)定和補(bǔ)償方法,并廣泛應(yīng)用于實(shí)際測量中。
這些算法均假設(shè)三軸磁傳感器測量時(shí)不包含誤差較大的采樣點(diǎn)或異常值。然而,在實(shí)際測量中,該假設(shè)不成立,并會(huì)對最終精度造成較大影響。M估計(jì)法是一種有效提高系統(tǒng)魯棒性的技術(shù)[15]。通過構(gòu)造權(quán)重代價(jià)函數(shù),獲得與采樣點(diǎn)殘差相對應(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)重,對數(shù)據(jù)進(jìn)行靈活處理,從而增強(qiáng)系統(tǒng)魯棒性,降低異常點(diǎn)對擬合的影響,最終提高地磁測量的精度。
地磁傳感器誤差主要由載體磁場誤差和自身誤差組成。
載體磁場誤差是由磁傳感器周圍的各種磁性材料造成的。載體磁場誤差是地磁場測量中所特有的一種誤差,可以分為硬磁誤差、軟磁誤差以及隨機(jī)磁場誤差共3類主要誤差。
三軸磁傳感器在實(shí)際使用過程中,存在安裝差異,即自身誤差。傳感器誤差來源主要有3類:三軸靈敏度不一致造成的靈敏度誤差、傳感器三軸未完全正交造成的非正交誤差以及傳感器的零位偏移誤差。自身誤差屬于機(jī)械誤差,比較固定,出廠后不易改變。
綜合考慮載體磁場誤差和自身誤差,參考文獻(xiàn)[17]建立誤差模型
Hm=CsiCn(CsHe+bh)+b0+ε
(1)
將式(1)簡化可得
Hm=CHe+b+ε
(2)
式中:C=CsiCnCs;b=CsiCnbh+b0。顯然,C是可逆矩陣,因此可以得到地磁場誤差校正模型
He=C-1(Hm-b-ε)
(3)
在求出校準(zhǔn)矩陣C-1和偏移量b的基礎(chǔ)上,當(dāng)獲得磁傳感器量測數(shù)據(jù)后,通過式(3)可以計(jì)算出當(dāng)?shù)氐恼鎸?shí)地磁強(qiáng)度,式中ε可以通過合理的硬件技術(shù)手段減小,因此計(jì)算過程中可以忽略。
在磁場強(qiáng)度恒定的區(qū)域內(nèi)對磁傳感器進(jìn)行任意旋轉(zhuǎn),測量出的磁場應(yīng)該是一個(gè)定值,即磁場矢量的軌跡應(yīng)該是一個(gè)正球體,可表達(dá)為
‖He‖2=Cconst
(4)
式中Cconst代表磁場強(qiáng)度的定值。結(jié)合式(2)(4)得
‖He‖2=HeTHe=
(Hm-b)T(C-1)TC-1(Hm-b)
(5)
進(jìn)一步轉(zhuǎn)化,得到
(6)
對式(6)進(jì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)式寫成向量形式,則有
Hiξ=1
(10)
至此,磁場量測值擬合的橢球面參數(shù)已經(jīng)獲得。下一步求解校準(zhǔn)矩陣C-1和偏移量b。將式(9)轉(zhuǎn)化為矩陣形式
HmTEHm+(2F)THm+G=0
(11)
建立式(11)與式(7)之間的聯(lián)系,對式(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)
對式(14)進(jìn)一步轉(zhuǎn)化可得
(15)
通過式(13)(15)得到校準(zhǔn)矩陣C-1和偏移矢量b,結(jié)合式(3)可求出實(shí)際地磁場數(shù)據(jù),完成誤差補(bǔ)償。
當(dāng)數(shù)據(jù)中包含離群點(diǎn)的時(shí)候,式(10)中ξ的估計(jì)將受到嚴(yán)重的影響,利用最小二乘法對磁場數(shù)據(jù)進(jìn)行橢球擬合時(shí),無法有效消除異常點(diǎn)帶來的影響,其精度無法得到保證。為此,本文引入M估計(jì)算法解決上述問題。魯棒M估計(jì)通過自適應(yīng)地為樣本分配不同的權(quán)值(為離群點(diǎn)分配接近于0的權(quán)值),消除離群點(diǎn)對模型參數(shù)估計(jì)結(jié)果的影響。因此ξ的M估計(jì)相當(dāng)于橢球面參數(shù)的優(yōu)化問題。
在數(shù)據(jù)采集過程中,會(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ù)。
對于式(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對各采樣點(diǎn)的權(quán)重會(huì)有一定的影響。通常k取為1.345,此時(shí)估計(jì)算法既是穩(wěn)健的,又有較高的估計(jì)效率。后文4.1小節(jié)實(shí)驗(yàn)也證明了該數(shù)值下本文算法的有效性。
M估計(jì)的估計(jì)方程寫成矩陣形式為
XTWXβ=XTWY
(24)
迭代公式為
(25)
式中:W是以Wi為對角元的權(quán)矩陣,i=1,2,…,n;X是解釋變量矩陣,X=[x1,x2,…,xn]T;Y是因變量向量,Y=[y1,y2,…,yn]T。
在磁場內(nèi)部旋轉(zhuǎn)磁傳感器,得到采樣數(shù)據(jù),將采樣數(shù)據(jù)代入地磁傳感器誤差模型中,然后設(shè)計(jì)M估計(jì)補(bǔ)充策略下的磁傳感器誤差補(bǔ)償算法,步驟如下。
步驟1 初始化。通過最小二乘法得到的橢球面參數(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ǔ)償。
在該橢球面上選取1 200個(gè)點(diǎn)并疊加高斯白噪聲ε作為磁傳感器量測值,其中ε=[ε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è)測量數(shù)據(jù),分布見圖1。
圖1 模擬的測量數(shù)據(jù)Fig.1 Simulated measurement data
對模擬出的測量數(shù)據(jù)分別用最小二乘法、遞推最小二乘法和M估計(jì)補(bǔ)充策略下的誤差補(bǔ)償算法進(jìn)行橢球擬合。3種算法擬合出的橢球面如圖2所示??梢钥闯?遞推最小二乘法擬合的橢球面受到異常值的影響最大,這是因?yàn)樵谶f推最小二乘法中,越晚采集到的數(shù)據(jù)所占權(quán)重越大,如果這些數(shù)據(jù)中有偏離值,則會(huì)嚴(yán)重影響擬合效果。
圖2 橢球擬合結(jié)果對比Fig.2 Comparison of results of ellipsoid fitting algorithm
將測量、最小二乘法校正后、遞推最小二乘法校正后、M估計(jì)法校正后的三軸磁場數(shù)據(jù)轉(zhuǎn)化為總磁場強(qiáng)度,結(jié)果如圖3所示。可以看出:未校正數(shù)據(jù)波動(dòng)較大;3種算法校正后的數(shù)據(jù)總體波動(dòng)幅度明顯降低;M估計(jì)法受到異常值的影響明顯小于最小二乘法和遞推最小二乘法。
圖3 補(bǔ)償前后總磁場強(qiáng)度對比Fig.3 Comparison of total magnetic field intensity before and after compensation
在仿真中,3種算法去除異常測量值后的數(shù)據(jù)統(tǒng)計(jì)特性如表1所示??梢钥闯?M估計(jì)法的校正效果明顯好于最小二乘法和遞推最小二乘法的,M估計(jì)法有效降低了異常測量值對誤差補(bǔ)償?shù)挠绊?提高了系統(tǒng)的魯棒性。
表1 仿真中3種算法的評估指標(biāo)
在南京某開闊地帶進(jìn)行地磁數(shù)據(jù)采集實(shí)驗(yàn),使用美國Honeywell公司的HMR2300磁力計(jì)作為三軸磁傳感器,如圖4所示。將磁力計(jì)分別繞x、y、z軸旋轉(zhuǎn)一周采集數(shù)據(jù)。對數(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ǔ)償算法的總場強(qiáng) Fig.6 Total field intensity of three algorithms for compensation
將校正前和3種算法校正后的三軸磁場數(shù)據(jù)轉(zhuǎn)化為總磁場強(qiáng)度,結(jié)果如圖6所示。可以看出,在第1 000個(gè)采樣點(diǎn)附近,磁場強(qiáng)度發(fā)生明顯波動(dòng),這可能是因?yàn)榇嬖谕饨绱艌龈蓴_,因此造成采樣出現(xiàn)異常值。本文認(rèn)為單位化殘差(Huber目標(biāo)函數(shù)中的k)大于1.345的為異常值。
測量數(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種算法去除異常測量值后的數(shù)據(jù)統(tǒng)計(jì)特性如表2所示。M估計(jì)法的指標(biāo)均優(yōu)于傳統(tǒng)最小二乘法和遞推最小二乘法的,這表明M估計(jì)法有效克服了最小二乘法和遞推最小二乘法對異常值敏感的缺點(diǎn),抑制了磁測噪聲,提高了系統(tǒng)魯棒性。
表2 實(shí)驗(yàn)中3種算法的評估指標(biāo)
地磁測量誤差校正是獲取高精度真實(shí)地磁數(shù)據(jù)的關(guān)鍵。但是,之前的誤差校正算法多假設(shè)數(shù)據(jù)采集過程中無異常點(diǎn)的出現(xiàn)。因此,本文提出了一種M估計(jì)補(bǔ)充策略下的三軸磁傳感器誤差的補(bǔ)償算法,能夠降低異常點(diǎn)對擬合的影響,提高地磁測量的精度。首先,分析了地磁測量中各種誤差源,建立了磁傳感器誤差模型;其次,利用最小二乘法擬合橢球的參數(shù);然后,在最小二乘法的基礎(chǔ)上引入了M估計(jì)法對橢球參數(shù)進(jìn)行迭代更新;最終,求取校準(zhǔn)矩陣和偏移矢量,完成誤差補(bǔ)償過程。仿真結(jié)果表明,M估計(jì)法優(yōu)于最小二乘法和遞推最小二乘法。實(shí)驗(yàn)結(jié)果表明,用M估計(jì)補(bǔ)償獲得的數(shù)據(jù)均方差比使用遞推最小二乘法的低78.12%,比使用最小二乘法的低47.74%。M估計(jì)在評估指標(biāo)方面均優(yōu)于傳統(tǒng)的最小二乘法和遞推最小二乘法,體現(xiàn)出較好的魯棒性,具有很好的工程應(yīng)用前景。