基于全球地磁異常場的自主導航仿真系統

田峰敏，卞雷祥，王宏波

(1 南京工程學院，南京 211167；2 南京理工大學，南京 210094)

0 引言

地磁場分為主磁場、異常場和變化場，其中異常場在時間上極為穩定，幾乎不隨時間變化，含有比地磁主磁場更豐富的細節信息。相比重力測量，地磁測量對載體的運動沒有限制[1]，測量設備價格適中，適合作為修正慣導誤差的信息源。

Goldenberg介紹了Goodrich公司在飛機上使用磁通門磁力計進行地磁導航的總體方案，并對地磁圖制作、載體磁補償等問題進行了分析[2]。Kato等提出了基于地磁圖和等深線圖的自主潛航器導航定位算法，并采用相關海域的地磁場數據和水深數據進行驗證性分析,得出采用地磁圖和等深線圖直接定位的效果與用來修正INS的定位效果差不多，且前者受水流影響更小的結論[3]。Canciani利用北美地磁異常數據庫和地磁調查飛機的實測數據，進行了地磁異常輔助慣性導航的實驗，實驗表明地磁圖分辨率越高、飛行高度越低、飛行速度越快，定位精度越高，在加州地區根據飛行高度和速度的不同，定位的均方根誤差從幾十米到幾百米不等[4]。穆華等用地磁場作為基準圖，提出了基于模標定的地磁測量誤差補償技術，并采用兩級濾波算法，底層濾波器估計慣導解算的位置誤差，然后將其作為上層濾波器的觀測輸入，用來估計慣導系統誤差。水下實驗表明，地磁測量誤差補償技術能夠將干擾磁場從 10 000 nT降低至100 nT，能夠將600 m的初始位置誤差降低至300 m[5]。劉穎等對無人機應用地磁匹配定位和地磁貫續導航進行研究，并設計了原理樣機，由于條件限制，并未進行飛行實驗，進行了陸面跑車實驗，平均定位誤差小于2個網格，湖面實驗的平均定位誤差小于1個網格[6]。呂志峰等設計了地磁匹配導航半實物仿真系統的總體方案，系統分析了地磁場數據庫仿真、載體干擾磁場補償、地磁場環境仿真和地磁匹配導航算法4種關鍵技術，并提出了相應的工程解決方案[7]。

由于缺乏全球性的地磁異常數據[2]，目前關于地磁導航研究主要集中在方案設計和短距離實驗驗證，對于采用全球地磁異常場的地磁導航還沒有文獻進行討論。

2009年美國NOAA下屬的全國地球物理資料中心對衛星測量、航海測量、航空測量的地磁數據進行整理，發布了分辨率為2′的全球地磁異常網格(global earth magnetic anomaly grid，EMAG2)，最新版本為2017年發布的EMAG2_v3[8]，EMAG2為在全球范圍內進行地磁導航提供了潛在的數據支持。

文中建立了以EMAG2為基準圖的慣性/地磁組合導航仿真系統，用來研究地磁輔助慣性導航的系統模型及算法，為半實物仿真和原理樣機的研制提供技術驗證。

1 系統結構

如圖1所示，組合導航仿真系統按功能結構共分為航跡生成、慣導、航行磁圖生成、地磁測量數據生成、地磁推算數據生成、濾波器6個模塊。

圖1 地磁輔助慣性組合導航仿真系統結構圖

1)“航跡生成模塊”根據載體的運動指令生成航跡點；

2)“慣導模塊”根據航跡點生成慣導運動參數及指示航跡；

3)“航行磁圖生成模塊”根據載體的運行高度，由基礎地磁異常圖延拓生成航行地磁異常圖；

4)“地磁測量數據生成模塊”把航跡點位置在航行地磁異常圖中插值生成地磁測量數據的異常場部分，把航跡點位置和時間戳代入國際地磁參考場(international geomagnetic reference field，IGRF)模型生成地磁測量數據主磁場部分，然后計算測量噪聲，三者之和為地磁測量數據；

5)“地磁推算數據生成模塊”把修正后位置在航行地磁異常圖中插值生成地磁異常推算數據，修正后位置和時間戳代入IGRF模型[9]生成主磁場推算數據，二者之和為地磁推算數據；

6)慣導模塊輸出的緯度、高度、速度、比力、姿態角信息傳送至狀態方程；把地磁推算數據和地磁測量數據傳送至觀測方程；由濾波器計算慣導誤差的修正值，并反饋至慣導模塊，得到載體位置的修正值。

2 方法

2.1 位場延拓

載體的運行高度與基礎地磁圖的制備高度可能不一致，解決這一問題的方法是位場延拓技術。事先根據航行計劃對基礎地磁數據進行向上延拓(upward continuation)或向下延拓(downward continuation)，使航行地磁異常圖的制備高度與載體的運行高度一致，延拓方法參見文獻[10]。

下面給出以EMAG2為基礎地磁異常圖的延拓試驗結果，圖2是制備高度海拔4 km處的EMAG2原始數據，圖3和圖4分別為經向上延拓和向下延拓得到的海拔4 000 m處的地磁異常與EMAG2原始數據之間的差異，都在±2×10-12nT之間。

圖2 海拔4 km處EMAG2原始數據三維圖

圖3 向上延拓海拔4 km處地磁異常的誤差

圖4 向下延拓海拔4 km處地磁異常的誤差

2.2 地磁測量數據和地磁推算數據生成

模擬的地磁測量值由IGRF模型計算的地磁主磁場值、由EMAG2插值生成的地磁異常值、測量誤差3部分組成。計算某一位置的地磁異常值，需要進行插值，目前常用的線性插值為C0連續，為了獲得更高的連續性，這里采用了自然鄰點插值算法[11]，由自然鄰點位置地磁異常的凸組合來表示插值點位置的地磁異常，獲得了在自然鄰點上是C0連續，在Delaunay三角形外接圓邊界上是C1連續，在定義域其它地方都是C∞連續的光滑性。

2.3 系統誤差模型及求解

在文獻[4]和[12-13]的基礎上，采用貫序濾波的方式實時估算慣導的位置誤差，系統的誤差模型由狀態傳遞方程和量測方程組成，選取15階狀態量，包括位置誤差、速度誤差、姿態誤差、加速計測量誤差、陀螺儀測量誤差。

1)狀態量

X=[δLδλδhδυeδυnδυuφeφnφu

?x?y?zεgxεgyεgz]

(1)

式中：[δLδλδh]為載體位置誤差；[δυeδυnδυu]為本地坐標系下速度誤差；[φeφnφu]為姿態誤差角，即數學本地系與真實本地系之間的夾角；[?x?y?z]為加速度計零偏；[εgxεgyεgz]為陀螺儀漂移[4,13]。

2)狀態方程

(2)

3)測量方程

(3)

4)離散化的系統狀態方程和量測方程

Xk=Φk,k-1Xk-1+Γk-1Wk-1

(4)

(5)

當采樣周期Δt較小時，式(2)中的F(t)和G(t)可以看作是非時變的，此時狀態轉移矩陣Φk,k-1、誤差傳遞矩陣Γk-1可由下式計算。

(6)

(7)

式中：E(Wk)=0,E(WkWj)=Qkδkj,Qk=Q(t)/Δt；E(Vk)=0,E(VkVj)=Rkδkj,Rk=R(t)/Δt。

5)濾波器

式(1)～式(7)推導了地磁輔助慣性導航系統的狀態方程和量測方程，其中量測方程為非線性方程，采用無跡卡爾曼濾波(unscented kalman filter, UKF)作為系統誤差模型的求解方法。

3 仿真驗證

測試條件：載體先以1g的加速度作20 s的加速，接著爬升至10°仰角并勻速飛行10 s，然后恢復平飛，再繼續加速10 s然后作2個90°轉彎，最后俯沖，運動持續392 s，飛行距離約139 km，飛行軌跡見圖5。采用了A、B、C三種不同精度等級的慣導，慣導及磁力計參數見表1，量測更新周期為1 s。仿真試驗基于EMAG2全球地磁異常數據，網格分辨率為2′(≈3.7 km)。表2是三組試驗結果的統計，圖5～圖6是仿真過程展示。

表2 地磁導航仿真結果

圖5 地磁修正后航跡與3組慣導指示航跡的對比

圖6 線性插值地磁導航和自然鄰點插值地磁導航對3組慣導誤差的跟蹤

從圖5和圖6可以看出，在純慣性導航條件下的位置誤差在200 s之后開始發散，最多達到22.5 km，采用地磁異常進行修正后，線性插值地磁導航和自然鄰點插值地磁導航都可以跟蹤慣導誤差，但后者位置誤差的最大值、均值、方差、均方根都優于前者。從表2可以看出，根據慣導精度的不同，采用自然鄰點插值地磁導航后，位置誤差最大值可控制在0.3～1.7個網格以內，位置誤差均方根可控制在0.1～0.7個網格以內。兩種插值濾波方法單點濾波耗時都小于0.1 s，滿足1 s的量測更新周期。

4 結論

設計了基于全球地磁異常場的地磁/慣性組合導航仿真系統，經過仿真試驗驗證，得到兩個結論：

1)基于全球地磁異常場的地磁輔助慣性導航方法，可以有效降低累積誤差，提高定位精度；

2)在地磁輔助慣性導航濾波器中，采用自然鄰點插值逼近任意位置地磁異常曲面比采用傳統線性插值方法能獲得更好的濾波效果，且運算耗時沒有明顯增加，適合作為地磁異常輔助慣性導航的濾波算法。

EMAG2全球地磁異常場由海測、航測、衛測多種數據加工處理而成，官方并沒有給出數據的具體精度，根據常理推斷在小范圍內其數據精度無法與實測數據相比，但作為大范圍的平均值，其精度還是可信的。為滿足自主導航的需要，這個范圍怎么取取多大，有待進一步的理論研究和實物驗證。