1 前言

瞬变电磁法（Transient Electromagnetic Method，TEM）是一种建立在电磁感应原理基础上的时间域人工源电磁探测方法，具有观测纯二次场、灵敏度高、施工效率高等特点，在矿产资源探测、工程地质勘察、地下水调查、地下未爆物识别等领域得到了广泛的应用。瞬变电磁正演是研究瞬变电磁响应规律的有效途径以及指导探测施工的理论依据，也是反演和资料处理解释的基础。

自上世纪70年代开始，国内外学者逐步开展瞬变电磁法一维正演。一维正演大多先在频率域中计算层状介质的电磁响应，然后将频域响应转换到时间域。从频率域到时间域响应的转换可采用的方法有：G-S逆拉氏变换方法、延迟谱方法、线性数字滤波方法、余弦变换法等。瞬变电磁一维正演问题较为成熟，有诸多优秀的研究著作可供参考，国外代表性著作有Kaufmann主编的《频率域和时间域电磁探深》以及Nabighian主编的《探测地球物理-电磁法》等，国内学者主编的书籍有李貅的《瞬变电磁测深的理论与应用》以及牛之琏的《时间域电磁法原理》等。

随着计算机技术的发展，瞬变电磁法由一维正演向三维正演发展。针对复杂地质结构，一维正演适用性差，利用数值模拟方法进行瞬变电磁三维正演更加符合实际情况。从上世纪八十年代开始，国内外研究学者开始在瞬变电磁三维正演方面做了大量的研究工作。瞬变电磁三维正演包含空间域离散和时间域离散两个关键计算问题。空间域离散常用的方法主要有：积分方程法、有限差分法、有限单元法和有限体积法；时间域离散的方法主要包括频时转换算法、时间步迭代算法（隐式向后差分、中心差分和显式向前差分）、及基于Krylov子空间技术的近似方法。随着观测环境的复杂化以及探测精度要求的提高，瞬变电磁法面临越来越多的挑战，进一步改进瞬变电磁法三维数值模拟的精度和效率，对提升瞬变电磁数据处理效果和解释水平具有重要的理论和现实意义。

2 三维数值模拟方法

电磁场方程组的描述形式分为积分法和微分法。常用的瞬变电磁三维数值模拟方法中，积分方程法和有限体积法属于积分法，有限差分法、有限单元法则为微分法。本节将简短的介绍这几种三维数值模拟方法，并通过国内外学者的工作，探究TEM三维正演模拟的发展方向。

2.1 积分方程法

积分方程法（Integral Equation Method, IEM）是最早实现三维电磁场模拟的数值计算方法，它只需对小体积的异常体进行离散化，不必计算整个区域的场。1985年SanFilpo和Hohmann利用积分方程法对矩形回线源求解了均匀半空间中三维地质体的瞬变电磁响应；1995年Xiong等基于积分方程法求解了层状大地中含多个三维棱柱体异常目标的瞬变电磁响应；2006年Zhdanov等采用积分方程法求解了非均匀背景下三维复杂模型的瞬变电磁响应；2010年Cox等利用积分方程法实现了对航空电磁数据的三维正演模拟。

积分方程法适用于简单模型的电磁响应计算，对内存需求较少，在数值计算的早期，比微分方程法具有更大的优越性。然而，面向复杂地电模型时，求解过程中需要计算大型密实矩阵和格林函数，计算成本激增，尽管后人采用近似算法提高了积分方程法的计算效率，也降低了计算结果的精度。

2.2 有限差分法

有限差分法(Finite Difference Method，FDM)原理简单且易于实现，在电磁正演计算中应用较为广泛。它将求解域划分为规则的Yee网格，用有限个网格节点代替连续的求解域。Yee网格单元将磁场分量置于各面的中心，将电场分量置于各棱边的中间，这样的网格方式使场分量在突变面上的连续性条件得到自然满足，有利于复杂结构的计算。有限差分法的优势是不需要求解大型方程，因此其求解速度快且能够较好地实现并行运算。

1993年Wang和Hohmann提出了一种基于Du Fort-Frankel时域有限差分格式的三维瞬变电磁场数值模拟算法；1999年日本学者Endo和Noguchi Wang以Wang和Hohmann的算法为基础，利用流体力学中坐标变换的方法将物理域转换至计算域，解决了带地形模型的回线源瞬变电磁三维正演问题；2004年Commer和Newman对Wang和Hohmann的算法进行了改进，直接求解电场的泊松方程以得到电流关断前的电场，实现了接地长导线源激发瞬变电磁场的三维数值模拟；2013年孙怀凤等对Wang和Hohmann提出的方法进行改进，将矩形回线源的电流密度直接代入麦克斯韦方程中，实现了矩形回线源瞬变电磁三维时域有限差分正演研究，并在激发电流源的计算中考虑关断时间；同年，邱稚鹏等利用非正交网格的时域有限差分方法对带地形的瞬变电磁三维问题进行了仿真，给出了地形对瞬变电磁场的影响规律；2014年李展辉和黄清华将卷积完全匹配层（CPML）边界条件应用于瞬变电磁场三维有限差分正演中，对CPML边界的吸收性能进行了分析；2016年常江浩等将有限差分法应用于矿井全空间瞬变电磁场的计算，对煤矿典型富水区的瞬变电磁响应进行了模拟和分析；2017年长安大学赵越等基于交错网格的时域有限差分方法对三维起伏地形条件下航空瞬变电磁进行正演模拟。

对于发射源加载问题，目前主要3种方法：1）利用发射电流关断后一个非常小时刻在均匀半空间产生的解析解作为电磁场初始值；2）取关断前建立的稳定场作为场的初始值；3）将源电流密度加入麦克斯韦方程，一般采用梯形电流波形，也可设置为其它任意波形。边界条件是影响晚期场精度的主要因素，有效的吸收边界不仅能提高解的精度，而且能减少网格数量。传统的有限差分模拟一般采用狄利克雷边界，随着研究的深入，目前的有限差分已经较多采用了吸收边界条件，大大提高了晚期解的精度。

随着研究的不断深入，有限差分法仍面临着很多问题：1）在时间迭代上采用的是显式差分格式，计算量大；2）只能采用结构化网格，限制了该方法的发展；3）在截断边界处多采用狄里克雷（Dirichlet）边界条件，影响晚期精度和计算速度。

2.3 有限元法

有限元法（Finite Element Method，FEM）是一种求解偏微分方程边值问题近似解的数值方法，它将求解域分解为有限个单元，在每个单元假设一个合适的近似解，推导求解域内边值问题满足的线性方程组，从而得到问题的解。有限元法既可使用规则矩形或六面体网格对计算区域进行剖分，又可使用不规则的三角形或四面体网格对计算区域进行剖分，非常适用于复杂模型的数值模拟。

有限元法于1956年由Turner和Clough首次提出，目的是解决弹性力学问题；1971年Coggon利用能量最小化原理推导了电磁变分方程，率先将有限元方法引入到地球物理领域中；1980年Kuo和Cho首先采用三维时域矢量有限元法对低阻矿体在低阻屏蔽层下的瞬变电磁响应进行了模拟。之后，国内外地球物理学者对于三维时域有限元有了大量的研究。2008年Börner等提出了一种高效求解瞬变电磁场的方法，该方法首先利用有限元法求解少量频率域电磁场，然后通过余弦变换转化为时间域电磁场响应；针对电各向异性介质计算问题，2010年Um等基于矢量有限元法研究了电各向异性条件下水平接地线源的海洋瞬变电磁响应特征，之后在2012年将该方法应用于垂直线源的瞬变电磁响应特征研究；为了提高大型方程组求解速度，2015年Fu等基于iChol分解和预处理共轭梯度（PCG），采用一种并行迭代算法求解大型稀疏线性方程组，提高了计算效率；2016年殷长春等提出了自适应非结构化网格的有限元法模拟三维频域航空电磁响应；针对复杂形状回线源问题，2018年李建慧等采用非结构化网格将实际回线源剖分为电偶极子，实现了时域有限元计算；2019年邱长凯等人采用块状有理Krylov子空间算法加快了多源瞬变电磁响应的求解，并将其成功应用于海洋瞬变电磁响应计算；同年，Zeng等采用时域有限元法对瞬变电磁全发射电流波形的影响进行了研究。

有限元法虽然可以用非结构化网格离散求解域，有利于带地形或不规则异常体的正演，但仍面临许多问题：

1）网格剖分：有限元法的计算精度和效率很大程度上取决于网格剖分的密度，当网格较密时，计算精度高，但计算效率低，因此如何在保证计算精度的条件下提高计算速度是目前主要考虑的问题。

2）时间离散：有限元法求解时域解的方式主要包括三种：一种是先求频域解，再经过时频转换技术得到时域解，计算精度受到时频转换技术的影响；第二种是采用差分法直接在时间域求解，计算稳定性受到差分的限制；第三种是基于有理Krylov子空间算法，计算精度和效率取决于模型电导率反差及子空间正交基的求解。选择合适的时间方式，兼顾计算精度的效率对瞬变电磁响应的快速有效正演至关重要。

3）方程组的求解：方程组的求解不仅决定了正演模拟结果的准确性，在采用向后差分或有理Krylov子空间算法时更决定了正演模拟的计算效率。目前常用的方程组求解方法主要是直接法，该方法计算精度较高，但是不适用于大型、超大型矩阵方程组的求解。而迭代法可用于巨型方程组的求解，但其精度及计算速度受到预条件的影响。由此看来，随着计算模型复杂化、多元化，大型线性方程组的求解对整个数值计算也十分重要。

2.4 有限体积法

有限体积法（Finite Volume Method，FVM）是一种基于物理量守恒的空间离散方法。该方法将计算区域离散为一系列规则或不规则的控制体积，在控制体积的离散网格上对待解的控制方程进行积分，之后对积分式进行离散化处理得到离散方程组，通过求解该离散方程组从而得到数值解。该方法适用于任意类型的单元网格，对复杂边界形状的区域亦可满足边界问题的需要。与有限差分法相比，有限体积法精度不仅取决于对导数处理的精度，还取决于积分的精度，而且积分过程能够保证物理量的守恒特性。与有限元法一样，有限体积法也适用于结构化或非结构化网格。有限体积法在瞬变电磁三维正演中已有较多的应用。2004年Haber等采用有限体积法对瞬变电磁场进行计算，采用一阶向后Euler法对时间进行离散；2018年周建美等采用自然边界条件将电磁场控制方程转化为弱形式，采用结构化交错网格进行有限体积空间离散，实现了双轴各向异性地层回线源瞬变电磁三维正演。

  从空间离散角度，有限体积法可视作有限元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律（即插值函数），并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑在网格点之间如何变化。有限体积法只寻求结点值，这与有限差分法类似，但有限体积法在寻求控制体积的积分时，又必须假定值在网格点之间的分布，这又与有限单元法类似。有限体积法的优点是可以用于不规则网格，且适用于并行计算，其关键问题在于提升精度困难，一般最高只能达到二阶精度。