瞬变电磁(Transient Electromagnetic, TEM)或时间域电磁(Time Domain Electromagnetic, TDEM)方法基于探测目标的电阻率差异, 依靠电磁感应(Electromagnetic Induction, EMI)信号实现探测，是一种低成本、高效率、施工简单、精度较高的地球物理探测手段，已经在金属矿产资源勘探、海洋天然气水合物探测、环境与水文地质调查、地下水评估与监测、矿山采空区及灾害水源探测、地下未爆弹药(Unexploded Ordnance, UXO)探查、地质填图等应用领域发挥了重要的作用(Zhdanov, 2010; 朴化荣, 1990; Buselli et al., 1992).随着经济社会的快速发展，越来越多的勘探需要在复杂地质条件下进行，这就需要更加合理的三维模拟方法.瞬变电磁勘探时，基于概化模型的正演响应规律与特征分析贯穿整个方案的设计、可行性分析、数据采集方式、资料反演与解释，这一工作主要依靠数值正演来实现.对于复杂地质条件下的瞬变电磁勘探，三维正演模拟显得尤为重要(孙怀凤等, 2013; 李貅, 2002).

复杂模型的瞬变电磁三维正演计算主要依靠数值方法完成，常用的有有限体积法(Finite Volume, FV)、体积分方程法(Volume Integral Equation, VIE)、谱Lanczos分解法(Spectral Lanczos Decomposition Method, SLDM)、有限单元法(Finite Element Method, FEM)、以及时域有限差分法(Finite Difference Time Domain, FDTD)，并且主要的工作首先由国外学者完成，我国的瞬变电磁三维正演工作很早就有人注意到并着手开始研究，但真正趋于实用和模型复杂化是近年才开始的(李建慧等, 2013; 薛国强等, 2008).自2012年至今，这一方面的研究在国内外均属于热点(余翔等, 2017; 孙怀凤等, 2013; 李建慧等, 2013; 许洋铖等, 2012; 邱稚鹏等, 2013; 周建美等，2018; Li et al., 2017; Li and Huang, 2015; Cai et al., 2014).Commer和Newman针对长偏移距瞬变电磁的三维正反演问题展开研究，开发了并行三维正反演程序(Commer and Newman, 2004, 2006; Newman and Commer, 2005)；Oristaglio和Hohmann(1984)、闫述等(2002)使用DuFort-Frankel有限差分法研究了二维地电模型的电磁场响应规律；Wang和Hohmann (1993)通过引入虚拟位移电流项构建显示差分格式解决了瞬变电磁FDTD三维正演问题；孙怀凤(2013)使用非均匀网格有限差分正演模拟了隧道的含水构造.

FDTD方法直接在时间域求解Maxwell方程组来获得瞬变电磁响应(葛德彪和闫玉波, 2005)，不经过余弦变换或Gaver-Stohfest变换之类的频率-时间转换问题，能够准确反映三维地质结构在瞬变电磁场激励下的响应规律.然而，时域有限差分采用Yee结构化网格结构，必须采用正六面体的网格形式进行剖分，即使发展了非均匀网格剖分方案，但由于差分方程需要使用计算节点在三个自由度上与相邻采样点距离作为计算参数，长宽比过大的网格易引入较大的计算误差，如图 1所示，长宽比过大易导致产生“扁盒子”状的网格，如果要求解E_z，则需要其周围的两个H_x和H_y，而两个H_x之间的距离远小于两个H_y之间的距离，对E_z的求解引入较大误差.一般情况下，时域有限差分的网格要求相邻网格尺寸比例小于1.2 : 1，最大与最小网格尺寸的比例限制在20 : 1(余文华等, 2005).瞬变电磁三维FDTD正演网格剖分的最主要矛盾是最小网格尺寸问题，如果目标较小，就必须使用小的网格来建立异常体，而最小网格尺寸又与时间采样步长密切相关，当最小网格尺寸设置较小时，为了保证拥有足够的模型计算区域，必须使用较多的网格数量来进行剖分建模，同时时间网格也会受到最小网格尺寸的影响；当采用的最小网格尺寸较大时，能够使用较少的网格数目剖分和较大的时间步长，但对于异常目标的细节描述就会大大降低(王长清和祝西里, 1994).

图 1

Fig. 1

图 1 FDTD网格长宽比示意图 Fig. 1 The aspect ratio of FDTD grid

针对上述问题，基于非均匀网格剖分的研究基础，引入多尺度网格剖分方法，即首先使用较大尺寸的网格进行第一次剖分(称为“粗剖分”)，然后在希望加密的区域进行二次剖分(称为“细剖分”)，在模型计算过程中，至少同时存在粗-细两套网格，其中：粗网格用来建立背景导电率模型或描述构造级别的较大目标区块，细网格用于建立关心区域的细节目标模型，例如起伏地形变化、金属矿脉、套管井壁等使用大网格时不能很好地描述，实现粗网格-细网格混合使用.尽管细网格包含在粗网格内部，但其具有Yee网格的全部属性，因而可以在网格中设置不同的电性参数模拟不同形状的目标.上述网格改进产生了3个需要解决的问题：①电磁场分别在粗网格与细网格区域的计算迭代问题：由于细网格的时间步长与粗网格存在差异，因而需要控制方程必须很好地描述细网格内部电磁场的扩散传播过程.②粗细网格的边界条件问题：电磁场的传播过程中电磁场既需要由粗网格传入细网格，也需要由细网格传入粗网格，因而在计算区域产生了新的内部边界条件.③粗细网格电阻率赋值问题：由于粗网格和细网格分别需要进行各自的迭代计算，细网格区域的电阻率是根据模拟物体的形态设置的，那么粗网格与细网格的电阻率取值必然相关，以适应粗细网格的电磁场变化.

本研究基于Maxwell方程组推导了细网格内电场和磁场的迭代公式，基于泰勒展开推导了粗-细网格的边界条件，使粗、细网格在时间域和空间域得到统一.首先，采用均匀半空间中包含三维低阻目标的经典模型进行精度验证，计算结果分别与积分方程法、Wang的FDTD方法以及作者原先开发的FDTD方法进行了对比验证；其次，建立三维垂直接触带复杂模型，将多分辨网格的计算结果与三维矢量有限元和有限体积法对比；之后，使用“L”型异常模型进行计算效率对比，发现多分辨网格算法能在满足精度要求的前提下显著提高运算效率.

本文研究的多分辨网格方法是指首先使用粗网格剖分，建立较大区域和较少网格数目的FDTD三维正演模型，然后在模型需要关注的部位进行局部细化剖分，形成细网格，使整个模型计算区域中至少包含两套网格用于描述尺寸差异较大目标的方法.如图 2所示，该方法能够在细化区域多次剖分，形成多套计算网格，通过大网格描述地质背景和较大体积的目标，通过小网格描述相对细部地质结构(例如较大的地形起伏、狭长型矿脉、钢筋干扰源等).这样处理能够最大限度的减小内存和计算时间，提高计算效率.对于单个的Yee晶胞单元而言，不管是粗网格还是细网格，都符合标准FDTD的Yee晶胞格式.

图 2

Fig. 2

图 2 Yee晶胞与多分辨网格Yee晶胞 Fig. 2 Yee cell and multi-resolution grid Yee cell

在均匀、有耗、无磁性、无源的各向同性介质中，Maxwell方程组为(葛德彪和闫玉波, 2005)：

(1)

式中：E是电场强度，B是磁通量密度，H是磁场强度，ε是介电常数，σ表示电导率.

在Maxwell方程组中对位移电流进行特殊处理，引入虚拟位移电流，并结合各向同性介质中的本构关系，将Maxwell方程组变换到直角坐标系(孙怀凤, 2013).直角坐标系下的瞬变电磁场方程的有限差分离散使用Yee晶胞格式，对部分网格进行局部细化构成多分辨网格.

对于均匀部分和中心细化网格的剖分，电磁场六个分量的迭代公式已有详细的解答(孙怀凤等, 2013)，而多分辨网格只是在空间和时间上对其细化，在细网格内部电磁场的值需要重新定义，需将电场分量、磁场分量、电导率和迭代时间等参数替换为细网格参数.对一个粗网格进行n_f次剖分，可以获得内部细网格区域电场分量的表达式如下：

(2)

(3)

(4)

(5)

(6)

(7)

原有FDTD三维正演的边界条件已经在文献(孙怀凤等, 2013)中充分讨论，然而，采用多分辨网格方法后，增加了一类新的边界条件，即粗-细网格边界条件.这一边界条件决定了电磁场迭代计算过程中粗网格数据与细网格数据的相互交换，对多分辨网格而言尤为重要.

如图 3所示，以二维TM网格为例，图中左侧网格为粗网格(以实线表示)，右侧为细网格(以虚线表示).图中黑色圆点代表电场在粗网格的空间采样，红色圆点表示电场在细网格的空间采样，绿色圆点表示电场采样点既在粗网上又在细网格上；而黑色箭头代表粗网格上的磁场采样，红色箭头表示细网格上的磁场采样.以电场采样为例，对于重合的(绿色)采样点，在边界条件处理上按照电磁场扩散传播的时间顺序迭代和更新节点值即可，对于仅在细网格采样而不在粗网格采样的(红色)节点，则需要单独处理.图 3中给出细网格采样节点E₁₁以及周边的粗网格采样点E₁、E₂、E₃、E₄和E₅，根据泰勒展开公式，可以获得E₁₁节点值为：

(8)

图 3

Fig. 3

图 3 多分辨网格的边界处理 Fig. 3 Boundary treatment of multi-resolution grid

其中，n_t代表细网格从n时刻到n+1时刻迭代分为n_t个小的时间步.

在空间域，由于整个计算区域存在至少两套网格，需首先对粗网格进行划分，然后对选定的粗网格进行二次剖分建立细网格区域.由于迭代过程具有一定的独立性，所以粗网格和细网格需要分别独立编号.如图 4所示，空间域的计算可以分为3部分，即粗网格计算部分、内部边界条件交换部分和细网格计算部分.

图 4

Fig. 4

图 4 时间域与空间域统一流程图 Fig. 4 Unified flow chart of time-domain and space-domain

在时间域，粗网格按照自己的迭代要求自n时刻计算到n+1时刻，然后将n时刻的电磁场值通过内部边界条件处理传递给细网格作为其外部边界条件，细网格按照粗网格类似的做法完成内部迭代过程到达n+1时刻，细网格再将边界处的电磁场值传递处理替换粗网格的值，完成一次的时间与空间统一.详细过程如图 4所示.

对于多分辨网格方法，细网格剖分形成的区域仍然符合FDTD的通用格式，因而可以设置不同的电阻率参数，如图 5所示的粗细网格示意图，左侧为一个粗网格，右侧为细网格区域，将粗网格边分为多份，图中每种灰度色块代表一个电阻率取值，可以看出，细网格部分的电阻率取值也是可以任意设置的，这样的好处是我们可以在粗网格内部剖分成的细网格设置电阻率异常体，形成更加细化的边界形态.

图 5

Fig. 5

图 5 粗-细网格及电阻率取值示意图 Fig. 5 Diagram of coarse and fine grid and the resistivity distribution

对于右侧的细网格部分，其所在的粗网格单元仍然需要参与粗网格迭代计算，因而也需要赋以合理的电阻值，本文采用细网格体平均电阻率值的方法计算获得粗网格电阻率值.

为验证提出的多分辨网格算法及程序的可靠性，选用Newman等(1986)建立的均匀半空间包含低阻三维异常经典模型，分别与FDTD和积分方程法进行对比.多分辨网格具体的剖分方式如图 6所示，传统FDTD方法剖分网格采用最小尺寸为10 m的网格，多分辨网格方法首先使用30 m的最小网格尺寸剖分，然后在异常区域细化剖分，如图 6所示给出的局部网格剖分示意，其中灰色区域代表模型中的低阻体.

图 6

Fig. 6

图 6 多分辨网格方法剖分模型局部示意图 Fig. 6 Localization diagram of multi-resolution grid method

如图 7所示，对比四种算法结果可以发现，多分辨网格方法与积分方程、传统FDTD方法的结果基本一致.事实上，本方法所使用的算法和程序依据作者早期开发的计算方法，并已经在文献(孙怀凤等, 2013)中证明了其算法有效性，增加多分辨网格的方法对原方法影响较小，因此，多分辨网格算法满足精度要求.

图 7

Fig. 7

图 7 多分辨网格算法精度对比图 Fig. 7 Computing precision comparison chart of multi-resolution grid method

本文采用Li等(2017)的三维垂直接触带的复杂模型进行计算精度的验证，该模型最早由Commer和Newman(2004)提出.具体模型如图 8所示，使用100 m×100 m回线源发射，四个接收点的位置分别为(0, 50, 0)、(0, 150, 0)、(0, 450, 0)、(0, 1050, 0).将本文使用多分辨网格剖分的计算结果与三维矢量有限元和三维有限体积法进行比较，结果如图 8所示.对比三种算法结果可以发现，多分辨网格方法和三维矢量有限元法、三维有限体积法结果基本一致，在感应电压变号时相对误差较大.

图 8

Fig. 8

图 8 三维模型及计算结果对比(Li et al., 2017; Zhou et al., 2018) Fig. 8 3D model and the comparison results (Li et al., 2017; Zhou et al., 2018)

多分辨网格方法能够针对模型中包含小异常目标的情况实现较好的网格剖分和优化，而三维正演的计算时间主要取决于网格数量的多少和最小网格的尺寸.为了对比多分辨网格的计算效率，建立隧道超前探测模型，使用3 m×3 m的小发射线圈，设计的模型如图 9所示.在隧道掌子面前方存在电阻率为1 Ωm的“L”型低阻异常，围岩背景电阻率为1000 Ωm.如表 1所示，模型1采用传统FDTD网格剖分，模型2和模型3首先采用粗网格剖分，然后再在粗网格内部进行细化剖分形成“L”型低阻异常，图 10给出了这3组模型的细部示意，其中灰色区域代表低阻异常体.

图 9

Fig. 9

图 9 模型主视图(a)和俯视图(b) Fig. 9 Front view (a) and top view (b) of the model

图 10

Fig. 10

图 10 三种网格剖分示意图 Fig. 10 Diagram of three meshing grid

表 1 Table 1

表 1 三组模型参数对比 Table 1 Comparison of three sets of model parameters 序号 背景电阻率/Ωm 异常体电阻率/Ωm 粗网格数量 粗网格最小尺寸/m 细网格尺寸/m 1 300×300×300 10 无 2 1000 1 150×150×150 20 10 3 150×150×150 40 10

表 1 三组模型参数对比 Table 1 Comparison of three sets of model parameters

从图 11a可以看出三种类型的衰减电压基本一致，表明多分辨网格对于异常体的剖分满足精度要求.图 11b显示三种类型计算效率的对比，运算时间随着细网格数的增加而变短.综上，多分辨网格剖分异常体不仅能满足瞬变电磁的精度要求，且能极大提高运算效率.

图 11

Fig. 11

图 11 多分辨网格不同剖分个数计算结果对比图 (a)计算精度对比; (b)计算效率对比. Fig. 11 Calculation results comparison figure of different multi-grid number (a) Calculation accuracy comparison; (b) Calculation efficiency comparison.

本文在已有的瞬变电磁三维FDTD正演基础上，提出了多分辨网格的方法，并给出了内部边界条件的处理方法和粗-细网格的时间-空间迭代统一策略.采用多分辨网格方法能够保证计算精度的前提下显著减少计算时间.

通过本文算例以及进行瞬变电磁时域有限差分三维正演的经验可知，限制计算时间的因素主要是显示格式对时间步长的严格要求.采用多分辨网格后的直接益处是：①减少了粗网格的数量；②增大了粗网格中最小网格的尺寸.一般情况下，粗网格的数量远大于细网格的数量.事实上，不管是传统方法还是多分辨网格方法，单次时域有限差分迭代计算时间与网格数量密切相关，这是由迭代算法决定的，降低了粗网格的数量自然会大幅降低单次迭代的时间，进而降低总体计算时间.另外，增大了最小网格尺寸，无形中相当于放宽了时间条件的限制，并且时间步长的放大倍数与最小网格尺寸的放大倍数成正比，例如，将最小网格尺寸由10 m放大到40 m，则相应的时间步长也会放大4倍，相当于大幅减少了整体迭代计算次数，对于降低整体计算时间非常显著.

说明  作者在文献(孙怀凤等, 2013)开发的瞬变电磁三维时域有限差分正演计算程序(TEM3dFDTD)代码已经向学术和科研机构免费开放，如有需要请与作者联系.本文开发的多分辨网格方法计算代码也将在论文投稿完成和代码审核之后更新到TEM3dFDTD之中，供研究者参考使用.