基于聯合估計迭代算法的三軸磁力儀標定補償方法

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

(空軍工程大學 信息與導航學院， 陜西 西安 710077)

0 引言

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

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

針對地磁傳感器的誤差補償問題，國內外學者進行了諸多研究。Renaudin等[8]在磁場域對各向異性磁阻傳感器的校正問題進行了研究，其詳細分析了傳感器自身儀表誤差、軟磁以及硬磁誤差的來源和產生機理，建模時將這些誤差都予以充分考慮；文獻[9-10]對地磁傳感器測量誤差和軟、硬磁誤差進行一體化建模，前者采用橢球擬合校準法進行補償，后者考慮到觀測向量和數據矩陣均存在誤差，因此用總體最小二乘進行了補償。對于以上矢量建模方式，仍然不能精確地表示載體內磁場的真實狀況。渦流磁場及相關的低頻磁場是由于載體的姿態變換所產生的，渦流磁場不依賴于載體的鐵磁性，而只依賴于其導電性能[11]；文獻[12] 采用多物理仿真軟件Comsol Multiphysics建立了載體渦流干擾場的仿真模型，并進行定性分析；文獻[13]提出了一體化建模的思想，并對旋轉速度等特性做了研究。

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

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

1 誤差補償模型

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

B=CH+b+ε，

(1)

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

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

(2)

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

變換得地磁場真實值的表達式為

(3)

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

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

(4)

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

Q=QARA.

(5)

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

(6)

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

各參數矩陣表達式為

當各個系統參量明確時，即能對地磁場測量值實現補償，因此設參數向量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為觀測噪聲。

則該問題的狀態方程和觀測方程可以表示為

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

(9)

(10)

fk(Xk)=Xk,

(11)

(12)

(13)

式中：設系統噪聲ηk為均值為0，協方差矩陣為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初值設置中，X0取經驗值，P0可以取較大，由于數據量充足，經過k次求解最后得到Xk即為所需參量值。

2.2 非線性最小二乘算法

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

(18)

則地磁場真實值的標量值即為

(19)

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

(20)

經過k次測量，對差值求和可得總差值表達式為

(21)

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

2.3 聯合估計迭代算法

根據上文，地磁場補償模型為

(22)

(23)

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

(24)

(25)

(26)

3 數值仿真

利用仿真數據對上述算法進行驗證，根據(2)式所述地磁場觀測值模型，設地磁場真實值為50 000 nT， 為了模擬載體姿態變換，設置球坐標系下轉動變量θ和γ，當地磁場真實值3軸分量如下產生：

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

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

(27)

(28)

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

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

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

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

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

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

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

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

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

表1 預設參數與不同算法解算值

注：XP為預設參數；XE為EKF算法解算值；XN為非線性最小二乘算法解算值;XJ為聯合估計迭代算法解算值。

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

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

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

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

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

表2 不同算法補償結果(仿真數據)

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

算法總體誤差均值/nT總體誤差標準差/nT補償前355.6542.9EKF算法26.380.0非線性最小二乘算法69.887.9聯合估計迭代算法21.061.6

4 實驗設計

4.1 實驗系統

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

圖8 地磁場測量系統Fig.8 Geomagnetic field measuring system

圖9 地磁場數據采集系統Fig.9 Geomagnetic data acquisition system

4.2 實驗步驟

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

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

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

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

5 實驗結果

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

經過3種算法補償，補償效果圖如圖10～圖12所示。

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

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

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

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

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

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

表3 不同算法補償結果(實驗數據)

6 結果分析

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

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

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

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

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

7 結論

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