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

doi:10.11720/wtyht.2020.1416

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

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

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

The calculation method of full tensor geomagnetic gradient based on IGRF model

ZHONG Yang, GUAN Yan-Wu, SHI Jia-Qiang, XIAO Feng

College of Geo-Exploration Science and Technology,Jilin University,Changchun 130026,China

责任编辑: 王萌

收稿日期: 2019-09-2 修回日期: 2019-10-14 网络出版日期: 2020-06-20

基金资助:

国家重点研发计划“航空磁场测量技术系统研制”. 2017YFC0602000

Received: 2019-09-2 Revised: 2019-10-14 Online: 2020-06-20

作者简介 About authors

钟炀(1995-),男,吉林长春人,硕士研究生,主要研究方向为地球探测与信息技术。Email:zhongyang18@mails.jlu.edu.cn 。

摘要

国际地磁参考场(IGRF)是一种描述地球主磁场的国际通用模型。目前利用该模型能够计算任意点位的地磁七要素,但随着航空全张量磁测技术的发展,对地球主磁场的全张量磁梯度数据有着迫切的需求。笔者梳理了IGRF模型的计算原理并进一步推导出球谐展开的全张量磁梯度表达式,实现了任意给定点位地磁场七要素和全张量的计算,并用地磁七要素与美国国家海洋和大气管理局(NOAA)网站计算数据进行对照,结果准确可靠。绘制了某地区地磁场的全张量磁梯度等值线图,结果满足Laplace方程,为航空全张量磁梯度测量中选择学习飞行工区和飞行高度提供了理论依据。

关键词： IGRF ; 全张量磁梯度 ; 球谐分析 ; 勒让德多项式

Abstract

The international geomagnetic reference field (IGRF) is a general international model for describing the earth’s main magnetic field. At present,this model can be used to calculate the seven elements of geomagnetic field at any point. However,with the development of aeronautical full tensor magnetic measurement technology,there is an urgent need for full tensor geomagnetic gradient data. In this paper,the calculation principle of the IGRF model is summarized and the expression of the full tensor geomagnetic gradient with spherical harmonic expansion is derived. The calculation of the seven elements of geomagnetic field and the full tensor geomagnetic gradient at any given point is realized. Comparing with the calculated data from the website of the National Oceanic and Atmospheric Administration of the United States (NOAA),the results are accurate and reliable. The contour map of the full tensor geomagnetic field in a region is drawn, and the results were in accordance with the Laplace equation. It provides the theoretical basis for the selection of learning flight working area and flight height in the aeromagnetic survey.

Keywords： IGRF ; FTMG ; spheric harmonic analysis ; legendre polynomials

PDF (3448KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

钟炀, 管彦武, 石甲强, 肖锋. 利用IGRF模型计算全张量地磁梯度. 物探与化探[J], 2020, 44(3): 582-590 doi:10.11720/wtyht.2020.1416

ZHONG Yang, GUAN Yan-Wu, SHI Jia-Qiang, XIAO Feng. The calculation method of full tensor geomagnetic gradient based on IGRF model. Geophysical and Geochemical Exploration[J], 2020, 44(3): 582-590 doi:10.11720/wtyht.2020.1416

0 引言

为了描述地磁场时空非线性变化特征,自1965年起国际地磁与高空大气物理学协会(IAGA)根据地磁场的球谐分析模式,建立了全球磁场模型,并每隔5年给出一个国际地磁参考场(IGRF)模型。IGRF的实质是一套球谐模型的高斯系数,即 $g_{n}^{m}$ 和 $h_{n}^{m}$ 。迄今为止,利用这些系数可以计算1900~2020年之间的地磁场矢量。总览与国际地磁参考场相关研究的发展历程,国外学者起步较早,将IGRF模型综合应用于航空磁测、地质构造分析等许多领域,并建立起多个高精度局部地磁场模型以满足地质勘探及其他相关科研工作的需要。而国内相关研究工作主要集中在分析并提高IGRF模型的精度,研究我国及周边国家利用该模型的可行性等方面。

国内外诸多学者针对IGRF模型的研究,早期集中在分析模型精度及验证模型适用范围等问题上。Barraclough等^[1]利用大量地磁台站数据推导出地磁场长期变化的球谐模型,从模型中导出了磁极、地磁极和偏心磁极的位置及其变化速率,并比较了该模型与IGRF模型的差别,指出在模型中加入加速变化系数对地磁场模型的拟合度提高具有显著的作用。安振昌^{[2,3,4,5,6,7,8,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)

\frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial V}{\partial r}) + \frac{1}{r^{2} \sin θ} \frac{\partial}{\partial θ} (\sin θ \frac{\partial V}{\partial θ}) + \frac{1}{r^{2} si n^{2} θ} \frac{\partial^{2} V}{\partial ϕ^{2}} = 0,

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

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

(2)

\begin{array}{l} V (r, θ, ϕ) = \overset{\infty}{\sum_{n = 0}} (A r^{n} + \frac{B}{r^{n + 1}}) \cdot \\ \overset{\infty}{\sum_{m = 0}} (a_{nm} \cos mϕ + b_{nm} \sin mϕ) P_{n}^{m} (\cos θ), \end{array}

其中:A、B是与球谐函数边界条件有关的系数;a_nm、b_nm为与球谐函数展开阶数有关的系数; $P_{n}^{m}$ (cosθ)为n阶m次的规格化缔和勒让德多项式。由式(2)可见,磁位可分为两部分,一部分含 $\frac{1}{r^{n + 1}}$ 因子,另一部分含rⁿ因子。根据地磁场高斯理论,按照场源位置,地磁场分为内源场和外源场(内源场即磁场源位于地球内部,外源场磁场源位于地球外部)。地磁场可进一步表示为^[24]:

(3)

\begin{array}{r} V & = V_{i} + V_{e} \\ = a \overset{k}{\sum_{n = 1}} \overset{n}{\sum_{m = 0}} [{(\frac{a}{r})}^{n + 1} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) P_{n}^{m} (\cos θ) + \\ {(\frac{r}{a})}^{n} (j_{n}^{m} \cos mϕ + k_{n}^{m} \sin mϕ) P_{n}^{m} (\cos θ)], \end{array}

其中:V_i、V_e分别代表地磁场中内源场和外源场的磁位,a是参考球的半径(a=6 371.2 km);k代表球谐级数展开的最大阶数, $g_{n}^{m}$ 和 $h_{n}^{m}$ 是高斯球谐系数,由国际地磁学与高空物理学会(IAGA)每5年发布一次; $j_{n}^{m}$ 、 $k_{n}^{m}$ 为外源场的高斯系数。地球的基本磁场是指地磁场中稳定缓慢变化的部分,高斯分析得出该部分磁场强度占总场强度的99%以上,来源于地球内部,即内源场。只考虑地磁场中的内源场,则磁位可表示为:

(4)

\begin{array}{l} V (r, θ, ϕ) = a \overset{k}{\sum_{n = 1}} {(\frac{a}{r})}^{n + 1} \overset{n}{\sum_{m = 0}} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) \cdot \\ P_{n}^{m} (\cos θ) 。 \end{array}

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

(5)

\begin{array}{l} \begin{array}{r} U_{r} & = - \frac{\partial V}{\partial r} \\ = \overset{k}{\sum_{n = 1}} (n + 1) {(\frac{a}{r})}^{n + 2} \overset{n}{\sum_{m = 0}} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) \cdot \\ P_{n}^{m} (\cos θ), \end{array} \end{array}

(6)

\begin{array}{r} U_{θ} & = - \frac{\partial V}{\partial θ} \\ = - a \overset{k}{\sum_{n = 1}} {(\frac{a}{r})}^{n + 1} \overset{n}{\sum_{m = 0}} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) \cdot \\ \frac{\partial P_{n}^{m} (\cos θ)}{\partial θ}, \end{array}

(7)

\begin{array}{r} U_{ϕ} & = - \frac{\partial V}{\partial ϕ} \\ = - a \overset{k}{\sum_{n = 1}} {(\frac{a}{r})}^{n + 1} \overset{n}{\sum_{m = 0}} m (- g_{n}^{m} \sin mϕ + h_{n}^{m} \cos mϕ) \cdot \\ P_{n}^{m} (\cos θ) 。 \end{array}

为计算全张量磁梯度,保持原坐标系不变,在式(5)~(7)中,分别对U_r、U_θ和U_ϕ求r、θ、ϕ的偏导数,得到磁位的二阶偏导数:

(8)

\begin{array}{l} \begin{array}{r} U_{θθ} & = - \frac{\partial^{2} V}{\partial θ^{2}} \\ = - a \overset{k}{\sum_{n = 1}} {(\frac{a}{r})}^{n + 1} \overset{n}{\sum_{m = 0}} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) \cdot \\ \frac{\partial^{2} P_{n}^{m} (\cos θ)}{\partial θ^{2}}, \end{array} \end{array}

(9)

\begin{array}{r} U_{ϕθ} & = - \frac{\partial^{2} V}{\partial θ \partial φ} \\ = - a \overset{k}{\sum_{n = 1}} {(\frac{a}{r})}^{n + 1} \overset{n}{\sum_{m = 0}} m (- g_{n}^{m} \sin mϕ + h_{n}^{m} \cos mϕ) \cdot \\ \frac{\partial P_{n}^{m} (\cos θ)}{\partial θ}, \end{array}

(10)

\begin{array}{r} U_{rθ} & = - \frac{\partial^{2} V}{\partial r \partial θ} \\ = \overset{k}{\sum_{n = 1}} (n + 1) {(\frac{a}{r})}^{n + 2} \overset{n}{\sum_{m = 0}} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) \cdot \\ \frac{\partial P_{n}^{m} (\cos θ)}{\partial θ}, \end{array}

(11)

\begin{array}{r} U_{ϕϕ} & = - \frac{\partial^{2} V}{\partial φ^{2}} \\ = a \overset{k}{\sum_{n = 1}} {(\frac{a}{r})}^{n + 1} \overset{n}{\sum_{m = 0}} m^{2} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) \cdot \\ P_{n}^{m} (\cos θ), \end{array}

(12)

\begin{array}{r} U_{rϕ} & = - \frac{\partial^{2} V}{\partial r \partial φ} \\ = \overset{k}{\sum_{n = 1}} (n + 1) {(\frac{a}{r})}^{n + 2} \cdot \\ \overset{n}{\sum_{m = 0}} m (- g_{n}^{m} \sin mϕ + h_{n}^{m} \cos mϕ) P_{n}^{m} (\cos θ), \end{array}

(13)

\begin{array}{r} U_{rr} & = - \frac{\partial^{2} V}{\partial r^{2}} \\ = - \frac{1}{r} \overset{k}{\sum_{n = 1}} (n + 2) (n + 1) {(\frac{a}{r})}^{n + 2} \cdot \\ \overset{n}{\sum_{m = 0}} (g_{n}^{m} \cos mϕ + h_{n}^{m} \sin mϕ) P_{n}^{m} (\cos θ) 。 \end{array}

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

(14)

\begin{array}{l} \{\begin{array}{l} B_{x} = - \frac{1}{r} U_{θ} \\ B_{y} = \frac{1}{r \sin θ} U_{ϕ} \\ B_{z} = U_{r} \end{array} \end{array}

(15)

\{\begin{array}{l} B_{xx} = \frac{1}{r} U_{r} + \frac{1}{r^{2}} U_{θθ} \\ B_{xy} = \frac{\cos θ}{r^{2} si n^{2} θ} U_{ϕ} - \frac{1}{r^{2} \sin θ} U_{ϕθ} \\ B_{xz} = \frac{1}{r^{2}} U_{θ} ⁃ \frac{1}{r} U_{rθ} \\ B_{yy} = \frac{1}{r} U_{r} + \frac{1}{r^{2} \tan θ} U_{θ} + \frac{1}{r^{2} si n^{2} θ} U_{ϕϕ} \\ B_{yz} = \frac{1}{r \sin θ} U_{rϕ} - \frac{1}{r^{2} \sin θ} U_{ϕ} \\ B_{zz} = U_{rr} \end{array}

其中:B_x、B_y、B_z分别表示局部“北东上”直角坐标系下地磁场三分量,B_xx、B_xy、B_xz、B_yy、B_yz、B_zz表示局部“北东上”直角坐标系下地磁场的全张量的6个分量。

然后将局部“北东上”直角坐标系xyz转换到符合右手定则且更加常用的局部“北东下”直角坐标系x₁y₁z₁,转换关系为:

(16)

\begin{array}{l} [\begin{array}{l} x_{1} \\ y_{1} \\ z_{1} \end{array}] = [\begin{array}{l} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & - 1 \end{array}] [\begin{array}{l} x_{1} \\ y_{1} \\ z_{1} \end{array}] \\ [\begin{array}{l} B_{x 1} \\ B_{y 1} \\ B_{z 1} \end{array}] = [\begin{array}{l} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & - 1 \end{array}] [\begin{array}{l} B_{x} \\ B_{y} \\ B_{z} \end{array}] \\ [\begin{array}{l} B_{xx 1} & B_{xy 1} & B_{xz 1} \\ B_{yx 1} & B_{yy 1} & B_{yz 1} \\ B_{zx 1} & B_{zy 1} & B_{zz 1} \end{array}] = [\begin{array}{l} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & - 1 \end{array}] [\begin{array}{l} B_{xx} & B_{xy} & B_{xz} \\ B_{yx} & B_{yy} & B_{yz} \\ B_{zx} & B_{zy} & B_{zz} \end{array}] \cdot \\ [\begin{array}{l} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & - 1 \end{array}] \end{array}

再将z₁指向地心的x₁y₁z₁坐标系转换到z₂沿旋转椭球面法线向下的局部直角坐标系x₂y₂z₂,即相当于x₁y₁z₁坐标系绕y₁旋转角度d。d为地心余纬度θ与地理余纬度θ'的差值。图1中P点表示测点。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 地理坐标系与地心坐标系

Fig.1 Geographic coordinate system and geocentric coordinate system

由图1可知:

(17)

\begin{array}{l} d = θ - θ', \end{array}

(18)

r = \{h [h + 2 (a^{2} si n^{2} θ + b^{2} co s^{2} {θ)}^{\frac{1}{2}}] \frac{}{} + \frac{a^{4} si n^{2} θ + b^{4} co s^{2} θ}{a^{2} si n^{2} θ + b^{2} co s^{2} θ}\},

(19)

\cos d = \frac{h + \sqrt[]{a^{2} si n^{2} θ + b^{2} co s^{2} θ}}{r},

(20)

\sin d = \frac{(a^{2} - b^{2}) \sin θ \cos θ}{r \sqrt[]{a^{2} si n^{2} θ + b^{2} co s^{2} θ}},

(21)

b = a \cdot (1 - f),

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

由式(17)可知:

(22)

\begin{array}{l} \sin θ = \sin θ' \cos d + \cos θ' \sin d, \end{array}

(23)

\cos θ = \cos θ' \cos d - \sin θ' \sin d,

(24)

\tan θ = \frac{\sin θ}{\cos θ} 。

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

再由图1所示坐标系x₁y₁z₁与x₂y₂z₂的旋转关系可知:

(25)

[\begin{array}{l} x_{2} \\ y_{2} \\ z_{2} \end{array}] = [\begin{array}{l} \cos d & 0 & \sin d \\ 0 & 1 & 0 \\ - \sin d & 0 & \cos d \end{array}] [\begin{array}{l} x_{1} \\ y_{1} \\ z_{1} \end{array}],

(26)

[\begin{array}{l} B_{x 2} \\ B_{y 2} \\ B_{z 2} \end{array}] = [\begin{array}{l} \cos d & 0 & \sin d \\ 0 & 1 & 0 \\ - \sin d & 0 & \cos d \end{array}] [\begin{array}{l} B_{x 1} \\ B_{y 1} \\ B_{z 1} \end{array}],

(27)

\begin{array}{l} [\begin{array}{l} B_{xx 2} & B_{xy 2} & B_{xz 2} \\ B_{yx 2} & B_{yy 2} & B_{yz 2} \\ B_{zx 2} & B_{zy 2} & B_{zz 2} \end{array}] = [\begin{array}{l} \cos d & 0 & \sin d \\ 0 & 1 & 0 \\ - \sin d & 0 & \cos d \end{array}] \cdot \\ [\begin{array}{l} B_{xx 1} & B_{xy 1} & B_{xz 1} \\ B_{yx 1} & B_{yy 1} & B_{yz 1} \\ B_{zx 1} & B_{zy 1} & B_{zz 1} \end{array}] [\begin{array}{l} \cos d & 0 & - \sin d \\ 0 & 1 & 0 \\ \sin d & 0 & \cos d \end{array}] 。 \end{array}

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

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

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

数学上,勒让德函数是指以下勒让德方程的解:

(28)

(1 - y^{2}) \frac{d^{2} P_{n} (v)}{d v^{2}} - 2 v \frac{d P_{n} (v)}{d v} + n (n + 1) P_{n} (v) = 0 。

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

(29)

P_{n} (v) = \frac{1}{2^{n} n!} {(\frac{d}{d v})}^{n} (v^{2} {- 1)}^{n},

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

表1 0~3阶勒让德多项式

Table 1 Legendre polynomials of order 0~3

n	P_n(cosθ)	P_n(v)
0	1	1
1	cosθ	μ
2	(3cos2θ+1)/4	(3v²-1)/2
3	(5cos3θ+3cosθ)/8	(5v³-3v)/2

新窗口打开| 下载CSV

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

(30)

\begin{array}{l} (1 - v^{2}) \frac{d^{2} P_{n, m} (v)}{d v^{2}} - 2 v \frac{d P_{n, m} (v)}{d v} + \\ [n (n + 1) - \frac{m^{2}}{1 - v^{2}}] P_{n, m} (v) = 0 。 \end{array}

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

(31)

P_{n, m} (v) = (1 - v^{2})^{m / 2} \frac{d^{m}}{d v^{m}} (P_{n} (v)) 。

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

(32)

P^{n, m} (\cos θ) = \frac{2^{n} (n - m)!}{(2 n)!} P_{n, m} (\cos θ),

其中:Pⁿ^,^m(cosθ)为高斯规格化伴随勒让德多项式,P_n_,_m(cosθ)为伴随勒让德多项式。

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

(33)

P_{n}^{m} (\cos θ) = \{\begin{array}{l} P_{n, m} (\cos θ), & m = 0 \\ {[\frac{2 (n - m)!}{(n + m)!}]}^{1 / 2} P_{n, m} (\cos θ), & m > 0 \end{array}

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

(34)

P_{n}^{m} (\cos θ) = S_{n, m} P^{n, m} (\cos θ), X

其中: $P_{n}^{m}$ (cosθ)表示施密特半规格化伴随勒让德多项式,Pⁿ^,^m(cosθ)表示高斯规格化伴随勒让德多项式,S_n_,_m定义如式(35)所示:

(35)

S_{n, m} = \{\begin{array}{l} \frac{(2 n)!}{2^{n} n!}, & m = 0 \\ {[\frac{2 (n - m)!}{(n + m)!}]}^{1 / 2} \frac{(2 n)!}{2^{n} (n - m)!}, & m > 0 \end{array}

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

表2 0~3阶伴随勒让德多项式

Table 2 Order 0-3 adjoint Legendre polynomials

n	m	伴随勒让德多项式 P_n_,_m(cosθ)	高斯规格化伴随勒让德多项式 Pⁿ^,^m(cosθ)	施密特拟规格化伴随勒让德多项式 $P_{n}^{m}$ (cosθ)
0	0	1	1	1
1	0	cosθ	cosθ	cosθ
1	1	sinθ	sinθ	sinθ
2	0	$\frac{1}{2}$ (3cos²θ-1)	$\frac{1}{6}$ (3cos²θ-1)	$\frac{1}{2}$ (3cos²θ-1)
2	1	3cosθsinθ	$\frac{1}{2}$ cosθsinθ	$\sqrt[]{3}$ cosθsinθ
2	2	3sin²θ	$\frac{1}{2}$ sin²θ	$\frac{\sqrt[]{3}}{2}$ sin²θ
3	0	$\frac{1}{2}$ cosθ(5cos²θ-3)	$\frac{1}{30}$ cosθ(5cos²θ-3)	$\frac{1}{2}$ cosθ(5cos²θ-3)
3	1	$\frac{3}{2}$ sinθ(5cos²θ-1)	$\frac{1}{30}$ sinθ(5cos²θ-1)	$\sqrt[]{\frac{3}{8}}$ sinθ(5cos²θ-1)
3	2	15cosθsin²θ	$\frac{1}{6}$ cosθsin²θ	$\sqrt[]{\frac{15}{4}}$ cosθsin²θ
3	3	15sin³θ	$\frac{1}{6}$ sin³θ	$\sqrt[]{\frac{5}{8}}$ sin³θ

新窗口打开| 下载CSV

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

(36)

\begin{array}{l} P_{0}^{0} = 1 \\ P_{1}^{1} = \sin θ \\ P_{n}^{n} = \sqrt[]{1 - \frac{1}{2 n}} \sin θ P_{n - 1}^{n - 1} \\ P_{n}^{m} = \frac{2 n - 1}{\sqrt[]{n^{2} - m^{2}}} \cos θ P_{n - 1}^{m} - \sqrt[]{\frac{{(n - 1)}^{2} - m^{2}}{(n^{2} - m^{2})}} P_{n - 2}^{m} \end{array}

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

(37)

\begin{array}{l} \begin{array}{l} \frac{\partial P_{0}^{0}}{\partial θ} = 0 \\ \frac{\partial P_{1}^{1}}{\partial θ} = \cos θ \\ \frac{\partial P_{n}^{n}}{\partial θ} = \sqrt[]{1 - \frac{1}{2 n}} (\sin θ \frac{\partial P_{n - 1}^{n - 1}}{\partial θ} + \cos θ P_{n - 1}^{n - 1}) \\ \frac{\partial P_{n}^{m}}{\partial θ} = \frac{2 n - 1}{\sqrt[]{n^{2} - m^{2}}} (\cos θ \frac{\partial P_{n - 1}^{m}}{\partial θ} - \sin θ P_{n - 1}^{m}) - \end{array} \\ \sqrt[]{\frac{{(n - 1)}^{2} - m^{2}}{n^{2} - m^{2}}} \frac{\partial P_{n - 2}^{m}}{\partial θ} \end{array}

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

(38)

\begin{array}{l} \frac{\partial^{2} P_{0}^{0}}{\partial θ^{2}} = 0 \\ \frac{\partial^{2} P_{1}^{1}}{\partial θ^{2}} = \sin θ \\ \frac{\partial^{2} P_{n}^{n}}{\partial θ^{2}} = \sqrt[]{1 - \frac{1}{2 n}} (- \sin θ P_{n - 1}^{n - 1} + 2 \cos θ \frac{\partial P_{n - 1}^{n - 1}}{\partial θ} + \\ \sin θ \frac{\partial^{2} P_{n - 1}^{n - 1}}{\partial θ^{2}}) \\ \frac{\partial^{2} P_{n}^{m}}{\partial θ^{2}} = \frac{2 n - 1}{\sqrt[]{n^{2} - m^{2}}} (- 2 \sin θ \frac{\partial P_{n - 1}^{m}}{\partial θ} + \\ \cos θ \frac{\partial^{2} P_{n - 1}^{m}}{\partial θ^{2}} - \cos θ P_{n - 1}^{m} -) \cdot \\ \sqrt[]{\frac{{(n - 1)}^{2} - m^{2}}{n^{2} - m^{2}}} \frac{\partial^{2} P_{n - 2}^{m}}{\partial θ^{2}} \end{array}

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日进行计算。

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

图2

新窗口打开| 下载原图ZIP| 生成PPT

图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°。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 地磁场偏角D(a)及倾角I(b)等值线

Fig.3 Contour map of geomagnetic field deviation(a) and dip angle(b)

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

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 地磁场三分量(a)B_x(b)B_y(c)B_z等值线

Fig.4 Three-component contour map of the geomagnetic field

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 全张量地磁梯度等值线

Fig.5 Full-tensor geomagnetic gradient contour map

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

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 主对角线上3个张量分量之和等值线

Fig.6 Contour map of the sum of the three tensor components on the main diagonal

如图6所示,全张量磁梯度的3个分量B_xx、B_yy、B_zz之和的变化范围为:-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

城市	来源	B_x/nT	B_y/nT	B_z/nT	B_h/nT	B/nT	D/°	I/°
成都市	本文计算	33866.2	-1213.9	37855.4	33887.9	50807.6	-2.0528	48.1652
成都市	NOAA	33866.2	-1213.9	37855.4	33887.9	50807.7	-2.0529	48.1652
自贡市	本文计算	34644.5	-1268.7	36051.7	34667.6	50015.7	-2.0972	46.1212
自贡市	NOAA	34644.5	-1268.7	36051.7	34667.7	50015.8	-2.0973	46.1212
泸州市	本文计算	34909.9	-1334.7	35359.9	34935.3	49707.1	-2.1895	45.346
泸州市	NOAA	34909.9	-1334.7	35359.9	34935.4	49707.2	-2.1895	45.346
德阳市	本文计算	33580.9	-1266.4	38443.8	33604.7	51060.8	-2.1597	48.8424
德阳市	NOAA	33580.9	-1266.4	38443.8	33604.8	51060.8	-2.1597	48.8424
最大绝对误差		0	0	0	0.1	0.1	0.0001	0

新窗口打开| 下载CSV

4 结论

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

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

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Barraclough D

, Harwood J

, Leaton B

, et al.

A definitive model of the geomagnetic field and its secular variation for 1965—I.Derivation of model and comparison with the IGRF

[J]. Geophysical Journal International, 1978,55(1):111-121.

DOI:10.1111/j.1365-246X.1978.tb04750.x URL [本文引用: 1]

[2]

安振昌, 徐元芳.

球谐分析与国际地磁参考场

[J]. 地震地磁观测与研究, 1982,3(5):22-26.