基于时域有限差分法的电性源感应极化效应三维数值模拟

摘要：时域电性源电磁探测法是一种可以有效快速探测矿产资源的方法。极化效应会导致电性源电磁响应快速衰减，发生符号反转现象。本文采用Cole-Cole模型描述极化效应，利用整数阶有理逼近算法实现任意分数阶Cole-Cole模型有理化；采用Yee氏网格对仿真区域进行剖分，基于时域有限差分（finite-difference time-domain，FDTD）法实现电性源感应极化效应三维数值模拟。本文对三种典型模型（均匀半空间模型、极化半空间模型和三维极化体模型）的电磁响应进行数值模拟，结果表明：均匀半空间模型的电磁响应与解析解基本吻合，并且相对误差小于10%，证明了三维数值模拟方法的正确性；极化半空间模型与三维极化体模型的电磁响应均在晚期出现了负响应，与极化理论结果相符合。

关键词：电性源；Cole-Cole模型；有理逼近；时域有限差分；极化效应

doi：10.13278/j.cnki.jjuese.20230125

中图分类号：P631

文献标志码：A

Supported by the National Natural Science Foundation of China （42030104）

Three-Dimensional Numerical Simulation of Induction-Polarization Effect of Electrical Sources Based on Finite-Difference Time-Domain Method

Ji Yanju， Deng Changwei， Wang Yuhang， Liu Hang， Wu Qiong

College of Instrument Science and Electrical Engineering， Jilin University， Changchun 130026， China

Abstract： The time-domain electrical source electromagnetic detection method is an effective and fast method for detecting mineral resources. The polarization effect can lead to rapid attenuation of the electromagnetic response and even a symbol reversal phenomenon.

In this paper，

the Cole-Cole model is used to describe the polarization effect， and arbitrary fractional Cole-Cole models are rationalized by using the integer order rational approximation algorithm. Yee’s grid is used to divide the simulation area， and the three-dimensional numerical simulation of induction-polarization effect of electrical sources is realized based on the finite difference time domain （FDTD） method. The electromagnetic responses of three typical models， namely uniform half-space model， polarized half-space model and three-dimensional polarized body model， are numerically simulated. The results show that the electromagnetic response of" uniform half-space model is basically consistent with the analytical solution， and the relative error is less than 10%， which proves the feasibility of the three-dimensional numerical simulation method. The electromagnetic responses of" polarized half-space model and three-dimensional polarized body model both have negative responses in the late stage， which is consistent with the polarization theoretical results.

Key words： electrical source; Cole-Cole model; rational approximation; finite-difference time-domain; induced polarization

0 引言

我国是矿产资源大国，同时也是世界上最大的矿产消费国之一。我国矿产资源总量大、种类丰富，但开采情况却不容乐观，尤其是铁、锰、铜等大宗矿产供应明显不足。资源供需矛盾的日益突出一定程度上阻碍了社会经济的发展[13]。

电性源时域电磁法是一种能够有效快速探测矿产资源的方法[46]。该方法的操作步骤如下：首先，利用电性源（接地长导线）发射一次脉冲场信号（一次场），在一次场的激励下，地下低阻介质中产生感应电流进而激起时变的电磁场（即二次场）；然后，使用接收线圈接收返回的二次场信号，通过分析二次场中包含的丰富的地下介质电性信息获得地下深部构造。

极化效应是一种重要的电化学效应，主要存在于侵染矿、硫化矿、金属矿等介质中。其原理是当有外加电场激励时，介质中正负电荷会发生定向运动，在界面两侧逐渐积累形成双电层；将电场撤掉后，介质中的正负电荷重新恢复到原始状态，形成了与外加电场方向相反的位移电流。极化率是金属矿的重要物性参数，基于极化效应可以实现对金属矿等介质的探测。国内外学者提出许多数学模型用于描述地下介质的极化特性，如Debye模型[7]、Dias模型[8]、Cole-Cole模型[9]等，其中最为典型的就是Cole-Cole模型；因此，本文在Cole-Cole模型基础上进行探究。

时域电磁数值模拟方法有很多种，例如时域有限差分（finite-difference time-domain， FDTD）法、时域有限体积法、时域积分方程法等。在这些算法中，FDTD算法具有简明直观的特点，非常适合宽频计算，被广泛应用于三维时域电磁数值模拟[1012]。

本文针对电性源激励下极化介质的感应极化效应进行三维数值模拟，使用Cole-Cole模型描述地下极化介质的极化特性，采用整数阶有理逼近方法对任意分数阶的Cole-Cole模型电导率表达式进行有理化逼近，获得Cole-Cole电导率时域表达方程。采用Yee氏网格对仿真区域进行剖分，将Cole-Cole电导率时域表达式代入麦克斯韦方程组中，推导电性源初始场公式与电磁迭代公式，并利用时域有限差分法计算电性源激励下极化介质的电磁响应。对均匀半空间模型、极化半空间模型与三维极化体模型的电磁响应进行数值模拟，以验证数值模拟算法的可行性，为实际野外探测的实测数据解释奠定基础。

1 感应极化效应建模

1.1 Cole-Cole模型频域表达

极化效应与电容的充放电效应类似，表现为大地电导率的频散特性。Cole-Cole模型电导率频域表达式为

σ（ω）=σSymboleB@1-η1+（1－η）（iωτ）c。（1）

式中：ω为角频率；σSymboleB@为无穷频率极限下的电导率；η为极化率；τ为时间常数；c为频散系数。η与c取值范围均为［0，1］。当c取值为1时，Cole-Cole模型简化为Debye模型[13]。

1.2 Cole-Cole模型有理函数逼近

极化介质中的电导率是时变函数，为对时域电磁响应进行数值模拟，需要将分数阶Cole-Cole模型电导率的时域表达方程代入欧姆定律中进行卷积计算；因此需要对Cole-Cole模型表达式进行频时转换。然而，与整数阶系统不同，分数阶系统在本质上是无穷维系统，对其直接进行计算非常困难；因此采用整数阶有理函数逼近法对任意分数阶Cole-Cole模型进行有理化处理[14]，得到三维空间坐标r位置处Cole-Cole模型电导率时域表达式为：

σ（r，t）=σSymboleB@（r）δ（t）－σ^（r，t）u（t）；σ^（r，t）=ησSymboleB@（r）·∑Nn=1rne－pnt。（2）

式中：t为时间；u（t）为阶跃函数；rn、pn为Cole-Cole模型的有理函数逼近系数；N为拟合阶次。

为验证分数阶Cole-Cole模型的逼近效果，对参数为c=0.6、η=0.5、τ=0.01 s的模型进行逼近，并将N依次设置为7、8、9、10。不同拟合阶次下拟合电导率与实际电导率的幅值与相角伯德图如图1所示，可以看到不同拟合阶次均取得了较好的拟合效果。

为直观观察拟合误差与拟合阶次的关系，分别计算幅值与相位有理逼近的误差，结果如图2所示。从图2中可以看出，随着拟合阶次的逐渐提高，拟合准确度也逐渐提升。但是较大的拟合阶次会导致e指数项的增多，进而提高了计算量、降低了计算速度。因此，在保证较高拟合精度的前提下应尽量选择相对低阶的拟合电导率进行后续的迭代操作。

2 三维时域电磁响应数值模拟

2.1 仿真区域Yee氏网格剖分

为了对电磁响应问题进行求解，Yee提出了时域有限差分方法，通过对电场磁场交替采样，利用变微分为差分的方法构建迭代方程，实现对电磁迭代过程的表征并得到计算结果[1516]。本文采用对仿真区域进行Yee氏网格剖分。

如图3所示，ex、ey与ez分别x、y与z方向的电场分量，分布在Yee元胞的棱中心点；hx、hy与hz分别为x、y与z方向的磁场分量，分布在Yee元胞的面中心点；二者在空间中交叉分布，彼此环绕。其排布规律与取样方式符合法拉第感应定律与安培

环路定律，能够恰当描述电磁波在空间的传播过程，且便于麦克斯韦方程的差分计算。此外，电场与磁场在时间上交替取值，彼此相差二分之一的时间步长。电场信号在整数时刻t采样，磁场信号在半整数时刻t+1/2采样。通过对麦克斯韦旋度方程组进行差分离散，构成显示差分格式，从而迭代出各个时刻的电磁场情况。

2.2 控制方程选取

麦克斯韦方程组全面总结了电磁场的变化规律，是宏观电磁场理论的基础，它表征着变化的电场与变化的磁场之间相互激励的密切联系。利用麦克斯韦方程，并带入初始值与边界条件，就可以对自然界中的电磁传播现象进行数值模拟。

对麦克斯韦方程组的时域表达式进行拉氏变换，获得频域表达式为：

式中：E为电场强度；μ为磁导率；H为磁场强度；σ为电导率；ε为介电常数；B为磁感应强度；D为电位移。

考虑极化效应的影响，将Cole-Cole模型频域表达式（式（1））代入到式（4）中，全电流定律公式可变为

对式（7）进行拉式逆变换后得到

因为极化介质中电导率为时变函数，所以欧姆定律由乘积形式改变为卷积形式：

J（t）=∫t－

式中，J为电流密度。式（8）与式（3）、式（5）的时域微分形式共同构成控制方程组。

2.3 电性源阶跃电场表达

对计算区域进行Yee氏网格剖分后，进行电磁响应迭代计算，迭代的基本思想为：利用前一时刻的电场值与前半整数时刻的磁场值推导下一时刻的电场值；利用上一步求出的电场值与前半整数时刻的磁场值求解下半整数时刻的磁场值。重复上述过程，完成全部时间段的电磁响应迭代计算：

Etl+1=acoeffEtl+bcoeffHtl+1/2；

Htl+3/2=ccoeffEtl+1+dcoeffHtl+1/2。（10）

式中：acoeff、bcoeff、ccoeff、dcoeff为相关系数；l为时刻序数。

由式（10）表征的迭代过程可以发现，需要将接地长导线的地下电场响应（初始场）作为计算初始条件带入到差分迭代中才能保证迭代的进行。接地长导线在均匀半空间下的阶跃电场表达式为[17]：

Ex=－I4π∫L-Lz^∫SymboleB@0rTE1e－u1zλu0J0（λR）dλdx′+

Λ（x，R2）－Λ（x，R1）；（11）

Ey=Λ（y，R1）－Λ（y，R2）；（12）

Λ（A1，A2）=

IA14πA2∫SymboleB@0u1y^1rTM1－z^u0rTE1e－u1zJ1（λA2）dλ。（13）

其中：

R=［（x－x′）2+y2］1/2；R1=［（x－L）2+y2］1/2；R2=［（x+L）2+y2］1/2。（14）

式中：I为长导线电流；2L为长导线长度；z^=iωμ；rTM1和rTE1为折射系数；u0、u1、y^1为rTE1与rTM1引入的相关系数，详见文献[14]；λ为波长；x′为源横坐标（沿长导线方向）。

为保证电磁迭代的稳定性，初始时刻的计算公式应表示为[1819]

t0=1.13μσh21。（15）

式中，h1为顶层第一个网格在z方向的高度。时间步长公式为

Δtn=α（μσmint6）12Δmin。（16）

式中：α的取值范围为0.1～0.2；σmin为最小电导率；Δmin为最小网格步长。

从式（16）可以看出，时间步长是一个随着时间推移逐渐变大的量，表明迭代后期计算的点之间的间隔越来越大。这与电磁波的衰减特性相符，在早期衰减速度快，而晚期衰减速度慢[20]。

2.4 初始场脉冲响应计算

因为采用电场与磁场的脉冲响应作为初始场，所以需要计算整数时刻的脉冲电场响应与半整数时刻的脉冲磁场响应[20]。本文通过电磁感应定律计算磁场脉冲响应：

SymbolQC@×E=－Bt。（17）

展开式（17）中的电场旋度，并利用中心差分近似：

ΔfΔww=f（w+Δw/2）－f（w－Δw/2）Δw。（18）

式中：f为关于变量w的函数；Δf为函数f的变化量；Δw为差分步长。

利用式（18）可以获得t0+Δt0/2时刻的脉冲磁场响应。以x方向脉冲磁场响应为例，表达式为

Bxt=Ey（i，j，k+1/2）－Ey（i，j，k－1/2）Δz－Ez（i，j+1/2，k）－Ez（i，j－1/2，k）Δy。（19）

式中：Bx为x方向的磁感应强度；Ey、Ez分别为y、z方向的电场强度；i、j、k分别为x、y、z方向的坐标序数；Δy、Δz分别为y、z方向的最小步长。

通过对全电流公式进行处理，可以获得电场脉冲响应。根据波数计算，当ωlt;105Hz时，位移电流远小于传导电流，可以忽略[17]，此时全电流定律为

SymbolQC@×B=μσE。（20）

同理，对磁场的旋度进行展开，两侧同时关于时间求导，以x方向为例，电场的脉冲响应公式为

Ext=Bz（i，j+1/2，k）t－Bz（i，j－1/2，k）tμσΔy－By（i，j，k+1/2）t－By（i，j，k－1/2）tμσΔz 。（21）

式中：By、Bz分别为y、z方向的磁感应强度；Ex为x方向的电场强度。

2.5 电磁场迭代公式推导

对麦克斯韦方程组中全电流定律与电磁感应定律进行处理，将电场与磁场在三维直角坐标系沿x、y、z方向展开，得到：

Hzy－Hyz=εExt+σEx；Hxz－Hzx=εEyt+σEy；Hyx－Hxy=εEzt+σEz。（22）

Ezy－Eyz=－μHxt；Exz－Ezx=－μHyt；Eyx－Exy=－μHzt。 （23）

式中：Hx、Hy、Hz分别为x、y、z方向的磁场强度；*表示卷积操作。

选取观察点（x，y，z），即（i+1/2，j，k）点在时间t=l+1/2的时间节点进行观察，以差商形式表示导数。通过全电流定律求解电场迭代方程，将式（9）带入式（22）中，并进行数学变换，整理得到电场迭代公式，以Ex为例（其余方向电场信号同理可推导），其表达式为

E（l+1）x（i+12，j，k）=2ΘH（l+12）z（i+12，j+12，k）Δyj－H（l+12）z（i+12，j－12，k）Δyj－H（l+12）y（i+12，j，k+12）Δzk+ H（l+12）y（i+12，j，k－12）Δzk+2∑Nn=1σ^（l+12）n（i+12，j，k）ξln（i+12，j，k）Θ+ΨElx（i+12，j，k）Θ。（24）

其中：

ξln（r）=

∑li=1Δt（i－1）［epnt（i－1）E（r，t（i－1））+epnt（i）E（r，t（i））］2； （25）

Θ=σSymboleB@（i+12，j，k）－

∑Nn=1ησSymboleB@（i+12，j，k）rnΔt（l）4+2εΔt（l+1）；（26）

Ψ=σSymboleB@（i+12，j，k）－

∑Nn=1ησSymboleB@（i+12，j，k）rnΔt（l）4+

2β（l+1）（i+12，j，k）－2εΔt（l+1）；（27）

β（l+1）（i+12，j，k）=－∑Nn=1σ^（12Δt（l））n（i+12，j，k）Δt（l）4； （28）

Δyj=Δyj－1+Δyj2；（29）

Δzk=Δzk－1+Δzk2。（30）

对电磁感应定律同样进行差分离散，可获得各磁场分量的计算公式；选取r=（i，j+1/2，k+1/2）和t=l时刻为磁场节点，以Hx为例，磁场迭代公式为

H（l+32）x（i，j+12，k+12）=

H（l+12）x（i，j+12，k+12）－Δt（l+1）μ（i，j+12，k+12）·

E（l+1）z（i，j+1，k+12）－E（l+1）z（i，j，k+12）Δyj－E（l+1）y（i，j+12，k+1）－E（l+1）y（i，j+12，k）Δzk。（31）

按照上述方法，还可以推导出Ey、Ez以及Hy的迭代公式。在推导Hz的迭代公式时，由SymbolQC@·H=0进行计算[21]。通过对上述公式进行逐步求解，就可以获得整个计算区域的电磁场变化情况，实现对感应极化效应的三维数值模拟。

3 典型模型电磁响应数值模拟

按照前面所述方法与流程，本文对均匀半空间模型、极化半空间模型与三维极化体模型的电磁响应进行数值模拟，并对模拟结果进行分析。

3.1 均匀半空间模型

均匀半空间模型模拟无极化均匀大地，用于验证数值模拟方法的可行性。本文设置计算区域x、y、z方向网格维度为201×201×75，最小网格步长为10 m，电导率取0.05 S/m。电性源的发射电流为50 A，长度为1 000 m。接收位置坐标为（0 m， 200 m， 0 m）。由于时域电性源电磁响应只在地面处存在解析解，因此接收位置设置在地面高度为0 m处，以便能与解析解进行对比[17]。接地导线垂直磁场时间导数的解析表达式为[4]

hzt=

Ids2πσμ0yρ53erf（θρ）－2π1/2θρ（3+2θ2ρ2）e－θ2ρ2。（32）

式中：ρ为径向距离（ρ2=x2+y2）；θ=μ0σ/4t。

均匀半空间模型的电磁响应数值模拟结果如图4所示。由图4可以看出，均匀半空间模型电磁响应的三维数值模拟结果与解析解几乎重合。为了对数值模拟结果进行量化评估，计算了数值模拟结果与解析解的相对误差：

e=E2－E1E2×100%。（33）

式中：E1为数值模拟结果；E2为解析解。

均匀半空间模型电磁响应数值模拟

相对误差如图5所示，相对误差均在10%以内，证明了本文提出的三维数值模拟方法具有可行性。

3.2 极化半空间模型

极化半空间模型即计算区域为均匀极化介质。设置地下介质电导率为0.05 S/m，极化率为0.3，时间常数为0.05 s，频散系数为0.5。利用本文的数值模拟方法计算该模型的电磁响应。

图6为极化半空间模型的电磁响应数值模拟结果。由图6可以看出，极化半空间模型早期的电磁响应曲线与电导率为0.05 S/m的均匀半空间模型电磁响应的解析解基本重合，但是极化效应导致衰减曲线快速衰减，在4.5 ms发生符号反转现象，负响应最大绝对值为0.10 μV。

极化半空间模型电磁响应曲线反号现象的发生证明该区域为极化介质，为本文数值模拟方法的正确性提供了依据。

3.3 三维极化体模型

对三维地下极化体的电磁响应进行模拟，设置大地模型尺寸为8 010 m×8 010 m×2 175 m，围岩电导率为0.01 S/m；极化体大小为800 m×800 m×800 m，极化体中心位置坐标为（0 m，0 m，-450 m），极化率为0.8，电导率为0.05 S/m；时间常数与频散系数分别设置为0.001 s与0.5。

为了分析接收位置对极化效应的影响，本文设置不同的接收位置x=0 m，y=200、300、400、450、600 m，z=0 m，电磁响应如图7所示。从图7中可以看出：随着偏移距离的增加，接收位置由三维极化体上方区域（y=200 m）向三维极化体边缘位置（y=400 m）移动，电磁响应反号出现的时间延迟；当接收位置超过极化体所在区域时（y=450 m），电磁响应反号现象出现时间发生明显延迟；当接收位置距离极化体较远时（y=600 m），电磁响应不出现负响应。

图8为三维极化体模型电磁响应切片图。其中，可以清楚观察到极化体的位置（红色矩形框区域）。对比图8a、b的电磁响应情况，可以观察到极化体中的电磁响应在晚期发生了符号反转，证明该区域产生了极化效应，成功探测到异常体。三维极化体模型电磁响应数值模拟的结果中出现符号反转现象，进一步证明了本文提出的数值模拟算法的正确性。

4 结论

1）本文提出了一种基于时域有限差分法的电性源感应极化效应三维数值模拟方法，实现了电性源电磁响应的三维数值模拟计算。

2）本文对均匀半空间模型进行数值模拟，电磁响应与解析解基本重合，且相对误差小于10%，验证了算法的可行性。

3）极化半空间模型和三维极化体模型的电磁响应均出现符号反转现象，与极化理论相符。

综上，本文实现了极化介质Cole-Cole模型的时域电性源感应极化效应三维数值模拟，为实际野外探测的实测数据处理奠定了基础。本文仅模拟了三种典型模型的电磁响应，但实际地质情况会更为复杂，该方法可能存在一定局限性。可继续优化有限差分算法、调整参数，如减小时间步长、增加网格精度等，提高算法准确性和计算效率。

参考文献（References）：

[1] 余丽秀，孙亚光.矿产资源综合利用与地质调查关系探讨[J].矿产保护与利用，2008（5）：1013.

Yu Lixiu， Sun Yaguang. Comprehensive Utilization of Mineral Resources and Geological Survey Relationships[J]. Mineral Protection and Utilization， 2008（5）：1013.

[2] 杨沈生.我国矿产资源现状与应对措施[J].科技创新导报，2010（12）：70.

Yang Shensheng. The Current Status and Response Measures of Mineral Resources in China[J]. Science and Technology Innovation Guidance， 2010（12）：70.

[3] 滕吉文.中国地球物理学研究面临的机遇、发展空间和时代的挑战[J].地球物理学进展，2007，22（4）：11011112.

Teng Jiwen. Opportunities， Development Space， and the Challenges of the Times in Chinese Earth Physics Research[J]. Earth Physics Progress， 2007， 22（4）： 11011112.

[4] 纳比吉安.勘查地球物理电磁法：第1卷：理论[M]. 赵经祥，王艳君，译. 北京：地质出版社，1992.

Misac N. Electromagnetic Methods in Applied Geophysics： Volume 1： Theory [M]. Translated by Zhao Jingxiang， Wang Yanjun. Beijing： Geological Publishing House， 1992.

[5] 雷松达，王显祥，刘遂明.复杂地质条件下浅海区电性源瞬变电磁法三维响应特征[J].吉林大学学报（地球科学版），2024，54（3）：10161030.

Lei Songda， Wang Xianxiang， Liu Suiming. Three-Dimensional Response Characteristics of Transient Electromagnetic Method with Electrical Sources in Shallow Sea Areas Under Complex Geological Conditions[J]. Journal of Jilin University （Earth Science Edition）， 2024， 54（3）： 10161030.

[6] 齐彦福，李貅，孙乃泉，等.电性源短偏移距瞬变电磁地形影响特征分析 [J].吉林大学学报（地球科学版），2022，52（1）：247260.

Qi Yanfu， Li Xiu， Sun Naiquan， et al. Analysis of Influence Characteristics of Topography on Grounded-Source Short-Offset Transient Electromagnetic Responses[J]. Journal of Jilin University （Earth Science Edition）， 2022， 52（1）： 247260.

[7] Debye P， Falkenhagen H. Dispersion of the Conductivity and Dielectric Constants of Strong Electrolytes[J]. Physikalische Zeitschrift， 1928， 29： 121132.

[8] Dias C A. Analytical Model for a Polarizable Medium at Radio and Lower Frequencies[J]. Journal of Geophysical Research Atmospheres， 1972， 77（26）： 49454956.

[9] Cole K S， Cole R H. Dispersion and Absorption in Dielectrics： I： Alternating Current Characteristics[J]. The Journal of Chemical Physics， 1941， 9（4）： 341351.

[10] 葛德彪，闫玉波.电磁波时域有限差分方法[M]. 西安：西安电子科技大学出版社，2002.

Ge Debiao， Yan Yubo. Electromagnetic Wave Time Domain Finite Difference Method[M]. Xi’an： Xi’an University of Electronic Science and Technology Press， 2002.

[11] 赵友超，张军，范涛，等. 地井瞬变电磁三维响应特征分析与异常体快速定位方法研究[J]. 物探与化探，2022，46（2）：383391.

Zhao Youchao， Zhang Jun， Fan Tao， et al. Analysis of 3D Ground Borehole TEM Response Characteristics and Rapid Positioning Method for Anomalous Bodies[J]. Geophysical and Geochemical Exploration， 2022， 46（2）： 383391.

[12] 周钟航，张莹莹.山峰对电性源地面瞬变电磁响应的影响及校正方法[J]. 物探与化探，2023，47（5）：12361249.

Zhou Zhonghang， Zhang Yingying. Correction of the Influence of Mountains on Grounded Source Transient Electromagnetic Responses[J]. Geophysical and Geochemical Exploration， 2023， 47（5）：12361249.

[13] 张鑫崇，殷长春.基于Debye模型的时间域航空电磁激电效应正演模拟[J].吉林大学学报（地球科学版），2023，53（5）：15731581.

Zhang Xinchong， Yin Changchun. Forward Modeling of Airborne Electromagnetic Induced Polarization Effect in Time-Domain Based on Debye Model[J]. Journal of Jilin University （Earth Science Edition）， 2023， 53 （5）： 15731581.

[14] 左信，庄馨，许鋆，等. 基于线性规划的分数阶系统有理函数逼近方法[J].信息与控制，2014，43（6）：697704， 710.

Zuo Xin， Zhuang Xin， Xu Jun，et al. Active Function Approximation Method Based on a Linear Planning-Based Segmentation and Order System [J]. Information and Control， 2014， 43（6）：697704， 710.

[15] 孙凤杰，周启明.时域有限差分（FDTD）法发展综述[C]//2009系统仿真技术及其应用学术会议（CCSSTA'2009）论文集.合肥：[s. n.]，2009：5.

Sun Fengjie， Zhou Qiming. Summary of the Development of Time Domain Finite Difference[C]//2009 System Simulation Technology and Its Application Academic Conference Papers. Hefei： [s. n.]， 2009：5.

[16] Yee K S. Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media[J]. IEEE Transactions on Antennas amp; Propagation， 1966， 14（5）： 302307.

[17] 黎东升. 时域地空电性源的三维电磁数值模拟及噪声抑制方法研究[D]. 长春：吉林大学， 2016.

Li Dongsheng. Study on the Three-Dimensional Electromagnetic Value Simulation and Noise Suppression Method of Time Domain Ground-Air Electrical Sources[D]. Changchun： Jilin University， 2016.

[18] Wang T， Hohmann G W. A Finite-Difference， Time-Domain Solution for 3Dimensional Electromagnetic Modeling [J]. Geophysics， 1993， 58（1）： 797809.

[19] Commer M， Newman G. A Parallel Finite-Difference Approach for 3D Transient Electromagnetic Modeling with Galvanic Sources [J]. Geophysics， 2004， 69（10）： 11921202.

[20] 关珊珊. 基于GPU的三维有限差分直升机瞬变电磁响应并行计算[D]. 长春：吉林大学，2012.

Guan Shanshan. The GPU-Based Three-Dimensional Finite Differential Helicopter Transient Electromagnetic Response Parallel Calculation[D]. Changchun： Jilin University， 2012.

[21] 吴燕琪. 极化介质的三维时域电磁响应数值模拟与智能识别[D]. 长春：吉林大学， 2021.

Wu Yanqi. The Three-Dimensional Time Domain Electromagnetic Response of Polarization Medium of Electromagnetic Response and Intelligent Identification[D]. Changchun： Jilin University， 2021.