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

關(guān)鍵詞：電性源；Cole-Cole模型；有理逼近；時(shí)域有限差分；極化效應(yīng)

doi：10.13278/j.cnki.jjuese.20230125

中圖分類號(hào)：P631

文獻(xiàn)標(biāo)志碼：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 引言

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

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

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

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

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

1 感應(yīng)極化效應(yīng)建模

1.1 Cole-Cole模型頻域表達(dá)

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

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

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

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

極化介質(zhì)中的電導(dǎo)率是時(shí)變函數(shù)，為對(duì)時(shí)域電磁響應(yīng)進(jìn)行數(shù)值模擬，需要將分?jǐn)?shù)階Cole-Cole模型電導(dǎo)率的時(shí)域表達(dá)方程代入歐姆定律中進(jìn)行卷積計(jì)算；因此需要對(duì)Cole-Cole模型表達(dá)式進(jìn)行頻時(shí)轉(zhuǎn)換。然而，與整數(shù)階系統(tǒng)不同，分?jǐn)?shù)階系統(tǒng)在本質(zhì)上是無窮維系統(tǒng)，對(duì)其直接進(jìn)行計(jì)算非常困難；因此采用整數(shù)階有理函數(shù)逼近法對(duì)任意分?jǐn)?shù)階Cole-Cole模型進(jìn)行有理化處理[14]，得到三維空間坐標(biāo)r位置處Cole-Cole模型電導(dǎo)率時(shí)域表達(dá)式為：

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

式中：t為時(shí)間；u（t）為階躍函數(shù)；rn、pn為Cole-Cole模型的有理函數(shù)逼近系數(shù)；N為擬合階次。

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

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

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

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

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

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

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

2.2 控制方程選取

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

對(duì)麥克斯韋方程組的時(shí)域表達(dá)式進(jìn)行拉氏變換，獲得頻域表達(dá)式為：

式中：E為電場(chǎng)強(qiáng)度；μ為磁導(dǎo)率；H為磁場(chǎng)強(qiáng)度；σ為電導(dǎo)率；ε為介電常數(shù)；B為磁感應(yīng)強(qiáng)度；D為電位移。

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

對(duì)式（7）進(jìn)行拉式逆變換后得到

因?yàn)闃O化介質(zhì)中電導(dǎo)率為時(shí)變函數(shù)，所以歐姆定律由乘積形式改變?yōu)榫矸e形式：

J（t）=∫t－

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

2.3 電性源階躍電場(chǎng)表達(dá)

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

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

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

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

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

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

t0=1.13μσh21。（15）

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

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

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

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

2.4 初始場(chǎng)脈沖響應(yīng)計(jì)算

因?yàn)椴捎秒妶?chǎng)與磁場(chǎng)的脈沖響應(yīng)作為初始場(chǎng)，所以需要計(jì)算整數(shù)時(shí)刻的脈沖電場(chǎng)響應(yīng)與半整數(shù)時(shí)刻的脈沖磁場(chǎng)響應(yīng)[20]。本文通過電磁感應(yīng)定律計(jì)算磁場(chǎng)脈沖響應(yīng)：

SymbolQC@×E=－Bt。（17）

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

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

式中：f為關(guān)于變量w的函數(shù)；Δf為函數(shù)f的變化量；Δw為差分步長(zhǎng)。

利用式（18）可以獲得t0+Δt0/2時(shí)刻的脈沖磁場(chǎng)響應(yīng)。以x方向脈沖磁場(chǎng)響應(yīng)為例，表達(dá)式為

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方向的磁感應(yīng)強(qiáng)度；Ey、Ez分別為y、z方向的電場(chǎng)強(qiáng)度；i、j、k分別為x、y、z方向的坐標(biāo)序數(shù)；Δy、Δz分別為y、z方向的最小步長(zhǎng)。

通過對(duì)全電流公式進(jìn)行處理，可以獲得電場(chǎng)脈沖響應(yīng)。根據(jù)波數(shù)計(jì)算，當(dāng)ωlt;105Hz時(shí)，位移電流遠(yuǎn)小于傳導(dǎo)電流，可以忽略[17]，此時(shí)全電流定律為

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

同理，對(duì)磁場(chǎng)的旋度進(jìn)行展開，兩側(cè)同時(shí)關(guān)于時(shí)間求導(dǎo)，以x方向?yàn)槔?，電?chǎng)的脈沖響應(yīng)公式為

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方向的磁感應(yīng)強(qiáng)度；Ex為x方向的電場(chǎng)強(qiáng)度。

2.5 電磁場(chǎng)迭代公式推導(dǎo)

對(duì)麥克斯韋方程組中全電流定律與電磁感應(yīng)定律進(jìn)行處理，將電場(chǎng)與磁場(chǎng)在三維直角坐標(biāo)系沿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方向的磁場(chǎng)強(qiáng)度；*表示卷積操作。

選取觀察點(diǎn)（x，y，z），即（i+1/2，j，k）點(diǎn)在時(shí)間t=l+1/2的時(shí)間節(jié)點(diǎn)進(jìn)行觀察，以差商形式表示導(dǎo)數(shù)。通過全電流定律求解電場(chǎng)迭代方程，將式（9）帶入式（22）中，并進(jìn)行數(shù)學(xué)變換，整理得到電場(chǎng)迭代公式，以Ex為例（其余方向電場(chǎng)信號(hào)同理可推導(dǎo)），其表達(dá)式為

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）

對(duì)電磁感應(yīng)定律同樣進(jìn)行差分離散，可獲得各磁場(chǎng)分量的計(jì)算公式；選取r=（i，j+1/2，k+1/2）和t=l時(shí)刻為磁場(chǎng)節(jié)點(diǎn)，以Hx為例，磁場(chǎng)迭代公式為

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）

按照上述方法，還可以推導(dǎo)出Ey、Ez以及Hy的迭代公式。在推導(dǎo)Hz的迭代公式時(shí)，由SymbolQC@·H=0進(jìn)行計(jì)算[21]。通過對(duì)上述公式進(jìn)行逐步求解，就可以獲得整個(gè)計(jì)算區(qū)域的電磁場(chǎng)變化情況，實(shí)現(xiàn)對(duì)感應(yīng)極化效應(yīng)的三維數(shù)值模擬。

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

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

3.1 均勻半空間模型

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

hzt=

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

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

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

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

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

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

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

3.2 極化半空間模型

極化半空間模型即計(jì)算區(qū)域?yàn)榫鶆驑O化介質(zhì)。設(shè)置地下介質(zhì)電導(dǎo)率為0.05 S/m，極化率為0.3，時(shí)間常數(shù)為0.05 s，頻散系數(shù)為0.5。利用本文的數(shù)值模擬方法計(jì)算該模型的電磁響應(yīng)。

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

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

3.3 三維極化體模型

對(duì)三維地下極化體的電磁響應(yīng)進(jìn)行模擬，設(shè)置大地模型尺寸為8 010 m×8 010 m×2 175 m，圍巖電導(dǎo)率為0.01 S/m；極化體大小為800 m×800 m×800 m，極化體中心位置坐標(biāo)為（0 m，0 m，-450 m），極化率為0.8，電導(dǎo)率為0.05 S/m；時(shí)間常數(shù)與頻散系數(shù)分別設(shè)置為0.001 s與0.5。

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

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

4 結(jié)論

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

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

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

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

參考文獻(xiàn)（References）：

[1] 余麗秀，孫亞光.礦產(chǎn)資源綜合利用與地質(zhì)調(diào)查關(guān)系探討[J].礦產(chǎn)保護(hù)與利用，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] 楊沈生.我國礦產(chǎn)資源現(xiàn)狀與應(yīng)對(duì)措施[J].科技創(chuàng)新導(dǎo)報(bào)，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] 滕吉文.中國地球物理學(xué)研究面臨的機(jī)遇、發(fā)展空間和時(shí)代的挑戰(zhàn)[J].地球物理學(xué)進(jìn)展，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)祥，王艷君，譯. 北京：地質(zhì)出版社，1992.

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

[5] 雷松達(dá)，王顯祥，劉遂明.復(fù)雜地質(zhì)條件下淺海區(qū)電性源瞬變電磁法三維響應(yīng)特征[J].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），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].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），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] 葛德彪，閆玉波.電磁波時(shí)域有限差分方法[M]. 西安：西安電子科技大學(xué)出版社，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] 趙友超，張軍，范濤，等. 地井瞬變電磁三維響應(yīng)特征分析與異常體快速定位方法研究[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] 周鐘航，張瑩瑩.山峰對(duì)電性源地面瞬變電磁響應(yīng)的影響及校正方法[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] 張?chǎng)纬?，殷長(zhǎng)春.基于Debye模型的時(shí)間域航空電磁激電效應(yīng)正演模擬[J].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），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ī)劃的分?jǐn)?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] 孫鳳杰，周啟明.時(shí)域有限差分（FDTD）法發(fā)展綜述[C]//2009系統(tǒng)仿真技術(shù)及其應(yīng)用學(xué)術(shù)會(huì)議（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í)域地空電性源的三維電磁數(shù)值模擬及噪聲抑制方法研究[D]. 長(zhǎng)春：吉林大學(xué)， 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] 關(guān)珊珊. 基于GPU的三維有限差分直升機(jī)瞬變電磁響應(yīng)并行計(jì)算[D]. 長(zhǎng)春：吉林大學(xué)，2012.

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

[21] 吳燕琪. 極化介質(zhì)的三維時(shí)域電磁響應(yīng)數(shù)值模擬與智能識(shí)別[D]. 長(zhǎng)春：吉林大學(xué)， 2021.

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