瞬变电磁Crank-Nicolson FDTD三维正演

孙怀凤， 柳尚斌， 杨洋

1 山东大学 岩土与结构工程研究中心， 济南 250061 2 山东大学 地球电磁探测研究所， 济南 250061

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使用的Yee网格Fig.1 Yee grid

1 CN-FDTD方法

1.1 CN-FDTD时空采样

图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.

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的电场和磁场在n和n+1/2时刻交替采样不同，CN-FDTD的电场和磁场同时在整数时刻采样.

1.2 控制方程与差分离散

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

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

其中，H和E分别为磁场强度和电场强度，ε、σ和μ分别是介电常数、电导率和磁导率，x、y、z构成直角坐标系.

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

(2a)

(2b)

(2c)

(2d)

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

(3)

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

(4a)

(4b)

其中，η(η=x,y,z)为坐标方向，满足循环移位规律，即若η=x, 则η-1和η+1分别对应z和y；Dη为沿η轴方向的一阶中心差分算子，系数a2=Δt/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吸收边界

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

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

(7a)

(7b)

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

(8)

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

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

(9a)

(9b)

(9c)

(9d)

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

(10a)

(10b)

(10c)

(10d)

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

(11a)

(11b)

(11c)

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

(12)

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

(13a)

(13b)

(14a)

(14b)

(14c)

(14d)

(14e)

(14f)

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层，包裹在计算区域之外.

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

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

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.

表1 层状模型参数Table 1 Parameters of the layered models

图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.

图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.

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

图6 三维垂直接触带复杂模型计算结果对比Fig.6 The numerical results of three-dimensional complex model

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

3 CN-FDTD时间步放大对计算精度的影响

图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程序基础上完成.作者在研究过程中与长安大学李貅教授团队鲁凯亮博士进行了深入讨论，在此一并表示感谢.