地球物理学报  2017, Vol. 60 Issue (1): 369-382   PDF    
基于瞬时电流脉冲的三维时间域航空电磁全波形正演模拟
齐彦福 , 殷长春 , 刘云鹤 , 蔡晶     
吉林大学地球探测科学与技术学院, 长春 130021
摘要: 传统时间域航空电磁全波形正演模拟主要采用间接法(褶积算法)和直接法(时域有限差分方法等),然而褶积算法需要获得精确的电流二阶导数,这给发射电流数据采集工作带来极大挑战;时域有限差分方法受到网格和时间步长的严格限制,缺乏灵活性.为解决这些问题,本文采用时域有限元方法,通过直接改变每个时间道上的瞬时电流强度模拟任意发射波形的电磁响应.由于无需计算电流二阶导数,大大提高了正演结果的精度.利用基于非结构四面体网格的矢量有限元方法和后推欧拉技术对时间域电场扩散方程进行空间和时间离散,实现三维航空电磁时间域全波形的直接正演模拟.由此不仅可以模拟复杂的地电结构,而且基于后推欧拉法的无条件稳定性,可以更加灵活地选取时间步长,提高计算效率.通过与1D数值模拟结果进行对比验证了该方法的准确性.本文对三维柱状体模型上HELITEM MULTIPULSE和VTEM系统实际发射波形电磁响应进行模拟,并与褶积算法的结果进行比较,验证了本文算法模拟实际发射波形电磁响应的优越性.对复杂三维地质体模型上不同发射波形电磁响应进行模拟,验证了时间域有限元算法可有效处理复杂地下地质结构.
关键词: 航空电磁      时间域      全波形正演      矢量有限元      后推欧拉法     
3D time-domain airborne EM full-wave forward modeling based on instantaneous current pulse
QI Yan-Fu, YIN Chang-Chun, LIU Yun-He, CAI Jing     
College of Geo-exploration Sciences and Technology, Jilin University, Changchun 130021, China
Abstract: Traditional airborne EM full-wave forward modeling may be problematic when simulating complex transmitting waveforms and underground structures. The convolution algorithm requires high-precision of second order of derivative of transmitting current, which brings big challenges to AEM modeling, while the FDTD method is restrained by the mesh size and time step. To solve these problems, we adopt the edge finite-element method based on unstructured grids in combination with backward Euler scheme to perform 3D time-domain airborne EM (ATEM) modeling. This research will lay a foundation for the processing and interpretation of 3D full-wave time-domain airborne EM data. We start from time-domain electric field diffusion equation and discretize it with edge finite-element method based on unstructured grids. The vector basis functions can automatically satisfy the divergence-free property of electrical fields. The tetrahedral grids have big advantages in simulating complex geological structures. We further adopt Galerkin's method to obtain the finite-element governing equation. The backward Euler scheme is then used to discretize the governing equation to establish an unconditionally stable implicit equations system. We simulate the full-wave responses of arbitrary transmitting waveform by directly changing the instantaneous current for each time channel. To check the accuracy of our algorithm, we first construct a homogeneous half-space model for a center-loop AEM system. We calculate the responses for a half-sine and trapezoid transmitting wave that shows a good agreement with 1D solutions. Secondly, we simulate the EM responses for the transmitting wave used by HELITEM MULTIPULSE system over a 3D geological model with both convolution algorithm and FETD method. The modeling results of dBz/dt from the convolutional algorithm shows drastic oscillations for both on- and off-time. That's due to the low precision of second derivative of transmitting current. However, the results obtained from FETD are stable and smooth. Thirdly, to further prove the capability of FETD to model complex-waveforms, we simulate the responses of VTEM-system. It shows that the results obtained by FETD method is more stable. Fourthly, we model the full-wave responses for half-sine, trapezoidal-, and triangular transmitting waves with the same energy over a vertical plate model. The comparison of the responses for the three waveforms shows that trapezoidal wave has the biggest investigation depth. Finally, we simulate the full-wave responses for a dipping plates model, which verifies the capability of FETD method to handle complex geological structures. We have successfully implemented the 3D ATEM forward modeling with edge finite-element method based on tetrahedral grids. Different from the convolution algorithm, the full-wave responses of arbitrary transmitting waveform is simulated by directly changing the instantaneous current of every time channel. The use of vector basis functions and unstructured grids avoid the false solutions for complex structures. As a stiff integration, backward Euler scheme is unconditional stable that relaxes the condition on the time steps, resulting high-precision AEM responses..
Key words: Airborne EM      Time-domain      Full-wave modeling      Edge finite-element      Backward Euler scheme     
1 引言

航空电磁法(AEM)利用飞机作为移动平台搭载电磁勘探设备,可以对高山、沙漠、湖泊沼泽、森林等难以开展地面工作的地区进行快速、高效的电磁探测(Pemberton,1962).从20世纪40年代航空电磁系统在加拿大首飞成功(Fountain,1998),经过近70年的发展,AEM方法已经成为一种十分成熟的地球物理勘查手段,并被广泛应用于金属矿产、油气资源、地下水和地热勘查、地质填图以及环境工程领域(Dobrin,1976Palacky,1981Beamish,2002Wynn,2002Smith,2010).

根据发射和观测信号的不同,航空电磁系统可分为频率域和时间域系统.与前者相比,时间域航空电磁法(ATEM)具有更大的勘探深度.通过发射大功率脉冲一次磁场,在地下导电介质中产生感应电流(涡流),观测由其产生随时间变化的感应二次场,可以实现对深部目标体的精确探测.

为了提高时间域航空电磁系统的勘探能力,人们尝试对ATEM发射波形进行改进.Liu(1998) 讨论了发射电流波形对航空电磁响应的影响特征,针对不同的勘探目标给出了最佳波形参数的选取方法.陈曙东等(2012) 利用自由空间中回线作为目标体,推导出方波、梯形波、半正弦和三角波激励的响应,并对发射电流波形对瞬变电磁响应的影响进行研究.Chen和Hodges(2015) 开发了HELITEM MULTIPULSE系统,在半正弦波发射结束后发射一个小的梯形波,补充发射信号的高频成分,实现对浅部和深部目标体的同步高精度观测.殷长春等(2015a)以半正弦波和梯形波为例,系统地研究了航空电磁法对地下典型目标体的勘测能力,为获得对地下目标体最佳勘探能力提供参数组合.

在数值模拟方面,传统的时间域航空电磁三维正演主要采用频-时转换算法(殷长春等,2015b),即首先利用有限差分(Newman and Alumbaugh,1995Smith,1996)、积分方程(Wannamaker et al.,1894Zhdanov and Tartaras,2002)或有限元(Mogi,1996Badea et al.,2001)等方法计算0.001 Hz到108 Hz多个频率的频率域电磁响应,再采用Gaver-Stehfast逆拉氏变换(Knight and Raich,1982)或余弦变换(考夫曼和凯勒,1987陈乐寿和王光锷,1991)等算法将其转换成时间域响应.除了间接算法之外,时域有限差分(Wang and Homhmann,1993)作为一种直接算法也被广泛应用于时间域电磁正演计算当中,此方法采用时间递推,交替求解电场和磁场的策略,实现时间域电磁扩散的模拟.

为了获得了ATEM系统全波形正演响应,殷长春等(2013) 提出采用阶跃响应与电流的时间导数褶积代替脉冲响应与发射电流褶积,有效地客服了脉冲响应在t→0时产生的奇异性(Smith and Lee,2004).此类方法在计算dB/dt响应时需要使用发射电流的二阶导数,考虑到实际工作中采样率有限,难以获得电流二阶导数的精确值,严重影响褶积结果.在直接算法方面,许洋铖等(2012) 基于无源有限差分法,将发射电流离散为若干梯形脉冲元,并将其产生的磁场转换为初始场,在迭代过程中分时引入差分方程,实现全波形时间域航空电磁响应模拟.孙怀凤等(2013) 对传统的时域有限差分法进行改进,将矩形回线源直接加入到安培环路定理中,模拟了考虑关断时间的TEM三维正演响应.然而时域有限差分方法采用规则六面体网格,在模拟复杂地电结构及地形影响时只能采用“楼梯”近似,不仅缺乏灵活性,而且耗费大量网格,影响计算效率.而且时域有限差分法采用显示方程,必须满足Courant稳定性条件,严格限制了时间步长的选取.在源处理方面,要求必须将发射源放到网格棱边上,这使得发射源形状受到限制,从而无法适应实际工作需要.

Um等(2010) 采用基于四面体网格的矢量有限元法对时间域电场扩散方程进行空间离散,然后利用后推欧拉法对其进行时间离散,获得无条件稳定的隐式方程,实现了海洋可控源三维时间域电磁法的正演模拟.此方法不仅可以准确地模拟地下地质体的复杂几何形态、起伏地形和任意形状的发射源,而且放宽了对时间步长的限制,减少了时间迭代次数.本文将该方法成功应用于时间域航空电磁正演模拟当中.需特别指出的是,本文通过改变每个时间道上发射的瞬时电流强度,我们实现了对任意发射波形电磁响应的直接正演模拟.本文首先采用Galerkin方法推导有限元控制方程,然后着重介绍利用时间域有限元(FETD)方法模拟任意波形响应的基本原理.通过与1D结果进行对比,验证该方法的准确性.利用FETD方法模拟三维地质体模型上HELITEM MULTIPULSE和VTEM系统实际发射波形电磁响应,检验该方法模拟任意发射波形响应的有效性.最后模拟复杂地质体条件下半正弦波、梯形波和三角波的电磁响应,并对不同波形的响应特征和勘探能力进行分析.

2 时间域航空电磁3D正演理论

时间域电磁场扩散满足如下麦克斯韦方程组

(1)

(2)

(3)

(4)

其中e(r,t)h(r,t)分别为rt时刻的电磁场,js(r,t)是源电流密度,σ 为介质的电导率,ε为介电常数,μ为磁导率,本文只考虑电导率随空间变化的情况,介电常数和磁导率均假设为常数.考虑到位移电流对航空电磁响应的影响可以忽略(Yin and Hodges,2005),我们得到时间域电场扩散方程

(5)

采用非结构四面体网格对整个计算区域进行剖分,利用其灵活性可以较为精细地刻画复杂地质模型,并对源、接收点以及异常体所在位置进行局部网格加密.我们利用矢量有限元方法对时间域电磁扩散方程进行空间离散,四面体单元中任意位置的电场可以表示为

(6)

其中eei(t)是第e个单元第i条棱边上的电场值(图 1),nei(r)是矢量插值基函数(Nédélec,1980),由下式给出:

(7)

图 1 四面体单元棱边电场分布图 Fig. 1 Edged electrical fields in a tetrahedral element

其中lie是第e个单元第i条棱边的长度,Li1eLi2e是第i条棱边两个结点i1i2的标量插值基函数(Jin,2002).nei(r)自动满足无散条件:

(8)

而且矢量基函数的切向分量在单元界分面上连续,这一特性保证了切向电场的连续性条件.我们采用Galerkin方法在(5) 式基础上建立有限元方程,设单元加权余量为

(9)

其中

(10)

将所有单元的加权余量组合起来可以得到

(11)

其中M是单元总数.令总体加权余量R为0,可以得到电场的有限元方程(Um et al.,2010)

(12)

针对每一个单元,AeBeSe由下列各式给出

(13)

(14)

(15)

其中V是单元的体积.我们将发射线圈分解成首尾相连的若干段导线,每段导线近似为一个电偶极子(Jahandari and Farquharson,2014; Ansari and Farquharson,2014),每个电偶极子可以表示为

(16)

其中I(t)t时刻的电流强度,dl是电偶极子长度,υ表示电流方向,δ是脉冲函数,rs是电偶极子的位置.

本文通过直接改变每个时间道的瞬时电流强度,实现对任意电流波形响应的模拟.根据Um等(2010) ,对于简单的阶跃波形,可以采用一阶后推欧拉法对(12) 式进行时间域离散,可得

(17)

式中Δt是时间步长.针对较为复杂发射波形,需采用二阶后推欧拉法对(12) 式进行离散,可得

(18)

与一阶离散形式(17) 不同,(18) 式显示当前时刻的电场不仅与当前时刻的电流源有关,而且与前两个时刻电场有关.这一离散格式更加准确地刻画了电磁场的复杂扩散规律.因为利用后推欧拉法获得的隐式方程是无条件稳定的,所以放宽了对时间步长的限制,可以针对不同的时间段选取不同的Δt.在on-time阶段,针对电流变化剧烈的部分采用较小的Δt,而对于电流变化较缓的部分可以适当增大Δt;在off-time阶段,根据电磁场随时间的扩散规律,可以逐渐增大Δt,时间步长的变化不宜太剧烈,应综合考虑数值稳定性和计算效率.(17) 和(18) 式均可以简写为

(19)

其中F是大型稀疏系数矩阵,P是右端项.根据电磁波在空间中的几何和指数衰减规律,通过加大计算区域的扩边范围,使边界上的切向电场分量满足Dirichlet边界条件:

(20)

其中Г是计算区域的外边界.我们采用多波前直接求解器MUMPS(Amestoy et al.,2006)对线性方程组进行求解.从(17) 和(18) 式可以发现,当时间步长Δt不变时,系数矩阵F不发生变化,所以只需要一次分解,通过不断的回代右端项即可实现对多个时间道电磁响应的正演模拟.根据法拉第电磁感应定律((1)式))可以计算磁场响应.利用上述理论获得的时间域响应是总场信号B(t),其中包含大地感应的二次场Bs(t)以及发射线圈产生的直流场Bp(t),我们在总场中去除直流场,获得纯大地响应,

(21)

3 理论模型正演 3.1 精度验证

为了检验3D时间域有限元算法模拟任意发射波形全波形正演响应的准确性,本文设计一个电阻率为10 Ωm的均匀半空间模型(图 2).我们采用中心回线装置,发射线圈是边长为6.66 m的正八边形,飞行高度30 m.发射两种波形:1) 单位电流强度、脉冲宽度为4 ms的半正弦波(图 4a);2) 上升和下降沿时间均为0.2 ms、稳定电流时间为3.6 ms的梯形波(图 4b).整个计算区域被剖分成169923个四面体单元,共产生197055个未知数(图 3).计算中心区域为600 m×600 m×600 m,为满足Dirichlet边界条件,每个方向设置20 km的扩边.针对这两种波形,由于在0时刻没有供电,所以初始场均为0. 图 56展示了3D正演响应及与1D模型正演结果的相对百分误差.为了模拟供电时段(on-time)电磁场随发射电流的剧烈变化,我们采用较小的时间步长;为提高计算效率,断电后(off-time)根据电磁场的扩散规律我们逐渐扩大时间步长.

图 2 模型示意图 Fig. 2 Sketch map of ATEM system over a homogeneous half-space
图 3 四面体网格剖分 Fig. 3 Tetrahedral mesh at cross section y=0
图 4 发射波形 (a)半正弦波;(b)梯形波. Fig. 4 Transmitting waveform (a)Half-sine wave;(b)Trapezoidal wave.
图 5 半正弦波响应及误差分析 (a)和(c)分别为dBz/dtBz响应;(b)和(d)分别为dBz/dtBz相对百分误差. Fig. 5 EM responses for a half-sine transmitting wave and error analysis (a)and(c)are respectively dBz/dt and Bz; while(b)and(d)are percentage errors.

针对半正弦波,on-time时段采用同一时间步长1 μs,off-time时段将计算时间分为4段,第一段:Δt=1 μs,100个步长;第二段:Δt=5 μs,180个步长;第三段:Δt=25 μs,160个步长;第四段:Δt=125 μs,200个步长.图 5为半正弦波响应,从图中可以看出FETD结果与1D褶积结果吻合很好.除了dBz/dt在670 μs以及Bz响应在2.33 ms附近误差较大,其余时刻误差均在5%以内.在这两个时刻,电磁响应由负值变为正值,幅值很小.

对于梯形波,由于上升和下降沿附近电流变化剧烈,所以选用较小的时间步长0.2 μs;在电流稳定的中间部分时间步长扩大为5 μs.断电后将计算时间分成3段,分别采用时间步长2 μs,20 μs和200 μs.图 6显示dBz/dtBz均与1D结果吻合的非常好,只在200 μs、3.8 ms和4 ms电流拐点处出现较大误差.

图 6 梯形波响应及误差分析 (a)和(c)分别为dBz/dtBz响应;(b)和(d)分别为dBz/dtBz相对百分误差. Fig. 6 EM responses for a trapezoidal transmitting wave and error analysis (a)and(c)are respectively dBz/dt and Bz; while(b)and(d)are percentage errors.
3.2 HELITEM MULTIPULSE系统电磁响应模拟

为了检验利用FETD算法模拟复杂发射波形响应的有效性,本文设计了一个三维柱状体模型.低阻异常体的电阻率为1 Ωm,埋藏在100 Ωm的均匀半空间中.异常体上顶面距地表 50 m,异常体尺寸为100 m×100 m×200 m,发射和接收装置参数与图 2相同.发射波形采用HELITEM MULTIPULSE系统的理论电流波形(图 9).发射基频为30 Hz,半正弦波on-time为4 ms,off-time为10.5 ms,梯形波on-time为1 ms.共使用42万个四面体网格对整个计算区域进行剖分(图 8),中心计算区域为1000 m×1000 m×800 m,在x,y,z三个方向各有20 km的扩边,共产生49万个未知数.图 9中电流数据的采样间隔为8 μs.图 10展示了利用褶积和FETD方法计算的dBz/dtBz响应,其中实线表示正值,虚线表示负值.可以看出两种方法的Bz响应吻合的比较好,然而利用褶积算法计算的dBz/dt响应不稳定,在半正弦波on-time早期出现明显的错误,在其他时段出现剧烈震荡.这是由于褶积算法中,利用差分方法计算实际发射电流的二阶导数不准确,误差较大造成的.与其相比,FETD方法直接利用瞬时发射电流作为源项,获得的dBz/dt响应非常稳定,曲线十分光滑,这一结果有力地证明了利用FETD算法模拟任意发射波形响应的稳定性.相比之下,利用褶积算法获得的Bz响应曲线较为光滑,仅在发射脉冲的晚期道出现局部震荡,这是因为只需要计算电流的一阶导数,误差较小的缘故.

图 7 三维柱状体模型 Fig. 7 3D model with a column embedded in a homogeneous half-space
图 8 柱状体模型网格剖分 Fig. 8 Tetrahedral mesh for a columnar model at cross section y=0
图 9 HELITEM MULTIPULSE系统发射波形 Fig. 9 Transmitting waveform of HELITEM MULTIPULSE system
图 10 MULTIPULSE波形三维模型电磁响应 (a)dBz/dt响应;(b)Bz响应.实线表示正值,虚线表示负值. Fig. 10 AEM responses for MULTIPULSE transmitting waves (a)dBz/dt;(b)Bz. The solid lines are for the positive value,while the dashed lines are for negative values.
3.3 VTEM系统实际发射波形电磁响应模拟

为了进一步检验FETD算法模拟复杂发射波形电磁响应的能力,本文对VTEM系统实际发射波形(图 11)电磁响应进行模拟,地电模型参数与图 7相同.采用图 2中的发射和接收装置.实际发射电流采样间隔为8 μs.图 12展示了利用褶积方法和本文算法获得的结果,其中实线表示正值,虚线表示负值.分析电流数据发现实际发射电流并不光滑,而是在on-time段成锯齿状.由图 12可以看出,FETD方法和褶积算法得到的dBz/dt均出现了与之相对应的异常.然而在off-time段,没有供电电流的情况下,利用褶积方法获得的正演响应仍然不稳定(图 11a),出现剧烈震荡.相反,FETD的计算结果非常稳定,衰减曲线光滑.对比图 12a12b,我们同样发现利用褶积算法获得的Bz比dBz/dt响应稳定,曲线光滑.

图 11 VTEM系统发射波形 Fig. 11 Transmitting waveform of VTEM system
图 12 VTEM发射波形电磁响应 (a)dBz/dt响应;(b)Bz场.实线表示正值,虚线表示负值. Fig. 12 Electromagnetic responses for VTEM transmitting waveform (a)dBz/dt;(b)Bz. The solid lines are for the positive value,while the dashed lines are for negative values.
3.4 不同发射波形航空电磁响应及探测能力分析

为了比较不同发射波形的勘探能力,本文设计了一个100 m×200 m×200 m的直立板状体模型(图 13).异常体电阻率为1 Ωm,顶面埋深为100 m,均匀半空间电阻率为100 Ωm.发射和接收装置参数与图 2相同.为便于比较,我们假设所采用的半正弦波、梯形波和三角波形具有相同的发射能量(发射电流与时间轴构成的面积相同).图 14展示了在板状体中心正上方对不同波形的dBz/dtBz响应.从图中可以看出,三种波形在on-和off-time阶段对直立良导板均有明显反映,而且dBz/dtBz响应之间存在良好的微分/积分关系,这再次证明了FETD算法模拟任意发射波形航空电磁响应的有效性.图 15给出不同波形电磁响应的对比(为简单起见,我们仅以off-time响应为例).从图中可以看出:对于相同的激发能量,梯形波响应的幅值最大,其次分别是半正弦波和三角波响应,这一结果表明梯形波具有比半正弦波和三角波更强的探测能力.在相同的噪声水平条件下,梯形波可以实现更大深度的目标体探测.

图 13 三维直立板状体模型 Fig. 13 A vertical plate embedded in a homogeneous half-space
图 14 不同发射波形全波形电磁响应 (a)(c)和(e)分别是半正弦波、梯形波和三角波的dBz/dt响应;(b)(d)和(f)是相应的Bz响应. Fig. 14 Full-wave EM responses for different transmitting waveforms (a)(c)and(e)are dBz/dt respectively for a half-sine,trapezoidal or triangular wave; while(b),(d)and(f)are Bz responses.
图 15 不同发射波形off-time响应对比 (a)dBz/dt响应;(b)Bz响应. Fig. 15 Comparison of off-time EM responses for different transmitting waveforms (a)dBz/dt;(b)Bz.
3.5 复杂模型航空电磁响应模拟

为了检验FETD算法模拟复杂地质条件电磁响应的能力,我们设计了如图 16所示的倾斜板状体模型.两个电阻率为1 Ωm的倾斜薄板平行放置于电阻率为100 Ωm的均匀半空间中.良导板状体与z轴夹角为30°,模型大小如图 16.板状体顶部埋深为20 m,薄板间的水平距离为100 m.发射和接收参数与图 2相同.发射电流为图 4所示的半正弦波和梯形波.图 1819展示了位于板状体正上方(y=0) 测线上半正弦和梯形波全波形电磁响应.从图可以看出:两种波形的dBz/dtBz在on-和off-time时段均出现了双峰,对地下两个目标体均有较好的反映能力.从图 18(c,d)图 19(c,d)可以发现:随着时间的增加,off-time响应的双峰逐渐减弱,位置逐渐向良导板倾向方向移动,而且双峰的界限逐渐模糊,直至最后难以分辨.这是由于在脉冲关断早期高频成分丰富,对两个相邻的异常体有较强的分辨能力;随着时间向晚期道推移,电磁场穿透深度增加,异常中心向良导体倾向一侧移动;同时,由于晚期道低频成分占据主导地位,且低频电磁波波长较长,导致相邻异常体分辨能力降低,直至无法区分.与off-time电磁响应随时间单调衰减规律不同,图 18(a,b)19(a,b)给出的on-time响应的变化规律较为复杂.其中半正弦波激发的电磁响应由正向双峰异常变换为负向异常,而梯形波激发的电磁响应基本呈随时间衰减的趋势.然而,无论on-time还是off-time电磁响应均对地下目标体有良好的探测能力.

图 16 三维倾斜板状体模型 Fig. 16 Two dipping plates embedded in a homogeneous half-space
图 17 倾斜板状体模型网格剖分 Fig. 17 Tetrahedral meshes at y=0 for the dipping-plates model in Fig. 16
图 18 半正弦波航空电磁响应(模型见图 16) (a)和(b)分别是dBz/dtBz的on-time响应;(c)和(d)分别是dBz/dtBz的off-time响应. Fig. 18 AEM responses for the model in Fig. 16 for a half-sine transmitting wave (a)and(b)are on-time dBz/dt and Bz responses; while(c)and(d)are off-time responses.
图 19 梯形波航空电磁响应(模型见图 16) (a)和(b)分别是dBz/dtBz的on-time响应;(c)和(d)分别是dBz/dtBz的off-time响应. Fig. 19 AEM responses for the model in Fig. 16 for a trapezoid wave (a)and(b)are on-time dBz/dt and Bz responses; while(c)and(d)are off-time responses.
4 结论

本文采用时间域矢量有限元方法实现了3D时间域航空电磁全波形正演模拟,与褶积算法不同,我们只需将每个时间道上瞬时电流直接加入到源项中,即可成功实现对任意发射波形时间域航空电磁响应的正演模拟.由于矢量插值基函数可以自动满足电场的无散条件以及切向电场的连续性条件,从而避免了伪解的出现;非规则四面体网格可实现对复杂地电模型的精细剖分和只对源及异常体位置进行局部网格加密,从而减少了网格使用量.利用后推欧拉法获得无条件稳定的隐式方程,放宽了对时间步长的限制,大大提高了计算效率.对HELITEM MULTIPULSE系统发射波形和VTEM实际发射波形响应的模拟结果显示,由于采样率有限,不能获得精确的电流时间导数,褶积算法结果十分不稳定,出现剧烈震荡;而利用FETD算法可以获得稳定、光滑的电磁响应曲线,证实了FETD算法模拟复杂发射波形时间域航空电磁响应的有效性.通过对倾斜板状体模型进行模拟,检验了本文算法模拟复杂地质情况任意发射波形响应的能力.与半正弦和三角发射波形响应的对比结果显示梯形波具有更强的激励能力,在相同的噪声水平条件下可以实现对更大深度目标体的勘探.本文算法可以推广到其他时间域电磁勘查技术的正演模拟中,这将是我们未来的研究方向.

致谢

作者向厦门大学刘娜博士、山东大学孙怀凤博士以及吉林大学电磁“千人计划”研究团队成员在文章准备过程中提供的帮助表示感谢.我们特别向各位审稿人和编辑同志对本文提出的建设性意见表示感谢.

参考文献
Amestoy P R, Guermouche A, L'Excellent J Y, et al. 2006. Hybrid scheduling for the parallel solution of linear systems. Parallel Computing, 32(2): 136-156. DOI:10.1016/j.parco.2005.07.004
Ansari S, Farquharson C G. 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids. Geophysics, 79(4): E149-E165. DOI:10.1190/geo2013-0172.1
Badea E A, Everett M E, Newman G A, et al. 2001. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials. Geophysics, 66(3): 786-799. DOI:10.1190/1.1444968
Beamish D. 2002. An assessment of inversion methods for AEM data applied to environmental studies. J. Appl. Geophys., 51(2-4): 75-96. DOI:10.1016/S0926-9851(02)00213-6
Chen L S, Wang G E. Structural Electromagnetic Exploration (in Chinese).Wuhan: China University of Geosciences Press, 1991.
Chen S D, Lin J, Zhang S. 2012. Effect of transmitter current waveform on TEM response. Chinese J. Geophys. (in Chinese), 55(2): 709-716. DOI:10.6038/j.issn.0001-5733.2012.02.035
Chen T Y, Hodges G, Miles P. 2015. MULTIPULSE-high resolution and high power in one TDEM system. Exploration Geophysics, 46(1): 49-57. DOI:10.1071/EG14027
Dobrin M. Introduction to Geophysical Prospecting.(3rd ed). New York: McGraw-Hill, 1976.
Fountain D. 1998. Airborne electromagnetic systems-50 years of development. Exploration Geophysics, 29: 1-11. DOI:10.1071/EG998001
Knight J H, Raich A P. 1982. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method. Geophysics, 47(1): 47-50. DOI:10.1190/1.1441280
Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids. Geophysics, 79(6): E287-E302. DOI:10.1190/geo2013-0312.1
Jin J M. The Finite Element Method in Electromagnetics.(2nd ed). New York: John Wiley and Sons, 2002.
Kaufman A A, Keller G V. 1987. Frequency and Transient Soundings. Wang Y M Trans. Beijing:Geological Publishing House.
Liu G M. 1998. Effect of transmitter current waveform on airborne TEM response. Exploration Geophysics, 29(2): 35-41. DOI:10.1071/EG998035
Mogi T. 1996. Three-dimensional modeling of magnetotelluric data using finite element method. J. Appl. Geophys., 35(2-3): 185-189. DOI:10.1016/0926-9851(96)00020-1
Nédélec J C. 1980. Mixed finite elements in R3. Numer. Math., 35(3): 315-341. DOI:10.1007/BF01396415
Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophysical Prospecting, 43(8): 1021-1042. DOI:10.1111/gpr.1995.43.issue-8
Palacky G J. 1981. The airborne electromagnetic method as a tool of geological mapping. Geophysical Prospecting, 29(1): 60-88. DOI:10.1111/j.1365-2478.1981.tb01011.x
Pemberton R H. 1962. Airborne electromagnetics in review. Geophysics, 27(5): 691-713. DOI:10.1190/1.1439081
Smith J T. 1996. Conservative modeling of 3-D electromagnetic fields, Part I:properties and error analysis. Geophysics, 61(5): 1308-1318. DOI:10.1190/1.1444054
Smith R. 2010. Airborne electromagnetic methods:applications to minerals, water and hydrocarbon exploration. CSEG 2010 distinguished lecture: 7-10.
Smith R S, Lee T J. 2004. Asymptotic expansions for the calculation of the transient electromagnetic fields induced by a vertical magnetic dipole source above a conductive halfspace. Pure and Applied Geophysics, 161(2): 385-397. DOI:10.1007/s00024-003-2445-6
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 J. Geophys. (in Chinese), 56(3): 1049-1064. DOI:10.6038/cjg20130333
Um E S, Harris J M, Alumbaugh D L. 2010. 3D time-domain simulation of electromagnetic diffusion phenomena:A finite-element electric-field approach. Geophysics, 75(4): F115-F126. DOI:10.1190/1.3473694
Wang T, Hohmann G W. 1993. A finite-difference time-domain solution for three-dimensional electromagnetic modeling. Geophysics, 58(6): 797-809. DOI:10.1190/1.1443465
Wannamaker P E, Hohmann G W, SanFilipo W A. 1984. Electromagnetic modeling of three-dimensional bodied in layered earth's using integral equations. Geophysics, 49(1): 60-74. DOI:10.1190/1.1441562
Wynn J C. 2002. Evaluating groundwater in arid lands using airborne magnetic/EM methods:An example in the southwestern U. S. and northern Mexico. The Leading Edge, 21(1): 62-64. DOI:10.1190/1.1445851
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 J. Geophys. (in Chinese), 55(6): 2105-2114. DOI:10.6038/j.issn.0001-5733.2012.06.032
Yin C C, Hodge G. 2005. Influence of displacement currents on the response oh helicopter electromagnetic systems. Geophysics, 70(4): G95-G100. DOI:10.1190/1.1993710
Yin C C, Huang W, Ben F. 2013. The full-time electromagnetic modeling for time-domain airborne electromagnetic systems. Chinese J. Geophys. (in Chinese), 56(9): 3153-3162. DOI:10.6038/cig20130928
Yin C C, Ren X Y, Liu Y H, et al. 2015a. Exploration capability of airborne TEM systems for typical targets in the subsurface. Chinese J. Geophys. (in Chinese), 58(9): 3370-3379. DOI:10.6038/cjg20150929
Yin C C, Zhang B, Liu Y H, et al. 2015b. Review on airborne EM technology and developments. Chinese J. Geophys. (in Chinese), 58(8): 2637-2653. DOI:10.6038/cjg20150804
Zhdanov M S, Tartaras E. 2002. Three-dimensional inversion of multitransmitter electromagnetic data based on the localized quasi-linear approximation. Geophys. J. Int., 148(3): 506-519. DOI:10.1046/j.1365-246x.2002.01591.x
陈乐寿, 王光锷. 构造电法勘探.武汉: 中国地质大学出版社, 1991.
陈曙东, 林君, 张爽. 2012. 发射电流波形对瞬变电磁响应的影响. 地球物理学报, 55(2): 709–716. DOI:10.6038/j.issn.0001-5733.2012.02.035
考夫曼 A A, 凯勒 G V. 1987. 频率域和时间域电磁测深. 王建谋译. 北京:地质出版社.
孙怀凤, 李貅, 李术才, 等. 2013. 考虑关断时间的回线源激发TEM三维时域有限差分正演. 地球物理学报, 56(3): 1049–1064. DOI:10.6038/cjg20130333
许洋铖, 林君, 李肃义, 等. 2012. 全波形时间域航空电磁响应三维有限差分数值计算. 地球物理学报, 55(6): 2105–2114. DOI:10.6038/j.issn.0001-5733.2012.06.032
殷长春, 黄威, 贲放. 2013. 时间域航空电磁系统瞬变全时响应正演模拟. 地球物理学报, 56(9): 3153–3162. DOI:10.6038/cig20130928
殷长春, 任秀艳, 刘云鹤, 等. 2015a. 航空瞬变电磁法对地下典型目标体的探测能力研究. 地球物理学报, 58(9): 3370–3379. DOI:10.6038/cjg20150929
殷长春, 张博, 刘云鹤, 等. 2015b. 航空电磁勘查技术发展现状及展望. 地球物理学报, 58(8): 2637–2653. DOI:10.6038/cjg20150804