地球物理学报  2021, Vol. 64 Issue (1): 343-354   PDF    
瞬变电磁Crank-Nicolson FDTD三维正演
孙怀凤1,2, 柳尚斌1,2, 杨洋1,2     
1. 山东大学 岩土与结构工程研究中心, 济南 250061;
2. 山东大学 地球电磁探测研究所, 济南 250061
摘要:时域有限差分(FDTD)方法使用Yee网格剖分电磁场的空间采样, 通过时间步迭代实现电磁场数值模拟, 具有内存消耗低、计算简单等特点, 常用于瞬变电磁三维正演.然而, 常规FDTD方法的时间迭代步长Δt受Courant-Friedrich-Lewy(CFL)条件严格限制, 过多的迭代次数以及过密的采样往往导致计算速度慢、累积误差不断增大.本文提出一种不受CFL条件约束的无条件稳定隐式差分算法Crank-Nicolson FDTD(CN-FDTD)用于瞬变电磁三维正演.基于Crank-Nicolson差分方法对Maxwell方程组重新离散, 空间网格仍然采用Yee元胞, 时间步进采用在整时间步电场、磁场同时采样的策略, 建立无条件稳定FDTD格式, 突破CFL条件限制.与常规FDTD交替采样相比, CN-FDTD电场、磁场同时采样的策略构成的隐式差分格式, 需要求解大型稀疏矩阵方程组.通常, 瞬变电磁三维正演模型中产生的矩阵阶数往往较大, 需要占用大量内存和求解时间.为解决上述问题, 采用Crank-Nicolson-cycle-sweep-uniform(CNCSU-FDTD)方法近似求解CN-FDTD方程, 在保证求解精度的同时, 计算效率大幅提高.在边界条件处理上, 采用双线性变换推导了复频率参数完全匹配层(CFS-PML)吸收边界.采用均匀半空间模型、四类三层模型进行精度验证, 发现CN-FDTD三维正演结果与解析解、线性数字滤波解吻合较好.之后, 与接触带上的低阻复杂模型进行对比, 结果显示CN-FDTD正演结果与矢量有限元、有限体积法以及FDTD计算结果吻合较好.在此基础上, 研究了时间步放大对CN-FDTD计算精度的影响, 发现最大时间步放大到常规FDTD的3200倍时才会在晚期出现较明显的误差.在一台CPU为Intel Core i5-7300HQ的笔记本电脑单线程计算条件下, 模拟到关断后30 ms仅需要50 min.在进行并行化后, 将有望实现复杂模型分钟级的三维正演, 从而为三维反演提供可靠、快速的正演方法.
关键词: 瞬变电磁      三维正演      CN-FDTD      双线性变换      CFS-PML     
Crank-Nicolson FDTD 3D forward modeling for the transient electromagnetic field
SUN HuaiFeng1,2, LIU ShangBin1,2, YANG Yang1,2     
1. Geotechnical and Structural Engineering Research Center, Shandong University, Jinan 250061, China;
2. Laboratory of Earth Electromagnetic Exploration, Shandong University, Jinan 250061, China
Abstract: The time-domain finite-difference (FDTD) method defines the spatial position of the electromagnetic field using the Yee grid and can simulate the propagation of the electromagnetic field with low memory cost. It is often used in three-dimensional forward modeling of transient electromagnetics. However, the time-step size of the conventional FDTD is strictly limited by the Courant-Friedrich-Lewy (CFL) stability condition. Too many iterations and intensive sampling result in an increase in calculation time and cumulative error. In this study, an implicit difference algorithm, called Crank-Nicolson FDTD (CN-FDTD), is proposed for forward modeling of three-dimensional transient electromagnetics. This algorithm is no longer constrained by the CFL condition and is unconditionally stable. The space grid still adopts the Yee cell. The electric field and magnetic field are simultaneously sampled in integer time steps to establish an unconditionally stable Crank-Nicolson scheme. Compared to conventional FDTD, CN-FDTD introduces an implicit difference with the strategy of sampling in space and time simultaneously. It needs to solve large sparse matrix equations at each step and will consume a lot of physical memory and time. Crank-Nicolson-cycle-sweep-uniform (CNCSU-FDTD) is applied to solve this problem. This approximation method greatly improves the calculation efficiency while ensuring the accuracy of the solution. The bilinear transformation method is used to derive the complex-frequency-shifted perfectly matched layer (CFS-PML) absorption boundary condition. Tests are conducted on the uniform half-space and three-layered models to examine the accuracy of CN-FCTD. The forward modeling results in this approach are in good agreement with the analytical solution and the linear digital filter solutions. The accuracy of the algorithm is further tested using a complex low-resistance model with a vertical contact zone. The CN-FDTD result is in good agreement with the finite element, finite volume, and original FDTD results. Then the effect of time-step magnification on the calculation accuracy of CN-FDTD is studied. When the time step is enlarged to 3200 times of the original FDTD algorithm, the obvious error appears in the late stage. All tests are conducted on a laptop with Intel Core i5-7300HQ CPU. It takes only 50 minutes to finish forward modeling on a complex model to 30 ms after the current is turned off with single-thread. By parallelization, it is expected to realize 3D forward modeling in minutes and provide a reliable and fast forward modeling algorithm for inversion.
Keywords: Transient electromagnetics    Three-dimensional forward modeling    CN-FDTD    Bilinear transformation    CFS-PML    
0 引言

瞬变电磁法(TEM)利用不接地回线或接地长导线向大地发射一次场, 在电流关断间隙, 通过采集到的感应电动势分析地下电性分布, 在工程勘察、金属矿产资源勘探、环境与水文地质调查、采空区探测等领域得到广泛应用(Chang et al., 2019a; Chang et al., 2019b; 薛国强等, 2007, 2008).瞬变电磁三维正演是进行勘探可行性研究、方案设计、观测数据反演与解释的必备工具.国际上, 瞬变电磁或时间域电磁三维正演研究主要有五类方法, 即体积分方程法(VIE)(Cox et al., 2010; Xiong, 1992; Zhdanov et al., 2006)、有限单元法(FEM)(Sugeng, 1998; Li et al., 2017; Li et al., 2011)、时域有限差分法(FDTD)(孙怀凤等, 2013; 许洋铖等, 2012; Commer et al., 2004; Wang et al., 1993)、有限体积法(FV)(Zhou et al., 2018; Ren et al., 2018; Oldenburg et al., 2013)以及谱Lanczos分解法(Spectral Lanczos Decomposition Method, SLDM)(Druskin et al., 1999).瞬变电磁正演的过程就是给定地电模型通过数值方法求解Maxwell方程组, 获得随时间衰减的电磁响应曲线.尽管求解方程的过程与数学方法各不相同, 但瞬变电磁的基本方程是相同的.

在时域有限差分瞬变电磁正演方面:Oristaglio和Hohmann(1984), 闫述等(2002)采用Du Fort-Frankel有限差分方法研究了二维地电断面平行导线电流源产生的瞬变电磁场, 并分析了均匀半空间中包含二维异常体时的瞬变电磁场分布特征.Wang和Hohmann(1993)采用改进的Du Fort-Frankel方法首次给出了通过求解一阶Maxwell方程组进行电磁探测建模的三维时域有限差分算法.之后, 宋维琪和仝兆歧(2000)针对电偶源瞬变电磁进行了三维有限差分法正演计算.Commer和Newman(2004;2005)针对电性源长偏移瞬变电磁(LOTEM)建模问题建立了求解非因果场的三维时域有限差分并行算法.针对井下瞬变电磁探测的问题, 岳建华和杨海燕等(2008)进行了矿井瞬变电磁探测三维时域有限差分正演研究.孟庆鑫、潘和平(2012)采用三维时域有限差分研究了地井和井中瞬变电磁响应.

自2012年以来, 瞬变电磁时域有限差分三维正演成为国内的研究热点, 许洋铖和林君等(2012)的航空时间域电磁响应三维有限差分数值计算, 邱稚鹏和黄清华等(2013)的起伏地形下瞬变电磁场的三维建模方法, 李建慧和胡祥云等(2013)基于电场Helmholtz方程的回线源瞬变电磁三维正演方法, 孙怀凤和李貅等(2013)考虑关段时间的回线源瞬变电磁三维正演方法, 余翔和王绪本等(2017)的三维时域有限差分的CPML边界以及赵越和李貅等(2017)针对航空电磁开展的复杂模型三维正演研究.之后, 孙怀凤等(2018)还开发了针对瞬变电磁三维时域有限差分网格剖分的多尺度网格方法, 用于解决小目标建模问题并提高正演速度.

然而, 常规的时域有限差分方法为了保证数值计算的稳定, 时间步长Δt需要严格满足Courant-Friedrich-Lewy(CFL)限制条件, 这使得迭代计算过程中时间采样过密、计算时间过长.为了克服上述问题, Fijany等(1995)首次提出采用Crank-Nicolson差分格式的有限差分方法求解麦克斯韦方程.Crank-Nicolson方法是一种隐式差分方法, 具有无条件稳定的特点, 由于Δt不再受CFL的限制, 可以有效的减少迭代次数.然而, 该方法的缺点也非常明显, 即在每次迭代中都需要求解大型稀疏矩阵, 当矩阵阶数较大时会占用大量计算资源和时间, 这制约了该方法的推广和应用.直到近年来, CN-FDTD近似快速算法陆续出现, 才使得CN-FDTD方法得到逐步推广, 主要的近似求解算法有:CNDS-FDTD(Sun et al., 2003), CNAFS-FDTD(Sun et al., 2004), CNCSU-FDTD(Sun et al., 2006)等, 这些方法保留了CN-FDTD方法的无条件稳定性, 求解速度得到了巨大提升, 同时, 计算精度远高于ADI等隐式方法(Garcia et al., 2002).

本研究基于Crank-Nicolson差分方法对Maxwell方程组重新离散, 空间网格仍然采用Yee元胞, 时间步进采用整时间步电场、磁场同时采样, 建立无条件稳定的FDTD格式.在时间采样上, 与常规FDTD交替采样相比, CN-FDTD电场、磁场同时采样, 构成了隐式差分格式, 需要求解稀疏矩阵方程组.采用Crank-Nicolson-Cycle-Sweep-Uniform(CNCSU)近似求解方法, 在保证精度的同时, 计算效率大幅提高, 且内存占用小.

由于计算资源是有限的, 无法模拟电磁波在开域情况下的传播过程, 因此需要将计算模型在某处截断.常规的瞬变电磁数值模拟中, 边界条件多采用Dirichlet边界, 这是因为该边界条件理论简单、容易实现, 但如果想要边界反射对目标区域的扩散场影响足够小, 则需要将边界设置的足够远, 采用更大尺度的模型, 不可避免的会增大计算量.因此, 一些学者开始研究其他的边界条件来代替传统的Dirichlet边界, 如Mur吸收边界(Mur, 1981)、Liao吸收边界(Liao et al., 1984)和完全匹配层边界(Berenger, 1994; Berenger, 1996)等.由于复频率参数完全匹配层(CFS-PML)(Kuzuoglu et al., 1996)对低频波有较好的吸收效果, 李展辉和黄清华(2014)以及余翔和王绪本等(2017)将CFS-PML应用于瞬变电磁FDTD的数值模拟中.本文采用双线性变换方法将CFS-PML边界条件施加于无条件稳定的CN-FDTD方法, 并将算法用于模拟三维瞬变电磁场的传播.

1 CN-FDTD方法 1.1 CN-FDTD时空采样

CN-FDTD是基于克兰克-尼科尔森(Crank-Nicolson)差分格式提出的时域有限差分方法, 该方法是一种既能够满足时间步长放大又能够满足计算精度且无条件稳定的数值计算方法(Fijany et al., 1995; Sun et al., 2003; Yang et al., 2006), 在微波、毫米波以及光学领域等高频电磁波问题中已经得到应用.与常规FDTD一样, CN-FDTD仍然采用Yee元胞作为基本空间离散单元, 如图 1所示, 电场在网格棱边中心采样, 磁场在网格各面的中心采样.然而, 在时间采样上, CN-FDTD与常规FDTD存在较大差别, 图 2给出了CN-FDTD与常规FDTD时间步进采样对比图, 与常规FDTD的电场和磁场在nn+1/2时刻交替采样不同, CN-FDTD的电场和磁场同时在整数时刻采样.

图 1 CN-FDTD使用的Yee网格 Fig. 1 Yee grid
图 2 常规FDTD与CN-FDTD电场、磁场时间采样分布对比 (a)常规FDTD时间采样;(b) CN-FDTD时间采样. Fig. 2 Comparison of original FDTD and CN-FDTD sampling in time for electric and magnetic fields (a) Original FDTD; (b) CN-FDTD.
1.2 控制方程与差分离散

在均匀、各向同性、非磁性、无源媒质中, Maxwell方程组旋度方程可以写成如下分量形式(葛德彪和闫玉波, 2005):

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

其中, HE分别为磁场强度和电场强度, εσμ分别是介电常数、电导率和磁导率, xyz构成直角坐标系.

Crank-Nicolson差分策略在空间域使用中心差分, 以式(1a)为例, 用nn+1时刻的平均值代替式中出现的n+1/2时刻的场值, 可得:

(2a)

(2b)

(2c)

(2d)

将(2a—2d)代入(1a)并化简可得:

(3)

式中, 上角标nn+1分别代表着相邻的两个整数采样时刻, DyDz分别为沿着y轴、z轴方向的一阶中心差分算子, .

类似地, 可以得到其他离散方程, 经过进一步压缩与简化可以概化为

(4a)

(4b)

其中, η(η=x, y, z)为坐标方向, 满足循环移位规律, 即若η=x, 则η-1和η+1分别对应zyDη为沿η轴方向的一阶中心差分算子, 系数a2t/2μ.

将(4b)代入(4a)并代换n+1时刻的磁场分量, 可得:

(5)

式中, D2η为沿着η轴方向的二阶中心差分算子, 例如η=x时, 公式(5)中的二阶差分项可以表示为

(6)

从(5)式可以看出, 等号的左侧为n+1时刻的电场值, 为待求的未知量, 等号右侧全部为n时刻的已知场值, 通过求解3个方程组就能得到待求电场.然而, 等号左侧的三个电场分量相互耦合在一起, 求解每一时间步上的电场分量都需要求解一个大型稀疏矩阵, 这将会消耗大量的计算资源与时间, 计算效率较低(Feng et al., 2018; Sun et al., 2004).

1.3 CNCSU近似及CFS-PML吸收边界

为了提高计算效率, 采用Crank-Nicolson Cycle-Sweep-Uniform(CNCSU)方法人为引入高阶误差项对CN-FDTD方程进行近似, 该过程引入虚拟电场, 首先求解虚拟电场, 然后通过虚拟电场反求真实电场, 其中虚拟场与真实场的求解顺序为:①顺序求解虚拟场, 即Ex*Ey*Ez*; ②反序求解真实电场, 即EzEyEx.在这个循环求解过程中, 原来必须求解的3个大型稀疏矩阵被6个对角占优的三对角矩阵代替, 可以用Thomas算法(Thomas, 1949)快速求解, 从而显著提高求解效率.

事实上, 在CNCSU-FDTD近似公式的推导过程中, 发现迭代公式在CFS-PML介质内和普通介质内具有类似的形式, 因而可以推导更一般的CFS-PML介质内的迭代公式, 将相应的系数替换就可以得到CNCSU-FDTD在普通介质内的求解迭代公式.

为了获得复频率参数下的PML边界条件, 在均匀、各向同性、非磁性无源媒质中频率域麦克斯韦旋度方程写为(Chew et al., 1994)

(7a)

(7b)

其中S为坐标伸缩因子, 在CFS-PML介质中, Roden等(2000)S定义为

(8)

其中, κη为网格延拓因子, σ是PML介质中人为添加的电导率, αη是一个大于零的实数.

通常情况下, FDTD计算会将公式(7)、(8)转换到时间域来求解, 然而, 从频率域通过傅里叶反变换到时间域会出现卷积项, 为避免三维模型中的复杂卷积运算, 可将公式(7)、(8)进行双线性变换(Ramadan and Oztoprak, 2002), 将其转换到Z域, 整理得

(9a)

(9b)

(9c)

(9d)

其中, 系数下标中的e代表求解电场, h代表求解磁场.

考虑Z变换的时移特性, 可将(9)式变换为

(10a)

(10b)

(10c)

(10d)

将(10b)代入(10a)并替换n+1时刻的磁场, 在η=x, y, z坐标方向上可以表示为

(11a)

(11b)

(11c)

公式(11)可以简化表示为:

(12)

其中, , b=a1a2, I为单位矩阵, 为人为引入的高阶误差项,

合并整理并进行因式分解, 可得下述子步骤公式:

(13a)

(13b)

其中, WE*是为方便计算而引入的虚拟场, 无实际物理意义, 将式(14)展开可得顺序求解方程如下:

(14a)

(14b)

(14c)

(14d)

(14e)

(14f)

按照式(14a—14c)的求解顺序, 求出Ex*, Ey*Ez*三个虚拟场值, 再按照(14d—14f)的求解顺序求出n+1时刻的三个电场分量Ez, EyEx, 再将它们代入(10b)可以求出n+1时刻的磁场分量.

上述公式(14)在CFS-PML介质和普通介质CNCSU-FDTD迭代具有统一的格式, 当wη=1, 时, 公式(14)退化为普通介质的CNCSU-FDTD近似迭代方程.尽管CNCSU-FDTD需求解的矩阵数量由原来的3个变为6个, 但新矩阵均为对角占优的三对角矩阵, 可以采用Thomas算法快速求解, 所用的计算资源和时间相比求解3个大型稀疏矩阵要少很多, 在普通PC上能很快完成计算.

2 精度验证与复杂模型对比

为了验证CN-FDTD瞬变电磁三维正演方法的计算精度以及对复杂模型的适用性, 与均匀半空间模型解析解、三层模型线性数字滤波解进行了对比, 之后, 选择三维垂直接触带复杂模型并与三维矢量有限元、三维有限体积以及常规时域有限差分的计算结果进行对比.

2.1 简单模型对比与精度验证

本节计算模型统一使用中心回线装置, 发射回线为100 m×100 m, 激发电流为1 A, 发射波形采用文献(孙怀凤等, 2013)中开关函数处理后的梯形波, 上升沿与下降沿均为1 μs, 电流持续时间为30 ms, 观测时间为30 ms, 占空比为1:1, 相当于实际观测中双极性矩形波条件下瞬变电磁基频为8.33 Hz.如图 3a所示, 模型分为计算区域与PML边界区域, 计算区域采用20 m均匀网格剖分, 网格数为61×61×30, 计算区域模型尺寸为1220 m×1220 m×600 m, 其中, 空气层采用6层网格, 大地采用24层网格;PML边界区域采用15层, 包裹在计算区域之外.

图 3 均匀半空间模型CN-FDTD与解析解对比 (a)模型示意图; (b) CN-FDTD数值解与解析解的感应电动势衰减曲线对比. Fig. 3 Comparison of CN-FDTD numerical solution and analytical solution for Homogeneous half space model (a) Schematic diagram of the model; (b) Decay curve comparison of the CN-FDTD and analytic solutions.

均匀半空间模型大地电阻率为100 Ωm, 空气电阻率设为106 Ωm, 数值模拟结果与解析解(李貅, 2002; Anderson, 1979)对比时间范围为关断后1 μs~30 ms, 图 3b给出了三维CN-FDTD计算结果与阶跃波解析解的感应电动势衰减曲线对比.可以看到, 在关断后的早期二者差异较大, 这是由于三维计算中考虑了关断时间而不是阶跃关断引起的(孙怀凤等, 2013), 而随着关断效应的消失, 误差逐渐减小, 10 μs之后两者吻合较好.

为了进一步测试, 采用典型的三层模型, 并与线性数字滤波解(李貅, 2002; Anderson, 1979)进行对比, 层状模型的层厚和电阻率如表 1所示, 计算时所使用的其他参数与均匀半空间模型一致.图 4给出了四种层状模型的感应电动势衰减曲线的对比结果.从图中可以看出, 层状模型的对比结果与均匀半空间模型类似, 均为早期误差较大, 这是由于关断时间的影响, 随着时间的推移, 关断时间的影响逐渐减弱, 误差逐渐减小, 在晚期, 数值解与解析解曲线吻合较好, 二者几乎重合.

表 1 层状模型参数 Table 1 Parameters of the layered models
图 4 典型三层模型CN-FDTD数值解与线性数字滤波解对比 (a) A型;(b) H型;(c) K型;(d) Q型. Fig. 4 Comparison of digital filter solution and CN-FDTD for typical there-layered models (a) A type Model; (b) H type Model; (c) K type Model; (d) Q type Model.
2.2 复杂模型对比与精度验证

采用三维垂直接触带复杂模型(Commer et al., 2004)进行进一步验证.模型参数如图 5所示, 地下介质可以分为四部分, 表层厚50 m, 电阻率为10 Ωm, 左侧与右侧部分的电阻率分别为100 Ωm和300 Ωm, 二者界面上有一个形状复杂的接触带, 其电阻率为1 Ωm, 接触带长400 m.发射回线边长为100 m, 设置四个接收点, 以最左侧的源边为相对坐标起点, 四个接收点坐标分别为(0, 50, 0)、(0, 150, 0)、(0, 450, 0)和(0, 1050, 0).模型采用均匀网格进行剖分, 网格大小为25 m, 计算区域尺寸为1500 m×1500 m×750 m(网格数60×60×30), 其中, 空气层层数为6, 大地层数为24, CFS-PML介质层数为15.

图 5 三维垂直接触带复杂模型(Li et al., 2017) Fig. 5 Three-dimensional complex model with a vertical contact zone

图 6给出了CN-FDTD计算的感应电动势绝对值随时间的变化曲线, 在曲线的尖点处数据符号发生改变, 在图中给出了标示, 并同时给出了和矢量有限元法(Li et al., 2017)、有限体积法(Zhou et al., 2018)、以及多尺度网格FDTD(孙怀凤等, 2018)方法的对比结果, 发现仅在感应电动势变号处有明显的差异, 其他区域吻合较好.

图 6 三维垂直接触带复杂模型计算结果对比 Fig. 6 The numerical results of three-dimensional complex model
3 CN-FDTD时间步放大对计算精度的影响

CN-FDTD算法的时间步长Δt虽然不再受CFL条件的限制, 但当Δt过大时, CN-FDTD的计算精度会受到较大影响, 导致晚期误差增大.以三维垂直接触带模型为例, 设计不同大小的Δt进行模拟, 研究满足精度条件下的最大时间步.为了对Δt的大小进行量化, 定义CFLN为CN-FDTD最大时间步与常规FDTD最大时间步的比值, 即CFLN=ΔtmaxCN-FDTDtCFLFDTD, 其中ΔtCFLFDTD是常规FDTD受CFL限制所能达到的最大Δt.图 6中与矢量有限元、有限体积法和常规有限差分法进行对比的模型使用的CFLN为200, 后续研究以此为参考结果.图 7给出了CFLN分别为400、800、1600以及3200时, 感应电动势的衰减曲线对比, 从图中可以看出, 当CFLN≤1600时的结果均与CFLN=200时结果吻合较好, 即使CFLN≤1600时, CN-FDTD算法依然有着较高的稳定性, 但当CFLN=3200时, 在晚期出现明显偏差.如图 7所示, 使用相同大小的符号, 当CFLN≤1600时, 所有的符号几乎重合, 而CFLN=3200时, 计算结果已经出现偏离, 尽管这一偏离非常小, 但在双对数坐标下晚期瞬变电磁响应的影响还是非常大的.

图 7 CN-FDTD在不同CFLNs下结果的对比 Fig. 7 Comparison of CN-FDTD numerical results under different CFLNs

对比试验在一台CPU为Intel Core i5-7300 HQ的笔记本电脑上完成, 且所有的对比试验均使用单线程计算, 未启用任何并行策略.表 2给出了不同CFLN时的计算时间和内存开销, 可以看出, 随着CFLN的增加, 迭代次数和计算时间总体减小, 但由于模拟过程考虑了发射波形, 激励源采用梯形波, 上升沿以及下降沿非常陡峭, 因此在计算过程中, 为了保证结果稳定性, 上述阶段的Δt采用逐渐增大和逐渐减小的操作(见图 8).经过测试, 在梯形波持续时间的前半段以及关断后Δt的增长速率都各有一个阈值, 不能增长过快, 如果超过这个阈值, 模拟结果与解析解会产生较大的误差, 这使得即使CFLN成倍数增长, 而迭代次数以及计算时间并不是同比例的减小.从图 8可以看出, 由于关断后Δt增长较慢, 当CFLN≥800, 计算结束时, Δt还没有增长到预设的大小, 因此, 对整体迭代步数贡献主要集中在on-time阶段.如果不考虑关断时间, 而是采用(Wang and Hohmann, 1993)给出的设定初始场的方法施加激励源, 则计算时间会大幅减小.表 2显示CFLN=1600时, 计算时间减少到了50 min, 后续如果考虑FDTD的天然并行性, 上述模型的计算时间有望减少到3 min以内, 这将为后续的三维反演工作带来便利条件.而且, 针对所设计的模型所需的内存消耗仅为192 Mb, 在普通PC就可以轻松完成.即使使用较多的网格模拟非常精细的地下目标时, 也能够保持较小的内存开销.

表 2 不同CFLNs下CN-FDTD算法内存开销及计算时间 Table 2 Memory overhead and computing time of CN-FDTD algorithm under different CFLNs
图 8 CN-FDTD在不同CFLNs下Δt的变化对比 Fig. 8 Changes of Δt in the CN-FDTD under different CFLNs
4 结论

常规FDTD进行瞬变电磁三维正演时, 迭代时间步长受CFL稳定性条件限制, 造成迭代次数过多、计算时间过长、晚期累积误差增大.本研究提出CN-FDTD隐式差分方法用于瞬变电磁三维正演, 突破了CFL条件的限制, 显著减少迭代步数, 通过CNCSU-FDTD近似算法实现方程组快速求解, 通过双线性变换方法施加CFS-PML吸收边界条件, 降低计算时间成本.通过与均匀半空间解析解、层状模型的线性数字滤波解, 以及三维垂直接触带复杂模型的矢量有限元、有限体积法、常规FDTD解进行了对比与精度验证, 结果均吻合较好.本研究给出的方法内存消耗非常小, 可以使用大量网格模拟非常复杂的模型, 计算过程可以在普通PC上快速完成.

致谢  研究在孙怀凤等(2013)开发的tem3dfdtd程序基础上完成.作者在研究过程中与长安大学李貅教授团队鲁凯亮博士进行了深入讨论, 在此一并表示感谢.
References
Adhidjaja J, Hohmann G. 1989. A finite-difference algorithm for the transient electromagnetic response of a 3-dimensional body. Geophysical Journal International, 98(2): 233-242.
Anderson W L. 1979. Numerical-integration of related hankel transforms of orders 0 and 1 by adaptive digital filtering. Geophysics, 44(7): 1287-1305.
Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic-waves. Journal of Computational Physics, 114(2): 185-200.
Berenger J P. 1996. Perfectly matched layer for the fdtd solution of wave-structure interaction problems. IEEE Transactions on Antennas and Propagation, 44(1): 110-117.
Chang J H, Xue G Q, Malekian R. 2019. A comparison of surface-to-coal mine roadway tem and surface tem responses to water-enriched bodies. IEEE Access, 7: 167320-167328.
Chang J H, Yu J C, Li J J, et al. 2019. Diffusion law of whole-space transient electromagnetic field generated by the underground magnetic source and its application. IEEE Access, 7: 63415-63425.
Chew W C, Weedon W H. 1994. A 3D perfectly matched medium from modified maxwells equations with stretched coordinates. Microwave and Optical Technology Letters, 7(13): 599-604.
Commer M, Newman G. 2004. A parallel finite-difference approach for 3d transient electromagnetic modeling with galvanic sources. Geophysics, 69(5): 1192-1202.
Cox L H, Wilson G A, Zhdanov M S. 2010. 3D inversion of airborne electromagnetic data using a moving footprint. Exploration Geophysics, 41(4): 250-259.
Druskin V, Knizhnerman L, Lee P. 1999. New spectral lanczos decomposition method for induction modeling in arbitrary 3-d geometry. Geophysics, 64(3): 701-706.
Feng N X, Zhang Y X, Sun Q T, et al. 2018. An accurate 3-D cfs-pml based crank-nicolson fdtd method and its applications in low-frequency subsurface sensing. IEEE Transactions on Antennas and Propagation, 66(6): 2967-2975.
Fijany A, Jensen M A, Rahmatsamii Y, et al. 1995. A massively parallel computation strategy for fdtd:time and space parallelism applied to electromagnetics problems. IEEE Transactions on Antennas and Propagation, 43(12): 1441-1449.
Garcia S G, Lee T W, Hagness S C. 2002. On the accuracy of the adi-fdtd method. IEEE Antennas and Wireless Propagation Letters, 1: 31-34.
Ge D B, Yan Y B. 2005. Finite Difference Time Domain Method of Electromagnetic Wave (in Chinese). Xi' an: Xidian University Press.
Kuzuoglu M, Mittra R. 1996. Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers. IEEE Microwave and Guided Wave Letters, 6(12): 447-449.
Li J H, Farquharson C G, Hu X Y. 2017. 3D vector finite-element electromagnetic forward modeling for large loop sources using a total-field algorithm and unstructured tetrahedral grids. Geophysics, 82(1): E1-E16.
Li J H, Zhu Z Q, Liu S C, et al. 2011. 3D numerical simulation for the transient electromagnetic field excited by the central loop based on the vector finite-element method. Journal of Geophysics and Engineering, 8(4): 560-567.
Li J H, Hu X Y, Zeng S H, et al. 2013. Three-dimensional forward calculation for loop source transient electromagnetic method based on electric field Helmholtz equation. Chinese Journal of Geophysics (in Chinese), 56(12): 4256-4267.
Li X. 2002. Theory and Application of Transient Electromagnetic Sounding (in Chinese). Xi'an: Shanxi Science and Technology Press.
Li S C, Sun H F, Lu X S, et al. 2014. Three-dimensional modeling of transient electromagnetic responses of water-bearing structures in front of a tunnel face. Journal of Environmental and Engineering Geophysics, 19(1): 13-32.
Liao Z, Huang K, Yang B, et al. 1984. A transmitting boundary for transient wave analyses. Science in China Series A-Mathematics, Physics, Astronomy & Technological Science, 27(0253-5831): 1063.
Li Z H, Huang Q H. 2014. Application of the complex frequency shifted perfectly matched layer absorbing boundary conditions in transient electromagnetic method modeling. Chinese Journal of Geophysics (in Chinese), 57(4): 1292-1299.
Meng Q X, Pan H P. Numerical simulation analysis of surface-hole TEM responses. Chinese Journal of Geophysics (in Chinese), 55(3): 1046-1053.
Mur G. 1981. Absorbing boundary-conditions for the finite-difference approximation of the time-domain electromagnetic-field equations. IEEE Transactions on Electromagnetic Compatibility, 23(4): 377-382.
Newman G, Commer M. 2005. New advances in three dimensional transient electromagnetic inversion. Geophysical Journal International, 160(1): 5-32.
Oldenburg D W, Haber E, Shekhtman R. 2013. Three dimensional inversion of multisource time domain electromagnetic data. Geophysics, 78(1): E47-E57.
Oristaglio M, Hohmann G. 1984. Diffusion of electromagnetic-fields into a two-dimensional earth-a finite-difference approach. Geophysics, 49(7): 870-894.
Qiu Z P, Li Z H, Li D Z, et al. 2013. Non-orthogonal-Grid-based three dimensional modeling of transient electromagnetic field with topography. Chinese Journal of Geophysics (in Chinese), 56(12): 4245-4255.
Ramadan O, Oztoprak A Y. 2002. DSP techniques for implementation of perfectly matched layer for truncating FDTD domains. Electronics Letters, 38(5): 211-212.
Ren X Y, Yin C C, Macnae J, et al. 2018. 3D time-domain airborne electromagnetic inversion based on secondary field finite-volume method. Geophysics, 83(4): 219-228.
Roden J A, Gedney S D. 2000. Convolution PML (CPML):an efficient FDTD implementation of the cfs-pml for arbitrary media. Microwave and Optical Technology Letters, 27(5): 334-339.
Song W Q, Tong Z Q. 2000. Forward finite differential calculation for 3-D transient electromagnetic field. Oil Geophysical Prospecting, 35(6): 751-756.
Sugeng F. 1998. Modeling the 3D TDEM response using the 3D full-domain finite-element method based on the hexahedral edge-element technique. Exploration Geophysics, 29(4): 615-619.
Sun G L, Trueman C W. 2003. Unconditionally stable crank-nicolson scheme for solving two-dimensional maxwell's equations. Electronics Letters, 39(7): 595-597.
Sun G L, Trueman C W. 2004. Unconditionally-stable fdtd method based on crank-nicolson scheme for solving three-dimensional maxwell equations. Electronics Letters, 40(10): 589-590.
Sun G L, Trueman C W. 2006. Efficient implementations of the crank-nicolson scheme for the finite-difference time-domain method. IEEE Transactions on Microwave Theory and Techniques, 54(5): 2275-2284.
Sun H F, Cheng M, Wu Q L, et al. 2018. A multi-scale grid scheme in three-dimensional transient electromagnetic modeling using FDTD. Chinese Journal of Geophysics (in Chinese), 61(12): 5096-5104.
Sun H F, Li X, Li S C, et al. 2013. Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time. Chinese Journal of Geophysics (in Chinese), 56(3): 1049-1064.
Thomas L H. 1949. Elliptic problems in linear differential equations over a network. Watson scientific computing laboratory Report, Columbia University, New York.
Wang T, Hohmann G W. 1993. A finite-difference, time-domain solution for three-dimensional electromagnetic modeling. Geophysics, 58(6): 797-809.
Xiong Z H. 1992. Electromagnetic modeling of 3-D structures by the method of system iteration using integral-equations. Geophysics, 57(12): 1556-1561.
Xu Y C, Lin J, Li S Y, et al. 2012. Calculation of full-waveform airborne electromagnetic response with three-dimension finite difference solution in time-domain. Chinese Journal of Geophysics (in Chinese), 55(6): 2105-2114.
Xue G Q, Li X. 2008. The technology of TEM tunnel prediction imaging. Chinese Journal of Geophysics (in Chinese), 51(3): 894-900.
Xue G Q, Li X, Di Q Y. 2007. The progress of TEM in theory and application. Progress in geophysics (in Chinese), 22(4): 1195-1200.
Yang H Y, Yue J H. 2008. Response characteristics of the 3D whole space TEM disturbed by roadway. Journal of Jilin University (Earth Science Edition), 38(1): 129-134.
Yang Y, Chen R S, Yun E. 2006. The unconditionally stable crank-nicolson fdtd method for three-dimensional maxwell's equations. Microwave and Optical Technology Letters, 48(8): 1619-1622.
Yan S, Chen M S, Fu J M. 2002. Direct time-domain numerical analysis of transient electromagnetic fields. Chinese Journal of Geophysics (in Chinese), 45(2): 275-284.
Yu x, Wang X B, Li X J, et al. 2017. Three-dimensional finite difference forward modeling of the transient electromagnetic method in the time domain. Chinese Journal of Geophysics (in Chinese), 60(2): 810-819.
Zhao Y, Li X, Wang Y P, et al. 2017. Characteristics of terrain effect for 3-D ATEM. Chinese Journal of Geophysics (in Chinese), 60(1): 383-402.
Zhdanov M S, Lee S K, Yoshioka K. 2006. Integral equation method for 3d modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity. Geophysics, 71(6): G333-G345.
Zhou J M, Liu W T, Li X, et al. 2018. Research on the 3d mimeticfinite volume method for loop-source tem response in biaxial anisotropic formation. Chinese Journal of Geophysics (in Chinese), 61(1): 368-378.
葛德彪, 闫玉波. 2005. 电磁波时域有限差分方法. 西安: 西安电子科技大学出版社.
李建慧, 胡祥云, 曾思红, 等. 2013. 基于电场Helmholtz方程的回线源瞬变电磁法三维正演. 地球物理学报, 56(12): 4256-4267.
李貅. 2002. 瞬变电磁测深的理论与应用. 西安: 陕西科技技术出版社.
李展辉, 黄清华. 2014. 复频率参数完全匹配层吸收边界在瞬变电磁法正演中的应用. 地球物理学报, 57(04): 1292-1299.
孟庆鑫, 潘和平. 2012. 地-井瞬变电磁响应特征数值模拟分析. 地球物理学报, 55(3): 1046-1053.
邱稚鹏, 李展辉, 李墩柱, 等. 2013. 基于非正交网格的带地形三维瞬变电磁场模拟. 地球物理学报, 56(12): 4245-4255.
宋维琪, 仝兆歧. 2000. 3D瞬变电磁场的有限差分正演计算. 石油地球物理勘探, 35(6): 751-756.
孙怀凤, 程铭, 吴启龙, 等. 2018. 瞬变电磁三维FDTD正演多分辨网格方法. 地球物理学报, 61(12): 5096-5104.
孙怀凤, 李貅, 李术才, 等. 2013. 考虑关断时间的回线源激发TEM三维时域有限差分正演. 地球物理学报, 56(3): 1049-1064.
许洋铖, 林君, 李肃义, 等. 2012. 全波形时间域航空电磁响应三维有限差分数值计算. 地球物理学报, 55(6): 2105-2114.
薛国强, 李貅. 2008. 瞬变电磁隧道超前预报成像技术. 地球物理学报, 51(3): 894-900.
薛国强, 李貅, 底青云. 2007. 瞬变电磁法理论与应用研究进展. 地球物理学进展, 22(4): 1195-1200.
闫述, 陈明生, 傅君眉. 2002. 瞬变电磁场的直接时域数值分析. 地球物理学报, 45(2): 275-284.
杨海燕, 岳建华. 2008. 巷道影响下三维全空间瞬变电磁法响应特征. 吉林大学学报(地球科学版), 38(1): 129-134.
余翔, 王绪本, 李新均, 等. 2017. 时域瞬变电磁法三维有限差分正演技术研究. 地球物理学报, 60(2): 810-819.
赵越, 李貅, 王祎鹏, 等. 2017. 三维起伏地形条件下航空瞬变电磁响应特征研究. 地球物理学报, 60(1): 383-402.