0 引言

地磁场的变化会在地球内感应出地电场。在获取地磁观测数据后，可根据电磁感应原理，计算相应的感应地电场。感应地电场由地磁及地下介质电性结构决定，实际计算需解决的问题是，准确的大地电性结构参数的确定。但高精度大地电性数据较为缺乏，难以满足地磁感应地电场计算需求。

磁暴期间，地电场观测台站能记录到地电场的剧烈变化，即地电暴。同一区域地电暴记录与地磁场变化率无论在形态、幅度还是周期成分上均呈现出较强的相关性（张素琴等，2010），事实上，二者是同源的，由分布于磁层—电离层的空间电流体系产生。地电场记录中通常包含大地电场、自然电场和干扰，一般难以将地磁感应地电场成分从记录中分离出来。磁暴时，地磁感应地电场成分增强。当磁暴规模较大时，地电场呈现剧烈的高频起伏变化形态，原有正常日变周期被破坏（张学民等，2006）。此时，其他场源对地电场记录形态影响相对较小，地电场记录中以磁暴感应成分为主。理论上，应可由此类强磁暴记录，基于电磁场理论，通过对地电暴波形地拟合，反演大地电性结构参数。

在地下介质电性结构研究方面，目前较成熟的大地电磁测深方法，通过观测变化的电、磁场水平分量，将电、磁场信号转换成视电阻率曲线，反演求得各地层的电阻率和厚度值（中国地震局监测预报司，2007）。该方法对视电阻率曲线进行反演，不直接拟合感应地电场变化形态。

本研究利用磁暴期间电、磁场观测数据，基于一维水平层状电性介质模型，对地电暴波形进行拟合，反演地下介质电性结构，并将反演得到的电性参数用于同一场地不同磁暴感应地电场的计算，以验证模型和方法的合理性和有效性。

1 地磁感应地电场计算模型

感应地电场波形反演首先需要确定正演问题，即已知入射地磁场和地下介质电性结构，计算感应地电场。

大地电磁感应麦克斯韦方程为

$ \nabla \times \overrightarrow E = - \frac{{\overrightarrow {\partial B} }}{{\partial t}} $

(1)

对于特定地电台站，磁暴可近似为均匀无限的平面入射电磁波，方程（1）可写为

$ \frac{{\partial {E_x}}}{{\partial Z}} = \frac{{\partial {B_y}}}{{\partial t}} $

(2)

$ \frac{{\partial {E_y}}}{{\partial Z}} = \frac{{\partial {B_x}}}{{\partial t}} $

(3)

相对于很小的地电台站观测极距而言，地下介质的横向几何形态及参数变化一般不大，台站地下电性介质可视为分层水平均匀结构。在此条件下，频域内感应地电场计算公式（Lundby S et al, 1985；刘同同等，2011）为

$ {E_x}\left(\omega \right) = \frac{1}{{{\mu _0}}}{B_y}\left(\omega \right){Z_1}\left(\omega \right) $

(4)

${E_y}\left(\omega \right) = \frac{1}{{{\mu _0}}}{B_x}\left(\omega \right){Z_1}\left(\omega \right) $

(5)

式中，μ₀为真空磁导率，E_x和B_x分别为NS向地电场和地磁场，E_y和B_y分别为EW向地电场和地磁场，Z₁为地面波阻抗，ω为圆频率。

第m层顶面波阻抗Z_m可按以下递推公式确定

$ {Z_m} = {Z_{0m}}\frac{{1 + \frac{{{Z_{m + 1}} - {Z_{0m}}}}{{{Z_{m + 1}} + {Z_{0m}}}}{{\rm{e}}^{ - 2{k_m}{h_m}}}}}{{1 - \frac{{{Z_{m + 1}} - {Z_{0m}}}}{{{Z_{m + 1}} + {Z_{0m}}}}{{\rm{e}}^{ - 2{k_m}{h_m}}}}} $

(6)

式中，h_m为第m层厚度，${k_m} = \sqrt {j\omega {\mu _0}{\sigma _m}} $（σ_m为第m层电导率），Z_0m = (jωμ₀)/k_m。

底层波阻抗Z_n按下式计算

$ {Z_n} = \sqrt {\frac{{j\omega {\mu _0}}}{{{\sigma _n}}}} $

(7)

据此，已知地磁场和地下介质分层电性结构，可根据式（7），从底层开始，递推求出地面波阻抗Z₁，基于式（4）和式（5），采用快速傅里叶变换（FFT）计算地磁感应地电场频域值，通过快速傅里叶逆变换（IFFT）得到地电场时程。

2 反演流程

已获得磁暴和地电暴观测数据，地电暴波形拟合反演问题可描述为，根据已知磁暴记录，按式（4）、式（5）反演地下介质各层电导率和厚度，使地电暴计算波形与观测波形之间误差最小。

地电暴波形反演流程如下：将磁暴水平分量时程通过FFT变换到频域，由大地电性结构参数的初始值计算地面波阻抗，与磁暴频域值相乘，按公式（4）、式（5）得到理论地电暴频域序列后，经IFFT得到理论地电暴水平分量时程，将地电暴实际观测波形作为目标值，进行非线性最小二乘拟合，反演大地电性结构参数。具体反演流程见图 1。

与其他地球物理反演问题一样，地下电性结构参数具有多解性和不唯一性。在反演过程中，应尽量依据已有区域地质和地球物理资料及电性介质参数的可能取值范围，对反演参数施加必要约束，使反演结果尽量反映区域实际情况。

3 拟合反演算例

选取文昌地震台记录的2015年8月15日—16日（K = 7）和2017年8月31日—9月1日（K = 7）2个磁暴事件作为算例，检验拟合反演效果。磁暴NS分量和地电暴EW测向观测记录见图 2。其中，前一个事件用于拟合反演，后一个事件用于结果检验。

2015年8月15日—16日磁暴事件中拟合地电暴波形与观测记录对比见图 3。图中，细实线为观测数据，虚线为拟合波形。从图 3可见，拟合的地电暴波形与观测记录一致性较好。图 3中的拟合波形是在假设本场地为一维水平层状电性结构的前提下，按照图 1所示地电场波形拟合反演流程，反演电性结构参数后，得到的理论地电暴波形。拟合波形与观测数据一致性较好，表明本研究场地一维水平层状电性结构模型应用的合理性。

根据以上波形拟合反演结果得到文昌台地下电性结构参数，见表 1。

根据表 1给出的大地电性结构参数，利用2017年8月31日—9月1日磁暴观测数据，计算理论感应地电场。调整基线后，与实际地电暴观测记录进行比较，结果见图 4。从图 4可见，除个别点的幅值外，理论地电暴波形与观测记录相符合，说明反演得到的一维水平层状电性介质模型及参数，可有效应用于本场地感应地电场的计算。

4 结论

（1）基于一维水平层状电性介质模型，对文昌台磁暴感应地电场波形进行拟合反演，结果表明，地电暴拟合波形与观测记录一致性较好。利用反演确定的地下电性结构参数和其他磁暴事件记录，计算得到的理论感应地电场与观测结果相符合，验证一维水平层状模型在本研究场地应用的合理性和反演方法的有效性。

（2）本研究通过直接拟合磁暴感应地电场记录，反演地下介质电性结构。在视电阻率曲线反演方法外，提供一种新的地下电性结构反演途径，反演结果可直接应用于水平层状电性介质地磁感应地电场计算。

（3）在地电暴反演的基础上，可根据地磁记录，计算其他时刻空间电流体系感应地电场，将有助于定量分析并理解地电场观测数据中不同场源成分的影响。