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

link管理

链接快照平台

输入网页链接，自动生成快照
标签化管理网页链接

相关文章推荐

谦逊的开心果 · Gibson Global Search ...· 3 月前 ·

安静的吐司 · 《天降神运完整版》高清免费在线播放 - ...· 11 月前 ·

深情的火车 · ORM查询-Group/Order/Havi ...· 1 年前 ·

健壮的冲锋衣 · Day 24 - Canvas 常用套件 ...· 1 年前 ·

光明磊落的高山 · [mysql] Include ...· 1 年前 ·

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

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

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： spheric harmonic analysis legendre polynomials

本文引用格式

钟炀, 管彦武, 石甲强, 肖锋. 利用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

为了描述地磁场时空非线性变化特征,自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模型计算地球主磁场的全张量磁梯度有着迫切的需求,故笔者梳理了地磁场七要素及张量的计算流程,推导出计算任意给定点位全张量地磁梯度的公式,绘制出某地区地磁场全张量等值线图,并验证了算法的正确性。

目前,球谐函数法被广泛应用于计算地磁场参考模型中。在该方法中,认为地球是均匀球体,在球坐标系下地磁场源区外任意一点坐标( 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 点表示测点。

(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 轴沿旋转椭球面法线向下的局部直角坐标系(北东下)的地磁场三分量及全张量公式。

勒让德微分方程(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 所示。

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	(3 v ² -1) / 2
3	(5cos3 θ +3cos θ ) / 8	(5 v ³ -3 v ) / 2

新窗口打开 | 下载CSV

常微分方程式(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 所示。

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}$ (3 cos ² θ -1)	$\frac{1}{6}$ (3 cos ² θ -1)	$\frac{1}{2}$ (3 cos ² θ -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

(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}

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

为进一步验证算法的正确性,选取计算工区内4个随机点位:①四川省成都市(30.67N,104.07E),②自贡市(29.35N,104.78E),③泸州市(28.87N,105.43E),④德阳市(31.13N,104.38E)),利用美国国家海洋和大气管理局(NOAA)官方网站(www.ngdc.noaa.gov)内全球地磁场数据,对比两种方式计算获得地磁场七要素误差见表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 In this paper, the seventh generation IGRF was introduced.Based on the IGRF1995, the grid values were calculated, the geomagnetic charts were drawn for the everygeomagnetic component（D, I, H, F, X, Y, Z） over China and adjacent areas at epoch1995.0. The ellipticity of the Earth was accounted, the effect of ellipticity was computedand analysed. 本文利用东亚及其邻近地区50个地磁台站1955—1977年地磁三要素（D、H、Z）的实测年平均值算出的非偶极场值X、Y、Z,选用“按距离加权最小二乘曲面拟合”方法,计算、绘制、分析了东亚地区地磁非偶极场及其长期变化的时空分布特征,并与国际参考场（IGRF）模式进行了粗略的对比。所得的主要结果是:1.1955.0—1965.0年代,东亚地区Z分量的增长速度在迅速减小;1965.0年代以后,Z值在逐年下降。2.从整体上看,东亚地区地磁非偶极场虽然形态与IGRF给出的结果基本吻合,但长期变化有明显差异。1965.0年代的IGRF长期变化模式与实测结果偏差甚大。为了使IGRF模式更好地拟合地磁长期变化处于转折时期的复杂地区的情况,看来时间间隔应更短些,例如5年为好,采用现在通用的方法是不适宜的。人们常用球谐分析方法来表示地磁场的全球分布(即全球模型,Definite Geomagnetic Re-ference Field,以下简称DGRF),用多项式表示地磁场的区域分布(即区域模型).本文以1980年中国地区地磁场区域模型(Chinese Geomagnetic Reference Field,以下简称CGRF)和DGRF1980为例,计算出地磁场区域模型与全球模型的差值(Dif=CGRF1980-DGRF1980),用电子计算机绘制出差值的等值线图,并进行了讨论.从差值等值线图可以看出:对于中国大部分地区,地磁场的区域模型与全球模型之间的差值比较小,但在西北和东北边境地区,其差值明显增大. Spherical harmonic analytical method is usually used to express the globaldistribution of the geomagnetic field (global model), and polynomial method isemployed to express the regional distribution of the geomagnetic field (regionalmodel). In this paper, on the basis of the regional models of the geomagneticfield in China at epoch 1980.0 (CGRF 1980) and the global models of the geo-magnetic field at epoch 1980.0 (DGRF 1980), the differences between CGRF1980 and DGRF 1980 were calculated (diff. =CGRF 1980-DGRF 1980), andcontour maps of the differences were constructed for total intensity (F), verti-cal intehsity (Z), horizontal intensity (H), declination (D) and inclination (I),with some analyses and discussions made. These maps show that the differencesbetween CGRF 1980 and DGRF 1980 are comparatively small in most parts ofChina, but obviously grow larger in border areas of Northwest and Northeast China. AbstractVertical cutoff rigidities derived from the International Geomagnetic Reference Fields (IGRF) are normally used to compute the radiation dose at a specific location and to organize the radiation dose measurements acquired at aircraft altitudes. This paper presents some of the usually ignored limits on the accuracy of the vertical cutoff rigidity models and describes some of the computational artifacts present in these models. It is noted that recent aircraft surveys of the radiation dose experienced along specific flight paths is sufficiently precise that the secular variation of the geomagnetic field is observable.]]>

A definitive model of the geomagnetic field and its secular variation for 1965—I.Derivation of model and comparison with the 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在电离层仍然具有很高的精度. ...

Ionosphere geomagnetic field: Comparison of IGRF model prediction and satellite measurements 1991—2010

Revision of the high-degree Gauss coefficients in the IGRF 1945—1955 models by using natural orthogonal component analysis

... 为适用于高精度地质构造研究工作,徐文耀 ^{[

12

]} 指出在1945~1955年IGRF模型中的高阶高斯系数( n >7)与其他时期的模型相比,存在较大偏差,并通过自然正交分量分析法对1945~1955年IGRF模型中的高次高斯系数进行修正,使这些模型的高斯系数具有较好的平滑时变特性.Maus等 ^{[

13

]} 美国地球物理数据中心(NGDC)的科学家们计算出16~720阶球谐系数,成功建立岩石圈磁场模型NGDC-720.张素琴等 ^{[

14

]} 利用中国部分地磁台站1990~2003 年的年均值资料, 研究了地磁台站年均值与IGRF模型计算值的一致程度,结果显示虽然存在差异,但是各地磁要素年均值与IGRF模型值的差值的标准差均低于IGRF模型的误差水平,故一致性较好. ...

A near-surface geomagnetic field model to spherical harmonic degree 720

Considerations and proposal for a best utilization of IGRF over areas including a geomagnetic observatory

... 在IGRF模型的应用方面, Molina(1987)等 ^{[

15

]} 针对意大利等欧洲中部国家,建立了基于IGRF模型的意大利地磁参考场(ITGRF)和欧洲中部地磁参考场(EGRF).任国泰 ^{[

16

]} 利用IGRF模型结合地面测量资料和地磁台资料分析了东亚大陆磁场.安振昌 ^{[

17

]} 比较和讨论了中国地区地磁场区域模型与全球模型,得出地磁场区域模型,更详细准确地表示我国地磁场的分布.Smart等 ^{[

18

,

19

]} 运用磁场模型与运动轨迹学原理,利用IGRF模型导出飞机飞行时截止刚度(VCR),进而通过线性拟合确定其飞行轨迹.Davis ^{[

20

]} 利用Matlab编程实现了IGRF模型关于场源外任意位置三分量地磁场的计算. ...

Magnetospheric models and trajectory computations

The limitations of using vertical cutoff rigidities determined from the IGRF magnetic field models for computing aircraft radiation dose

Mathematical modeling of earth’s magnetic field

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

24

]} 如式(2)所示: ...

... 其中: A 、 B 是与球谐函数边界条件有关的系数; a _nm 、 b _nm 为与球谐函数展开阶数有关的系数;

P_{n}^{m}

( cos θ )为 n 阶 m 次的规格化缔和勒让德多项式.由式(2)可见,磁位可分为两部分,一部分含

\frac{1}{r^{n + 1}}

因子,另一部分含 r ⁿ 因子.根据地磁场高斯理论,按照场源位置,地磁场分为内源场和外源场(内源场即磁场源位于地球内部,外源场磁场源位于地球外部).地磁场可进一步表示为 ^{[

24

]} : ...

Simulation of gravity gradients: a comparison study

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

25

]} : ...

Handbook of mathematical functions with formulas,graphs,and mathematical tables

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

26

]} .勒让德多项式的微分形式,可用罗德里格斯公式表示: ...