0 引言

地磁导航需要高精度的地磁模型以绘制地磁基准图（周能兵等，2018）。国际地磁参考场（international geomagnetic reference field，IGRF）是国际公认的地磁场模型，是广泛收集卫星、台站、航空和船载磁力数据的研究成果。当前最新版本是2015年发布的第12版，即IGRF-12（Thébault et al，2015），以球谐函数形式表达，具有195个系数。输入经度、余纬和地心距，可获得该点的磁场分量。根据获得的磁场分量数据可以绘制等值线，得到地磁基准图。

等值线的生成包括离散点网格化、等值点生成及等值线绘制等步骤（孙桂茹等，2000）。其中网格化方法包括矩形网格法和三角形网格法。前者易产生二义性（廖国忠等，2014）。后者在离散点分布不均匀时，多采用delaunay三角形网格法（Wan et al，2013；姚家振，2015），从包含所有离散点的原始三角网格出发，逐步迭代划分出分配合理的三角形网格，适用性好。如：苗润忠（2004）在分析有限元网格化算法后，提出在四边形网格形心上引入虚节点的三角形网格生成算法，取得较好的效果。若离散点均匀分布，则多采用基于构型表的三角形网格法（郑盛贵等，2004；侯凤贞等，2008）绘制等值线，首先将离散点分成均布的三角形网格，后依照已构建的构型表来绘制等值线。以上方法在构建网格时占用内存大，不利于IGRF模型的计算。

本文采用C#（Griffiths，2012；Albahari，2017），编程实现地磁基准图软件，以给定目标区域的经纬高数值为输入，解算IGRF模型获取地磁矢量，并计算地磁场要素；针对由IGRF模型生成的均匀分布的离散点集，提出具有虚节点的移动矩形网格算法，利用该方法，仅需1个矩形网格即可绘制给定区域的地磁等值线，具有效率高的特点。

1 地磁基准图磁场数据计算方法

地磁基准图的绘制需待测区域内的磁场要素值。将给定的经纬高区域均匀采样后转化为经度、余纬和地心距，解算IGRF磁场模型，得到并储存各点的磁场要素值。

1.1 IGRF模型

用IGRF模型描述的地磁场的磁位势函数为

$ \begin{array}{*{20}{l}} {V = {R_{\rm{e}}}\sum\limits_{n = 1}^\infty {\sum\limits_{m = 0}^n {{{\left( {\frac{{{R_{\rm{e}}}}}{r}} \right)}^{n + 1}}} } (g_n^m{\rm{cos}}(m\lambda ) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} h_n^m{\rm{sin}}(m\lambda ))P_n^m({\rm{cos}}\theta )} \end{array} $

(1)

其中，g_n^m、h_n^m为高斯系数；λ为格林尼治天文台起算的东经，范围为[0，360°)；θ为余纬，定义为90°- L，L为纬度，θ的范围为[0，180°]；r为地心距（坐标轴定义见图 1）；P_n^m为n次m阶的勒让德多项式；采用WGS-84参考椭球模型，赤道平面半径R_e= 6 378 137 m，1/f = 298.257。

图 1所示坐标系O-XYZ是将地球近似为球体时的北东地坐标系。根据位场转换原理，将磁位势分别向北、向东与向地求导，可获得其在3个方向上的磁感应强度，形式如下

$ \left\{ {\begin{array}{*{20}{l}} {{B_x} = - \frac{1}{r}\frac{{\partial V}}{{\partial \theta }} = \sum\limits_{n = 1}^N {\sum\limits_{m = 0}^n {{{\left( {\frac{{{R_{\rm{e}}}}}{r}} \right)}^{n + 2}}} } (g_n^m{\rm{cos}}(m\lambda ) + h_n^m{\rm{sin}}(m\lambda ))\frac{{{\rm{d}}P_n^m({\rm{cos}}\theta )}}{{{\rm{d}}\theta }}}\\ {{B_y} = - \frac{1}{{r{\rm{sin}}\theta }}\frac{{\partial V}}{{\partial \theta }} = \sum\limits_{n = 1}^N {\sum\limits_{m = 0}^n {\frac{m}{{{\rm{sin}}\theta }}} } {{\left( {\frac{{{R_{\rm{e}}}}}{r}} \right)}^{n + 2}}(g_n^m{\rm{cos}}(m\lambda ) - h_n^m{\rm{sin}}(m\lambda ))P_n^m({\rm{cos}}\theta )}\\ {{B_z} = - \frac{{\partial V}}{{\partial \theta }} = - \sum\limits_{n = 1}^N {\sum\limits_{m = 0}^n {(n + 1)} } {{\left( {\frac{{{R_{\rm{e}}}}}{r}} \right)}^{n + 2}}(g_n^m{\rm{cos}}(m\lambda ) + h_n^m{\rm{sin}}(m\lambda ))\frac{{{\rm{d}}P_n^m({\rm{cos}}\theta )}}{{{\rm{d}}\theta }}} \end{array}} \right. $

(2)

地磁场的总磁场强度T、水平强度H、磁倾角I和磁偏角D等的定义和计算参见Campbell（2003）的文章，文中不再赘述。

1.2 坐标系转换

通常目标区域的位置信息在地理坐标系下给定，即经度、纬度和高度；IGRF模型建立在地固坐标系下，解算它得到的是地磁场在该坐标系下的解。此外，地球是近似的椭球体，因此1.1节中的北东地坐标系还需修正地球扁率带来的数据差异。因此，要获得给定位置的地磁要素，需先将其坐标从地理系转换到地固系，之后解算IGRF模型，最后修正误差。坐标转换具体描述如下。

假设待测点的地理坐标为(λ，L，h)，其在地固坐标系中的坐标为

$ \left\{ {\begin{array}{*{20}{l}} {x = ({R_{\rm{N}}} + h){\rm{cos}}L{\rm{cos}}\lambda }\\ {y = ({R_{\rm{N}}} + h){\rm{cos}}L{\rm{sin}}\lambda }\\ {z = [{R_{\rm{N}}}{{(1 - f)}^2} + h]{\rm{sin}}L} \end{array}} \right. $

(3)

其中，R_N = R_e(1 + f sin² L)，是当地的卯酉圈半径。

利用所在地的直角坐标(x，y，z)计算经度、余纬和地心距，得

$ \left\{ {\begin{array}{*{20}{l}} {r = \sqrt {{x^2} + {y^2} + {z^2}} }\\ {y = {\rm{atan}} (y/x)}\\ {\theta = {\rm{acos}} (z/r)} \end{array}} \right. $

(4)

将式（4）的解输入IGRF模型，可得到北东地坐标系下的地磁要素，与实际的北东地坐标系结果仍有差别，需修正数据差异。

由式（3）—（4），在地心坐标系中的纬度L′可以表示为

$ {\rm{tan}}{L^\prime } = {\rm{tan}}(90 - \theta ) = \frac{{{R_{\rm{N}}}{{(1 - f)}^2}}}{{\sqrt {{x^2} + {y^2}} }}{\rm{sin}}L\frac{{{R_{\rm{N}}}{{(1 - f)}^2}}}{{{R_{\rm{N}}} + h}}{\rm{tan}}L $

(5)

根据图 2，将北东地坐标系绕y轴旋转δ（δ = L - L′），可与实际的北东地坐标系重合。北半球旋转公式如下

$ \left[ {\begin{array}{*{20}{c}} {{{\bar B}_x}}\\ {{{\bar B}_y}}\\ {{{\bar B}_z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\rm{cos}}\delta }&0&{{\rm{sin}}\delta }\\ 0&1&0\\ { - {\rm{sin}}\delta }&0&{{\rm{cos}}\delta } \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{B_x}}\\ {{B_y}}\\ {{B_z}} \end{array}} \right] $

(6)

将式（6）矩阵转置，可实现地磁矢量在南半球的旋转。

2 等高线算法设计

将上述获得的地磁要素统称为高程，记作h，可采用具有虚节点的移动矩形网格算法绘制等值线图。也就是，在矩形的形心引入虚节点，将其等分为4个三角形的网格，移动矩形网格并绘制等值线，直至覆盖整个区域，即获得完整的等值线图。

在矩形网格内绘制等值线前，需计算各等值点的高程h_i，取高程的最大值h_M和最小值h_m，则可得

$ {h_i} = {h_{\rm{m}}} + i \times ({h_{\rm{M}}} - {h_{\rm{m}}})/n $

(7)

其中，n为等高线数量，可由操作者给定；i = 1，…，n。等值线绘制步骤如下。

（1）构建初始矩形网格。取给定区域左上角的点(x₁，y₁)，与相邻的3个点(x₂，y₂)、(x₃，y₃)、(x₄，y₄)构成网格，其各自高程为h₁、h₂、h₃、h₄。如图 3所示，在矩形的形心引入虚节点，将其等分为4个三角形，按顺时针顺序标为①、②、③、④。通过式（8）—（10）求得虚节点的坐标和高程。

$ {{x_5} = ({x_1} + {x_2} + {x_3} + {x_4})/4} $

(8)

$ {{y_5} = ({y_1} + {y_2} + {y_3} + {y_4})/4} $

(9)

$ {{h_5} = ({h_1} + {h_2} + {h_3} + {h_4})/4} $

(10)

（2）判断三角形网格内等值点的存在性。将三角形①各顶点的高程排序，记作h₁₁、h₁₂、h₁₃，下标中的数字第1个为三角形序号，第2个表示高程大小（h₁₁≥h₁₂≥h₁₃），各自坐标为(x₁，y₁)、(x₂，y₂)、(x₃，y₃)。对比高程和目标值h_i（i = 1，…，n），若h₁₁≤h_i≤h₁₃，则该网格中存在等值点；否则，移动至下一网格，重复步骤（2）。若4个三角形网格均不存在等值点，执行步骤（5）。

（3）连接等值点。假定三角形①中存在等值点，先对边h₁₁h₁₃进行插值并储存坐标，公式如下

$ {x = {x_3} + \frac{{{x_1} - {x_3}}}{{{h_{11}} - {h_{13}}}}({h_{11}} - {h_i})} $

(11)

$ {y = {y_3} + \frac{{{y_1} - {y_3}}}{{{y_{11}} - {y_{13}}}}({h_{11}} - {h_i})} $

(12)

其中，h_i是该网格内的等值点。将h_i和h₁₂进行比较，若h_i≥h₁₂，则在边h₁₁h₁₂进行插值并储存坐标；否则，在边h₁₂h₁₃进行插值并储存坐标。连接对应各个h_i (i =1，…，n)的等值点，获得在该三角形网格内的等值线。

（4）重复步骤（2）、（3），直至该矩形网格内的等值线绘制完成。

（5）移动矩形网格，使之不重复地覆盖给定区域。重复步骤（2）、（3）、（4），获取等值线图。

3 实验验证与磁图分析 3.1 实验验证

（1）验证本软件解算的地磁模型的正确性。由于某极轨卫星的磁强计只能测得所在位置的磁场总强度。故以IGRF模型计算了2018年7月该卫星各位置的磁场强度，和实测值比较，结果见表 1，其中西经、北纬为正。由表 1可见，本软件获得的地磁场幅值的误差小于200 nT，与IGRF模型误差的数量级相当，验证了该程序的正确性。

（2）验证本软件等值线生成功能的正确性。Thébault等（2015）公开了IGRF-12的高斯系数和构建该模型所使用的观测数据来源。此外，如图 4(a)所示，该论文给出地心距为R_e的地磁场幅值T的地磁基准图供导航使用，需要注意的是其纬度并非线性变化。设定等值线数量为8，地心距为R_e，生成2015年地磁场幅值T的地磁基准图，见图 4(b)。

由图 4可见，二者等值线多数与纬线平行，并存在4个闭合等值线圈；二者在30°S处幅值最小，之后向南北两侧增大，至60°N（S）处达到峰值。若忽略因等值线数量不同导致的差异，二图无明显区别。因此，验证了利用本软件可正确生成地磁基准图。

3.2 磁图分析

等值线数量n = 8，地心距为R_e，生成2015年地球主磁场在地理坐标系北东地方向下的磁场分量（B_x，B_y，B_z）地磁基准图（B_x、B_y、B_z单位为nT），见图 5、图 6、图 7。由等值线分布可知，B_x和B_z分量等值线多与纬线平行，而B_y分量等值线多呈闭合状态。由各分量幅值变化趋势可知，B_x分量幅值在赤道附近最大，达35 000 nT，向南北减小，最小达-10 000 nT；B_y分量幅值在东西半球呈不同变化趋势，在东半球上，由南极处近-14 000 nT向北增加至近1 300 nT，在西半球上，由13 000 nT向北减小至1 500 nT；B_z分量幅值由南向北迅速变大，最小值不足-50 000 nT，最大值可达46 000 nT。结合图 4(a)和图 5—图 7可知，地磁场幅值T在低纬度地区受B_x分量影响大，随着纬度增加，B_z分量权重逐渐上升，因此，其幅值在高纬度地区取峰值，与B_x分量呈相反的变化规律。

4 结束语

研究基于IGRF模型的地磁要素计算方法，提出具有虚节点的移动矩形网格等值线绘制方法。采用C#，实现给定待测区域地理信息及自动绘制该区域地磁基准图的功能。经实验，该软件可正确解算IGRF模型，获得目标区域的地磁要素值；可设定等值线数量，生成的地磁基准图不存在等值线交叉等错误，成图效果好。利用本软件可生成地球主磁场在任意区域任意高度的地磁基准图，满足对定位精度要求不高的飞行器的使用要求。目前，本软件对存在空白或分布不均的数据成图尚无法实现，且对实测数据的追踪能力不佳，仍需进一步完善。