文章快速检索    
  地震地磁观测与研究  2020, Vol. 41 Issue (3): 98-104  DOI: 10.3969/j.issn.1003-3246.2020.03.013
0

引用本文  

陈贵芳, 郁丰, 王润. 基于IGRF的地磁基准图技术[J]. 地震地磁观测与研究, 2020, 41(3): 98-104. DOI: 10.3969/j.issn.1003-3246.2020.03.013.
CHEN Guifang, YU Feng, WANG Run. Geomagnetic reference map technology based on IGRF[J]. Seismological and Geomagnetic Observation and Research, 2020, 41(3): 98-104. DOI: 10.3969/j.issn.1003-3246.2020.03.013.

基金项目

国家自然科学基金(项目编号:61673212)

作者简介

陈贵芳(1995-), 男, 浙江湖州人, 硕士研究生, 主要研究方向:地磁导航。E-mail:cgf56579@163.com

文章历史

本文收到日期:2019-06-11
基于IGRF的地磁基准图技术
陈贵芳 1,2, 郁丰 1,2, 王润 1,2     
1. 中国南京 210016 空间光电探测与感知工信部重点实验室;
2. 中国江苏 210016 南京航空航天大学航天学院
摘要:地磁导航需要获得目标区域的地磁基准图。研究了基于国际地磁参考场的地磁场要素计算和等值线绘制方法。以给定目标区域的经纬高为输入,计算高阶地磁球谐函数获取地磁矢量并计算地磁场要素,修正了因地球扁率带来的数据差异。采用具有虚节点的移动矩形网格算法来绘制地磁参数等值线,生成目标区域的地磁基准图。采用C#实现了地磁基准图软件,测试表明,获得的地磁等高线数据满足使用要求,等值线精度较高。
关键词地磁导航    IGRF    地磁场    等值线    虚节点    
Geomagnetic reference map technology based on IGRF
CHEN Guifang 1,2, YU Feng 1,2, WANG Run 1,2     
1. Key Laboratory of Space Photoelectric Detection and Perception of Ministry of Industry and Information Technology, Nanjing 210016, China;
2. Academy of Astronautics, Nanjing University of Aeronautics and Astronautics, Jiangsu Province 210016, China
Abstract: It is necessary to obtain a geomagnetic reference map of the target area to perform geomagnetic navigation. Based on the international geomagnetic reference field (IGRF), this paper studied the methods to calculate the geomagnetic elements and draw the geomagnetic elements contour maps. The geomagnetic data is obtained for given geographical parameters by solving the IGRF. Then use this data to draw the geomagnetic elements contour map by using the method of moving quadrilateral within virtual nodes. The program is prepared by C# and the experimental results show that this program can calculate at any given point and draw its contour map.
Key words: geomagnetic navigation    IGRF    geomagnetic field    contour map    virtual nodes    
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,2012Albahari,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)

其中,gnmhnm为高斯系数;λ为格林尼治天文台起算的东经,范围为[0,360°);θ为余纬,定义为90°- LL为纬度,θ的范围为[0,180°];r为地心距(坐标轴定义见图 1);Pnmnm阶的勒让德多项式;采用WGS-84参考椭球模型,赤道平面半径Re= 6 378 137 m,1/f = 298.257。

图 1 坐标系统定义 Fig.1 The definition of the coordinate system

图 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模型,最后修正误差。坐标转换具体描述如下。

假设待测点的地理坐标为(λLh),其在地固坐标系中的坐标为

$ \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)

其中,RN = Re(1 + f sin2 L),是当地的卯酉圈半径。

利用所在地的直角坐标(xyz)计算经度、余纬和地心距,得

$ \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′),可与实际的北东地坐标系重合。北半球旋转公式如下

图 2 地心坐标和地理坐标系转换关系 Fig.2 The relationship between the geographic coordinate system and the geocentric coordinate system
$ \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个三角形的网格,移动矩形网格并绘制等值线,直至覆盖整个区域,即获得完整的等值线图。

在矩形网格内绘制等值线前,需计算各等值点的高程hi,取高程的最大值hM和最小值hm,则可得

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

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

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

图 3 有虚节点的矩形网格及其三角形网格 (a)矩形网格;(b)三角形网格 Fig.3 A rectangular grid and a triangle grid with a virtual node
$ {{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)判断三角形网格内等值点的存在性。将三角形①各顶点的高程排序,记作h11h12h13,下标中的数字第1个为三角形序号,第2个表示高程大小(h11h12h13),各自坐标为(x1y1)、(x2y2)、(x3y3)。对比高程和目标值hii = 1,…,n),若h11hih13,则该网格中存在等值点;否则,移动至下一网格,重复步骤(2)。若4个三角形网格均不存在等值点,执行步骤(5)。

(3)连接等值点。假定三角形①中存在等值点,先对边h11h13进行插值并储存坐标,公式如下

$ {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)

其中,hi是该网格内的等值点。将hih12进行比较,若hih12,则在边h11h12进行插值并储存坐标;否则,在边h12h13进行插值并储存坐标。连接对应各个hi (i =1,…,n)的等值点,获得在该三角形网格内的等值线。

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

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

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

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

表 1 地磁数据标准值和计算值统计 Table 1 The calculated and the actual geomagnetic field magnitude

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

图 4 2015年地磁总强度的地磁基准图 (a)软件生成的地磁基准图;(b)IGRF-12提供的地磁基准图 Fig.4 Contour map of T in 2015

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

3.2 磁图分析

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

图 5 Bx地磁基准 Fig.5 Contour map of Bx
图 6 By地磁基准 Fig.6 Contour map of By
图 7 Bz地磁基准 Fig.7 Contour map of Bz
4 结束语

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

参考文献
侯凤贞, 庄建军, 宁新宝. 2008. 基于构型表的等值线绘制算法及程序实现[J]. 南京大学学报(自然科学), 44(4): 371-378.
廖国忠, 张伟, 梁生贤, 吴文贤, 李富. 2014. 基于规则三角网的等值线追踪与填充算法的实现和应用[J]. 物探化探计算技术, 36(1): 120-123.
苗润忠. 2004. 光滑的等值线生成算法[J]. 长春理工大学学报, 27(1): 16-18.
孙桂茹, 马亮, 路登平, 赵国瑞, 郝嘉林. 2000. 等值线生成与图形填充算法[J]. 天津大学学报, 33(6): 816-818.
姚家振.基于Delaunay三角网的等值线生成算法的研究[D].安徽大学, 2015.
郑盛贵, 颜七笙, 黄临平. 2004. 基于点的三角形构网算法及等值线自动生成方法[J]. 计算机与现代化, 11061106(5): 7-9.
周能兵, 王亚斌, 王强. 2018. 地磁导航技术研究进展综述[J]. 导航定位学报, 6(2): 15-19.
Albahari J, Albahari B. 2017. C# 7.0 in a nutshell:the definitive reference[M]. O'Reilly Media, Inc.
Campbell W H. 2003. Introduction to geomagnetic fields[M]. Cambridge University Press.
Griffiths I. 2012. Programming C# 5.0:Building Windows 8, Web, and Desktop Applications for the. NET 4.5 Framework[M]. O'Reilly Media, Inc.
Thébault E, Finlay C C, Beggan C D, et al. 2015. International geomagnetic reference field:the 12th generation[J]. Earth, Planets and Space, 67: 79.
Wan M, Lim C, Zhang J, et al. Reconstructing patient-specific cardiac models from contours via delaunay triangulation and graph-cuts[C]//2013 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE, 2013: 2 976-2 979.