地磁日變影響下的地磁匹配算法

解偉男，李清華，屈楨深，奚伯齊

（哈爾濱工業大學空間控制與慣性技術研究中心，哈爾濱150001）

地磁場為近地空間提供了一個天然的導航坐標系，許多動物能夠利用地磁信息進行導航，如海龜[1]、鳥類[2]等。地磁匹配作為地磁導航中一類主要方式，由于其不需要對地球局部地磁場建模，而是在地磁數據庫中進行線圖匹配，因此特別適用于低空的導航定位。此外，地磁導航具有無源、無輻射、誤差不隨時間累積等優點，與慣性導航組合可以實現兩種導航方式的優勢互補，達到長時間高精度自主導航的目的，具有重要的軍事意義[3-5]。

地磁日變是地磁場的一個基本組成部分[6-7]，其起源于大氣電離層中一直存在的電流體系，在地磁導航中相當于地磁測量時的干擾。由于地磁日變的平均變化幅度約為幾十納特，因此會對地磁匹配導航定位精度和可靠性產生嚴重的影響[8]。然而地磁日變并未形成較為統一的數學模型，不同地點的地磁日變曲線有較大差異，同一地點的地磁日變曲線隨著季節不同也不盡相同，在地磁匹配過程中難以進行有效的補償或剔除，因此只有少數學者考慮了地磁匹配中的地磁日變影響。文獻[9]指出了地磁日變是地磁匹配中的一個主要噪聲源，并通過數值仿真比較了不同的相關準則對噪聲的適用性。文獻[7]采用FMI方法擬合地磁日變場，在地磁測量值中減去地磁日變場得到實時圖，該方法不需選取基線值，減小了人為誤差，所提取的日變更客觀、真實，然而該方法需以擬合當天及前后共三天的數據為計算基礎，在實際應用中實現困難。

本文在文獻[10]的基礎上，提出一種地磁日變影響下的地磁匹配算法。該算法根據地磁日變的特性，在每一段地磁匹配過程中，把地磁日變看成一個恒定的磁場偏差，在此基礎上構建地磁日變影響下的均方差相關性準則，引入匹配曲線的參數化模型，將地磁匹配問題轉化為非線性方程組的求解，最終通過快速迭代實現地磁日變影響下的匹配定位。該算法可以有效地消除地磁日變對匹配精度的影響，并具有較高的實時性。

1 地磁場的組成和地磁日變

地磁場根據其性質不同可以劃分為穩定磁場和變化磁場兩部分：

其中：T表示地磁總磁場；0T表示穩定磁場；TΔ表示變化磁場。變化磁場相對于總磁場來說很弱，在特殊的磁暴條件下也只有總磁場的2%～4%，而多數情況僅有總磁場的1%以下。

穩定磁場根據其性質不同可以劃分為主磁場和異常場兩部分：

其中：mT表示主磁場；aT表示異常場。主磁場主要是由地核內電流的對流形成，其強度約為穩定磁場的94%；異常場主要是由地殼淺部具有磁性的巖石或礦石引起的，其強度約為穩定磁場的6%。

變化磁場根據其變化特點又可分為平靜變化和擾動變化。平靜變化起源于電離層中一直存在著的穩定電流體系，其特點是變化連續，持續存在，并具有一定周期特性，平靜變化包括太陽靜日變化Sq和太陰日變化 L。擾動變化起源于太陽噴射出來的帶電粒子流所產生的多種短暫的電流體系的變化而產生的，其特點是缺乏長期連續性，出現時間不規則，變化形態復雜，擾動變化包括磁暴、亞磁暴、地磁脈動等。

地磁匹配導航依據不同地點的地磁場強度差異進行匹配定位。雖然變化磁場強度只占地磁總場強度很小一部分，但是變化磁場強度相對于不同地點的地磁變化梯度來說卻很巨大，會嚴重影響匹配定位精度或產生誤匹配，因此在匹配定位時必須予以考慮。擾動變化是一種復雜、不規則的短期磁干擾，應避免在擾動變化條件下進行地磁導航，因此本文不考慮擾動變化，只考慮平靜變化中的太陽靜日變化Sq和太陰日變化L。

太陽靜日變化Sq簡稱靜日變化或日變，變化周期為1個太陽日（24 h），其幅度約為幾十nT。太陰日變化L變化周期為1個太陰日（約為25 h），其幅度很微弱，僅有1～2nT，很難將太陰日變化從太陽靜日變化中分離出來，因此通常將太陰日變化包括在太陽靜日變化之內。

太陽靜日變化Sq具有如下特點：

1）Sq場隨著地方時變化而周期改變；

2）Sq場白天變化劇烈，夜晚變化平緩；

3）Sq場隨著緯度變化而改變；

4）Sq場隨著季節變化而改變，夏季變化幅度較大，冬季變化幅度較??；

5）Sq場的不同分量在赤道兩側對稱分布或反對稱分布。

2 地磁日變影響下的地磁匹配算法

地磁日變是地磁場的一個基本組成部分。地磁日變并未形成較為統一的數學模型，不同地點的地磁日變曲線有較大的差異，同一地點的地磁日變曲線隨著季節不同也不盡相同，因此如何實現地磁日變影響下的地磁匹配定位具有重要的意義。

2.1 問題描述

地磁匹配可以表示為平面上兩條曲線的相關性計算，如圖1所示的地磁匹配示意圖：曲線Rf表示慣導系統輸出的參考軌跡；曲線M為匹配算法求得的匹配軌跡；曲線Rf′與匹配軌跡M 平行，參考軌跡Rf與曲線Rf′之間的夾角為α，即航向誤差角。從圖中可以看出參考軌跡與匹配軌跡間存在著平移和旋轉，對應著慣導系統的初始位置誤差和初始航向誤差。令點為參考軌跡上的第i個采樣點，其對應的匹配軌跡點為(u ,v )T，其中 a 、u表示點的經度，b、v表示點的緯度。

圖1 地磁匹配示意圖Fig.1 Schematic diagram of geomagnetic matching

2.2 性能指標

若采用均方差準則（MSD）計算曲線的相關性，需計算匹配軌跡曲線上各點所對應的地球磁場特征值與真實磁場測量值之差的平方和的均值。然而如前文所述，若存在地磁日變時，該相關性計算方法并不能真實反映曲線的相關性，需要對地磁日變進行補償。根據地磁日變特性可知，變化磁場TΔ是時間的函數。在不考慮擾動變化的條件下，相對于地磁匹配過程中的每一段匹配周期而言，地磁日變變化緩慢。因此，在每一段軌跡的匹配過程中，地磁日變相當于在匹配局部區域疊加了一個近似恒定的磁場偏差，由此可得地磁日變影響下的MSD相關性準則：

考慮到參考軌跡在匹配軌跡附近，即可將 I(ui, vi)泰勒展開：

將泰勒展開式(4)帶入式(3)，并忽略高階小項 O2，可以得到如式(5)的相關性約束：

在不引起混淆的情況下，可將相關性約束(5)寫成如下簡化形式：

其中，

與文獻[10]相比，相關性約束(6)中不僅存在 ui和vi兩組未知變量，而且存在未知的地磁日變補償場δ，因此式(6)為地磁日變影響下的相關性約束。只要找到相應的 ui、vi、δ使相關性約束(6)取極小值即可求得地磁日變影響下的匹配結果。

為了實現匹配結果的快速求解，建立了匹配軌跡的參數化模型，并將該模型引入相關性約束，從而將遍歷搜索方法轉化為快速迭代算法。如圖1所示，當慣導系統存在初始位置誤差和初始航向誤差時，匹配曲線M的參數化模型為

將式(8)帶入相關性約束(6)，可以得到地磁日變影響下的均方差約束指標函數：

因此地磁匹配問題可以轉化為尋求Δx、Δy、α和δ，使指標函數(9)取極小值。

2.3 匹配求解

為使指標函數(9)取極小值，可分別將指標函數對自變量xΔ、yΔ、α、δ求一階偏導數，并令其為零，即：

將均方差約束指標函數(9)帶入式(10)可得：

因此地磁匹配可以轉化為由xΔ、yΔ、α、δ組成的非線性方程組(15)的求解問題，該非線性方程組可通過Broyden迭代法求解。匹配結果可以通過式(8)計算得到。

3 迭代計算

地磁日變影響下的地磁匹配算法的實現步驟如下：

第1步：將參考軌跡離散化成N個點，對應離散點的坐標為 ( ai, bi)T， i = 1 ,2,…,N ；讀取載體在離散點(a,b )T時磁傳感器實時測量得到的地磁特征值

第2步：根據參考軌跡離散化的N個位置，從預先存儲的地磁數據庫中讀取該位置的地磁特征值，并根據地磁數據庫信息計算該位置的地磁特征值梯度信息 Ix,i和 Iy,i。

第3步：初始化迭代變量

第4步：根據式(17)～(19)解算迭代參數G、F和H：

其中，

第5步：計算位置誤差的增量、航向誤差的增量和地磁日變誤差的增量δM：

第6步：更新位置誤差、航向誤差和地磁日變誤差M：

第7步：判斷是否滿足終止迭代條件，若滿足則停止迭代并跳到第10步，否則跳到第8步。

終止迭代的條件有兩個，滿足任意一個即終止迭代：1）迭代次數達到預設的最大迭代次數；2）迭代誤差增量δM 的2范數小于設定值，即

其中，0ε＞為預先給定的迭代最小誤差。

第8步：根據更新后的M 計算參數K和δK：

第9步：根據公式(37)和公式(38)更新迭代變量H和G，然后跳到第5步。

第10步：根據迭代計算得到位置誤差、航向誤差和地磁日變誤差xΔ、yΔ、α、δ，將所得到結果代入匹配曲線的參數化模型(8)即得到匹配軌跡。

根據上述匹配步驟可知，迭代計算過程中地磁特征值梯度信息Ix,i和Iy,i并不改變，因此可以將匹配區域地磁特征值梯度信息事先存入載體計算機中，匹配計算前直接在載體計算機中檢索得到，也可以通過檢索載體計算機中的地磁特征值信息后計算得到。

需要指出的是，地磁日變影響下的地磁匹配算法假設參考軌跡位于匹配軌跡附近，在忽略高階小項2O的基礎上，以一階差分的形式近似相關性約束。然而當慣導系統的初始誤差較大時，忽略的高階小項2O會嚴重影響算法的精度。為了克服上述問題，可以采用多次迭代的方式提高匹配算法的精度。在多次迭代時，可以將當前的迭代結果作為下一次迭代的參考軌跡。

4 仿真實驗

為了驗證上述算法的有效性，本節對算法進行仿真驗證。根據地球磁場異常網格（EMAG2）繪制某局部區域地磁異常圖作為背景磁圖，其中地球磁場異常網格間距為2′，分辨率為0.01 nT，主要仿真參數設定如表1所示。此外設定載體初始經度為127.00°，初始緯度為45.05°，以50 m/s的速度向北運動，慣性導航系統的解算頻率為100 Hz，地磁場參數的測量頻率為0.025 Hz。初始誤差假設為：經度誤差0.014°，緯度誤差0.014°，航向誤差3°。圖2給出了真實軌跡和慣導系統輸出的參考軌跡。

表1 仿真參數Tab.1 Main parameters of the simulation

圖2 真實軌跡和參考軌跡Fig.2 The real track and the reference track

仿真時首先假設不存在地磁日變，圖3和圖4分別給出了仿真結果的經度誤差曲線和緯度誤差曲線，其中算法 1為文獻[10]中二次迭代的基于迭代計算地磁輪廓線匹配算法，算法2為本文所提出的地磁日變影響下的地磁匹配算法，算法3為本文所提出的二次迭代的地磁日變影響下的地磁匹配算法。為了驗證本文所提出算法對地磁日變的有效性，分別在磁場測量數據中加入10 nT和50 nT地磁日變。圖5～6和圖7～8分別給出不同地磁日變影響下的誤差曲線。表2給出了不同仿真條件下的匹配結果。根據匹配結果可以看出：在沒有日變影響時，本文所提出的算法與文獻[10]算法的匹配精度相當；當地磁日變為10 nT時，本文所提出算法的匹配誤差僅為文獻[10]算法匹配誤差的27.3%；當地磁日變為50 nT時，文獻[10]所提出的算法發生了誤匹配，而本文所提出的算法由于在相關性準則中考慮了地磁日變的影響，依然能夠保證較高的匹配定位精度。

圖3 無地磁日變時匹配經度誤差曲線Fig.3 Longitude errors of matching results without geomagnetic diurnal variation

圖4 無地磁日變時匹配緯度誤差曲線Fig.4 Latitude errors of matching results without geomagnetic diurnal variation

圖5 10nT地磁日變時匹配經度誤差曲線Fig.5 Longitude errors of matching results with geomagnetic diurn al variation10 nT

圖6 10 nT地磁日變時匹配緯度誤差曲線Fig.6 Latitude errors of matching results with geomagnetic diurnal variation10nT

圖7 50nT地磁日變時匹配經度誤差曲線Fig.7 Longitude errors of matching results with geomagnetic diurnal variation50nT

圖8 50nT地磁日變時匹配緯度誤差曲線Fig.8 Latitude errors of matching results with geomagnetic diurnal variation50nT

由匹配結果可以看出，地磁日變影響下的地磁匹配算法可以有效補償地磁日變對匹配精度的影響，實現地磁日變影響下的高精度地磁匹配。此外，該算法采用多次迭代，匹配精度可以進一步提高。利用Matlab中的tic和toc語句計算匹配算法的運行時間，采用二次迭代的地磁日變影響下的地磁匹配算法在仿真計算機中總耗時為 16 ms。因此所提出的算法具有較高的實時性，可以保證地磁匹配導航的在線計算。

表2 匹配結果Tab.2 Matching results

5 結 論

本文提出了一種地磁日變影響下的地磁匹配新方法。該方法通過建立地磁日變影響下的相關性準則函數，補償地磁日變對匹配結果的影響。仿真結果表明：在無地磁日變時，所提出算法的匹配誤差與文獻[10]相當；當地磁日變為10 nT時，所提出算法的最大匹配誤差為文獻[10]的27.3%；當地磁日變為50 nT時，文獻[10]發生誤匹配，而本文所提出算法仍然可以實現高精度的匹配。此外，所提出的算法在Matlab仿真中的總耗時為16 ms。因此提出的地磁日變影響下的地磁匹配算法可以有效地補償地磁日變對匹配結果的影響，實現高精度的實時地磁匹配。

（

）：

[1] Brothers J R,Lohmann K J.Evidence for geomagnetic imprinting and magnetic navigation in the natal homing of sea turtles[J].Current Biology,2015,25(3):392-396.

[2] Qi X K,Ye D X,Sun Y Z,et al.Simulations to true animals’long-distance geomagnetic navigation[J].IEEE Transactions on Mag-netics,2017,53(1):1-8.

[3] Guo C F,Cai H,van der Heijden G H M.Feature extraction and geo-magnetic matching[J].The Journal of Navigation,2013,66(6):799-811.

[4] Wang P,Wu Y X,Hu X P,et al.Geomagnetic aided navigation suitability evaluation based on principal component analysis[C]//Inter-national Conference on Industrial Control and Electronics Engineering.2012:324-329.

[5] Xie W N,Qu Z S,Li Q H.A fast algorithm of the geomagnetic correlation matching based on MSD[C]//The3rd International Conference on Control,Automation and Systems Engineering.2013:59-62.

[6] Nizam M,Abdullah E J.Nanosatellite navigation with the WMM2005geomagnetic field model[J].Turkish Journal of Engineering&Envi-ronmental Sciences,2006,30:43-55.

[7] 鄧翠婷.地磁匹配中的地磁日變效應與修正[D].南京:南京航空航天大學,2013.Deng C T.Research on affection and modification of magnetic diurnal variation for geomagnetic matching[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2013.

[8] 黃朝艷,鄧翠婷,趙華.地磁匹配中的地磁日變效應[J].測繪工程,2013,22(4):17-21.Huang C Y,Deng C T,Zhao H.The affection of magnetic diurnal variation in geomagnetic matching[J].Engineering of Surveying and Mapping,2013,22(4):17-21.

[9] 喬玉坤,王仕成,張金生,等.基于相關分析的飛行器地磁匹配算法仿真研究[J].地球物理學進展,2009,24(2):721-727.Qiao Y K,Wang S C,Zhang J S,et al.Simulation of geomagnetic matching algorithm based on correlation analysis for flying objects[J].Progress in Geophysics,2009,24(2):721-727.

[10] 解偉男,李清華,奚伯齊,等.基于迭代計算的地磁輪廓線匹配算法[J].中國慣性技術學報,2015,23(5):631-635.Xie W N,Li Q H,Xi B Q,et al.Geomagnetic contour matching algorithm based on iterative method[J].Journal of Chinese Inertial Technology,2015,23(5):631-635.