田峰敏，卞雷祥，王宏波

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

0 引言

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

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

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

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

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

1 系統(tǒng)結(jié)構(gòu)

如圖1所示，組合導(dǎo)航仿真系統(tǒng)按功能結(jié)構(gòu)共分為航跡生成、慣導(dǎo)、航行磁圖生成、地磁測量數(shù)據(jù)生成、地磁推算數(shù)據(jù)生成、濾波器6個模塊。

圖1 地磁輔助慣性組合導(dǎo)航仿真系統(tǒng)結(jié)構(gòu)圖

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

2)“慣導(dǎo)模塊”根據(jù)航跡點生成慣導(dǎo)運動參數(shù)及指示航跡；

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

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

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

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

2 方法

2.1 位場延拓

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

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

圖2 海拔4 km處EMAG2原始數(shù)據(jù)三維圖

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

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

2.2 地磁測量數(shù)據(jù)和地磁推算數(shù)據(jù)生成

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

2.3 系統(tǒng)誤差模型及求解

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

1)狀態(tài)量

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

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

(1)

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

2)狀態(tài)方程

(2)

3)測量方程

(3)

4)離散化的系統(tǒng)狀態(tài)方程和量測方程

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

(4)

(5)

當(dāng)采樣周期Δt較小時，式(2)中的F(t)和G(t)可以看作是非時變的，此時狀態(tài)轉(zhuǎn)移矩陣Φ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)推導(dǎo)了地磁輔助慣性導(dǎo)航系統(tǒng)的狀態(tài)方程和量測方程，其中量測方程為非線性方程，采用無跡卡爾曼濾波(unscented kalman filter, UKF)作為系統(tǒng)誤差模型的求解方法。

3 仿真驗證

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

表2 地磁導(dǎo)航仿真結(jié)果

圖5 地磁修正后航跡與3組慣導(dǎo)指示航跡的對比

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

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

4 結(jié)論

設(shè)計了基于全球地磁異常場的地磁/慣性組合導(dǎo)航仿真系統(tǒng)，經(jīng)過仿真試驗驗證，得到兩個結(jié)論：

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

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

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