基于地磁台站的电离层场源研究

Study of Ionospheric Sources Based on Geomagnetic Stations

摘要: 自90年代以来，随着地磁观测卫星的不断增加，地磁测深方法迎来了快速发展的阶段。地磁测深的激发场源为磁层环电流和电离层日变(Sq)电流。其中，磁层环电流产生的感应磁场变化周期为几天到几个月，能够提供地幔转换带下部到下地幔(500~1600 km)深度的电导率信息；电离层Sq电流产生的感应磁场变化周期主要为4小时到1天，能够提供上地幔到转换带上部(100~500 km)的电性结构。为了研究电离层电流产生的电磁响应，本文根据球谐分析方法，通过已有的分钟采样的地磁台站数据还原电离层Sq电流，然后根据不同的地球模型，通过积分方程法模拟其在地表产生的电磁场。

Abstract: Since the 1990s, with the continuous increase of geomagnetic observation satellites, the method of geomagnetic sounding has entered a stage of rapid development. The excitation sources for geomagnetic sounding are the magnetic layer ring currents and the ionospheric daily variation (Sq) currents. Among them, the induced magnetic field changes caused by the magnetic layer ring currents have a period of several days to several months, providing information about the conductivity of the lower part to the lower mantle (500~1600 km) in the mantle transition zone. The induced magnetic field changes caused by the ionospheric Sq currents have a period mainly ranging from 4 hours to 1 day, providing information about the electrical structure from the upper mantle to the upper part of the transition zone (100~500 km). In order to study the electromagnetic response generated by ionospheric currents, this paper, based on the spherical harmonic analysis method, reconstructs the ionospheric Sq currents using existing minute-sampled geomagnetic station data. Subsequently, based on different Earth models, it simulates the electromagnetic field generated by these currents on the Earth’s surface using the integral equation method.

1. 前言

受太阳活动的影响，地球外部磁层和电离层存在时变电流体系，这些电流产生时变磁场，由于电磁感应，这些时变磁场会在导电的地球内部产生感应电磁场。地表或者卫星观测的感应电磁场数据包含了丰富的频谱成分，根据电磁波的趋肤效应，不同频率的数据携带了地球内部不同深度的电导率信息 ‎[1] ‎[2] 。因此，通过深部电磁感应成像研究地幔电导率对进一步研究地球深部的热结构和流变性特征、探讨地幔动力学演化过程、解释重大地质现象等具有重要的科学意义。

地磁测深的激发场源为磁层环电流和电离层日变(Sq)电流。其中，磁层环电流产生的感应磁场变化周期为几天到几个月，能够提供地幔转换带下部到下地幔(500~1600 km)深度的电导率信息；电离层Sq电流产生的感应磁场变化周期主要为4小时到1天，能够提供上地幔到转换带上部(100~500 km)的电性结构 ‎[3] ‎[4] ‎[5] ‎[6] ‎[7] 。然而，目前的地磁测深研究主要考虑空间结构比较简单、周期较长的磁层环电流，很少考虑空间结构更为复杂、周期较短的电离层Sq电流。

由于长周期感应电场幅值微弱观测困难，地磁测深主要使用磁场数据，主要数据来源为全球地磁台站数据和地磁卫星数据。全球地磁台站数据观测时间长达几年至几十年，经过数据处理可以得到周期长达几天到一百多天的响应函数进行地球深部的电性结构成像，但是台站数据分布稀疏，在南半球和海洋区域缺少台站覆盖。地磁卫星数据具有覆盖全球特别是海洋区域的优点，可以提高南半球和海洋区域地幔电导率成像的分辨率。然而，卫星数据存在时空混叠，由于缺乏可靠的数据处理和正反演等手段，目前地磁测深三维反演仍以分布稀疏的地磁台站数据为主，没有或很少使用覆盖全球的地磁卫星数据，导致已有的全球三维电导率模型在南半球和海洋区域分辨率低且不同模型之间差异较大。本文针对目前已有的全球地磁数据，建立不同的全球和区域一维及三维电导率模型，并结合积分方程法计算不同的电磁响应。

2. 数据挑选

围绕地球的电磁变化涵盖了一个广泛的频率范围，受到各种内部和外部物理过程的影响。在地磁台站记录的信号中，除了外部电流引起的磁场变化，还包括地球主磁场的长期变化、海水运动产生的场以及磁层的贡献。这些场的存在使得在电离层Sq场源时面临着挑战。本文使用的数据来自World Data Center (WDC)，包括了从1997年到2019年的地磁台站时间序列数据.通过以下步骤进行数据处理：

首先，本文选择了台站数据的时间段，遵循了ISGI (International Service of Geomagnetic Indices) 2018年的定义，将3小时内的aa指数低于13 nT作为磁静时期的标准。这个标准确保了所选取的数据处于相对较为静态的磁场环境中，降低了外部因素的影响。

另外，考虑到研究的电离层电流来自极区电离层，限制了AE指数(极区指数)小于200 nT。AE指数是描述极区磁活动程度的指标，限制其值是为了确保所选取的数据受极地电流系统的影响较小，同时增加对极地地区的数据约束 ‎[8] 。

Figure 1. Global distribution of geomagnetic stations

图1. 全球地磁台站分布

3. 方法

3.1. 估计场源

首先通过频域中的麦克思韦方程组来表述。

$\nabla \times \stackrel{\to }{B}={\mu }_0\sigma \stackrel{\to }{E}+{\mu }_0\stackrel{\to }{J^{ext}}$ (1)

$\nabla \times \stackrel{\to }{E}=i\omega \stackrel{\to }{B}$ (2)

其中， $\text{<mover accent="true"> B <mo> → </mo> </mover>}\equiv \text{<mover accent="true"> B <mo> → </mo> </mover>}\left(\text{r},\text{ω}\right)$ 和 $\text{<mover accent="true"> E <mo> → </mo> </mover>}\equiv \text{<mover accent="true"> E <mo> → </mo> </mover>}\left(\text{r},\text{ω}\right)$ 分别是磁场和电场的傅立叶变换； $\text{<mover accent="true"> r <mo> → </mo> </mover>}=\left(\text{r},\text{θ},\varphi \right)$ ，其中r，θ，ϕ分别表示距离地心的距离、极角和经度；而 $\omega =\frac{2\pi }{T}$ 表示角频率，其中T为周期。 $\stackrel{\to }{J^{ext}}\equiv \stackrel{\to }{J^{ext}}\left(r,\omega \right)$ 是外部源电流密度

的傅立叶变换。σ(r)是介质中的电导率分布，μ₀是自由空间的磁导率。这个表述忽略了位移电流，因为在考虑的小时和更长周期范围内可以忽略它们。采用了傅立叶变换为。

$f(t)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }f}\left(\omega \right){\text{e}}^{-i\omega t}$ (3)

在导电地球的上方的无源区域，但在电流 $\stackrel{\to }{J^{ext}}$ 所包围的区域之下，方程(1)简化为 $\nabla \times \stackrel{\to }{B}=0$ 。因此， $\stackrel{\to }{B}$ 是一个势场，并且可以被写成标量磁势V的梯度，即：

$\stackrel{\to }{B}=-\nabla V$ (4)

由于 $\stackrel{\to }{B}$ 是旋度零场，即 $\nabla \cdot \stackrel{\to }{B}=0$ ， 满足拉普拉斯方程 ${\nabla }^2\cdot V=0$ ，可以表示为外部(感应)部分和内部(感应)部分之和，

$V=V^{ext}+V^{\mathrm{int}}$ (5)

其中，V可按球谐函数展开成级数：

$V={\displaystyle \underset{n,m}{\sum }a}\left[{\epsilon }_n^m\left(\omega \right){\left(\frac{r}{a}\right)}^n+{\iota }_n^m\left(\omega \right){\left(\frac{a}{r}\right)}^{n+1}\right]Y_n^m\left(\vartheta ,\phi \right)$ (6)

在这里，a = 6371.2 km是地球的平均半径， ${\epsilon }_n^m\left(\omega \right)$ 和 ${\iota }_n^m\left(\omega \right)$ 分别是势场感应部分和感生部分的球谐系数， $Y_n^m$ 是n次和m阶的球谐函数，其中

$\underset{\text{n},\text{m}}{{\displaystyle \sum }}={\displaystyle \underset{\text{m}=\text{P}-1}{\overset{\text{P}+1}{\sum }}{\displaystyle \underset{\text{n}=\text{m}}{\overset{\text{m}+3}{\sum }}}}$ (7)

$Y_n^m\left(\vartheta ,\phi \right)=P_n^{\left|m\right|}\left(\cos \vartheta \right){\text{e}}^{im\varphi }$ (8)

P为电离层阶数，代表的周期为T = 86,400/P，因此，由公式(4)~(8)可知，地球表面上方或上方的磁场r ≥ a可写成

$B_r\left(\stackrel{\to }{r},\omega ,\sigma \right)=-{\displaystyle \underset{n,m}{\sum }\left[{\epsilon }_n^m\left(\omega \right)n{\left(\frac{r}{a}\right)}^{n-1}-\left(n+1\right){\iota }_n^m\left(\omega \right){\left(\frac{a}{r}\right)}^{n+2}\right]}{Y}_n^m\left(\vartheta ,\phi \right)$ (9)

$B_H\left(\stackrel{\to }{r},\omega ,\sigma \right)=-{\displaystyle \underset{n,m}{\sum }\left[{\epsilon }_n^m\left(\omega \right){\left(\frac{r}{a}\right)}^{n-1}+{\iota }_n^m\left(\omega \right){\left(\frac{a}{r}\right)}^{n+2}\right]}{\nabla }_{\perp }{Y}_n^m\left(\vartheta ,\phi \right)$ (10)

$J^{ext}=\frac{\delta \left(r-b\right)}{{\mu }_0}{\displaystyle \underset{n,m}{\sum }\frac{2n+1}{n+1}}{\epsilon }_n^m\left(\omega \right){\left(\frac{b}{a}\right)}^{n-1}{\nabla }_{\perp }{Y}_n^m\left(\vartheta ,\phi \right)$ (11)

这里 $B_H=\left(B_{\theta },B_{\varphi }\right)$ 并且 ${\nabla }_{\perp }$ 表示梯度的角度部分，即

${\nabla }_{\perp }=\stackrel{\to }{{\text{e}}_{\vartheta }}\frac{\partial }{\partial \vartheta }+\stackrel{\to }{{\text{e}}_{\varphi }}\frac{1}{\sin \vartheta }\frac{\partial }{\partial \varphi }$ (12)

$\stackrel{\to }{{\text{e}}_{\vartheta }}$ 和 $\stackrel{\to }{{\text{e}}_{\varphi }}$ 为球面坐标系的单位切向量。

Figure 2. Distribution of surface magnetic field at P = 1 order

图2. P = 1阶地表磁场分布

根据公式(9)-(10)可以计算出不同阶数的电离层电流与其产生的磁场，其地表台站分布见图1。由图2所示，当P=1阶时，外部磁场的大小在10nT左右，并且主要分布在中低纬度地区，这与该频率下电离层电流产生的磁场的主要特征一致。

3.2. 电离层电磁电流响应

Figure 3. Global 1D conductivity model used in forward modeling

图3. 正演使用的全球电导率1维模型

Figure 4. Global 3D conductivity model used in forward modeling

图4. 正演使用的全球电导率3维模型

Figure 5. Distribution of total surface magnetic field in P = 1 Order 1D earth model forward modeling

图5. P = 1阶1D地球模型地表磁场总场分布

Figure 6. Distribution of total surface magnetic field in P = 1 Order 3D earth model forward modeling

图6. P = 1阶3D地球模型地表磁场总场分布

本文仅展示当阶数P = 1, 2时根据不同地球介质计算的地表电磁场，不同的地区介质见图3，图4。使用的方法是Mikhail Kruglyakov开发的开源程序积分方程法，详情可见 ‎[2] 。

根据不同介质的地表电磁场响应，P = 1阶下不同地球介质感应的电磁场几乎一致，这与电磁场的趋肤效应有关，当周期越大时，电磁场穿透的深度越大，受地表电导率的影响就越小。当P的阶数逐渐增大时，不同地表电导率感应出来的电磁场差异逐渐明显，并且表现出强烈的海洋–陆地电导率差距。图5、图6分别为当P = 1阶时，利用积分方程法计算的地表电磁总场分布，可以看出，当P = 1阶时，由于电磁场趋肤深度的原因，不同地表结构对于电磁场的影响不大，磁场均在10 nT左右。当P = 2阶时，图7和图8的磁场总场大小在5 nT左右变化，且受地表电导率的影响严重，在电导率变化严重的地区，磁场的变化通常也较大。

Figure 7. Distribution of total surface magnetic field in P = 2 Order 1D earth model forward modeling

图7. P = 2阶1D地球模型地表磁场总场分布

Figure 8. Distribution of total surface magnetic field in P = 2 Order 3D earth model forward modeling

图8. P = 2阶3D地球模型地表磁场总场分布

4. 结论

球谐分析作为能够将地磁场分离成内外部场贡献的一种经典方法，广泛应用于地磁测深方法中，用于估计地球电导率。本文通过地磁台站的分钟采样的时间序列磁场数据，基于球谐分析理论，找到产生其磁场的电流的最小二乘解。最后通过积分方程法实现了不同阶数的以及不同介质的电磁场正演模拟，并且验证了正演的程序，为接下来的地球电导率反演提供了一定基础。