地磁日变影响下的地磁匹配算法

解伟男，李清华，屈桢深，奚伯齐

（哈尔滨工业大学空间控制与惯性技术研究中心，哈尔滨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.