孫歡， 楊賓峰， 管樺， 王潤

(空軍工程大學(xué) 信息與導(dǎo)航學(xué)院， 陜西 西安 710077)

0 引言

近年來，地磁導(dǎo)航發(fā)展迅猛，其隱蔽性強、自主性高、抗干擾能力強的特點，使得地磁導(dǎo)航在水下、電磁環(huán)境惡劣的空域，有著明顯優(yōu)勢和廣闊的發(fā)展前景[1-3]。

地磁傳感器的觀測精度將直接影響導(dǎo)航精度，但在實際測量過程中，地磁信號和各類誤差信號相互耦合，組成了地磁傳感器輸出，嚴(yán)重影響了地磁場測量精度，抑制了磁力儀潛在性能的發(fā)揮。因此對于地磁傳感器進行標(biāo)定和補償，對提高地磁匹配導(dǎo)航的精度有著重要意義[4-7]。

針對地磁傳感器的誤差補償問題，國內(nèi)外學(xué)者進行了諸多研究。Renaudin等[8]在磁場域?qū)Ω飨虍愋源抛鑲鞲衅鞯男Ｕ龁栴}進行了研究，其詳細(xì)分析了傳感器自身儀表誤差、軟磁以及硬磁誤差的來源和產(chǎn)生機理，建模時將這些誤差都予以充分考慮；文獻(xiàn)[9-10]對地磁傳感器測量誤差和軟、硬磁誤差進行一體化建模，前者采用橢球擬合校準(zhǔn)法進行補償，后者考慮到觀測向量和數(shù)據(jù)矩陣均存在誤差，因此用總體最小二乘進行了補償。對于以上矢量建模方式，仍然不能精確地表示載體內(nèi)磁場的真實狀況。渦流磁場及相關(guān)的低頻磁場是由于載體的姿態(tài)變換所產(chǎn)生的，渦流磁場不依賴于載體的鐵磁性，而只依賴于其導(dǎo)電性能[11]；文獻(xiàn)[12] 采用多物理仿真軟件Comsol Multiphysics建立了載體渦流干擾場的仿真模型，并進行定性分析；文獻(xiàn)[13]提出了一體化建模的思想，并對旋轉(zhuǎn)速度等特性做了研究。

對于磁場補償算法，目前主流的算法有航磁補償法、兩步估計法、橢圓擬合法及無跡卡爾曼濾波(UKF)算法。航磁補償法[11,14]是基于Tolles-Lawson方程，將地磁場測量值的三分量均投影在地磁場方向，但忽略了地磁場真實值和測量值方向之間的誤差，因此只適用于干擾場較小的情況；在兩步估計算法[5,10]中，由于中間變量的存在，各個參量的相關(guān)性可能出現(xiàn)系數(shù)矩陣奇異，從而無法正確解算；橢圓擬合法[15-17]要求測量值要擬合出一個橢圓，從而建模方式比較簡單，難以適應(yīng)復(fù)雜的測量環(huán)境；UKF算法[18-19]能夠?qū)崿F(xiàn)高精度的實時補償，但對初值敏感度高，選取不當(dāng)可能造成濾波發(fā)散，得不到正確結(jié)果。

通過分析，以往研究對載體內(nèi)渦流磁場及雙噪聲研究不足，缺少一體化模型和高效的補償算法。本文針對此問題進行一體化建模，提出了基于聯(lián)合估計的迭代算法，并與擴展卡爾曼濾波(EKF)算法和非線性最小二乘算法進行性能比較。該方法具有良好的仿真和實驗性能，可用于提高三軸磁力儀的測量精度。

1 誤差補償模型

根據(jù)以往研究[20]，地磁傳感器誤差模型可以表示為

B=CH+b+ε，

(1)

式中：B為地磁傳感器測量值；H為地磁場真實值；C=CMCNOCSFCSI，CM為地磁傳感器安裝于載體過程中傳感器和載體3個軸完全不重合所產(chǎn)生的誤差矩陣，CNO為傳感器內(nèi)部三軸不完全正交所引起的誤差矩陣，CSF為制作工藝所產(chǎn)生的刻度因子誤差，CSI為載體內(nèi)軟磁誤差系數(shù)矩陣；b=CMCNOCSFCSIBh+w+bs，Bh為載體內(nèi)硬磁材料所產(chǎn)生的干擾，w為傳感器內(nèi)部零刻度漂移誤差，bs為傳感器的剩磁誤差；ε為測量過程中的觀測噪聲，以往研究表明，ε呈高斯分布。

對于高速、頻繁作姿態(tài)變換的載體，其內(nèi)部會產(chǎn)生渦流，進而產(chǎn)生渦流磁場，影響傳感器的精確測量。研究表明，渦流磁場強度與在體內(nèi)部渦流呈線性比例關(guān)系，而渦流強度又與單位時間磁通量的變化呈線型比例關(guān)系，將比例系數(shù)矩陣整合，考慮渦流磁場在內(nèi)，可得模型[21]：

(2)

式中：k代表第k次求解，k=1,2,…，n，n為數(shù)據(jù)總量；P為整合后的渦流磁場與單位時間內(nèi)磁通量變化的比例系數(shù)矩陣。

變換得地磁場真實值的表達(dá)式為

(3)

式中：Q=C-1，在(3)式中，參數(shù)空間維數(shù)Θ為

dimΘ=elements(Q)+elements(b)+
elements(P)=21.

(4)

由上文分析可知，Q在本文所設(shè)環(huán)境下為非奇異矩陣，即儀表誤差和載體干擾磁場誤差必然存在(即使很小)，因此對參數(shù)矩陣Q進行QR分解如下：

Q=QARA.

(5)

經(jīng)過QR[22]分解，將Q參數(shù)矩陣分解為正交矩陣QA和上三角矩陣RA，將(5)式代入(3)式并取模，得

(6)

由于QA、RA是唯一的，因此在保證參數(shù)矩陣的唯一性基礎(chǔ)上，將參數(shù)空間由21維降到了18維，得模型如下：

各參數(shù)矩陣表達(dá)式為

當(dāng)各個系統(tǒng)參量明確時，即能對地磁場測量值實現(xiàn)補償，因此設(shè)參數(shù)向量Xk為待求量：

Xk=[vecs(RA)T,bT,vecs(P)T]T=
[q11,q12,q13,q22,q23,q33,b1,b2,b3,p11,
p12,p13,p21,p22,p23,p31,p32,p33]T.

(7)

2 補償算法

2.1 EKF算法

(8)

式中：RAk為RA在第k次的值；νk為觀測噪聲。

則該問題的狀態(tài)方程和觀測方程可以表示為

Xk=fk(Xk-1)+ηk,

(9)

(10)

fk(Xk)=Xk,

(11)

(12)

(13)

式中：設(shè)系統(tǒng)噪聲ηk為均值為0，協(xié)方差矩陣為Q的高斯噪聲。

由(13)式，觀測噪聲的均值μk和方差R可表示為

μk=E[vk]=-tr(Σ),

(14)

(15)

Σ={εεT}.

(16)

k|k-1=fk(k-1),


k=k|k-1+Gk(Zk-hk(k|k-1)),
Pk=(I-GkJk)Pk|k-1,

(17)

在EKF初值設(shè)置中，X0取經(jīng)驗值，P0可以取較大，由于數(shù)據(jù)量充足，經(jīng)過k次求解最后得到Xk即為所需參量值。

2.2 非線性最小二乘算法

在非線性最小二乘算法處理時，由于可以經(jīng)過多次測量取平均值，因此可以將系統(tǒng)噪聲ηk與觀測噪聲εk忽略掉，根據(jù)(6)式，地磁場補償模型變換為

(18)

則地磁場真實值的標(biāo)量值即為

(19)

測量值與真實值的差值可以表示為

(20)

經(jīng)過k次測量，對差值求和可得總差值表達(dá)式為

(21)

根據(jù)非線性最小二乘算法原理，當(dāng)E取得最小值時，X取得最優(yōu)解，即X=argmin(E)，進而得到補償模型中的各個系統(tǒng)參量。

2.3 聯(lián)合估計迭代算法

根據(jù)上文，地磁場補償模型為

(22)

(23)

這同樣是關(guān)于X和w的二次優(yōu)化問題，因此存在閉合解，如(25)式所示。

(24)

(25)

(26)

3 數(shù)值仿真

利用仿真數(shù)據(jù)對上述算法進行驗證，根據(jù)(2)式所述地磁場觀測值模型，設(shè)地磁場真實值為50 000 nT， 為了模擬載體姿態(tài)變換，設(shè)置球坐標(biāo)系下轉(zhuǎn)動變量θ和γ，當(dāng)?shù)卮艌稣鎸嵵?軸分量如下產(chǎn)生：

地磁場3個軸的分量值將會隨θ和γ的變化而變化，從而模擬載體姿態(tài)變化。

設(shè)觀測過程中系統(tǒng)噪聲、觀測噪聲服從均值為0 nT，標(biāo)準(zhǔn)差分別為10 nT和50 nT的高斯分布，且相互統(tǒng)計獨立，為了更好地說明補償效果，本文采用誤差均值E(n)和標(biāo)準(zhǔn)差δ進行分析：

(27)

(28)

式中：Bi為每次的測量總值；Bm為測量值均值；B0為該地磁場的真實值。

在仿真中，設(shè)置不同轉(zhuǎn)動頻率，從而改變載體內(nèi)渦流磁場的強弱，得地磁場測量值如圖1所示。

圖1 不同轉(zhuǎn)動頻率下地磁場測量值Fig.1 Geomagnetic field measurements at different rotational frequencies

由圖1可知，當(dāng)載體瞬時轉(zhuǎn)動頻率較大時，其內(nèi)部的渦流磁場將嚴(yán)重影響地磁傳感器的精確測量。提取仿真數(shù)據(jù)中50 Hz轉(zhuǎn)動下的地磁場觀測數(shù)據(jù)，用本文所提聯(lián)合估計迭代算法分別對有無渦流磁場項的模型進行補償，驗證含渦流磁場一體化模型的有效性，結(jié)果如圖2所示。

圖2 不同模型補償結(jié)果Fig.2 Compensated results of different algorithms

圖2結(jié)果表明，含渦流磁場項的一體化模型有效濾除了由渦流磁場引起的波動干擾，使補償結(jié)果更加精確穩(wěn)定。對以上仿真數(shù)據(jù)，分別采用EKF算法、非線性最小二乘算法和聯(lián)合估計迭代算法進行處理，得到各參量解算值，而后用這些值進行誤差補償，比較各個算法實用性。經(jīng)過解算，各系統(tǒng)參量的預(yù)設(shè)值和解算值如表1所示。

由表1可知，3種算法均能對系統(tǒng)參量實現(xiàn)補償，其中，非線性最小二乘算法由于未考慮噪聲影響，因此補償效果次于其他兩種算法。3種算法補償情況如圖3～圖7、表2所示。

如圖3所示，EKF算法補償時，任取參數(shù)中的兩個值分析其收斂情況。由圖3可看出，在經(jīng)過大約400組數(shù)據(jù)的運算，各參數(shù)收斂并趨于穩(wěn)定，因此在補償初始階段，整體補償效果不理想。EKF算法整體補償效果如圖4所示。非線性最小二乘算法由于未進行降噪處理，補償結(jié)果有較大誤差(見圖5)。聯(lián)合估計迭代算法無論是地磁場分量還是總量均能夠?qū)崿F(xiàn)較好的補償(見圖6、圖7)，各算法補償結(jié)果數(shù)據(jù)如表2所示。

仿真結(jié)果表明：非線性最小二乘算法補償結(jié)果誤差較大，另外兩種算法均可實現(xiàn)精度較高的補償，但EKF算法補償精度對數(shù)據(jù)量的依賴較大(由圖1可知，在本文仿真中，當(dāng)數(shù)據(jù)量少于400組時，參數(shù)不能收斂，無法做到高精度的補償)；聯(lián)合估計迭代算法相比于其他兩種算法，補償精度最高，抑制比可以達(dá)到94.08%.

表1 預(yù)設(shè)參數(shù)與不同算法解算值

注：XP為預(yù)設(shè)參數(shù)；XE為EKF算法解算值；XN為非線性最小二乘算法解算值;XJ為聯(lián)合估計迭代算法解算值。

圖3 EKF算法部分參數(shù)收斂情況Fig.3 Partial parameter convergence of EKF algorithm

圖4 EKF算法補償效果(仿真數(shù)據(jù))Fig.4 Compensation effect of EKF algorithm (simulation data)

圖5 非線性最小二乘算法補償效果(仿真數(shù)據(jù))Fig.5 Compensation effect of nonlinear least square algorithm (simulation data)

圖6 聯(lián)合估計迭代算法整體補償效果(仿真數(shù)據(jù))Fig.6 Compensation effect of joint estimation iterative algorithm (simulation data)

圖7 聯(lián)合估計迭代算法各分量補償效果(仿真數(shù)據(jù))Fig.7 Compensation effect of joint estimation iterativealgorithm (simulation data)

表2 不同算法補償結(jié)果(仿真數(shù)據(jù))

Tab.2 Compensated results of different algorithms(simulation data)

算法總體誤差均值/nT總體誤差標(biāo)準(zhǔn)差/nT補償前355.6542.9EKF算法26.380.0非線性最小二乘算法69.887.9聯(lián)合估計迭代算法21.061.6

4 實驗設(shè)計

4.1 實驗系統(tǒng)

在實驗設(shè)計中，采用英國Bartington公司的Mag-690-FL100(測量精度0.1 nT)三軸磁力儀進行補償，將其置于1∶48飛行器縮比模型上方模擬真實情況下處于飛行載體內(nèi)部的情況。為模擬飛行載體姿態(tài)變化，選用3FHT30C無磁轉(zhuǎn)臺提供姿態(tài)信息，該轉(zhuǎn)臺剩磁為7 nT，測量系統(tǒng)采用12 V直流電源進行供電以隔絕市電電壓不穩(wěn)產(chǎn)生的不良影響。數(shù)據(jù)采集軟件用National Instruments，實驗系統(tǒng)如圖8、圖9所示，后期用數(shù)值分析Matlab軟件進行數(shù)據(jù)處理。采用DM050質(zhì)子磁力儀提供地磁場真實值。

圖8 地磁場測量系統(tǒng)Fig.8 Geomagnetic field measuring system

圖9 地磁場數(shù)據(jù)采集系統(tǒng)Fig.9 Geomagnetic data acquisition system

4.2 實驗步驟

1)實驗場所選在遠(yuǎn)離強磁干擾的空曠地點，利用水平儀對無磁轉(zhuǎn)臺進行調(diào)整，將3軸磁力儀固定好之后，在下方放置飛行器縮比模型，用以仿真飛行載體內(nèi)部產(chǎn)生的軟硬磁誤差及其他干擾。

2)連接好設(shè)備，接通直流電源，將數(shù)據(jù)采集軟件采集頻率設(shè)為50 Hz，進行連續(xù)數(shù)據(jù)采集，沿轉(zhuǎn)軸快速旋轉(zhuǎn)無磁轉(zhuǎn)臺進行全方位姿態(tài)信息測量。模擬飛行載體的姿態(tài)轉(zhuǎn)動，快速轉(zhuǎn)動時，瞬時速度較大，因此所采集數(shù)據(jù)中含有渦流磁場。采集所得數(shù)據(jù)用Matlab進行后期處理。

3)轉(zhuǎn)動無磁轉(zhuǎn)臺，將地磁傳感器3軸分別對準(zhǔn)東、北、天，所得數(shù)據(jù)即為導(dǎo)航坐標(biāo)系下地磁場值的各軸分量；無磁轉(zhuǎn)臺轉(zhuǎn)動過程中，每一個姿態(tài)對應(yīng)一個三維角度值，分別為內(nèi)框、中框和外框的轉(zhuǎn)動角度(顯示于數(shù)顯箱，記錄于采集數(shù)據(jù)中)，即地磁傳感器3軸偏離初始值的角度，分別記為α、β、γ，由此可以解算得到每一個姿態(tài)下地磁場3軸分量的真實值[22]：

式中：Bxk,Byk,Bzk為第k個姿態(tài)時地磁場各軸分量真實值；Bx0,By0,Bz0為初始狀態(tài)地磁場各軸分量真實值。

5 實驗結(jié)果

在實驗中，首先采用質(zhì)子磁力儀測得實驗場所的地磁場值為46 157 nT，對于實驗數(shù)據(jù)依次采用EKF算法、非線性最小二乘算法和聯(lián)合估計算法進行補償，其中，EKF算法初值設(shè)置為

經(jīng)過3種算法補償，補償效果圖如圖10～圖12所示。

圖10 EKF算法補償效果(實驗數(shù)據(jù))Fig.10 Compensation effect of EKF algorithm (experimental data)

圖11 非線性最小二乘算法補償效果(實驗數(shù)據(jù))Fig.11 Compensation effect of nonlinear least square algorithm (experimental data)

圖12 聯(lián)合估計迭代算法補償效果(實驗數(shù)據(jù))Fig.12 Compensation effect with joint estimation iterative algorithm (experimental data)

由圖10～圖12可以看出：實驗結(jié)果與仿真結(jié)果基本一致；實驗采用的1 000組數(shù)據(jù)，因為初值設(shè)為0向量，EKF算法經(jīng)過約300組數(shù)據(jù)后趨于收斂，在前期補償效果次于其他兩種算法；非線性最小二乘算法補償結(jié)果依然有較大噪聲干擾；聯(lián)合估計迭代算法法收斂速度快，且補償效果良好。

各算法補償結(jié)果均值和方差如表3所示。

由表3可知，聯(lián)合估計迭代算法相較于其他兩種算法優(yōu)勢比較明顯，可以將補償誤差均值由142.4 nT縮小至12.5 nT，標(biāo)準(zhǔn)差由170.0 nT縮小至37.2 nT，抑制比達(dá)到84.51%. 考慮到無磁轉(zhuǎn)臺剩磁誤差達(dá)到7.0 nT，且實驗環(huán)境周邊其他干擾，因此本次實驗結(jié)果對本文所提理論進行了驗證，具有實際意義。

表3 不同算法補償結(jié)果(實驗數(shù)據(jù))

6 結(jié)果分析

1)EKF算法原理簡單，能夠在較短的時間內(nèi)實現(xiàn)較高精度的補償。在模擬仿真中，EKF算法經(jīng)過約500組數(shù)據(jù)達(dá)到收斂，所解算參數(shù)可用于補償，而其他兩種算法皆用2 000組數(shù)據(jù)進行參數(shù)解算，體現(xiàn)出EKF算法在時效性方面的優(yōu)越性，但卡爾曼濾波算法對初值的依賴度很高。在本文實驗中，由于對系統(tǒng)參量的先驗知識不足，因此將初值取0向量，協(xié)方差矩陣取得比較大，由于數(shù)據(jù)量大，以此后期也能達(dá)到較好的補償效果。對于1階EKF算法補償精度不夠的問題，因為在本文所研究模型為弱非線性，因此其補償精度在可接受范圍之內(nèi)。

2)在實驗過程中，快速轉(zhuǎn)動無磁轉(zhuǎn)臺使得飛行器縮比模型內(nèi)部產(chǎn)生渦流磁場，對3軸磁力儀的測量產(chǎn)生影響，體現(xiàn)在測量值有較大的波動，含渦流磁場的一體化模型在建模過程中考慮到渦流磁場信息，因此相比于不含渦流磁場的模型，能夠較好地濾除波動，使補償結(jié)果更加穩(wěn)定。

3)對于觀測系統(tǒng)系統(tǒng)不穩(wěn)定的問題，在數(shù)學(xué)上體現(xiàn)為各系統(tǒng)參量的變化，在該實驗中，選取新的實驗場所、調(diào)平無磁轉(zhuǎn)臺及進行全姿態(tài)測量過程中都會使得系統(tǒng)參量發(fā)生微小改變，勢必將產(chǎn)生系統(tǒng)誤差。對于該誤差，本文將其簡化為零均值高斯噪聲，取得了較好的補償效果，為了更進一步提高補償精度，應(yīng)該對該誤差進行更深入的研究。

4)在采用非線性最小二乘算法處理時，對于噪聲采取多組數(shù)據(jù)求均值的方法進行處理，但求均值的本質(zhì)依然是含噪的。仿真和實驗結(jié)果表明，該算法機理對于本問題的局限性導(dǎo)致補償結(jié)果相比于其他兩種算法誤差偏大。

5)質(zhì)子磁力儀測量誤差很小，因此將所測地磁場值作為實驗場所真實值。本文提出含渦流磁場和雙噪聲的一體化模型，含有信息種類更全，模型更完整，補償中更有優(yōu)勢；在聯(lián)合估計迭代算法中將觀測系統(tǒng)噪聲和測量噪聲進行聯(lián)合估計處理，補償效果表明該方法比非線性最小二乘算法有著更高的補償精度，比EKF算法補償結(jié)果收斂更快，證明了該補償模型和補償方法的開發(fā)潛力。

7 結(jié)論

地磁場測量值誤差補償向來是地磁導(dǎo)航中一個關(guān)鍵問題。本文針對載體中所含渦流磁場和雙噪聲的問題，提出了雙噪聲聯(lián)合估計迭代算法對Mag-690-FL100地磁傳感器進行標(biāo)定，并與EKF算法、非線性最小二乘算法進行比較。仿真和實驗結(jié)果表明，聯(lián)合估計迭代算法補償效果更好，該方法原理簡單，對數(shù)據(jù)量依賴小，實用性強，有利于提高地磁導(dǎo)航系統(tǒng)的精度。