2021, Vol. 64
2. 煤炭资源高效开采与洁净利用国家重点实验室(煤炭科学研究总院), 北京 100013;
3. 北京市煤矿安全工程技术研究中心, 北京 100013;
4. 中国科学院大学, 北京 100049;
5. 中国科学院地质与地球物理研究所 页岩气与地质工程重点实验室, 北京 100029
, WANG GuangJie4,5
, LI HongJie1,2,3, LIAN YuGuang1,2,3, LI Wen1,2,3, QIU Hao1,2,3, MOU Yi1,2,3
2. State Key Laboratory of Coal Mining and Clean Utilization(China Coal Research Institute), Beijing 100013, China;
3. Beijing Coal Mine Safety Engineering Technology Research Center, Beijing 100013, China;
4. University of Chinese Academy of Sciences, Beijing 100049, China;
5. Key Laboratory of Shale Gas and Geoengineering, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
瞬变电磁法又称时间域电磁法，是近年来发展迅速的一种地球物理勘探方法.它具有可在近区测量、旁侧效应影响小、对低阻体敏感等优点(朴化荣，1990)，被广泛地应用于水文地质调查、金属矿勘探、工程勘查等领域.目前，实际工作中瞬变电磁法的反演解释仍以一维方法为主(薛国强等，2008)，分辨率较低，在某些地质条件复杂的地区难以取得理想效果.Di等(2021)耦合弹性波和电磁波数据，解译了华南陆块深部-浅部的地质结构.为更好地指导瞬变电磁法的资料处理、反演解释、进一步改善勘探效果，开发速度快、精度高、对复杂地电模型适用性好的三维正演技术已成必由之路.
目前常用的三维数值模拟方法主要有：积分方程法(Sanfilipo and Hohmann, 1985；Zhdanov et al., 2006)、有限差分法(Wang and Hohmann, 1993；Commer and Newman, 2004；孙怀凤, 2018, 2021)和有限元法(殷长春等, 1994, 2013；李瑞雪等, 2016a, 2016b；余翔等，2017；周建美等，2018；Zeng et al., 2019)，其中有限元法具有对复杂地质模型的适用性强、形成的总体系数矩阵是稀疏且对称、易于并行等优点，成为电磁法领域解决正演问题应用最为广泛的方法(Jin，2014).Kuo和Cho(1980)利用有限元法对低阻覆盖层下的低阻矿体的瞬变电磁法响应进行了三维数值模拟，得到的曲线形状和数量级与野外实测数据基本一致，证明了该方法的可行性.此后，Pridmore等(1981)也对有限元求解三维瞬变电磁法正演问题进行了探索性研究.Sugeng(1998)采用异常场和背景场分离的思路，使用六面体矢量单元对回线源瞬变电磁法进行了三维正演.Jin等(1999)基于兰乔斯谱分解法(SLDM)，分别运用节点有限元法和矢量有限元法对时间域和频率域的三维电磁场进行了数值模拟.Ralph-Uwe等(2008)提出了一种快速的基于频域有限元法离散、Krylov子空间投影技术建立的模型简化方法求得频域响应，然后通过余弦变换求得时间域电磁场响应.Um等(2010)提出了用于电性源瞬变电磁法的时域矢量有限元三维正演算法，其模型剖分采用四面体网格，时间离散使用隐式后退欧拉法，边界条件选用Dirichlet齐次边界条件，模拟了海洋时间域电磁场.李建慧(2011, 2013)，刘亚军(2019)利用矢量有限元法对大回线源瞬变电磁法进行了三维正演，其算法从电场双旋度方程出发，模型剖分采用六面体单元，边界条件为电磁场切向分量为零的齐次边界条件，先通过异常场/背景场法求得频域响应，然后通过Gaver-Stehfest变换得到时间域响应.张博等(2016)采用矢量有限元法对带有起伏地形的三维航空瞬变电磁法进行了数值模拟，其思路与李建慧的基本相同，不同之处在于为了更好地模拟地形起伏采用了非结构化的Delaunay四面体进行网格剖分，时频变换则采用余弦变换进行.李贺(2016)也基于时间域电场双旋度方程和电磁场切向分量为零的齐次边界条件，对典型地电模型的瞬变电磁法响应进了三维正演，模型剖分采用六面体矢量单元，时间离散则采用向后差分方式.李建慧(2017)基于齐次边界条件，采用将回线源视为多个电偶源的策略和非结构化四面体矢量单元正演了复杂形态回线源激发的瞬变电磁场分布.
综上所述，针对瞬变电磁法的三维有限元正演，目前常采用齐次边界条件作为模型的截断边界条件，这虽然简化了有限元弱形式方程的推导和后续的单元分析，但为了满足该条件，需要构建较大尺寸的模型，否则会产生由边界引起的非物理反射，降低正演结果的精度.而较大尺寸的模型意味着更多的网格单元，这会降低正演问题的求解速度.吸收边界条件能使边界的非物理反射最小化，因此可以缩小模型的尺寸、减少单元数目、加快求速解速度，或者在模型大小相同时提高正演结果的精度(Jin，2014).鉴于此，本文在前人研究的基础上引入吸收边界条件，以加快瞬变电磁法三维有限元正演的求解速度.
1 基于吸收边界条件的时间域矢量有限元法 1.1 矢量势控制方程推导考虑场源的情况下，时间域中微分形式的麦克斯韦方程组如下
|
(1) |
|
(2) |
|
(3) |
|
(4) |
式中, H为磁场强度，E为电场强度，D为电位移矢量，B为磁感应强度，Js为场源电流密度，J为场源激发的传导电流密度，ρ表示净电荷体密度，t为时间.在导电介质中，
此外，还有三个辅助的介质特性方程如下
|
(5) |
式中μ、ε、σ分别为磁导率、介电常数和电导率.针对地学问题，常见矿物的磁导率与真空磁导率十分接近，故μ=μ0.
由(2)、(4)式可假设
|
(6) |
|
(7) |
式中，A称为矢量势，φ称为标量势.
在引入库伦规范
根据(3)、(5)和(7)式可得
|
故
|
(8) |
将(6)、(8)式和介质特性方程代入(1)式，可得
|
(9) |
与频率域电磁法相比，瞬变电磁法的吸收边界条件较为困难.针对这一难题，本文的解决思路为：对于瞬变电磁法，由模型边界引起的非物理反射主要影响晚期数据，而晚期数据的频段集中在低频，针对这一特性，借鉴频率域吸收边界条件的思想，构建对低频吸收效果好的时间域吸收边界条件.
由于模型规模远大于发射线圈尺寸，外表面处的电磁波可以近似为以发射线圈为中心的球面波.借鉴Webb和Kanellopoulos(1989)的研究成果，采用形如下式的对球面波吸收较好的一阶吸收边界条件来减小模型规模，提升正演速度.
|
(10) |
式中n为模型外表面的单位法向量，r为发射线圈中心点到外表面上任意一点的距离，ε′=(ε+|ε*|)/2，ε*=ε-iσ/2πf，为复介电常数，与介电常数、电导率和指定的主吸收频率f相关，f可以在n~n×10 Hz中选取.
1.3 弱形式推导根据Galerkin加权余量法，结合矢量恒等式和散度定律，基于(9)式控制方程和(10)式边界条件的电磁场问题的弱形式方程为
|
(11) |
式中，N为矢量基函数，Ω为模型区域，Γ为其外表面.
1.4 单元分析矢量单元(又称棱边单元)很好地解决了节点单元由于未强加散度条件或强加法向连续性边界条件造成的“伪解”问题以及与结构相关的(如导体或介质的边缘或锐角)奇异性问题(Jin，2014)，因此本文采用非结构化的一阶四面体矢量单元进行单元分析，相比六面体单元不仅能更好地模拟地形起伏或倾斜界面等复杂地电模型，还能在场源和异常体附近细分网格，提高结果精度.一阶四面体矢量单元的节点与棱边间的编码规则见图 1.
|
图 1 四面体矢量单元的节点和棱边示意图 Fig. 1 Sketch map of nodes and edges in a tetrahedral vector unit |
单元e中，点(x, y, z)处的矢量势A在某一时刻t的值可由六条棱边的矢量基函数表示为(Jin，2014)
|
(12) |
式中，上标e表示棱边i等编码均为基于单元的局部编码.其中矢量基函数
|
(13) |
式中的下标i1和i2为棱边i的对应节点的编号，lie为棱边长度，Lei1和Lei2为两个节点对应的标量基函数.
矢量基函数具有如下性质：
① 散度为零，即
|
(14) |
② Nie在其对应的棱边i上的切向分量为1，在其他5条棱边上的切向分量为0，在不包含棱边i的面上的切向分量也为0.
将(13)式代入(11)式，即可对弱形式方程进行空间离散.首先对(11)式左端体积分的第一部分进行离散，可得
|
(15) |
式中{Ae}为6×1的列向量，其元素为Aie(t)，[Ee]为6×6的对称矩阵.
同理，(11)式左端体积分中第二、三部分离散后可表示为[Fe] 
|
矩阵[Ee]、[Fe]的元素Eije、Fije的表达式详见金建铭的相关著作(Jin，2014).
若单元e位于外表面Γ上则需要进一步计算(11)式左端的面积分部分：
|
(16) |
同样，[Ke]、[Me]为6×6的对称矩阵，由基函数性质的②可知，仅当Nie、Nje对应的棱边均位于法向n对应平面时，Kije和Mije才不为零，否则Nie×n或Nje×n至少有一个为零.非零的Kije一般难以得到解析解，Mije的表达式也较为复杂，因此笔者参照Rathod等(2004)的做法，首先进行仿射变换将积分区域变换为标准三角形，然后采用增参变换变为边长为2的标准正方形，最后用高斯-勒让德积分法(Gauss-Legendre Quadrature)求数值解.
对于(11)式的右端的计算，由于线电流的奇异性需将回线源看做多个较短导线元的组合，并在源附近采用较密的网格剖分.Ansari和Farquharson(2014)、李建慧等(2016)、李建慧等(2017)以该方法处理接地长导线源，取得了很好的效果.下面以图 2中的Ⅰ边为例，详述发射电流加载过程.
|
图 2 发射回线示意图 Fig. 2 Sketch map of transmitter loop |
设发射源位于z=0的平面上，四边中通过的电流方向如图 2所示.网格剖分后，设Ⅰ边分布于若干个单元中，且总是与这些单元的某一个棱边重合，其源电流密度为
|
(17) |
式中，I(t)为随时间变化的电流，指定不同的I(t)即可实现不同发射电流波形的加载，δ为冲激函数(又称狄拉克函数)，dl为导线元的长度，(0, 0，0)为其中心坐标，ex表示x方向的单位矢量.
令{Qe}为6×1的列向量，其元素为
|
(18) |
根据基函数的性质，只有导线元所在棱边的Nje·Jse≠0，其余棱边均为零，对{Qe}中的非零元素，将(17)式代入(18)式，其值为
|
(19) |
令[Re]=[Ee]+[Ke]、[Pe]=[Fe]+[Me]，将上述单元矩阵后组装起来后，最终得到形如下式的总体矩阵方程：
|
(20) |
(20) 式中未知数A(t)的求解可采用时间逐步积分法，它包括以中心差分法为代表的时间步长受算法稳定性限制的显式算法和以Newmark法为代表的时间步长只受精度限制的隐式算法(王长清和祝西里，2011).显示算法不需要计算逆矩阵，但为满足稳定性条件，时间步长必须很小，选择显式算法有极大的计算量.根据瞬变电磁场前期变化剧烈、后期变化则相对平缓的特点，本文选择时间步长可逐渐增大的Newmark法进行求解.
该方法原理如下，设t时刻的A(t)已经求得，则t+Δt时刻
|
(21) |
|
(22) |
式中α和δ是按积分精度和稳定性要求决定的参数.当参数满足δ＞0.5以及
将上述两式代入(20)式，得到从A(t)、
|
(23) |
由上式可见，递推时需对与时间步长相关的
针对系数矩阵为对称稀疏矩阵的特点，对(23)式的求解，采用针对稀疏矩阵的Intel MKL PARDISO求解器直接求解.
2 算法验证 2.1 均匀半空间均匀半空间模型的大小为6 km×6 km×6 km，其中空气层厚度2 km，电导率10-12 S·m-1，岩石厚度4 km，电导率0.01 S·m-1，发射线框为边长100 m的正方形，水平置于岩石和空气交界面(z=0)的中心，接收线圈面积100 m2，发射电流20 A，电流下降缘波形设为指数下降(如图 3所示)，关断时间10-5s.
|
图 3 电流关断波形 Fig. 3 Waveform of the turn-off current |
模型的网格剖分采用约束Delaunay算法，共剖分3764个单元，起始时间步长设为10-10 s，每5次将步长放大100.1倍，共计算336步，在CPU为i5-8250U，内存8GB的笔记本电脑上，计算时间为164 s.计算完成后通过插值提取了按对数等间距分布的10-5 s至0.01 s共31个时间门的数据，并将其与解析解进行对比(如图 4所示)，除一个时间门的相对误差稍大于3%外其余时间门的相对误差的绝对值均小于3%，31门数据的均方相对误差为1.37%，两者吻合良好，证明了本文算法正确.
|
图 4 均匀半空间数值解及其相对误差(σ=0.01 S·m-1, 6 km×6 km×6 km) Fig. 4 Numerical solution and relative errors of homogeneous half-space (σ=0.01 S·m-1, 6 km×6 km×6 km) |
其他设置不变的情况下，将吸收边界条件替换为齐次边界条件(即强制式(11)中的面积分部分为零)，此时的数值解与吸收边界条件相比，早期的精度基本一致，但随着时间的增加，电磁场逐渐向边界扩散，误差也逐渐增大，最大相对误差接近-10%，在晚期的精度显著低于吸收边界条件.
图 5为均匀半空间模型尺寸为10 km×10 km×10 km(其中空气层厚3 km)时，采用齐次边界条件得到的数值解与解析解的对比，计算时间为231 s.31门数据的均方相对误差为1.38%，与采用吸收边界条件、模型大小为6 km×6 km×6 km时的精度相近.
|
图 5 均匀半空间数值解及其相对误差(σ=0.01 S·m-1, 10 km×10 km×10 km) Fig. 5 Numerical solution and relative errors of homogeneous half-space (σ=0.01 S·m-1, 10 km×10 km×10 km) |
通过以上对比可见，采用吸收边界条件确实可以缩小模型尺寸、加快正演速度，提高正演的精度.
2.2 H型地电断面将H型地电断面模型、中心回线装置的数值解与CR1Dmod一维正演程序的结果进行对比，进一步验证本文算法的正确性.模型围岩的电导率为0.01 S·m-1，其中有一顶部埋深80 m、厚20 m、电导率0.1 S·m-1的低阻层，其余设置与均匀半空间模型相同.网格剖分时在发射线框的正下方、低阻层的顶面设置一个正方形以控制该层的网格剖分，改善求解精度(任政勇和汤井田，2009).由于低阻层及其附近需要剖分较小尺寸的网格，因此单元数较均匀半空间大幅增长，共剖分16528个单元.在10-4 s之前，时间步长的设置与均匀半空间模型相同，10-4 s之后，可以认为受低阻层的影响，电磁场的扩散速度变慢，因此每5次将步长放大100.2倍，共计算294步，计算时间为3720 s.
由图 6可见，有限元解和CR1Dmod解在早期偏离较大，在0.01 ms时相对误差接近25%，这是因为受发射电流“暂态效应”的影响，早期场值与关断时间为零的理想条件(CR1Dmod即采用此种假设)相比会有明显的增大(白登海和Maxwell，2001).此后，随着时间的增加，暂态效应的影响降低，有限元解和CR1Dmod解的相对误差也快速减小，在0.1 ms后，相对误差绝对值小于2.5%，进一步证明了本文算法的正确性.
|
图 6 H型地电断面的有限元解和CR1Dmod解及其相对误差 Fig. 6 FEM solution and CR1Dmod solution of H type geoelectric section and their relative errors |
如图 7所示，在电导率为0.01 S·m-1的围岩中有一规模为600 m×300 m×20 m、顶部埋深80 m、电导率0.1 S·m-1的低阻异常体，异常体中心在地表的投影为坐标原点，其余设置与均匀半空间模型相同.正演时，在异常体的地表投影基础上外扩200 m范围内(即x：-350至350，y：-500至500)，以50 m× 50 m的网度、移动中心回线的方式计算不同位置的瞬变电磁响应.根据回线位置的不同，模型剖分的单元数4635~5749个不等，时间步长的设置与均匀半空间模型相同，计算时间为961~1872 s.
|
图 7 水平低阻异常体模型示意图 Fig. 7 Sketchmap of horizontal low resistivity anomalous body model |
图 8为水平低阻异常体与前述的H型地电断面、均匀半空间模型在坐标原点O的瞬变电磁响应对比.前两者在0.4 ms之前的响应基本相同，与均匀半空间相比，在0.03~0.06 ms都出现了由于反射波干涉引起的电压低于均匀半空间的Undershoot现象，之后感应电压(V)逐渐高于均匀半空间.0.4 ms之后，水平低阻异常体的响应开始快速下降并最终趋于均匀半空间，而H型地电断面受低阻屏蔽效应的影响，感应电压始终高于均匀半空间.
|
图 8 O点不同模型的瞬变电磁响应 Fig. 8 TEM response of different models at point O |
图 9为具有代表性时间门的瞬变电磁响应切片.电流完全关断的时刻(0.01 ms)感应电压随位置有幅度不大的无规律波动，推测是由发射框移动引起的网格剖分的差异造成的.0.05 ms时刻，由于Undershoot现象出现了一片低电压区域，其基本呈椭圆形，与异常体的地表投影相比范围略有扩大，中心区域的电压相比正常值下降约16%.0.4 ms时刻，受异常体影响出现的高电压区域，范围相比异常体的投影有所扩大，与正常值相比，测点在距异常体边缘约50 m处感应电压的增幅大于10%，在异常体边缘电压的增幅大于25%，进入异常体投影后电压开始迅速增长并在达到某一值后趋于平稳，电压增幅最大约为160%.10 ms时刻，随着波场的扩散，异常体的影响减弱，电压趋于均匀.
|
图 9 水平低阻异常体的瞬变电磁响应时间切片 Fig. 9 Time slice of TEM response of horizontal low resistivity anomalous body model |
下面尝试对Undershoot现象做出解释.根据正演结果，抽取了发射框中心位于O点时，均匀半空间和水平低阻异常体模型在0.05 ms、0.06 ms时刻XZ平面处的电流密度，将两者相减得到前后电流密度变化ΔJ，结果见图 10，图示以外的部分尚未被感应涡流的扩散波及或电流密度变化极小.与均匀半空间相比，受低阻异常体的影响ΔJ主要集中于异常体(红色虚线所示)内，上部的ΔJ减小.以10 m×10 m的网格剖分图 11中的区域，假设在每个网格内ΔJ是均匀的，并将ΔJ形成的电流环视为磁偶极子，则其在O点引起的磁感应强度变化可按下式计算
|
(24) |
|
图 10 XZ平面0.05~0.06 ms电流密度变化 (a) 均匀半空间；(b) 水平低阻异常体. Fig. 10 Change of current density in XZ plane from 0.05 ms to 0.06 ms (a) Homogeneous half-space; (b) Horizontal low resistivity anomalous body model. |
|
图 11 倾斜低阻异常体模型示意图 Fig. 11 Sketch map of sloping low resistivity anomalous body model |
式中ds为网格面积，(y, z)为网格中心坐标.均匀半空间的计算结果为-3.99×10-7T，而低阻异常体为-3.19×10-7T.上面的估算是十分粗略的，因为很多时候并不满足源点和场点之间的距离远大于电流环半径这一磁偶极子的假设条件，但仍在一定程度上说明之所以出现Undershoot现象是因为在某些时段低阻异常体对ΔJ的集中不足以补偿其对上部ΔJ减少的影响.
3.2 倾斜低阻异常体在前节模型的基础上，将地质异常体倾角调整为30°，维持异常体顶面的中心O′至原点O的距离仍为80 m.正演时，以50 m×50 m的网度计算x：-250至250、y：-400至400矩形区域内的瞬变电磁响应.根据回线位置的不同，模型剖分的单元数5349~8792个不等，时间步长的设置与均匀半空间模型相同，计算时间为1117~2616 s.
从图 12可见，倾斜和水平低阻异常体在O点的瞬变电磁响应基本相同，只是由于倾斜低阻异常体的垂向等效厚度增大，在0.05~0.2 ms之间感应电压有7%~18%的提升.
|
图 12 O点水平和倾斜低阻异常体模型的瞬变电磁响应 Fig. 12 TEM response of models of horizontal and sloping low resistivity anomalous body at point O |
y=0线的多测道剖面见图 13，图中只显示了部分具有代表性的时间门.0.05 ms时刻，感应电压随x坐标的增长呈先降后升的波浪式起伏；0.1~0.8 ms之间，感应电压呈不对称的单波峰形态，波峰偏向于异常体埋深较浅的一侧并随时间的增加向中心移动且振幅逐渐减小；1.6 ms时，感应电压在异常体投影范围内有一定的升高但开始趋于均匀.
|
图 13 y=0线瞬变电磁多测道剖面 Fig. 13 TEM multichannel profile of y=0 line |
图 14为具有代表性时间门的瞬变电磁响应切片.0.01 ms时刻，在异常体埋深较浅的一侧出现了一片低电压条带.0.05 ms时刻，在异常体投影范围内，埋深较浅的一侧呈现高电压，较深一侧呈现低电压，且高压升幅大于低压降幅.0.4 ms时刻，感应电压的变化与水平低阻异常体的特征基本相同，但是电压最大增幅更大(约180%)且波峰偏向于异常体埋深较浅一侧.10 ms时刻，电压的分布基本均匀.
|
图 14 倾斜低阻异常体的瞬变电磁响应时间切片 Fig. 14 Time slice of TEM response of sloping low resistivity anomalous body model |
(1) 从时间域麦克斯韦方程组出发，推导了基于矢量势的微分控制方程，在此基础上结合一阶球面波吸收边界条件推导了弱形式方程，采用一阶四面体矢量单元进行单元分析、带“数值阻尼”的Newmark法进行时间离散、并将回线源视为多个电偶源的策略，实现了任意发射电流波形的瞬变电磁法响应的时间域快速求解.
(2) 通过有限元数值解与均匀半空间模型的解析解、H型地电断面模型的CR1Dmod解的对比，证明了本文算法的正确性.均匀半空间模型采用吸收边界条件和齐次边界条件的结果对比则表明：吸收边界条件确实可以提高瞬变电磁法三维正演的精度或者缩小模型尺寸、加快计算速度.
(3) 根据正演结果，其他条件相同时，水平有限体积的低阻异常体和H型地电断面的瞬变电磁响应前期基本相同，但在晚期前者会快速趋于均匀半空间的响应.对比均匀半空间和水平低阻异常体模型在0.05 ms、0.06 ms时刻电流密度的变化，表明感应电压的Undershoot现象之所以出现，是因为低阻异常体对电流的集中不足以补偿其引起的上部电流减少的影响.对于倾斜低阻异常体，由于垂向上等效厚度增大，中心处的感应电压会高于相同条件的水平低阻异常体，在沿倾向的多测道剖面上感应电压的波峰会向异常体埋深较浅的一侧偏移.
致谢 中国地质大学(武汉)胡祥云教授、中国科学院地质与地球物理研究所安志国副研究员对本文提出了有益的建议，在此表示感谢.此外，对审稿人提出宝贵的修改意见和编辑们的辛苦工作表示感谢.
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 |
Bai D H, Maxwell Meju. 2001. The effect of two types of turn-off current on TEM responses and the correction techniques. Seismology and Geology, 23(2): 118-124. |
Borner R U, Ernst O G, Spitzer K. 2008. Fast 3D simulation of transient electromagnetic fields by model reduction in the frequency domain using Krylov subspace projection. Geophysics Journal International, 173(3): 766-780. DOI:10.1111/j.1365-246X.2008.03750.x |
Commer M, Newman G. 2004. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources. Geophysics, 69(5): 1192-1202. DOI:10.1190/1.1801936 |
Di Q Y, Tian F, Suo Y, et al. 2021. Linkage of deep lithospheric structures to intraplate earthquakes: A perspective from multi-source and multi-scale geophysical data in the South China Block. Earth-Science Reviews, 214(1): 103504. DOI:10.1016/j.earscirev.2021.103504 |
Jin J, Zunoubi M R, Donepudi K C, et al. 1999. Frequency-domain and time-domain finite-element solution of Maxwell's equations using spectral Lanczos decomposition method. Computer Methods in Applied Mechanics and Engineering, 169: 279-296. DOI:10.1016/S0045-7825(98)00158-3 |
Jin J M. 2014. The Finite Element Method in Electromagnetics. Hoboken, New Jersey: Wiley.
|
Kuo J T, Cho D H. 1980. Transient time-domain electromagnetics. Geophysics, 45(2): 271-291. DOI:10.1190/1.1441082 |
Li H. 2016. Three-dimensional transient electromagnetic forward modeling in the direct time-domain by vector finite element[M.S. thesis]. Xi'an: Chang'an University.
|
Li J H. 2011. 3D numerical simulation for transient electromagnetic field excited by large source loop based on vector finite element method[Ph. D. thesis]. Changsha: Central South University.
|
Li J H, Farquharson C G, Hu X Y, et al. 2016. A vector finite element solver of three-dimensional modelling for a long grounded wire source based on total electric field. Chinese J. Geophys. (in Chinese), 59(4): 1521-1534. DOI:10.6038/cjg20160432 |
Li J H, Hu X Y, Chen B, et al. 2017. 3D electromagnetic modeling with vector finite element for a complex-shaped loop source. OGP, 52(6): 1324-1332. DOI:10.13810/j.cnki.issn.1000-7210.2017.06.024 |
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 J. Geophys. (in Chinese), 56(12): 4256-4267. DOI:10.6038/cjg20131228 |
Li R X, Wang H, Xi Z Z, et al. 2016. The 3D transient electromagnetic forward modeling of volcanogenic massive sulfide ore deposits. Chinese J. Geophys. (in Chinese), 59(12): 4505-4512. DOI:10.6038/cjg20161213 |
Li R X, Wang H, Xi Z Z, et al. 2016. Fast 3D forward modeling of transient electromagnetic. Journal of Central South University (Science and Technology), 47(10): 3477-3482. DOI:10.11817/j.issn.1672-7207.2016.10.026 |
Liu Y J, Hu X Y, Peng R H, et al. 2019. 3D forward modeling and analysis of the loop-source transient electromagnetic method based on the finite-volume method for an arbitrarily anisotropic medium. Chinese J. Geophys. (in Chinese), 62(05): 1954-1968. DOI:10.6038/cjg2019M0532 |
Piao H R. 1990. Principles of Electromagnetic Sounding. Beijing: Geological Publishing House.
|
Pridmore D F, Hohmann G W. 1981. An investigation of finite-element modeling for electrical and electromagnetic data in three dimensions. Geophysics, 46(5): 1009-1024. |
Rathod H T, Nagaraja K V, Venkatesudu B, et al. 2004. Gauss Legendre quadrature over a triangle. Journal of the Indian Institute of Science, 84(5): 183-188. |
Ren Z Y, Tang J T. 2009. Finite element modeling of 3-D DC resistivity using locally refined unstructured meshes. Chinese J. Geophys. (in Chinese), 52(10): 2627-2634. DOI:10.3969/j.issn.0001-5733.2009.10.023 |
Sanfilipo W A, Hohmann G W. 1985. Integral equation solution for the transient electromagnetic response of a three-dimensional body in a conductive half-space. Geophysics, 50(5): 798-809. DOI:10.1190/1.1441954 |
Sihong Z, Hu X Y, Li J H, et al. 2019. Effects of full transmitting-current waveforms on transient electromagnetics: Insights from modeling the Albany graphite deposit. Geophysics, 84(4): 255-268. DOI:10.1190/geo2018-0573.1 |
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 H F, Cheng M, Wu Q L, et al. 2018. A multi-scale grid scheme in three-dimensional transient electromagnetic modeling using FDTD. Chinese J. Geophys. (in Chinese), 61(12): 5096-5104. DOI:10.6038/cjg2018L0659 |
Sun H F, Liu S B, Yang Y. 2021. Crank-Nicolson FDTD 3D forward modeling for the transient electromagnetic field. Chinese J. Geophys. (in Chinese), 64(01): 343-354. DOI:10.6038/cjg2021O0229 |
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): 115-126. DOI:10.1190/1.3473694 |
Wang C Q, Zhu X L. 2011. Transient Electromagnetic Field-Theory and Calculation. Beijing: Peking University Press.
|
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 |
Wang X C. 2003. Finite element method. Beijing: Tsinghua University Press.
|
Webb J P, Kanellopoulos V N. 1989. Absorbing boundary conditions for the finite element solution of the vector wave equation. Microwave Opt. Tech. Lett., 2(10): 370-372. DOI:10.1002/mop.4650021010 |
Xue G Q, Li X, Di Q Y. 2008. Research Progress in TEM Forward Modeling and Inversion Calculation. Progress in Geophysics, 23(04): 1165-1172. |
Yin C C, Huang W, Pen 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/cjg20130928 |
Yin C C, Liu B. 1994. Three dimensional forward modeling of transient electromagnetic method and characteristics of IP effect. Chinese J. Geophys. (in Chinese), 37(Supp. 1): 486-492. |
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 J. Geophys. (in Chinese), 60(02): 810-819. DOI:10.6038/cjg20170231 |
Yun L, Wang X B, Wang Y. 2013. Numerical modeling of the 2D time-domain transient electromagnetic secondary field of the line source of the current excitation. Applied Geophysics, 10(2): 134-144. DOI:10.1007/s11770-013-0376-2 |
Zhang B, Yin C C, Liu Y H, et al. 2016. 3D modeling on topographic effect for frequency/time domain airborne EM systems. Chinese J. Geophys. (in Chinese), 59(4): 1506-1520. DOI:10.6038/cjg20160431 |
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): 333-345. DOI:10.1190/1.2358403 |
Zhou J M, Liu W T, Li X, et al. 2018. Research on the 3D mimetic finite volume method for loop-source TEM response in biaxial anisotropic formation. Chinese J. Geophys. (in Chinese), 61(01): 368-378. DOI:10.6038/cjg2018K0598 |
Zunoubi M R, Jin J, Donepudi K C, et al. 1999. A spectral Lanczos decomposition method for solving 3-D low-frequency electromagnetic diffusion by the finite-element method. IEEE Transactions on Antennas and Propagation, 47(2): 242-248. DOI:10.6038/cjg2018K0598 |
白登海, Maxwell Meju. 2001. 瞬变电磁法中两种关断电流对响应函数的影响及其应对策略. 地震地质, 23(2): 118-124. |
李贺. 2016. 直接时间域矢量有限元瞬变电磁三维正演模拟[硕士论文]. 西安: 长安大学.
|
李建慧. 2011. 基于矢量有限单元法的大回线源瞬变电磁法三维数值模拟[博士论文]. 长沙: 中南大学.
|
李建慧, Farquharson C G, 胡祥云, 等. 2016. 基于电场总场矢量有限元法的接地长导线源三维正演. 地球物理学报, 59(4): 1521-1534. DOI:10.6038/cjg20160432 |
李建慧, 胡祥云, 陈斌, 等. 2017. 复杂形态回线源激发电磁场的矢量有限元解. 石油地球物理勘探, 52(6): 1324-1332. DOI:10.13810/j.cnki.issn.1000-7210.2017.06.024 |
李建慧, 胡祥云, 曾思红, 等. 2013. 基于电场Helmholtz方程的回线源瞬变电磁法三维正演. 地球物理学报, (12): 4256-4267. DOI:10.6038/cjg20131228 |
李瑞雪, 王鹤, 席振铢, 等. 2016. 深海热液硫化物矿体3D瞬变电磁正演. 地球物理学报, 59(12): 4505-4512. DOI:10.6038/cjg20161213 |
李瑞雪, 王鹤, 席振铢, 等. 2016. 瞬变电磁快速三维正演. 中南大学学报(自然科学版), 47(10): 3477-3482. DOI:10.11817/j.issn.1672-7207.2016.10.026 |
朴化荣. 1990. 电磁测深法原理. 北京: 地质出版社.
|
任政勇, 汤井田. 2009. 基于局部加密非结构化网格的三维电阻率法有限元数值模拟. 地球物理学报, 52(10): 2627-2634. DOI:10.3969/j.issn.0001-5733.2009.10.023 |
孙怀凤, 柳尚斌, 杨洋. 2021. 瞬变电磁Crank-Nicolson FDTD三维正演. 地球物理学报, 64(01): 343-354. DOI:10.6038/cjg2021O0229 |
孙怀凤, 程铭, 吴启龙, 等. 2018. 瞬变电磁三维FDTD正演多分辨网格方法. 地球物理学报, 61(12): 5096-5104. DOI:10.6038/cjg2018L0659 |
刘亚军, 胡祥云, 彭荣华, 等. 2019. 回线源瞬变电磁法有限体积三维任意各向异性正演及分析. 地球物理学报, 62(5): 1954-1968. DOI:10.6038/cjg2019M0532 |
王长清, 祝西里. 2011. 瞬变电磁场——理论和计算. 北京: 北京大学出版社.
|
王勖成. 2003. 有限单元法. 北京: 清华大学出版社.
|
薛国强, 李貅, 底青云. 2008. 瞬变电磁法正反演问题研究进展. 地球物理学进展, 23(4): 1165-1172. |
殷长春, 黄威, 贲放. 2013. 时间域航空电磁系统瞬变全时响应正演模拟. 地球物理学报, (9): 3153-3162. DOI:10.6038/cjg20130928 |
殷长春, 刘斌. 1994. 瞬变电磁法三维问题正演及激电效应特征研究. 地球物理学报, 37(增Ⅰ): 486-492. |
张博, 殷长春, 刘云鹤, 等. 2016. 起伏地表频域/时域航空电磁系统三维正演模拟研究. 地球物理学报, 59(4): 1506-520. DOI:10.6038/cjg20160431 |
余翔, 王绪本, 李新均, 等. 2017. 时域瞬变电磁法三维有限差分正演技术研究. 地球物理学报, 60(2): 810-819. DOI:10.6038/cjg20170231 |
周建美, 刘文韬, 李貅, 等. 2018. 双轴各向异性介质中回线源瞬变电磁三维拟态有限体积正演算法. 地球物理学报, 61(1): 368-378. DOI:10.6038/cjg2018K0598 |