利用IGRF模型计算全张量地磁梯度

钟炀，管彦武，石甲强，肖锋

(吉林大学 地球探测科学与技术学院，吉林 长春 130026)

0 引言

国内外诸多学者针对IGRF模型的研究，早期集中在分析模型精度及验证模型适用范围等问题上。Barraclough等[1]利用大量地磁台站数据推导出地磁场长期变化的球谐模型，从模型中导出了磁极、地磁极和偏心磁极的位置及其变化速率，并比较了该模型与IGRF模型的差别，指出在模型中加入加速变化系数对地磁场模型的拟合度提高具有显著的作用。安振昌[2-9]首次在国内较为详细地介绍了国际地磁参考场、发展历史及其主要内容，并利用第七代IGRF模型计算1900～2000年非偶极子磁场的全球变化和1995年中国及邻区地磁场。林云芳等[10]利用东南亚及其邻近地区的地磁台的实测资料与IGRF模式进行了对比，证明IGRF模型存在偏差，并指出应用时应该选取时间间隔短的IRGF模型。Matteo等[11]通过分析1991～2010年间电离层磁场强度的卫星数据，与IGRF模型计算的结果进行比较，结果表明即使在地磁场活跃的条件下，在一定的位置和高度范围内，IGRF在电离层仍然具有很高的精度。

为适用于高精度地质构造研究工作，徐文耀[12]指出在1945～1955年IGRF模型中的高阶高斯系数(n>7)与其他时期的模型相比，存在较大偏差，并通过自然正交分量分析法对1945～1955年IGRF模型中的高次高斯系数进行修正，使这些模型的高斯系数具有较好的平滑时变特性。Maus等[13]美国地球物理数据中心(NGDC)的科学家们计算出16～720阶球谐系数，成功建立岩石圈磁场模型NGDC-720。张素琴等[14]利用中国部分地磁台站1990～2003 年的年均值资料， 研究了地磁台站年均值与IGRF模型计算值的一致程度，结果显示虽然存在差异，但是各地磁要素年均值与IGRF模型值的差值的标准差均低于IGRF模型的误差水平，故一致性较好。

在IGRF模型的应用方面， Molina(1987)等[15]针对意大利等欧洲中部国家，建立了基于IGRF模型的意大利地磁参考场(ITGRF)和欧洲中部地磁参考场(EGRF)。任国泰[16]利用IGRF模型结合地面测量资料和地磁台资料分析了东亚大陆磁场。安振昌[17]比较和讨论了中国地区地磁场区域模型与全球模型，得出地磁场区域模型，更详细准确地表示我国地磁场的分布。Smart等[18-19]运用磁场模型与运动轨迹学原理，利用IGRF模型导出飞机飞行时截止刚度(VCR)，进而通过线性拟合确定其飞行轨迹。Davis[20]利用Matlab编程实现了IGRF模型关于场源外任意位置三分量地磁场的计算。

此外，冯春[21]、柴松均[22]及杨梦雨[23]等人，分别运用MATLAB、C及Java语言实现了任意给定点位的地磁七要素计算。但随着航空全张量磁梯度测量技术的发展，对利用IGRF模型计算地球主磁场的全张量磁梯度有着迫切的需求，故笔者梳理了地磁场七要素及张量的计算流程，推导出计算任意给定点位全张量地磁梯度的公式，绘制出某地区地磁场全张量等值线图，并验证了算法的正确性。

1 利用IGRF模型计算全张量地磁梯度

目前，球谐函数法被广泛应用于计算地磁场参考模型中。在该方法中，认为地球是均匀球体，在球坐标系下地磁场源区外任意一点坐标(r,θ,φ)磁位V的拉普拉斯方程形式如式(1)所示：

(1)

其中：r、θ和φ是地心球坐标系构成要素，r是从地心起算的径向距离，单位为km，θ是地心余纬度(θ=90°-L，L是地心纬度)，φ是从格林尼治子午线算起的地心经度。

利用分离变量法求解式(1)的Dirichlet问题(常微分方程第一类边界条件)，获得地磁场源区外任意一点坐标(r，θ，φ)的磁位V表达式[24]如式(2)所示：

(2)

(3)

(4)

在物理学中，磁场强度定义为磁位V的负梯度，在地心球坐标系中，对式(4)分别求r、θ和φ的偏导数，并取负得到磁场强度的三个分量Ur、Uθ、Uφ：

(5)

(6)

(7)

为计算全张量磁梯度，保持原坐标系不变，在式(5)～(7)中，分别对Ur、Uθ和Uφ求r、θ、φ的偏导数，得到磁位的二阶偏导数：

(8)

(9)

(10)

(11)

(12)

(13)

Balmino等对地心球坐标系下磁位的一阶及二阶偏导数(即式(5)～(13))进行坐标变换，将地心球坐标系转换为局部“北东上”直角坐标系，其变换方式如式(14)、(15)所示[25]：

(14)

(15)

其中：Bx、By、Bz分别表示局部“北东上”直角坐标系下地磁场三分量，Bxx、Bxy、Bxz、Byy、Byz、Bzz表示局部“北东上”直角坐标系下地磁场的全张量的6个分量。

然后将局部“北东上”直角坐标系xyz转换到符合右手定则且更加常用的局部“北东下”直角坐标系x1y1z1，转换关系为：

(16)

再将z1指向地心的x1y1z1坐标系转换到z2沿旋转椭球面法线向下的局部直角坐标系x2y2z2，即相当于x1y1z1坐标系绕y1旋转角度d。d为地心余纬度θ与地理余纬度θ′的差值。图1中P点表示测点。

图1 地理坐标系与地心坐标系Fig.1 Geographic coordinate system and geocentric coordinate system

由图1可知：

d=θ-θ′，

(17)

(18)

(19)

(20)

b=a·(1-f)，

(21)

其中：a为椭球长半轴，a=6 378.137 km，b为椭球短半轴，f为地球扁率，f=1/298.257 223 563，h为测点的大地高(椭球高度)。

由式(17)可知：

sinθ=sinθ′cosd+cosθ′sind，

(22)

cosθ=cosθ′cosd-sinθ′sind，

(23)

(24)

在实际应用中，通常测量得到测点的大地坐标，通过由式(22)～(24)可求得地心余纬度的三角函数值，并将其应用于前述公式。

再由图1所示坐标系x1y1z1与x2y2z2的旋转关系可知：

(25)

(26)

(27)

利用式(26)及(27)进行坐标变换，获得更普遍的z轴沿旋转椭球面法线向下的局部直角坐标系(北东下)的地磁场三分量及全张量公式。

2 伴随勒让德多项式及其递推关系

利用IGRF模型，为计算地磁场三分量及全张量磁梯度，还需要计算满足正交性的施密特规格化伴随勒让德多项式。

数学上，勒让德函数是指以下勒让德方程的解：

(28)

勒让德微分方程(28)的解，可以写成标准的幂级数形式，并且当方程满足|v|<1且n为非负整数时，方程的解将随n值的变化而变化并构成一组由正交多项式组成的多项式序列，即勒让德多项式[26]。勒让德多项式的微分形式，可用罗德里格斯公式表示：

(29)

其中：0～3阶勒让德多项式表达式如表1所示。

表1 0～3阶勒让德多项式Table 1 Legendre polynomials of order 0～3

在球坐标系下求解拉普拉斯方程，可得到如下常微分方程：

(30)

常微分方程式(32)的解序列称为伴随勒让德多项式，又称为缔合勒让德多项式、连带勒让德多项式或关联勒让德多项式。伴随勒让德多项式也可以通过对勒让德多项式求m次导数，经式(31)换算得到：

(31)

如式(4)所示，磁位的表达式是球谐函数的k阶展开，可表示为与展开阶数k有关的多项式之和。为使这些多项式在求和计算中的相对重要性更加接近，并方便数值计算，需要将伴随勒让德多项式与一个随阶次变化的因子相乘，即对伴随勒让德多项式规格化。伴随勒让德多项式的规格化主要有高斯规格化和施密特规格化两种方式[27]。采用高斯规格化伴随勒让德多项式的表达式如下：

(32)

其中：Pn,m(cosθ)为高斯规格化伴随勒让德多项式，Pn,m(cosθ)为伴随勒让德多项式。

利用施密特方程计算的规格化勒让德多项式表达式为：

(33)

两种规格化方法的关系见式(34)所示：

(34)

(35)

0～3阶规格化伴随勒让德多项式如表2所示。

表2 0～3阶伴随勒让德多项式Table 2 Order 0-3 adjoint Legendre polynomials

为了方便计算机进行编程计算，需根据施密特半规则化勒让德多项式(33)求取伴随勒让德多项式的递推关系。采用向前列递推法，施密特拟规格化伴随勒让德多项式的递推关系为：

(36)

地磁场三分量及全张量磁梯度表达式(26)、(27)中，还含有施密特拟规格化伴随勒让德多项式的微分形式。其递推关系如式(37)所示。

在全张量磁梯度表达式中的Bxx和Bzz分量，除含有勒让德多项式的一阶微分形式外，还含有二阶微分形式，其递推关系如式(38)所示。

3 实际地磁场模拟

实际应用中，需要选取地球主磁场梯度较小的区域和高度作为航磁或FTMG测量学习飞行的工区。同时，为验证上述算法在计算地磁场七要素及全张量地磁梯度的可行性，选取4°×4°的计算工区(经度变化范围为：103.305 6°～107.305 6°，纬度变化范围为：27.305 6°～31.305 6°)，网格大小为0.1°×0.1°。选取大地高1 km，选定时间为2019年4月7日进行计算。

(37)

(38)

图2为所选区域的总地磁场等值线图，图中等值线光滑，场值变化范围为48 693～51 200 nT。

图2 总地磁场B等值线Fig.2 Contour map of the total geomagnetic field

图3a为地磁偏角D等值线图，其变化范围-2.871 7°～-1.652 5°，图3b为地磁倾角I等值线图，其变化范围为42.748°～49.133 6°。

图4从左至右依次为地磁场在x、y、z轴的分量Bx、By、Bz。 由于目前还没有其他直接获得任意给定点位地球主磁场全张量磁梯度数据的方法，无法逐点对比图5中所选测区内全张量磁梯度场值。但根据地磁场三分量等值线图(图4)，其沿3个坐标轴方向的导数即为全张量地磁梯度值，可以大致判断6个张量分量符合地磁场梯度量级，且等值线形态合理。

图3 地磁场偏角D(a)及倾角I(b)等值线Fig.3 Contour map of geomagnetic field deviation(a) and dip angle(b)

为进一步验证图5中全张量地磁梯度计算值的正确性，将Bxx、Byy、Bzz相加如图6所示。

图4 地磁场三分量(a)Bx(b)By(c)Bz等值线Fig.4 Three-component contour map of the geomagnetic field

图5 全张量地磁梯度等值线Fig.5 Full-tensor geomagnetic gradient contour map

如图6所示，全张量磁梯度的3个分量Bxx、Byy、Bzz之和的变化范围为：-0.001 1～0.001 1 nT/km，可以认为利用该计算方法所获得全张量地磁梯度满足Laplace方程。

为进一步验证算法的正确性，选取计算工区内4个随机点位：①四川省成都市(30.67N，104.07E)，②自贡市(29.35N，104.78E)，③泸州市(28.87N，105.43E)，④德阳市(31.13N，104.38E))，利用美国国家海洋和大气管理局(NOAA)官方网站(www.ngdc.noaa.gov)内全球地磁场数据，对比两种方式计算获得地磁场七要素误差见表3。三分量矢量场在保留小数点一位的情况下误差很小，可忽略不计，从而验证了算法的正确性。

表3 测区内点位磁场值Table 3 Table of point position magnetic field values in the test area

图6 主对角线上3个张量分量之和等值线Fig.6 Contour map of the sum of the three tensor components on the main diagonal

4 结论

1)笔者梳理了利用IGRF模型计算地磁场的原理与步骤，利用4个点位的地磁七要素计算结果与美国国家海洋和大气管理局计算数据进行对比，结果准确可靠。

2)依据IGRF模型，推导出全张量地磁梯度的球谐展开形式，给出了施密特规格化伴随勒让德多项式二阶导数的递推公式。实验工区全张量地磁梯度的计算结果形态与量级合理，且满足Laplace方程，为航空全张量磁梯度测量中选择学习飞行工区和飞行高度提供了重要的依据。