基于地磁方向熵的高空地磁匹配航迹选取方法

邓小波，赵军瑞，宋 胜

(1.中国人民解放军61711部队，新疆 喀什 844000；2.中国电子科技集团公司第五十四研究所，河北 石家庄 050081)

0 引言

地磁场是地球本身固有的基本物理场[1]。在主地磁场基础上，不同位置的地形和组成物质不同，导致局部的异常场[2-4]，而且异常场与空间地理位置存在一定的对应关系。基于异常场匹配的地磁匹配导航是一种无源、自主的导航方式，具有良好的应用前景[5-8]。

在地磁导航算法的研究中，穆华等人[9]采用了扩展卡尔曼滤波的方法研究了惯性/地磁组合导航技术在船舶导航方面的应用，仿真结果证明惯性/地磁导航鲁棒性好，可实现长航时工作；杨功流等[10]研究了基于Sigma点采样的UKF方法构建惯性/地磁卡尔曼滤波器，相比于使用EKF方法避免了非线性系统的线性化误差；朱占龙等[11]研究了使用新息正交法的自适应滤波方法，仿真结果证明该方法提高了组合导航系统对野值的适应能力，提高了系统的可靠性；刘睿[12]研究了惯性/地磁/天文的组合导航方法，综合了多种导航方式的优点，构造了长航时高精度的组合导航系统。上述研究主要针对地面、水下以及低空域进行地磁匹配算法研究，然而高空高速飞行器由于地磁场环境和运动特性有所不同，需要对其导航方法进行针对性研究。

地磁图的准确性和包含可用信息量的大小是影响地磁匹配导航精度的重要因素。而不同位置地磁异常场特性不同，能用于匹配导航的信息量也不同。因此，进行地磁匹配导航前需要选择适合匹配的区域[13-16]。同时，穿越地磁区域的航迹方向也会影响到地磁匹配的效果，这也是进行载体航迹规划的重要依据[17-18]。

本文根据对高空地磁图特性的分析，利用高空大范围地磁梯度呈现趋势性的特点，提出了地磁方向熵的概念，由此给出了一种高空高速飞行器的地磁匹配轨迹选取方法。

1 高空地磁图匹配的特性分析

目前，高空地磁异常场测量在实现上面临较大困难，通常通过采集地面或低空地磁异常场数据，然后通过延拓方法得到高空地磁场的方案[19]。这使得高空地磁异常场的特性与地面或低空地磁异常场有所不同[20]。取1 000 m高度地磁图和30 000 m高度地磁图绘制三维曲面图，如图1所示。

(a) 1 000 m

(b) 30 000 m图1 同一区域不同高度的地磁场三维曲面图Fig.1 3D surface map of the geomagnetic field at different altitudes in the same area

通过对比图1不同高度的地磁三维图可以看出，高空地磁异常场的变化相比低空更加平滑。因此，高空地磁匹配算法需要适应高空地磁图的特点。地磁异常场的削弱意味着对地磁传感器测量精度的要求比地面地磁匹配要高。而高空地磁场这种平滑有明显趋势性变化的特性也会对地磁匹配方法提出新的要求和需要解决的问题，即飞行器飞入方向和飞行轨迹问题。高空地磁异常图的波动较小，起伏变化具有明显的与地理位置有关的趋势性，这使得飞行轨迹的航向和等值线变化趋势的角度会影响到匹配的性能。飞行轨迹与地磁等值线处于近平行状态时，匹配性能较好，而处于近垂直状态时匹配误差增大，且匹配概率大大降低。基于以上技术背景，提出地了磁梯度网络和地磁方向熵结合的地磁匹配航迹选择方法。

2 数字地磁图的地磁梯度向量场

然而，实测的地磁图一般都是离散的，以二维数组的形式存储在计算机中。二维数组可以以更加简单的方式计算地磁梯度，即使用地磁矩阵的数值梯度(Numerical Gradient)来计算数字地磁图沿着x方向和y方向的地磁场梯度。

矩阵的数值梯度的定义为二维函数在x方向和y方向的偏导，其值为：

(1)

对于矩阵这种离散数据，可以通过计算差分的方式来取代求偏导。例如对于单位间距的矩阵A，其内部的数据可以采用中心差分的方式计算梯度。内部梯度值G(:,j)为：

G(:,j)=0.5×(A(:,j+1)-A(:,j-1))，

(2)

式中，j∈(2～N-1)，对于矩阵边缘的数据，可以采用单边差分的方法计算梯度：

(3)

通过差分计算的方法可以得到离散的数字地磁矩阵关于x方向和y方向的梯度分量。其向量和为当前位置的地磁梯度向量，该向量的方向即为地磁场强度增大的方向。该向量的角度可利用x梯度与y梯度按照四象限反正切的方法计算出来：

(4)

式中，Y和X表示地磁场y方向和x方向地磁场梯度值;θ表示该向量按照北偏东为正方向的角度。将地磁梯度向量以箭头的形式，在地磁图内绘制出来，就称为向量场地磁图，或向量网格地磁图。以高空30 000 m地磁图的一个区域为例，其梯度向量场地磁图如图2所示，向量场地磁图中的网格数为100×100，横纵坐标分别表示网格的序号。因较为密集，图中给出了部分区域的放大。

图2 高空30 000 m区域的梯度向量场地磁图Fig.2 Gradient vector magnetic field map at 30 000 m altitude

可以看出，向量场中的向量箭头指向的是地磁场强度增大的反向，即由深蓝色等值线变化到黄色等值线，与右侧色标相对应。通过将式(1)中的四象限反正切公式改为：

(5)

由此画出的地磁梯度网格如图3所示。

图3 高空30 000 m区域的平行等值线向量场地磁图Fig.3 Parallel isoline vector magnetic field map at 30 000 m altitude

相比式(4)，式(5)求解出的θ结果相当于逆时针旋转90°，变为与地磁场强度梯度垂直、与地磁等值线平行的曲线；同样，可以从图3放大部分看出，由式(5)绘制出的向量场地磁图中的箭头与等值线方向平行。

随着高空地磁场随高度衰减，地磁等值线的变化趋势也逐渐趋于平缓。换言之，地磁场等值线就是高空地磁图中逐渐在大范围内趋于平行的直线。高空30 000 m的最大地磁信息熵区块中地磁图区块的地磁图等值线就符合这种特性，整个图内的等值线几乎都沿着一个方向呈现平行分布的趋势。因此，在该地磁区域内，地磁匹配方法应当适应这种特性。常用的MSD方法在航迹与等值线呈现过大角度，例如接近垂直于等值线时，平行于等值线方向的轨迹的惯导误差较大，但这些误差无法通过MSD方法匹配到真实位置附近，地磁匹配几乎失效。而当航迹平行于等值线方向时，惯导误差就能通过MSD方法匹配到，匹配的精度就能大大提高。

3 地磁方向熵

向量场地磁图表现了地磁数据中每个位置地磁场等值线的方向情况。实际上，为了衡量一个地磁图区域内整体性的地磁场等值线的复杂程度，可以构建出地磁方向熵的概念，其公式为：

(6)

该式的物理意义为：通过式(4)和式(5)可以计算出向量场中每个地磁梯度或等值线平行矢量的方向(该方向以北为起点，北偏东顺时针为正，范围为-180°～180°)。将这些矢量的角度排序后成为矢量角度数据集，将-180°～180°的范围划分为K个等跨度区间。矢量角度数据集落在对应K个等跨度区间内，每个区间的数量为ni。按照式(6)计算得出的Hdir就是对应地磁场区域的地磁方向熵。

从地磁方向熵的物理意义来看，地磁方向熵反映出地磁场内等值线的变化复杂程度。地磁方向熵越小，地磁矢量角的分布越集中，这意味着该区域内的地磁场等值线在大范围内呈现出平行直线的形状，其角度方向统一。地磁方向熵越大，说明地磁场等值线在大范围内呈现出较为丰富的形状变化，或大范围内地磁等值线形成闭环。对高空地磁图来说，其本身伴随着地磁等值线峰值的削弱，地磁等值线越平直，越适合降低地磁匹配算法的匹配误差。

4 仿真与分析

仿真实验在图2所示的地磁区域内进行。通过地磁方向熵计算方法，将地磁区域内的所有地磁矢量角划分为36个等跨度区间，计算得到地磁方向熵值为0.126 6。这个值意味着地磁等值线中的地磁矢量角高度集中于跨度区间[40°，50°]，即和等值线近似平行的角度。根据预测，航向角属于该区间内的轨迹的匹配效果更好。为了证明这个预测结果，设计3条航迹，分别为平行等值线、呈45°和垂直等值线来验证仿真结果。根据地磁方向上设定3条运动轨迹的飞入点都在(38.39°N,115.45°E)，飞行速度为3 400 m/s。地磁等值线的方向大约为北偏东45°，本文选择的3条运动轨迹由飞入点进入后的飞行方向与地磁等值线间的初始夹角分别选为：航迹① 0°(近平行状态)、航迹② 45°(中间状态)和航迹③ 90°(近垂直状态)；对应的航向分别为：航迹① 北偏东50°、航迹② 北偏东90°和航迹③ 北偏东130°。3条运动轨迹如图4所示。

图4 3条预设仿真轨迹Fig.4 Three preset simulation tracks

与3条运动轨迹相对应，设定的惯导系统初始误差、惯性器件误差和地磁误差如表1所示。

表1 仿真误差设置Tab.1 Simulation error setting

飞行器的仿真导航算法采用惯性/地磁匹配松组合，将MSD法计算得到的匹配位置输入松组合卡尔曼滤波器，用以修正惯导轨迹。

飞行器分别按照以上3条运动轨迹运动，运动过程中匹配定位结果分别如图5、图6和图7所示。

(a) 原始大小

(b) 部分放大图5 轨迹1的匹配效果图Fig.5 Matching result of track 1

(a) 原始大小

(b) 部分放大

(a) 原始大小

(b) 部分放大

图5、图6和图7的(a)图为匹配区域内的航迹的匹配定位图，(b)图为相应3幅图中黑色细实线框内小段轨迹的放大示意图。各幅图中，红色圆圈为轨迹真值位置，绿色实线为纯惯导输出位置，黑色星号点为地磁匹配定位结果，蓝色细实线为地磁匹配算法校正惯性导航后输出的轨迹。

图5中，轨迹①匹配效果良好，定位结果基本和真值重叠，匹配误差为530.19m(1σ)。轨迹②匹配结果接近真实位置，但其定位误差较大，匹配误差为2 051.11m(1σ)。轨迹③地磁匹配定位结果逐渐偏离真值，呈发散趋势，仿真中匹配定位误差达到4 162.34m(1σ)。

5 结束语

针对高空地磁图等值线呈现趋势性和全局随机性下降的现象，本文结合地磁信息熵与地磁方向熵2种量化指标，研究了高空地磁匹配区域飞行航迹选取方法。由仿真可见，在轨迹与等值线呈近平行状态时，其匹配定位误差较小，随着运动轨迹与等值线夹角增加，匹配定位误差加大，且在近垂直状态时其匹配结果可信度很差。数据处理结果说明了应用地磁方向熵选择飞行轨迹方向对提高地磁匹配定位精度具有重要作用和效果。