地球物理学报  2013, Vol. 56 Issue (12): 4256-4267   PDF    
基于电场Helmholtz方程的回线源瞬变电磁法三维正演
李建慧1 , 胡祥云1 , 曾思红1 , 路金阁1 , 霍光谱2 , 韩波1 , 彭荣华1     
1. 中国地质大学地球物理与空间信息学院, 武汉 430074;
2. 河南省地质调查院, 郑州 450001
摘要: 正演是电磁法勘探野外工作参数选取、室内资料处理与解释的基础,精确、稳定、高效的三维正演算法尤为重要.本文采取先求解拉普拉斯域电场、再由Gaver-Stehfest算法获得时间域磁场的思路,基于电场异常场Helmholtz方程实现了交错网格有限差分法和有限体积法对回线源瞬变电磁法的三维正演.通过对比低阻块状体的积分方程法、时域有限差分法、矢量有限单元法和SLDM法的数值解,验证了交错网格有限差分法和有限体积法的正确性.由于交错网格有限差分法、有限体积法和基于矩形块单元的矢量有限单元法将待求电场均定义在矩形块单元棱边上,因此三种数值算法可采用相同方法进行电场待求量编码、计算背景场和后处理.然而,与矢量有限单元法相比,交错网格有限差分法和有限体积法的系数矩阵更加稀疏,求解效率更高.通过对水平低阻板状体三维模型的数值模拟,我们发现本研究中交错网格有限差分法比有限体积法精度更高;再利用一维解析法求解相应三层层状地电模型的感应电动势,我们还发现两种数值算法和一维解析法计算的感应电动势等值线形状吻合程度高,只是数值范围略有差异.
关键词: Helmholtz方程      回线源      瞬变电磁法      有限差分法      有限体积法     
hree-dimensional forward calculation for loop source transient electromagnetic method based on electric field Helmholtz equation
LI Jian-Hui1, HU Xiang-Yun1, ZENG Si-Hong1, LU Jin-Ge1, HUO Guang-Pu2, HAN Bo1, PENG Rong-Hua1     
1. Institute of Geophysics and Geometics, China University of Geoscience, Wuhan 430074, China;
2. Henan Institute of Geological Survey, Zhengzhou 450001, China
Abstract: Forward calculation is the base of parameters selection in field work, data processing, and interpretation; particularly the accurate, stable and efficient three-dimensional forward solver is very important. In this paper, based on the strategy that solving the Laplace-domain electrical field firstly, then acquiring the time-domain magnetic field by Gaver-Stehfest algorithm, we employed staggered finite difference (SFD) and finite volume (FV) to discretize the electrical field Helmholtz equation for three-dimensional forward calculation of loop source transient electromagnetic method. Comparing with the numerical solutions of integral equation (IE), finite difference in time domain (FDTD), vector finite element method (VFEM) and spectral Lanczos decomposition method (SLDM) for low resistivity body in homogeneous half-space, we validated the SFD and FV algorithms. Because the unknown electrical field is arranged on the edge of brick for SFD, FV and VFEM based on brick element, so these three numerical algorithms could utilize the same method for coding the unknown electrical field, computing the background field and post processing. However, comparing with the VFEM, the coefficient matrixes produced by SFD and FV are sparser, so these two algorithms are more efficient. For low resistivity horizontal sheet model, we use SFD, FV and 1D analytical method to compute the induced electromagnetic force (EMF). The results show that SFD has higher accuracy than FV, and the contours of EMF for these three algorithms coincide with each other very well although the range of EMF has fine distinctions..
Key words: Helmholtz equation      Loop source      Transient electromagnetic method      Finite difference      Finite volume     
1 引言

瞬变电磁法已广泛应用于金属矿产勘探、煤田勘探、环境水文调查等领域[1-8].瞬变电磁法勘探流程中, 可行性分析、勘探方案设计和资料处理与解释三个环节均需要正演参与, 因此精确、稳定、高效的正演算法对电磁法勘探具有重要意义.

瞬变电磁法三维数值模拟主要有两种思路:第一种思路是时步法, 即对控制方程在时间域利用数值差分算法进行显式或隐式离散, 利用解析解或数值算法获得初始时刻电磁场后, 按照一定时间步长开始迭代, 每迭代一次均需利用数值算法, 比如有限差分法、有限单元法, 计算一次所有剖分网格的电磁场, 直到获得指定时刻的电磁场[9].Wang与Hohmann采用时步法对回线源激励的瞬变电磁场进行了三维数值模拟.他们利用基于交错网格的显式有限差分法, 对麦克斯韦方程组采用DuFort-Frankel法进行离散, 并以解析法计算的均匀半空间浅层网格电场、磁场脉冲响应作为初始值[10]. Commer与Newman沿用了Wang与Hohmann的数值模拟思路, 并做了如下改进:对接地长导线激发的电磁场进行了数值模拟, 使用有限差分法计算了复杂地电模型的电场值作为初始值, 并使用并行算法提高了计算效率[11].Um等学者继承了Commer和Newman的思路, 不同之处在于:以电场扩散方程作为控制方程代替了麦克斯韦方程组; 时间域离散时, 以隐式的向后欧拉法代替了显式的DuFort-Frankel法; 数值方法选用矢量有限单元法[12-13]. 2012年, 林君等将Wang和Hohmann的算法应用于全波形时间域航空电磁法的数值模拟, 为全波形三维反演和仪器设计奠定了基础[14].2013年, 李貅等对Wang与Hohmann算法做了重大改进, 主要体现在两方面:通过将矩形回线源电流密度加入麦克斯韦方程组的安培环路定理方程, 实现回线源瞬变电磁激发源加入; 在激发电流源的计算中考虑关断时间[15].另外, Druskin与Knizhnerman使用了频谱微分-差分技术用于模拟三维瞬变电磁场, 后来这种方法被称为spectral Lanczos decomposition method (SLDM)[16].

瞬变电磁法三维数值模拟的第二种思路为:先求得频率域电磁场, 再通过余弦变换或Gaver-Stehfest算法获得时间域电磁场.Newman等先利用积分方程法获得频率域回线源激发的电磁场, 再利用频率域电磁场的实部项由余弦变换获得了磁场脉冲响应[17-18].Sugeng和本文作者先后使用矢量有限单元法获得了回线源激发的瞬变电磁场响应, 不同的是Sugeng使用了逆傅里叶变换, 而本文作者使用的是Gaver-Stehfest算法[19-21].

本文三维正演采用先求解拉普拉斯域电场、再应用Gaver-Stehfest算法获得时间域磁场的思路, 数值方法选用基于交错网格的有限差分法和有限体积法以及基于矩形块网格的矢量有限单元法.交错网格有限差分法精度高、速度快, 目前该算法已广泛应用于大地电磁法、可控源电磁法的三维正、反演[22-26], 有限体积法也已应用于可控源电磁法的三维正演[27-28]; 而这两种数值算法却鲜见于瞬变电磁法的三维正演.

2 电场Helmholtz方程

本文中, 三种数值方法使用的控制方程均为拉普拉斯域电场异常场Helmholtz方程.在采用正谐时eiωt, 并不考虑位移电流的情况下, 拉普拉斯域电场异常场Helmholtz方程可写作[20-21]

(1)

式中为拉普拉斯域电场异常场; σ为地电模型的电导率, σb表示不存在电性异常体的均匀半空间的电导率; σa表示某一点的异常电导率, 即σa=σ-σb, 其值仅在异常体内不为0;上标b和a分别表示背景场和异常场; s(=iω)为拉普拉斯域变量, μ0为真空中磁导率.为拉普拉斯域电场背景场, 可采用一维解析法获得.

由于电磁场在导电介质中的自然衰减, 只要选取足够大的计算区域, 即可认为电场异常场在边界处的影响可忽略不计, 那么边界条件为

(2)

3 数值算法 3.1 网格剖分

当采用有限差分法和有限体积法时, 我们使用交错网格剖分计算区域.为了与基于矩形块单元的矢量有限元相对比(如图 1), 我们采用电场定义在棱边、磁场定义在表面的交错网格(如图 2).

图 1 矩形块单元示意图 Fig. 1 The sketch map of brick element
图 2 交错网格示意图 Fig. 2 The sketch map of staggered grid
3.2 交错网格有限差分法

交错网格有限差分法思路如下:根据旋度定义, 将(1)式中双旋度算子展开, xyz三个方向的电场定义分别在(i+1/2, j, k)、(i, j+1/2, k)和(i, j, k+1/2)处, 可得到如下三式[22]:

(3)

(4)

(5)

式中, Δxi、Δyj和Δzk分别是单元(i, j, k)在xyz三方向的棱边长度; Δx'i为单元(i-1/2, j, k)和(i+ 1/2, j, k)之间的棱边长度, Δy'j是单元(i, j-1/2, k)和(i, j+1/2, k)之间的棱边长度, Δz'i是单元(i, j, k-1/2)和(i, j, k+1/2)之间的棱边长度; 其中1/2表示棱边中点.

分别为xyz三方向的模型电导率, 分别为xyz三方向的背景电导率.以为例简述一下其定义方式, 定义于棱边(i+1/2, j, k), 那么就是共用该棱边的四个长方体电导率的加权平均, 权重为单个矩形块截面积与四个矩形块截面积之和的比值[22].

3.3 有限体积法

有限体积法思路如下:根据矢量恒等式

(6)

式中, Ω为单个剖分单元体积、Γ是该剖分单元的外表面, n为外表面单位法向分量.那么式(1)可改写为

(7)

对式(7)在xyz三个方向分别离散, 可得到有限体积法的方程组[27-28]:

(8)

(9)

(10)

电导率定义方式如同交错网格有限差分法.事实上, 有限差分法与有限体积法的离散方程可以相互转换, 分别用ΔxiΔy'jΔz'i、Δx'iΔyjΔz'i和Δx'iΔy'jΔzi乘以式(3)、(4)和(5)即可得到式(8)、(9)和(10).

3.4 矢量有限单元法

以式(1)为控制方程, 应用Galerkin法来推导有限元方程.设余量r

(11)

那么在单元e内有

式中N为各个方向矢量基函数组成的向量.将(11)式代入(12)式, 经过矢量运算, 得

(13)

其中

,详见参考文献[20-21].

4 数值算法后处理

由于交错网格有限差分法、有限体积法和矢量有限单元法将待求电场均定义在棱边上, 因此三种数值方法在电场待求量编码、计算背景场、后处理等环节具有相似性.对计算区域经三种数值算法离散后, 均要得到大型线性方程组, 本文采用对称超松弛预处理的双共轭梯度稳定法(SSOR-BICGSTAB)求解这些方程组[29-30], 进而获得拉普拉斯域电场异常场.

获得拉普拉斯域电场异常场后, 利用Gaver-Stehfest算法获得时间域电场异常场[31-32], 其中Gaver-Stehfest算法系数数目为12.之后, 再根据法拉第电磁感应定律获得磁场异常场, 其垂直分量为:

(14)

由一维解析解获得磁场时间域背景场后, 与磁场异常场相加获得磁场二次场, 进而得到感应电动势.

5 算法验证

本文共设置了3个地电模型, 模型1和模型2为低阻块状体模型, 模型3为水平低阻板状体模型.模型1和模型2的网格剖分方案见表 1.交错网格有限差分法和有限体积法中, 对系数矩阵按行压缩存储; 而矢量有限单元法中, 对半带宽系数矩阵按行压缩存储.对于模型1, 交错网格有限差分法(SFD)、有限体积法(FV)和矢量有限单元法(VFEM)的网格剖分方案一致、待求量相同, 但矢量有限单元法生成的系数矩阵非零元数目约为交错网格有限差分法和有限体积法系数矩阵非零元数目的5.3倍, 因此生成系数矩阵、求解大型方程组更耗时.对于模型2, 交错网格有限差分法和有限体积法的网格剖分比矢量有限单元法网格剖分更加精细, 但是系数矩阵非零元数目却较少.因此, 从减少计算量、提高效率的角度来看, 交错网格有限差分法和有限体积法比矢量有限元法更适合于简单地电模型; 而矢量有限单元法的优势在于能够精确描述复杂地质体、提高数值模拟精度.

表 1 VFEM、SFD和FV的网格剖分方案 Table 1 The meshing schemes for VFEM, SFD and FV

本文的计算环境为:Win7 64位操作系统, Inter I3双核四线程处理器、主频为2.4GHz, 4G内存.我们使用Matlab语言编制了矢量有限单元法、交错网格有限差分法和有限体积法的系数矩阵生成、计算背景场以及后处理程序, 并采用基于核数并行的Parfor循环加快计算速度; 而使用Fortran语言编制了求解方程组程序, 并采用基于线程数并行的应用程序编程接口OpenMP加快计算速度.

5.1 模型1

地电模型如图 3所示:一个电阻率为0.5Ωm块状体埋藏于电阻率为10Ωm的均匀半空间, 块状体顶部埋深30m、长100m、宽40m和高30m, 采用中心回线装置.由图 4可看出:SFD、FV与VFEM[20]结果吻合非常良好, 这是由于三种算法所使用的控制方程、时频转换方法、网格剖分以及求解大型方程组的方法都是相同的; 而与积分方程法(IE)[17]和时域有限差分法(FDTD)[10]的结果略有不同:在0.0001s至0.0016s段, SFD和FV与IE结果吻合较好, 而比FDTD结果略小; 在0.0016s至0.0064s段, SFD和FV与FDTD结果吻合较好, 而比IE结果略大.

图 3 模型1断面图 Fig. 3 The section map for Model
图 4 模型1的感应电动势曲线 Fig. 4 The curves of EMF for Model

Gaver-Stehfest算法中, 拉普拉斯域变量s的离散形式为m·ln2/t, 其中m为Gaver-Stehfest算法系数编号, 取值范围为1≤mM, M为Gaver-Stehfest算法系数数目, t为观测时刻.当Gaver-Stehfest算法系数数目为12时, 获得一个时刻的感应电动势, 需求解12个拉普拉斯域方程组.图 5图 6分别为采用交错网格有限差分法和有限体积法时, 模型1的SSOR-BICGSTAB法收敛曲线.从图中可看出:对于两个时刻, m值越小, 求解方程组收敛速度越慢; 观测时刻越大, 求解方程组收敛速度也越慢.这种现象类似于频率域电磁法三维数值模拟所遇到的情况:频率低时, 式(1)左端第二项和右端项均很小, 数值算法离散形成的方程组为病态方程组; 采用迭代法求解病态方程组时, 受计算机字长所限, 舍入误差随迭代次数的增加而累积、对数值解影响很大, 进而出现收敛速度慢、求解精度低的现象[25, 33-34].

图 5 采用交错网格有限差分法时, 模型1的SSOR-BICGSTAB法收敛曲线 Fig. 5 The convergence of the SSOR-BICGSTAB for Model 1 using SFD
图 6 采用有限体积法时, 模型1的SSOR-BICGSTAB法收敛曲线 Fig. 6 The convergence of the SSOR-BICGSTAB for Model 1 using FV

我们从场论角度进一步分析这一现象.采用数值算法求解(1)式时, 由于迭代法获得的数值解无法保证电流密度的散度为零、相当于引入了"虚拟源", 而采用散度校正技术可以消除这些"虚拟源"对电磁场的不利影响, 进而提高迭代法求解精度和速度[35-36].

瞬变电磁法三维数值模拟中, 当观测时刻较小时, "虚拟源"产生的电磁场与二次场相比而言较弱、对数值模拟精度的影响可忽略不计; 观测时刻较大时, 随着二次场的衰减, "虚拟源"产生的电磁场对数值模拟精度的影响逐渐增大, 数值解精度下降、甚至出现畸变.本文数值模拟时, 观测时刻内三种数值算法精度较高, 因此未再使用散度校正技术.

5.2 模型2

地电模型如图 7所示:一个电阻率为0.333Ωm块状体埋藏于电阻率为100Ωm的均匀半空间, 块状体顶部埋深80m、长100m、宽50m和高100m.坐标原点位于发射回线几何中心, 测点x坐标分别为0m和140m, y坐标均为0m.由图 8可看出:对于两个测点(0, 0), SFD、FV和VFEM[20]计算的感应电动势曲线介于FDTD[10]和SLDM[16]计算的曲线之间; 时域有限差分计算结果略小, 而SLDM算法计算结果略大.对于测点(140, 0), 0.0024s之前, SFD、VFEM、FDTD和SLDM法计算的感应电动势曲线吻合良好, FV计算的感应电动势略大; 0.0024s之后, 5种数值方法计算结果差别较大, 尤其是0.0032s时差别最大.

图 7 模型2断面图 Fig. 7 The section map for Model 2
图 8 模型2的感应电动势曲线 (a)测点(0, 0);(b)测点(140, 0). Fig. 8 The curves of EMF for Model 2 (a) Measure point (0, 0); (b) Measure point (140, 0).

对比图 5图 6图 9图 10, 我们可发现对于相同的m·ln2/t, 模型2比模型1的收敛速度慢.这主要是由以下原因引起:①模型2中电导率较小、式(1)左端第二项和右端项较小; ②模型2中, 网格剖分更加复杂、方程组阶数更大.

图 9 采用交错网格有限差分法时, 模型2的SSOR-BICGSTAB法收敛曲线 (a)0.0001s;(b)0.0064s. Fig. 9 The convergence of the SSOR-BICGSTAB for Model 2 using SFD (a)0.0001s;(b)0.0064s.
图 10 采用有限体积法时, 模型2的SSOR-BICGSTAB法收敛曲线 (a)0.0001s;(b)0.0064s. Fig. 10 The convergence of the SSOR-BICGSTAB for Model 2 using FV (a)0.0001s;(b)0.0064s.
5.3 模型3

模型3是一个水平低阻板状体模型.该板状体水平方向尺度为1880m×1880m、厚度为20m、埋深为100m;板状体电阻率为1Ωm, 均匀半空间电阻率为100Ωm.发射回线形状为正方形、边长为500m, 坐标原点位于发射回线几何中心.地表测点xy方向的分布范围均为[-640, 640], 一共1681个测点.本例中, 我们使用交错网格有限差分法、有限体积法和一维解析法[37]计算感应电动势.

图 11为测点(0, 0)和(220, 220)的感应电动势曲线.由于水平板状体可近似看作三层层状介质, 我们还计算了一维解析解.对于两个测点, 交错网格有限差分法和有限体积法计算的感应电动势非常吻合, 与0.0001~0.0008s之间4个时刻的感应电动势一维解析解差异较为明显; 而与0.0016~0.0128s之间4个时刻的感应电动势一维解析解吻合较好.

图 11 模型3的感应电动势曲线 (a)测点(0, 0);(b)测点(220, 220). Fig. 11 The curves of EMF for Model 3 (a) Measure point (0, 0); (b) Measure point (220, 220)

图 12图 13为交错网格有限差分法、有限体积法和一维解析法计算的不同时刻感应电动势等值线图.在较小时刻, 根据两种数值算法计算的感应电动势等值线图非常吻合; 从0.0032s开始, 随着观测时刻的增大, 这种吻合程度逐渐下降.对比两种数值算法计算的感应电动势, 我们发现有限体积法计算的晚期(尤其是0.0064s和0.0128s)感应电动势等值线图在矩形发射回线4个顶点附近不光滑、畸变程度大.因此, 在本研究中, 有限体积法比交错网格有限差分法的计算精度稍差.而一维解析法计算的感应电动势等值线图始终非常光滑、不存在畸变现象.

图 12 0.0001s时刻(a1, b1, c1)、0.0002s时刻(a2, b2, c2)、0.0004s时刻(a3, b3, c3)、0.0008s时刻(a4, b4, c4), 感应电动势等值线图 (a1, a2, a3, a4)有限差分法; (b1, b2, b3, b4)有限体积法; (c1, c2, c3, c4)一维解析解. Fig. 12 The contour maps of EMF at 0.0001s (a1, b1, c1), 0.0002s (a2, b2, c2), 0.0004s (a3, b3, c3) and 0.0008s (a4, b4, c4) (a1, a2, a3, a4) SFD; (b1, b2, b3, b4) FV; (c1, c2, c3, c4) 1D solution.
图 13 0.0016s时刻(a1, b1, c1)、0.0032s时刻(a2, b2, c2)、0.0064s时刻(a3, b3, c3)、0.0128s时刻(a4, b4, c4), 感应电动势等值线图 (a1, a2, a3, a4)有限差分法; (b1, b2, b3, b4)有限体积法; (c1, c2, c3, c4)一维解析解. Fig. 13 The contour maps of EMF at 0.0016s (a1, b1, c1), 0.0032s (a2, b2, c2), 0.0064s (a3, b3, c3) and 0.0128s (a4, b4, c4) (a1, a2, a3, a4) SFD; (b1, b2, b3, b4) FV; (c1, c2, c3, c4) 1D solution.

将两种数值算法与一维解析法计算的感应电动势作比较, 我们发现:3种方法计算的感应电动势等值线图形状非常吻合, 只是在感应电动势数值范围略有差别, 比如在0.0001s, 一维解析法计算的感应电动势数值范围更大, 又如在0.0002s、0.0004s和0.0008s, 两种数值算法计算的感应电动势数值范围更大.而在0.0128s, 两种数值算法计算的感应电动势仍存在感应电动势负值, 而一维解析法计算的感应电动势已不存在负值, 这主要与三维模型计算区域、算法精度有关.

观测时刻从小到大, 感应电动势等值线的形状、正值范围逐渐变化.从0.0001s到0.0008s, 感应电动势正值只集中于发射回线范围内; 正值等值线形状为发射回线形状, 而负值等值线形状总体为圆形、局部呈闭合准椭圆形.从0.0016s开始, 感应电动势正值逐渐分布于发射回线外部, 正值等值线逐渐变为圆形、闭合的准椭圆形负值等值线消失.至0.0128s时, 感应电动势正值范围已基本占满计算区域, 等值线形状均为圆形.感应电动势正值范围和等值线形状的变化可由瞬变电磁场的传播方式合理解释.

6 结论与讨论

本文实现了基于电场Helmholtz方程的交错网格有限差分法和有限体积法对回线源瞬变电磁法的三维正演, 并通过与积分方程法、时域有限差分法、矢量有限单元法和SLDM法的数值解对比, 验证了算法的正确性.交错网格有限差分法、有限体积法和基于矩形块网格的矢量有限单元法均将待求电场定义在棱边上, 因此三种数值算法在电场待求量编码、计算背景场、后处理等方面可采用相同的处理方式.通过模型1和模型2, 我们发现交错网格有限差分法和有限体积法生成的系数矩阵比矢量有限单元法的系数矩阵更加稀疏, 算法全部用时约为矢量有限单元法的三分之一.

我们通过水平低阻板状体模型比较了交错网格有限差分法和有限体积法的计算精度.结果表明:本例中, 有限体积法计算的晚期(尤其是0.0064s和0.0128s)感应电动势等值线在矩形发射回线4个顶点附近不光滑、畸变程度大, 算法精度较低.除此之外, 我们还利用一维解析法计算了相应三层层状地电模型的感应电动势; 对比两种数值算法与一维解析法后发现:三种算法计算的感应电动势等值线形状吻合程度高, 只是感应电动势数值范围略有差别, 进一步验证了本文数值算法的正确性.

通过与多种数值算法和一维解析解比较, 本文所给出的实例具有较高计算精度.但是, 为了能够获得更晚延时的精确瞬变电磁响应, 在未来研究中应考虑使用散度校正条件.另外, 鉴于交错网格有限差分法和有限体积法的精确性和高效性, 我们可以预见未来真正能够达到实用程度的瞬变电磁法三维反演也很有可能基于这两种数值算法.

致谢

感谢两位匿名审稿专家的专业评述与建议, 使得本文更加完善.

参考文献
[1] Flores C, Peralta-Ortega S A. Induced polarization with in-loop transient electromagnetic soundings: A case study of mineral discrimination at El Arco porphyry copper, Mexico. Journal of Applied Geophysics , 2009, 68(3): 423-436. DOI:10.1016/j.jappgeo.2009.03.009
[2] Yang D K, Oldenburg D W. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit. Geophysics , 2012, 77(2): B23-B34. DOI:10.1190/geo2011-0194.1
[3] Xue G Q, Qin K Z, Li X, et al. Discovery of a Large-scale Porphyry Molybdenum Deposit in Tibet through a Modified TEM Exploration Method. Journal of Environmental & Engineering Geophysics , 2012, 17(1): 19-25.
[4] 薛国强, 陈卫营, 周楠楠, 等. 接地源瞬变电磁短偏移深部探测技术. 地球物理学报 , 2013, 56(1): 255–261. Xue G Q, Chen W Y, Zhou N N, et al. Short-offset TEM technique with a grounded wire source for deep sounding. Chinese J. Geophys. (in Chinese) , 2013, 56(1): 255-261.
[5] Xue G Q, Bai C Y, Yan S, et al. Deep sounding TEM investigation method based on a modified fixed central-loop system. Journal of Applied Geophysics , 2012, 76: 23-32. DOI:10.1016/j.jappgeo.2011.10.007
[6] Ezersky M, Legchenko A, Al-Zoubi A, et al. TEM study of the geoelectrical structure and groundwater salinity of the Nahal Hever sinkhole site, Dead Sea shore, Israel. Journal of Applied Geophysics , 2009, 75(1): 99-112.
[7] Porsani J L, Bortolozo C A, Almeida E R, et al. TDEMsurvey in urban environmental for hydrogeological study at USP campus in So Paulo city, Brazil. Journal of Applied Geophysics , 2012, 76: 102-108. DOI:10.1016/j.jappgeo.2011.10.001
[8] Porsani J L, Almeida E R, Bortolozo C A, et al. TDEM survey in an area of seismicity induced by water wells in Paraná sedimentary basin, Northern So Paulo State, Brazil. Journal of Applied Geophysics , 2012, 82: 75-83. DOI:10.1016/j.jappgeo.2012.02.005
[9] 李建慧, 朱自强, 曾思红, 等. 瞬变电磁法正演计算进展. 地球物理学进展 , 2012, 27(4): 1393–1400. Li J H, Zhu Z Q, Zeng S H, et al. Progress of forward computation in transient electromagnetic method. Progress in Geophys. (in Chinese) , 2012, 27(4): 1393-1400.
[10] Wang T, Hohmann G W. A finite-difference, time-domain solution for three-dimensional electromagnetic modeling. Geophysics , 1993, 58(6): 797-809. DOI:10.1190/1.1443465
[11] Commer M, Newman G. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources. Geophysics , 2004, 69(5): 1192-1202. DOI:10.1190/1.1801936
[12] Um E S, Harris J M, Alumbaugh D L. 3D time-domain simulation of electromagnetic diffusion phenomena: A finite-element electric-field approach. Geophysics , 2010, 75(4).
[13] Um E S, Harris J M, Alumbaugh D L. An iterative finite element time-domain method for simulating three-dimensional electromagnetic diffusion in earth. Geophysical Journal International , 2012, 190(2): 871-886. DOI:10.1111/gji.2012.190.issue-2
[14] 许洋铖, 林君, 李肃义, 等. 全波形时间域航空电磁响应三维有限差分数值计算. 地球物理学报 , 2012, 55(6): 2105–2114. Xu Y C, Lin J, Li S Y, et al. Calculation of full-waveform airborne electromagnetic response with three-dimension finite-difference solution in time-domain. Chinese J. Geophys. (in Chinese) , 2012, 55(6): 2105-2114.
[15] 孙怀凤, 李貅, 李术才, 等. 考虑关断时间的回线源激发TEM三维时域有限差分正演. 地球物理学报 , 2012, 56(3): 1049–1064. Sun H F, Li X, Li S C, et al. Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time. Chinese J. Geophys. (in Chinese) , 2012, 56(3): 1049-1064.
[16] Druskin V, Knizhnerman L. Spectral differential-difference method for numeric solution of three-dimensional nonstationary problems of electric prospecting. Izvestiya: Earth Phys. , 1988, 24(8): 641-648.
[17] Newman G A, Hohmann G W, Anderson W L. Transient electromagnetic response of a three-dimensional body in a layered earth. Geophysics , 1986, 51(8): 1608-1627. DOI:10.1190/1.1442212
[18] Newman G A, Hohmann G W. Transient electromagnetic response of high-contrast prisms in a layered earth. Geophysics , 1988, 53(5): 691-706. DOI:10.1190/1.1442503
[19] Sugeng F. Modeling the 3D TDEM response using the 3D full-domain finite-element method based on the hexahedral edge-element technique. Exploration Geophysics , 1998, 29(4): 615-619. DOI:10.1071/EG998615
[20] Li J H, Zhu Z Q, Liu S C, et al. 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 , 2011, 8(4): 560-567. DOI:10.1088/1742-2132/8/4/008
[21] 李建慧, 朱自强, 鲁光银, 等. 回线源瞬变电磁法的三维正演研究. 地球物理学进展 , 2013, 28(2): 754–765. Li J H, Zhu Z Q, Lu G Y, et al. Study on three-dimensional forward of transient electromagnetic method excited by loop source. Progress in Geophys. (in Chinese) , 2013, 28(2): 754-765.
[22] Newman G A, Alumbaugh D L. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophysical Prospecting , 1995, 43(8): 1021-1042. DOI:10.1111/gpr.1995.43.issue-8
[23] 谭捍东, 余钦范, BookerJ, 等. 大地电磁法三维交错采样有限差分数值模拟. 地球物理学报 , 2003, 46(5): 705–711. Tan H D, Yu Q F, Booker J, et al. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method. Chinese J. Geophys. (in Chinese) , 2003, 46(5): 705-711.
[24] 陈辉, 邓居智, 谭捍东, 等. 大地电磁三维交错网格采样有限差分数值模拟中的散度校正方法研究. 地球物理学报 , 2011, 54(6): 1649–1659. Chen H, Deng J Z, Tan H D, et al. Study on divergence correction method in three-dimensional magnetotelluric modeling with staggered-grid finite difference method. Chinese J. Geophys. (in Chinese) , 2011, 54(6): 1649-1659.
[25] Sasaki Y, Meju M A. Useful characteristics of shallow and deep marine CSEM responses inferred from 3D finite-difference modeling. Geophysics , 2009, 74(5).
[26] Streich R. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy. Geophysics , 2009, 74(5).
[27] Weiss C J, Constable S. Mapping thin resistors and hydrocarbons with marine EM methods, Part Ⅱ-Modeling and analysis in 3D. Geophysics , 2006, 71(6).
[28] 杨波, 徐义贤, 何展翔, 等. 考虑海底地形的三维频率域可控源电磁响应有限体积法模拟. 地球物理学报 , 2012, 55(4): 1390–1399. Yang B, Xu Y X, He Z X, et al. 3D frequency-domain modeling of marine controlled source electromagnetic responses with topography using finite volume method. Chinese J. Geophys. (in Chinese) , 2012, 55(4): 1390-1399.
[29] Barrett R, Berry M, Chan T F, et al. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (2nd edition). Philadelphia, PA: SIAM , 1994: 19-20, 37-39.
[30] Saad Y. Iterative Methods for Sparse Linear Systems. Second Edition. Philadelphia, PA: SIAM , 2002: 205-229.
[31] Knight J H, Raiche A P. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method. Geophysics , 1982, 47(1): 47-50. DOI:10.1190/1.1441280
[32] Everett M E. Transient electromagnetic response of a loop source over a rough geological medium. Geophysical Journal International , 2009, 177(2): 421-429. DOI:10.1111/gji.2009.177.issue-2
[33] Mitsuhata Y, Uchida T. 3D magnetotelluric modeling using the T-Ω finite-element method. Geophysics , 2004, 69(1): 108-119. DOI:10.1190/1.1649380
[34] Nam M J, Kim H J, Song Y, et al. 3D magnetotelluric modelling including surface topography. Geophysical Prospecting , 2007, 55(2): 277-287. DOI:10.1111/gpr.2007.55.issue-2
[35] Smith J T. Conservative modeling of 3-D electromagnetic fields, Part Ⅱ: Biconjugate gradient solution and an accelerator. Geophysics , 1996, 61(5): 1319-1324. DOI:10.1190/1.1444055
[36] Farquharson C G, Miensopust M P. Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction. Journal of Applied Geophysics , 2011, 75(4): 699-710. DOI:10.1016/j.jappgeo.2011.09.025
[37] 李建慧, 朱自强, 刘树才, 等. 基于Gaver-Stehfest算法的矩形发射回线激发的瞬变电磁场. 石油地球物理勘探 , 2011, 46(3): 489–492. Li J H, Zhu Z Q, Liu S C, et al. Rectangular loop transient electromagnetic field expressed by Gaver-Stehfest algorithm. Oil Geophysical Prospecting (in Chinese) , 2011, 46(3): 489-492.