基于時域有限差分法的電性源感應極化效應三維數(shù)值模擬

摘要：時域電性源電磁探測法是一種可以有效快速探測礦產資源的方法。極化效應會導致電性源電磁響應快速衰減，發(fā)生符號反轉現(xiàn)象。本文采用Cole-Cole模型描述極化效應，利用整數(shù)階有理逼近算法實現(xiàn)任意分數(shù)階Cole-Cole模型有理化；采用Yee氏網(wǎng)格對仿真區(qū)域進行剖分，基于時域有限差分（finite-difference time-domain，F(xiàn)DTD）法實現(xiàn)電性源感應極化效應三維數(shù)值模擬。本文對三種典型模型（均勻半空間模型、極化半空間模型和三維極化體模型）的電磁響應進行數(shù)值模擬，結果表明：均勻半空間模型的電磁響應與解析解基本吻合，并且相對誤差小于10%，證明了三維數(shù)值模擬方法的正確性；極化半空間模型與三維極化體模型的電磁響應均在晚期出現(xiàn)了負響應，與極化理論結果相符合。

關鍵詞：電性源；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 引言

我國是礦產資源大國，同時也是世界上最大的礦產消費國之一。我國礦產資源總量大、種類豐富，但開采情況卻不容樂觀，尤其是鐵、錳、銅等大宗礦產供應明顯不足。資源供需矛盾的日益突出一定程度上阻礙了社會經(jīng)濟的發(fā)展[13]。

電性源時域電磁法是一種能夠有效快速探測礦產資源的方法[46]。該方法的操作步驟如下：首先，利用電性源（接地長導線）發(fā)射一次脈沖場信號（一次場），在一次場的激勵下，地下低阻介質中產生感應電流進而激起時變的電磁場（即二次場）；然后，使用接收線圈接收返回的二次場信號，通過分析二次場中包含的豐富的地下介質電性信息獲得地下深部構造。

極化效應是一種重要的電化學效應，主要存在于侵染礦、硫化礦、金屬礦等介質中。其原理是當有外加電場激勵時，介質中正負電荷會發(fā)生定向運動，在界面兩側逐漸積累形成雙電層；將電場撤掉后，介質中的正負電荷重新恢復到原始狀態(tài)，形成了與外加電場方向相反的位移電流。極化率是金屬礦的重要物性參數(shù)，基于極化效應可以實現(xiàn)對金屬礦等介質的探測。國內外學者提出許多數(shù)學模型用于描述地下介質的極化特性，如Debye模型[7]、Dias模型[8]、Cole-Cole模型[9]等，其中最為典型的就是Cole-Cole模型；因此，本文在Cole-Cole模型基礎上進行探究。

時域電磁數(shù)值模擬方法有很多種，例如時域有限差分（finite-difference time-domain， FDTD）法、時域有限體積法、時域積分方程法等。在這些算法中，F(xiàn)DTD算法具有簡明直觀的特點，非常適合寬頻計算，被廣泛應用于三維時域電磁數(shù)值模擬[1012]。

本文針對電性源激勵下極化介質的感應極化效應進行三維數(shù)值模擬，使用Cole-Cole模型描述地下極化介質的極化特性，采用整數(shù)階有理逼近方法對任意分數(shù)階的Cole-Cole模型電導率表達式進行有理化逼近，獲得Cole-Cole電導率時域表達方程。采用Yee氏網(wǎng)格對仿真區(qū)域進行剖分，將Cole-Cole電導率時域表達式代入麥克斯韋方程組中，推導電性源初始場公式與電磁迭代公式，并利用時域有限差分法計算電性源激勵下極化介質的電磁響應。對均勻半空間模型、極化半空間模型與三維極化體模型的電磁響應進行數(shù)值模擬，以驗證數(shù)值模擬算法的可行性，為實際野外探測的實測數(shù)據(jù)解釋奠定基礎。

1 感應極化效應建模

1.1 Cole-Cole模型頻域表達

極化效應與電容的充放電效應類似，表現(xiàn)為大地電導率的頻散特性。Cole-Cole模型電導率頻域表達式為

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

式中：ω為角頻率；σSymboleB@為無窮頻率極限下的電導率；η為極化率；τ為時間常數(shù)；c為頻散系數(shù)。η與c取值范圍均為［0，1］。當c取值為1時，Cole-Cole模型簡化為Debye模型[13]。

1.2 Cole-Cole模型有理函數(shù)逼近

極化介質中的電導率是時變函數(shù)，為對時域電磁響應進行數(shù)值模擬，需要將分數(shù)階Cole-Cole模型電導率的時域表達方程代入歐姆定律中進行卷積計算；因此需要對Cole-Cole模型表達式進行頻時轉換。然而，與整數(shù)階系統(tǒng)不同，分數(shù)階系統(tǒng)在本質上是無窮維系統(tǒng)，對其直接進行計算非常困難；因此采用整數(shù)階有理函數(shù)逼近法對任意分數(shù)階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）為階躍函數(shù)；rn、pn為Cole-Cole模型的有理函數(shù)逼近系數(shù)；N為擬合階次。

為驗證分數(shù)階Cole-Cole模型的逼近效果，對參數(shù)為c=0.6、η=0.5、τ=0.01 s的模型進行逼近，并將N依次設置為7、8、9、10。不同擬合階次下擬合電導率與實際電導率的幅值與相角伯德圖如圖1所示，可以看到不同擬合階次均取得了較好的擬合效果。

為直觀觀察擬合誤差與擬合階次的關系，分別計算幅值與相位有理逼近的誤差，結果如圖2所示。從圖2中可以看出，隨著擬合階次的逐漸提高，擬合準確度也逐漸提升。但是較大的擬合階次會導致e指數(shù)項的增多，進而提高了計算量、降低了計算速度。因此，在保證較高擬合精度的前提下應盡量選擇相對低階的擬合電導率進行后續(xù)的迭代操作。

2 三維時域電磁響應數(shù)值模擬

2.1 仿真區(qū)域Yee氏網(wǎng)格剖分

為了對電磁響應問題進行求解，Yee提出了時域有限差分方法，通過對電場磁場交替采樣，利用變微分為差分的方法構建迭代方程，實現(xiàn)對電磁迭代過程的表征并得到計算結果[1516]。本文采用對仿真區(qū)域進行Yee氏網(wǎng)格剖分。

如圖3所示，ex、ey與ez分別x、y與z方向的電場分量，分布在Yee元胞的棱中心點；hx、hy與hz分別為x、y與z方向的磁場分量，分布在Yee元胞的面中心點；二者在空間中交叉分布，彼此環(huán)繞。其排布規(guī)律與取樣方式符合法拉第感應定律與安培

環(huán)路定律，能夠恰當描述電磁波在空間的傳播過程，且便于麥克斯韋方程的差分計算。此外，電場與磁場在時間上交替取值，彼此相差二分之一的時間步長。電場信號在整數(shù)時刻t采樣，磁場信號在半整數(shù)時刻t+1/2采樣。通過對麥克斯韋旋度方程組進行差分離散，構成顯示差分格式，從而迭代出各個時刻的電磁場情況。

2.2 控制方程選取

麥克斯韋方程組全面總結了電磁場的變化規(guī)律，是宏觀電磁場理論的基礎，它表征著變化的電場與變化的磁場之間相互激勵的密切聯(lián)系。利用麥克斯韋方程，并帶入初始值與邊界條件，就可以對自然界中的電磁傳播現(xiàn)象進行數(shù)值模擬。

對麥克斯韋方程組的時域表達式進行拉氏變換，獲得頻域表達式為：

式中：E為電場強度；μ為磁導率；H為磁場強度；σ為電導率；ε為介電常數(shù)；B為磁感應強度；D為電位移。

考慮極化效應的影響，將Cole-Cole模型頻域表達式（式（1））代入到式（4）中，全電流定律公式可變?yōu)?/p>

對式（7）進行拉式逆變換后得到

因為極化介質中電導率為時變函數(shù)，所以歐姆定律由乘積形式改變?yōu)榫矸e形式：

J（t）=∫t－

式中，J為電流密度。式（8）與式（3）、式（5）的時域微分形式共同構成控制方程組。

2.3 電性源階躍電場表達

對計算區(qū)域進行Yee氏網(wǎng)格剖分后，進行電磁響應迭代計算，迭代的基本思想為：利用前一時刻的電場值與前半整數(shù)時刻的磁場值推導下一時刻的電場值；利用上一步求出的電場值與前半整數(shù)時刻的磁場值求解下半整數(shù)時刻的磁場值。重復上述過程，完成全部時間段的電磁響應迭代計算：

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

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

式中：acoeff、bcoeff、ccoeff、dcoeff為相關系數(shù)；l為時刻序數(shù)。

由式（10）表征的迭代過程可以發(fā)現(xiàn)，需要將接地長導線的地下電場響應（初始場）作為計算初始條件帶入到差分迭代中才能保證迭代的進行。接地長導線在均勻半空間下的階躍電場表達式為[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為折射系數(shù)；u0、u1、y^1為rTE1與rTM1引入的相關系數(shù)，詳見文獻[14]；λ為波長；x′為源橫坐標（沿長導線方向）。

為保證電磁迭代的穩(wěn)定性，初始時刻的計算公式應表示為[1819]

t0=1.13μσh21。（15）

式中，h1為頂層第一個網(wǎng)格在z方向的高度。時間步長公式為

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

式中：α的取值范圍為0.1～0.2；σmin為最小電導率；Δmin為最小網(wǎng)格步長。

從式（16）可以看出，時間步長是一個隨著時間推移逐漸變大的量，表明迭代后期計算的點之間的間隔越來越大。這與電磁波的衰減特性相符，在早期衰減速度快，而晚期衰減速度慢[20]。

2.4 初始場脈沖響應計算

因為采用電場與磁場的脈沖響應作為初始場，所以需要計算整數(shù)時刻的脈沖電場響應與半整數(shù)時刻的脈沖磁場響應[20]。本文通過電磁感應定律計算磁場脈沖響應：

SymbolQC@×E=－Bt。（17）

展開式（17）中的電場旋度，并利用中心差分近似：

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

式中：f為關于變量w的函數(shù)；Δf為函數(shù)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方向的坐標序數(shù)；Δy、Δz分別為y、z方向的最小步長。

通過對全電流公式進行處理，可以獲得電場脈沖響應。根據(jù)波數(shù)計算，當ω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的時間節(jié)點進行觀察，以差商形式表示導數(shù)。通過全電流定律求解電場迭代方程，將式（9）帶入式（22）中，并進行數(shù)學變換，整理得到電場迭代公式，以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時刻為磁場節(jié)點，以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]。通過對上述公式進行逐步求解，就可以獲得整個計算區(qū)域的電磁場變化情況，實現(xiàn)對感應極化效應的三維數(shù)值模擬。

3 典型模型電磁響應數(shù)值模擬

按照前面所述方法與流程，本文對均勻半空間模型、極化半空間模型與三維極化體模型的電磁響應進行數(shù)值模擬，并對模擬結果進行分析。

3.1 均勻半空間模型

均勻半空間模型模擬無極化均勻大地，用于驗證數(shù)值模擬方法的可行性。本文設置計算區(qū)域x、y、z方向網(wǎng)格維度為201×201×75，最小網(wǎng)格步長為10 m，電導率取0.05 S/m。電性源的發(fā)射電流為50 A，長度為1 000 m。接收位置坐標為（0 m， 200 m， 0 m）。由于時域電性源電磁響應只在地面處存在解析解，因此接收位置設置在地面高度為0 m處，以便能與解析解進行對比[17]。接地導線垂直磁場時間導數(shù)的解析表達式為[4]

hzt=

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

式中：ρ為徑向距離（ρ2=x2+y2）；θ=μ0σ/4t。

均勻半空間模型的電磁響應數(shù)值模擬結果如圖4所示。由圖4可以看出，均勻半空間模型電磁響應的三維數(shù)值模擬結果與解析解幾乎重合。為了對數(shù)值模擬結果進行量化評估，計算了數(shù)值模擬結果與解析解的相對誤差：

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

式中：E1為數(shù)值模擬結果；E2為解析解。

均勻半空間模型電磁響應數(shù)值模擬

相對誤差如圖5所示，相對誤差均在10%以內，證明了本文提出的三維數(shù)值模擬方法具有可行性。

3.2 極化半空間模型

極化半空間模型即計算區(qū)域為均勻極化介質。設置地下介質電導率為0.05 S/m，極化率為0.3，時間常數(shù)為0.05 s，頻散系數(shù)為0.5。利用本文的數(shù)值模擬方法計算該模型的電磁響應。

圖6為極化半空間模型的電磁響應數(shù)值模擬結果。由圖6可以看出，極化半空間模型早期的電磁響應曲線與電導率為0.05 S/m的均勻半空間模型電磁響應的解析解基本重合，但是極化效應導致衰減曲線快速衰減，在4.5 ms發(fā)生符號反轉現(xiàn)象，負響應最大絕對值為0.10 μV。

極化半空間模型電磁響應曲線反號現(xiàn)象的發(fā)生證明該區(qū)域為極化介質，為本文數(shù)值模擬方法的正確性提供了依據(jù)。

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；時間常數(shù)與頻散系數(shù)分別設置為0.001 s與0.5。

為了分析接收位置對極化效應的影響，本文設置不同的接收位置x=0 m，y=200、300、400、450、600 m，z=0 m，電磁響應如圖7所示。從圖7中可以看出：隨著偏移距離的增加，接收位置由三維極化體上方區(qū)域（y=200 m）向三維極化體邊緣位置（y=400 m）移動，電磁響應反號出現(xiàn)的時間延遲；當接收位置超過極化體所在區(qū)域時（y=450 m），電磁響應反號現(xiàn)象出現(xiàn)時間發(fā)生明顯延遲；當接收位置距離極化體較遠時（y=600 m），電磁響應不出現(xiàn)負響應。

圖8為三維極化體模型電磁響應切片圖。其中，可以清楚觀察到極化體的位置（紅色矩形框區(qū)域）。對比圖8a、b的電磁響應情況，可以觀察到極化體中的電磁響應在晚期發(fā)生了符號反轉，證明該區(qū)域產生了極化效應，成功探測到異常體。三維極化體模型電磁響應數(shù)值模擬的結果中出現(xiàn)符號反轉現(xiàn)象，進一步證明了本文提出的數(shù)值模擬算法的正確性。

4 結論

1）本文提出了一種基于時域有限差分法的電性源感應極化效應三維數(shù)值模擬方法，實現(xiàn)了電性源電磁響應的三維數(shù)值模擬計算。

2）本文對均勻半空間模型進行數(shù)值模擬，電磁響應與解析解基本重合，且相對誤差小于10%，驗證了算法的可行性。

3）極化半空間模型和三維極化體模型的電磁響應均出現(xiàn)符號反轉現(xiàn)象，與極化理論相符。

綜上，本文實現(xiàn)了極化介質Cole-Cole模型的時域電性源感應極化效應三維數(shù)值模擬，為實際野外探測的實測數(shù)據(jù)處理奠定了基礎。本文僅模擬了三種典型模型的電磁響應，但實際地質情況會更為復雜，該方法可能存在一定局限性。可繼續(xù)優(yōu)化有限差分算法、調整參數(shù)，如減小時間步長、增加網(wǎng)格精度等，提高算法準確性和計算效率。

參考文獻（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] 楊沈生.我國礦產資源現(xiàn)狀與應對措施[J].科技創(chuàng)新導報，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] 滕吉文.中國地球物理學研究面臨的機遇、發(fā)展空間和時代的挑戰(zhàn)[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]. 趙經(jīng)祥，王艷君，譯. 北京：地質出版社，1992.

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

[5] 雷松達，王顯祥，劉遂明.復雜地質條件下淺海區(qū)電性源瞬變電磁法三維響應特征[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] 左信，莊馨，許鋆，等. 基于線性規(guī)劃的分數(shù)階系統(tǒng)有理函數(shù)逼近方法[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）法發(fā)展綜述[C]//2009系統(tǒng)仿真技術及其應用學術會議（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] 黎東升. 時域地空電性源的三維電磁數(shù)值模擬及噪聲抑制方法研究[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] 吳燕琪. 極化介質的三維時域電磁響應數(shù)值模擬與智能識別[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.