基于EKF算法的三軸磁力儀一體化補償方法*

孫 歡，楊賓峰，吳 昊，李 馳，王 潤

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

0 引言

地磁場作為一個天然的物理場，具有全地域、全天候、全天時等特點，以地磁場為基礎的地磁導航，相比于傳統無線電導航和衛星導航，具有隱蔽性好，自主性高，抗干擾能力強的特點；相比于地形匹配導航和天文導航，地磁導航受天氣影響小，地磁場隨時間變化緩慢，因此導航定位更穩定[1-2]。但是在實際測量中，在諸多因素的影響下，磁場測量結果將會受到多方面的干擾。首先，三軸磁力儀制造及在安裝載體過程中存在工差；其次，載體內部設備結構將形成硬磁誤差和軟磁誤差；再者，對于頻繁做姿態變換的載體，比如戰斗機，機體姿態轉動會產生較強的渦流，進而產生渦流磁場，也會極大干擾磁力儀的測量[3-9]。

針對三軸磁力儀的標定和誤差補償問題，國內外學者進行了諸多研究，Geber-Egziabher提出two-step算法補償三軸磁力儀刻度因子誤差、零偏誤差和硬磁誤差[10]；Vasconcelos對三軸磁力儀誤差和安裝誤差進行一體化建模，并提出用橢球面旋轉、尺度變換、平移變換的方法進行估計，并設計了基于最大似然算法的自校正算法[11]；文獻[12-13]對三軸磁力儀的儀表誤差、安裝誤差和軟硬磁誤差分別進行建模，后者采用ANSYS軟件仿真建模，并采用最小二乘法進行補償；文獻[14]對三軸磁力儀本身誤差和軟硬磁誤差進行一體化建模，并考慮到觀測向量和數據矩陣均存在誤差，因此用總體最小二乘進行了補償；針對載體姿態變換的問題，文獻[15] 采用COMSOL Multiphysics仿真軟件建立了載體渦流干擾場的仿真模型，并進行定性分析；文獻[16-18]建立等效磁強計誤差模型，并采用兩步估計算法和圓約束非對準誤差估計算法實現了補償。

通過分析，以往的研究存在以下問題：1)模型中沒有考慮渦流磁場的影響，對于低速、姿態穩定的載體尚可補償，但對于高速、頻繁進行姿態變換的載體，其渦流磁場的影響將會降低補償精度；2)在以往研究所普遍采用的補償方法中，最小二乘算法解算依賴大量先采數據，對數據量敏感度高，且不能做到實時補償；橢圓約束法和橢球約束法只適用于以軟磁硬磁為核心的簡單建模方式，沒有考慮儀表、安裝、渦流所帶來的誤差。文中對三軸磁力儀的儀表誤差和安裝誤差進行了分析，結合硬磁誤差和軟磁誤差進行統一處理，考慮到飛行器載體進行頻繁姿態變換的現狀，引進渦流磁場進行一體化建模，并采用EKF(擴展卡爾曼濾波)算法進行補償，最后通過模擬仿真和實驗來檢驗建模方法和補償方法的實用性和有效性。

1 地磁場觀測值誤差分析

在磁力儀標定過程中，三軸磁力儀的三軸x′、y′、z′和真實三軸x、y、z的不完全重合，將會造成非正交誤差，從而影響磁力儀的測量精度。令z軸和z′軸重合，yoz和y′oz′共面， 記oy軸和oy′軸之間的夾角為β，ox軸和面x′oy′之間的夾角為γ，ox軸和面x′oz′之間的夾角為α，其對應關系如圖1所示。

圖1 三軸非正交誤差

則非正交誤差矩陣CNO可以表示為：

(1)

在磁力儀安裝過程中，原則上傳感器坐標系應該與載體坐標系嚴格一致，但是在安裝過程中，將不可避免的出現安裝誤差，從而造成安裝誤差，設磁力儀三軸與載體坐標系三軸的誤差角為εx、εy、εz，則安裝誤差CM為：

(2)

由于制作工藝精度的問題，三軸磁力儀三軸存在細微差異，由此產生刻度因子誤差，設三軸的刻度因子誤差分別為sx、sy、sz，則刻度因子誤差矩陣CSF為：

CSF=diag(sx,sy,sz)

(3)

制作工藝精度和安裝誤差同樣也會造成零刻度漂移誤差w和剩磁誤差bs，矩陣表達式分別為：

w=[δ1,δ2,δ3]T，bs=[bx1,by1,bz1]T

(4)

硬磁誤差是載體內部硬磁材料產生的磁場所造成的誤差，常見的硬磁材料有永磁體和高碳鋼，硬磁材料的磁場隨時間變化十分緩慢，在短時間內可以認為是時不變的，因此在3個軸上的投影為常值，記硬磁誤差Bh為：

Bh=[bh1,bh2,bh3]T

(5)

軟磁誤差是軟磁性材料被外部磁場激勵所產生的誤差，外部磁場包括地磁場和載體內部所產生的磁場，在載體內部磁場中，由于集成電路和磁場擾動所產生的干擾磁場遠遠小于硬磁磁場和地磁場，因此只考慮硬磁磁場Bh和地磁場H的激勵效應。軟磁材料被激勵產生的激勵磁場強度和外部激勵強度呈線性比例關系，比例系數CSI取決于軟磁材料的磁化特性，當外部激勵磁場方向發生變化時，所產生的激勵磁場方向也會隨之變化，因此軟磁誤差BSI可以表示為：

BSI=CSI(H+Bh)

(6)

綜上所述，地磁場的觀測值B可以表示為：

B=CMCNOCSFCSI(H+Bh)+w+bs+ε

(7)

式中：ε為零均值高斯地磁觀測噪聲，其方差為∑。將其化簡，則可以表示為：

B=CH+b+ε

(8)

式中：

C=CMCNOCSFCSI

(9)

b=CMCNOCSFCSIBh+w+bs

(10)

載體做姿態變換的過程中，由于自身做切割地球磁感線運動，根據法拉第電磁感應定律，將會在載體內部產生渦流，進而產生渦流磁場，對磁力儀的測量產生影響，渦流磁場的強度與感應渦流的強度成正比，而渦流的強度又與載體磁場變化率成正比，設該比例系數矩陣為P， 將其考慮在內，可得到地磁場測量值的數學模型為：

B=CH+b+P(ΔB/Δt)+ε

(11)

2 補償算法研究

2.1 EKF算法

在利用卡爾曼濾波進行誤差補償的時候，如果出現了二次及二次以上的高次項，整個過程呈非線性，就不能用線性方法進行處理。對于非線性系統，理論上很難推導嚴格的遞推濾波公式，因此常用近似的方法進行處理。非線性系統的線性化濾波就是利用近似的方法研究非線性問題的途徑之一。EKF(擴展卡爾曼濾波)就是一種典型的方法。采用EKF進行處理時，先將非線性的傳遞函數f(·)和觀測函數h(·)進行泰勒展開，并略去二階及以上的高次項，從而得到非線性系統的線性化模型，再用標準卡爾曼濾波進行線性化處理。

由系統方程

(12)

Zk=h(Xk,k)+Vk

(13)

(14)

令

(15)

(16)

則狀態方程為:

(17)

同理，對非線性函數h(·)進行處理，最終得觀測方程為：

Zk=HkXk+yk+Vk

(18)

利用卡爾曼基本方程可得:

(19)

(20)

(21)

(22)

Pk=[I-KkHk]Pk,k-1

(23)

其中，濾波初值和濾波誤差的協方差矩陣為：

(24)

(25)

2.2 EKF模型及校正原理

根據上文，地磁場觀測值模型如下：

B=CH+b+P(ΔB/Δt)+ε

(26)

由于在進行姿態變換的過程中，B的不規則度較大，而地磁場真實值H的模值‖H‖不隨姿態變換而變化，且以B作為自變量，H作為因變量的模型解算出的參數在實際問題中可以不經過變換直接使用，因此在研究中以‖H‖2作為研究對象，將上述模型進行變換處理，得地磁場校正模型為：

Zk=‖H‖2=(Q(B-b-P(ΔB/Δt))+ε)2=
hk(Xk)+vk

(27)

式中：Q=C-1，則變量X為：

(28)

hk(Xk)=(B-b-P(ΔB/Δt))TQTQ(B-b-P(ΔB/Δt))

(29)

可見校正模型Zk是關于變量X的二次系統，而并非常規線性系統，采用EKF進行降階處理。

vk=(2Q(B-b-P(ΔB/Δt)))Tε+‖ε‖2

(30)

觀測噪聲vk統計模型服從高斯分布，其均值和方差分別為：

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

(31)

Σ(Q(B-b-P(ΔB/Δt)))+2trΣ2)

(32)

Σ={εεT}

(33)

由于在此模型中，不含過程噪聲，只有觀測噪聲，因此EKF模型可以簡化為：

(34)

(35)

(36)

3 模擬仿真分析

在仿真實驗中，根據上文地磁場觀測值的模型

B=CH+b+P(ΔB/Δt)+ε

其中ε為均值為0、方差為30的高斯觀測噪聲。設地磁場真實值H總量為50 000 nT，為了得到各個方向的磁場分量，設參數p=(-180:0.3:510)×π/180，q=(250:0.5:1400)×π/180，則地磁場分量值為：

根據上文校正模型利用EKF算法進行校正， 為保證較好的權值收斂性，令所求參數初始值X0全為0，初始協方差矩陣P0為：

P0=107×I21×21

表1 參數估計值

各個參數均能收斂，由于參數較多，選取部分參數收斂過程，如圖2～圖7所示。

圖2 q11收斂過程

圖3 q12收斂過程

圖4 q13收斂過程

圖5 p31收斂過程

圖6 p32收斂過程

圖7 p33收斂過程

經過補償，從數據值隨機抽取10組數據進行比較，結果如表2所示。

由表2中數據可以看出，對地磁場觀測值進行補償之后，誤差降低兩個數量級，測量精度得到明顯提高。經過補償，數據的標準差明顯減小且趨于穩定，如圖8所示，表明用EKF算法補償有效且穩定。

采用不考慮渦流磁場的模型，同樣用EKF算法進行補償，將所得結果與文中所得結果進行比較，如圖9所示，可以看到，不考慮渦流磁場模型可以對地磁場實現補償，但效果欠佳，存在400 nT左右的誤差，并且會上下波動，考慮渦流磁場的一體化模型將誤差降至80 nT以下，從而實現更高精度的補償。

表2 補償前后結果比較

圖8 補償前后標準差比較

圖9 不同模型補償結果比較

4 實驗結果

在該實驗中，采用Bartington公司的Mag-03MSL70三軸磁通門傳感器進行補償，該磁通門傳感器能夠精確測量靜態或交變磁場，采用National Instruments軟件進行數據采集。為了模擬飛行載體在飛行過程中的不同姿態，該實驗采用3FHT30C三軸手動無磁轉臺來提供姿態變換信息。后期用Matlab進行數據處理。采用標準磁圖確定地磁場真實值。

選取一處磁場穩定的郊外地點作為試驗點，將磁通門傳感器置于無磁轉臺，快速轉動并進行全方位姿態信息測量。通過對比磁圖得到所處位置地磁場真實值總量為42 587 nT，將NI軟件所采數據套用文中所提出的模型，并采用EKF算法進行誤差補償，補償結果如圖11所示，在補償前，由于各種因素的影響，誤差幅值達到1 800 nT左右，在用不考慮渦流磁場的模型補償后，誤差幅值明顯減小，但峰值誤差仍能達到500 nT，并且補償結果不穩定，存在波動趨勢，在用考慮渦流磁場的一體化模型補償之后，補償結果明顯趨于穩定，峰值誤差減小到120 nT，證明了該補償模型和補償方法在實際應用中的有效性。

圖10 三軸磁通門傳感器和無磁轉臺

圖11 三軸磁通門傳感器校正結果

5 結論

文中采用EKF算法對地磁場觀測值進行補償，在建模過程中，在考慮磁力儀測量誤差和軟硬磁誤差的基礎上，結合飛行載體頻繁進行姿態變換的實際情況，加入渦流磁場誤差進行一體化建模，相比于以往研究，提出的方法具有以下優點：1)一體化建模方式對于進行頻繁姿態變換的載體，能夠實現更高精度的補償；2)卡爾曼濾波算法相比于最常用的最小二乘算法，無須事先采集大量數據進行解算，對數據量敏感度低，具有較強的魯棒性，并且可以實現即時補償；3)EKF相比于其他非線性算法，原理簡單，易于實現，降低了解算難度，有效縮短了補償時間。

在實驗中，相比于不考慮渦流磁場的模型，一體化建模的模型在有效減小誤差的基礎上，濾除了姿態變換所引起的渦流磁場對于補償結果的波動干擾，極大增強了補償結果的穩定性。補償精度的提升和穩定性的增加使得地磁導航在飛行載體的導航中有了更加廣闊的前景。考慮到EKF算法不適用于強非線性問題，下一步將針對于在其線性化過程中忽略的二次及更高項所引起的誤差進行研究。