基于改進Padé 逼近法的含激發極化效應瞬變電磁三維正演

關鍵詞： 激發極化效應，Padé 級數逼近，瞬變電磁法三維正演，Cole?Cole 模型中圖分類號：P631 文獻標識碼：A DOI：10. 13810/j. cnki. issn. 1000?7210. 20240265

3D forward modeling of transient electromagnetic method considering induced polarization effect based on improved Padé approximation method

YANG Wei1，CAI Panpan1，LI Qingkun1，LIU Lihuang2，HU Xiangyun2，LIU Yajun ^3，4_.

（1. Research Institute of Exploration and Development，SINOPEC Northwest Oilfield Company，Urumqi，Xinjiang 830011， China；2. School of Geophysics and Geomatics，China University of Geosciences，Wuhan，Hubei 430074，China；3. Badong Na? tional Observation and Research Station of Geohazards，China University of Geosciences，Wuhan，Hubei 430074，China；4. Hu? bei Key Laboratory of Resources and Eco?Environment Geology，Hubei Geological Bureau，Wuhan，Hubei ，China）

Abstract：The exploration targets such as metallic minerals，sulfides，and oil and gas resources often exhibit a strong induced polarization（IP）effect，which causes distortion of transient electromagnetic（TEM）signals. Neglecting the IP effect frequently leads to difficulties in fitting TEM data and even results in erroneous inter? pretations. This study employs the time?domain finite volume method（FVM）and approximates the Cole?Cole model by using a Padé series to achieve 3D forward modeling of the TEM method considering the IP effect. Traditional methods for selecting the central frequency in the Padé series rely on qualitative analysis，which makes the forward modeling results highly sensitive to changes in the frequency?dependent parameter. This sen sitivity limits the application of Padé approximation in regions with significant variability in the frequency ?depen? dent parameter. To address this limitation，this study quantitatively evaluates the errors introduced during the Padé approximation process and proposes an adaptive selection algorithm based on error information. Multiple numerical experiments demonstrate the accuracy and effectiveness of the proposed algorithm . Finally，the dif? ferences in TEM responses between polarized and non?polarized media models are analyzed，and the charging

and discharging processes of polarization anomaly bodies in TEM fields are simulated using an instance model. Keywords：induced polarization effect，Padé approximation，3D forward modeling of transient electromagnetic method，Cole?Cole model

楊威，蔡盼盼，李清坤等 . 基于改進 Padé 逼近法的含激發極化效應瞬變電磁三維正演［J］. 石油地球物理勘探， 2025，60（3）：794?805.

YANG Wei，CAI Panpan，LI Qingkun，et al. 3D forward modeling of transient electromagnetic method con? sidering induced polarization effect based on improved Padé approximation method［J］. Oil Geophysical Prospecting，2025，60（3）：794?805.

0 引言

瞬 變 電 磁 法 （Transient ElectromagneticMethod， TEM）是一種利用人工源激發時域電磁場的電磁探測方法，它通過記錄場源關斷后的電磁響應探測地下介質電阻率分布[1]。該方法具有探測成本低、勘探效果好等特點[2]，已經在金屬資源勘查[3]、采空區探測與定位[4-5]等方面獲得了廣泛應用。硫化物、金屬礦化物等極化介質的電導率通常會表現出頻散特性，與頻率呈復函數關系，瞬變電磁信號會因激發極化（Induced polarization， IP）效應發生畸變。在同一收發裝置下，IP 響應與 TEM 響應疊加在一起，會出現符號反轉現象，這一畸變現象會嚴重影響數據處理和解釋[6]。

電阻率頻散介質的瞬變電磁三維正演通常采用兩種策略： ① 頻譜法，即首先根據觀測時間范圍確定頻點范圍，然后計算異常體在頻域的電磁響應，再利用時頻轉換算法將其轉換到時域[7-9]。頻譜法能夠直接將復電阻率函數代入正演計算，無需進行其他處理，但正演精度依賴于頻點的選擇。近年來，有學者將離散化時域 Maxwell 方程組直接投影到 Krylov 子空間，選擇合適的頻點[10]，一定程度上提高了計算效率和精度[11]。 ② 時步法，首先采用隱式或顯式差分格式對時域 Maxwell 方程組進行時間域離散，求取初始場后，再根據時間步長逐步遞歸計算每個時間點的時域電磁場 [12-14]。由于電阻率頻散介質僅在頻域中有歐姆定律的解析式，相較于頻譜法，時步法在三維正演中的計算效率更高，但仍需要求解復電阻率在時域的表達式。

關于頻散介質歐姆定律的時域表達，一種可行的方法是直接對時域場進行卷積運算[15]。通常對復電阻率進行時頻變換，利用傅里葉逆變換或拉普拉斯逆變換將歐姆定律改寫成時域場關于時間的導數，再采用時間域差分法將其離散為時域場的卷積形式。這種方法適用于不存在分數次冪的復電阻率模型（如Debye 模型）。對于 Cole-Cole 模型，則會產生分數階導數。現有的離散方法只能離散時間域整數階導數，不能離散時間域分數階導數 [16]。為此，有學者相繼 提 出 了 Debye 近 似 法 [17]和 Padé 級 數 逼 近 法 [18]。前者對復電阻率進行一階泰勒展開，采用不同階次的Debye 模型近似 Cole-Cole 模型，這通常需要求解非線性最優化問題，不僅費時，所得結果也容易陷入局部最優值[14]。后者通常應用于隱式差分迭代格式，在某個中心頻率對 Cole-Cole 模型進行 Padé 級數逼近，再進行傅里葉逆變換，將歐姆定律改寫為時域場的卷積形式 [12]。在 Padé 級數逼近法中，中心頻率的選取對正演模擬誤差影響較大。基于此，Cai 等 [13]提出一種自適應選取中心頻率算法，能在不同的迭代時刻選取合理的中心頻率，提高了正演精度。

本文利用 Padé 級數逼近方法開發了一種基于時間域隱式迭代的瞬變電磁三維有限體積正演算法。關于 Padé 級數中心頻率選取的問題，本文在前人研究[13]的基礎上進一步提出一種方法，以量化Padé 逼近過程對正演所產生的誤差，并基于誤差信息分布實現中心頻率自適應選取，選取的中心頻率更合理，適用于更廣泛的激電參數范圍，允許在不同區域選取不同的中心頻率，使得該方法更加適用于數據的反演解釋。為分析激電效應產生的影響，本文分別對極化介質模型和非極化介質模型進行瞬變電磁法建模，分析了兩者間的特征差異。最后，以蒙古昌圖錫力錳銀鉛鋅多金屬礦床為正演模型，模擬分析了低充電率實例模型異常體的充放電過程。

1 含激電效應的瞬變電磁正演理論

本文采用時步法進行瞬變電磁場模擬。從Maxwell 方程組出發，對于給定的空間區域 和時間范圍[0，T]開展輸變電磁場模擬。時間域電磁場

的擴散滿足 Maxwell 微分方程[19]

式中： E 和 B 分別為電場強度和磁感應強度； μ 為磁導率； ε 為介電常數； J_r 為外部激發源的電流密度； J 表示傳導電流密度。

1. 1 歐姆定律的時域場卷積形式

考慮傳導電流密度項的情況下，歐姆定律可表示為

J=σ（ω）E

對于極化介質，上式中的電導率 σ（ω）=1/ρ ，是一個與頻率有關的函數，其中 ρ 為電阻率， ω=2πf 表示角頻率， ，f 表示頻率。 σ（ω） 可 以 用 Cole-Cole 模 型表 示[20]

式中： σ₀ 為直流電導率； η 為充電率，取值范圍為 0～ 1;τ 為時間常數； c 為頻率相關系數，取值范圍為 0～ 1。將式（4）代入式（3）并消去分母，得到

σ₀E+σ₀τ^c（iω）^cE=J+（1-η）τ^c（iω）^cJ

上式中，當 clt;1 時，直接對式（5）進行時頻轉換會產生分數階微分方程。由于求解這類分數階微分方程的數值解過于復雜，通常不會對其直接求解[21]。因此，本文利用 Padé 級數逼近式（5）中與頻率相關的 項[12，22]

式中： 表示 Padé 級數的相關系數； M，N 分別為 P，Q 的階數，其值取得越高，逼近效果越好，但所需求解的初始場數量也越多，相關系數的求解過程詳見 2. 1 節。令 M=N=K ，即上、下 Padé 級數取同階，將式（6）代入式（5）并消去分母，可以得到頻域歐姆定律逼近式

式中

式 （7）是 利 用 Padé 級 數 逼 近 后 得 到 的 Cole-

Cole 模型，它可以表示一定頻段內的頻率域復電阻率。但要想得到關于時域的歐姆定律，就要對其進行傅里葉逆變換

式中的導數項 ?^kf/?t^k 可以用高階后退歐拉法進行離 散[13]

式中：δt 表示時間步長 。將式（10）代入式（9），得到時域歐姆定律的離散形式

其中

1. 2 有限體積法正演

完成 Cole-Cole 模型的時域表達式求解后，首先對式（1）和式（2）在時間域進行離散[23]

式中 j 表示時間域離散的第 j 個時刻，例如 B^（j） 表示第 j 個時刻的磁場分布。再采用擬態有限體積法進行空間域離散。考慮式（11）、式（13）和式（14）的離散弱解形式[24-25]，采用空間交錯網格 [26]對其進行空間域離散。Cole-Cole 模型相關參數均定義在單元格中心，磁場 B 定義在單元格面心，電場 E 和電流密度 J 定義在單元棱邊。利用微分算子的標準交錯離散化，可以得到離散公式

式中：下標 e、f、r 分別代表與棱邊、面和源相關的參數； α 為時間步長 δt 的倒數； c 為旋度算符的離散形式； 和 分別表示有限體積離散后的棱邊平均值、介電常數、磁導率的倒數、式（12）

中 和 的平均矩陣； S_r^（j+1） 為 j+1 時刻離散后的源項。上式的詳細離散過程和表達式可參考文獻[27]。

將式（17）和式（15）代入式（16），可消去電流密度項和磁場項，得到

電場和電流密度可以通過以下迭代式得到

通過求解式（18）即可得到 j+1 時刻的電場分布，然后采用式（19）和式（20）計算得到磁場分布和電流密度分布。隨時間遞歸求解式（18）～式（20），便可實現含激發效應的瞬變電磁三維正演模擬。本文采用多波前分解法的 MUMPS 直接求解器[28]實現對式（18）大型稀疏系數矩陣的直接求解。

相較于常規瞬變電磁正演方法，含激電效應的瞬變電磁正演模擬需要考慮 Padé 級數逼近歐姆定律過程，而 Padé 級數的逼近效果直接影響正演精度。Padé 級數中心頻率的選擇是直接影響逼近效果的重要參數，一般認為瞬變電磁早期信號對應高頻信息而晚期信號對應低頻信息，所以 Padé 級數中心頻率應當隨瞬變電磁正演時間的不同而做出不同的選擇。本文提出一種新的自適應選取 Padé 級數中心頻率的方法。

2 Padé 級數中心頻率自適應選擇算法

本節先通過 Padé 系數的求取過程展示頻域下中心頻率與 Cole-Cole 模型逼近精度關系，再進一步量化這種逼近誤差對瞬變電磁正演響應的影響，最后基于誤差信息實現自適應中心頻率選取算法。

2. 1 中心頻率與逼近效果的關系

Padé 級數相關系數主要是通過泰勒展開式來求取，令 ，那么在 s=1 處的 Padé 級數和泰勒級數為

式中： ω_o 表示 Padé 級數展開的中心頻率； c_l 為 s=1 處泰勒展開式的相關系數。

值得注意的是，式（21）與式（6）的區別僅在于式（21）多了系數項 ω₀^c ，其中 ω₀ 稱為特征頻率，與中心頻率 f₀ 的關系為 f₀=ω₀/2π 。 參 照 Rekanos 等

[18]的推導過程，可以在 s=1 處（即 ω=ω_?0） 得到關于式（6）的 Padé 系數計算公式

式中 q_z 和 ρ_z 可以用泰勒展開的系數項求取，它們之間的關系為

顯然，正演結果依賴于式（6）的 Padé 逼近精度。與泰勒展開類似，當選取的特征頻率 ω_?0 與 ω 相差較大時，逼近精度往往不好，要想得到好的逼近效果就要選取合適的中心頻率 f₀ 。

下面通過一個例子來分析頻域的中心頻率選取對 Cole-Cole 模型逼近精度的影響。Cole-Cole 模型的參數設置為： σ₀=0.01 S?m^-1，η=0.2，τ= 。

圖 1 是中心頻率分別選取 10Hz 和 1000Hz 情況下的 Padé 逼近，其中取 Padé 級數的階數 K=5 。明顯看出，中心頻率不同，有效近似頻段也不同。當 K= 5 時，有效近似頻率約為中心頻率的兩個數量級，而且Padé 級數的階數越高，有效近似頻率范圍越廣。

雖然 Padé 級數越高階意味著逼近效果越好，但是不宜通過選取太高階的級數來減小逼近誤差以提高正演精度，其原因是：一是使用高階 Padé 級數會增加計算量，無論是 Padé 系數的求取還是時域下歐姆定律的卷積形式，增大階數都會增加計算量；二是高階 Padé 級數會提高初始條件的要求，這主要體現在式（10），所采用的 Padé 級數的階數越高，所需計算的初始場也越多。

圖1 中心頻率為 10，1000Hz 時電導率Padé 逼近效果

2. 2 逼近誤差的量化方法

在不改變中心頻率的時間范圍內，正演所用的電阻率模型可以認為是由單個中心頻率進行 Padé級數逼近后的電阻率模型，這顯然與真實的 Cole-Cole 電阻率模型存在差異。因此，想要在時域實現正演，存在兩個問題：逼近模型正演后的誤差能否忽略？怎樣在時域（而非頻域）選取恰當的中心頻率，使誤差最小？

對此，Marchant 等[12]認為觀測時間段內存在主頻段，只要選取的中心頻率足以覆蓋主頻段，就可以忽略逼近誤差對正演的影響，因而可以在整個觀測時間段內使用某個固定的中心頻率。但這一方法對Padé 級數的逼近階數要求較高，正演精度也比自適應選取中心頻率的算法低。Cai 等[13]提出可以針對不同的迭代時刻選取不同的中心頻率（圖 2）。本文借鑒這一思路，提出一種方法，能夠量化中心頻率對逼近電阻率模型正演帶來的誤差，為 Padé 級數自適應選取算法提供評價依據。

理論上，Padé 級數逼近產生的誤差與維度無關，而一維模型正演計算量遠遠小于三維模型，所以本文采用一維模型研究 Padé 級數過程。對于一維模型，采用某個中心頻率展開進行 Padé 級數逼近計算的正演響應，與直接采用 Cole-Cole 電阻率模型計算得到的正演響應，它們之間的相對誤差（RE）在一定程度上可以表征逼近模型的正演誤差。RE 的表達式為

圖2 中心頻率、時間步長隨時間變化的關系圖

式中 F₁（t）、F₂（t） 分別表示采用該中心頻率展開得到的正演響應和采用 Cole-Cole 電阻率模型得到的正演響應。根據上式可以得到該中心頻點的相對誤差與時間的函數關系。本文采用的一維正演算子是基于均勻層狀大地水平電偶極源推導的電場表達 式[29]

式中： G=Ids 為電偶源強度； J₁ 為一階貝塞爾函數；R=（x²+y²）^1/2 ； r_TE 為 TE 模式反射系數，它是一個與地下介質的電導率和地層厚度有關的函數； u₀ 表示波數域的波數函數； z 表示垂直坐標。將復電阻率代入式（25）可以得到不同頻點下水平電偶極源的正演電磁響應；再沿發射回線進行積分即可得到回線源的正演響應[30]；最后對其做余弦變換[31]，可得到矩形線圈的正演電磁響應。因 Cole-Cole 復電阻率模型有解析形式，后文將 Cole-Cole 模型的正演響應作為一維解析解。

通過式（24）可以確定某一中心頻率適用的時間范圍，即在該時間段內，采用該中心頻率所得到的逼近誤差是在可接受范圍內的。通過在對數域內均勻離散選擇擬定的中心頻率范圍，再利用式（24）求取每個元素的相對誤差函數，得到圖 3 所示階數 K= 3、5 下的相對誤差與采集時間、中心頻率之間的二維分布關系。

圖 3 不同階數下時域中心頻率相對誤差

圖 3 中出現一個低誤差區域（淺色區域），說明這個區間的正演誤差在可接受的范圍內，所對應的中心頻率是恰當的。利用該方法能夠量化逼近誤差對正演結果的影響，快速有效地評估三維模擬中所選擇的中心頻率是否合理。

Padé 級數的階數對于高誤差區域（深色區域）的影響較小，對低誤差區域的影響較大。對比圖 3a和圖 3b，可見在降低階數后，晚期的低誤差區間幾乎消失了。這意味著，如果要在晚期保持可靠的正演精度，就需要選擇更高的 Padé 級數階數。對于大部分復電阻率模型而言，選取 5 階 Padé 級數就能獲得良好的逼近效果，且計算量和初始條件也在可接受范圍內（后續模型計算中將展示誤差對比）。因此，本文中模型正演所使用的 Padé 級數均為 5。

2. 3 中心頻率選取策略

三維瞬變電磁正演選取中心頻率是在時間域進行的。為獲取可靠的中心頻率參數，本文基于相對誤差與采集時間、中心頻率之間的二維分布關系，采用最小誤差原則，提出了一種自適應中心頻率選取策略。對于圖 2 所示擬選中心頻率范圍內，遍歷計算等時間步長時段內的平均相對誤差，選出平均相對誤差最小的中心頻率，開展三維正演計算。該過程的詳細描述見圖 4 所示流程。

對于三維異常體模型，各個離散單元的頻率相關系數 c 存在一定差異。傳統的自適應算法選取的中心頻率對頻率相關系數的變化較為敏感。在頻率相關系數變化大的區域，通常使用自適應算法選取多組中心頻率，對不同網格單元采用不同的中心頻率組，但這可能導致后續的反演工作計算成本很高，本文提出的自適應算法則能有效避免這一問題。

圖4 自適應中心頻率選取策略流程

為驗證本文自適應算法選擇的中心頻率可適用于頻率相關系數變化大的區域，設置參數為： ，利用本文算法，得到了不同時間的中心頻率（圖 5a）。然后，將這些中心頻率值分別應用于 c=0.4，0.8 的模型，得到的磁場響應模擬結果見圖 5b、圖 5c。對比圖 5b 與圖 5c 可見，在一定 Ψ_c 值范圍內，中心頻率可以簡化為一個常數，不會對三維異常體的正演結果產生顯著影響。

圖5 三維異常體中心頻率選取算法驗證結果

3 理論模型正演分析

3. 1 精度驗證和特征分析

采用均勻半空間模型檢驗算法的準確性。對于該均勻半空間，均采用 Cole-Cole 模型， 設 σ₀= 0.01S?m^-1 ， τ=0.1s ，分別對 η=0.2，c=0.6 及η=0.1，c=0.4 兩種情況進行模擬分析。三維計算區域尺寸為 10km×10km×6km ，網格剖分總數為 49×49×47 。發射采用邊長為 50m 的矩形線圈，發射電流波形為矩形波，電流強度為 1 A，在回線中心進行信號接收。時間步長策略、擬態有限體積法相關參數等與文獻[27]保持一致，中心頻率選取采用前文所述的自適應選取算法。模擬計算采用一臺配置了兩個主頻為 2.4GHz 的 24 核 心Intel Xeon CPU 的服務器，內存為 512GB 。

圖 6 為該模型的三維數值解與一維解析解的對比。圖中實線變為虛線的位置稱為變號處，按照 Sei?del 等[32]的定義，變號處所對應的時間稱為變號時間。由圖可見，按照 Padé 逼近算法得到的三維正演響應與一維解析解吻合很好，除了變號附近誤差稍大，其余位置的相對誤差均小于 5% ，證明了本文算法的有效性和準確性。

圖 7 為均勻半空間模型取 η=0.2，c=0.6 時Padé 逼近階數分別選取 3 和 5 情況下的相對誤差曲線，對應的模擬時間分別為 720.6s，850.6s， 。由圖可見， K=5 時的計算精度更穩定，尤其在信號變號附近及 20ms 之后的晚期時間，穩定性更明顯，這

圖6 均勻半空間模型的電磁正演結果對比（上）及其與解析解的相對誤差（下）虛線和實線分別對應響應為正和負的情況。

圖 7 均勻半空間模型 η=0.2，c=0.6 時 Padé 級數不同階數相對誤差曲線一結論與圖 3 是一致的。

3. 2 極化/非極化異常體響應特征分析

設計圖 8 所示三維極化模型驗證本文算法對三維極化模型的有效性，并對極化體與非極化體模型（圖 8 中）正演響應進行對比分析。背景模型為非極化均勻半空間，異常體為一個低阻極化三維體。Cole-Cole 相 關 參 數 為 σ₀=0.1s/m ， η=0.5 ， τ= 0.2s，c=0.5 。發射源為 50m×50m 的中心回線，置于地面異常體正上方，二者的水平中心位置一致。

另外，為對比極化體模型與非極化體模型的電磁響應特征，對該模型充電率 η=0 的情況進行模擬，對地面脈沖響應、多測道和電磁響應曲線進行分析。

圖 8 三維極化體模型

圖 9 為極化與非極化介質異常體模型的正演響應信號對比。由圖 9a 和 9b 可見，兩種介質的響應早期 （0.05ms） ）在空間分布上均能清晰標示出目標體的邊界位置；隨時間增長（圖 9c、圖 9e），非極化介質模型的響應異常變得不明顯，尤其在 1ms 時刻，根據地面響應幾乎無法觀察到明顯的異常。極化介質的情況則不同，無論是在變號前還是變號后，響應信號的差異均較明顯，且瞬變電磁信號的衰減明顯慢于非極化介質。

圖 9 三維非極化模型（上）和極化模型（下）不同時刻電磁響應黑色虛線為異常體邊界。

為深入分析變號時間附近的信號空間分布特征，繪制 y=0 剖面不同時刻 （0.016～63.100ms） ）的磁場響應曲線（圖 10a）；分別選取異常體中心點（0， 0）、異常體邊界附近（ 80m ， 0）和 200m ， 0）三個位置的響應曲線（圖 10b 和 10c）進行對比。可見在接近變號時刻，異常體上方區域出現信號快速下降并變號（負響應）的現象，且信號衰減速度明顯低于圍巖上方區域。這種現象是由于極化電流產生的信號衰減速度遠低于傳導電流所致。在早期階段，由于傳導電流占主導地位，極化電流對總響應的貢獻較小；在變號時刻，兩者信號強度相近，導致總響應接近于零；變號后，極化電流的強度相對較高，因而地面觀測到的是負響應。隨著觀測時間從 2.51ms 延 長 到 63.1ms ，地 面 觀 測 到 負 響 應 的 范 圍 從[-55m ， 55m ]擴大至 [-135m ， 135m] ，同時，最大電 磁 響 應 幅 值 從 2.17×10^-10V?A^-1?m^-2 衰 減 至1.51×10^-12V?A^-1?m^-2 。

3. 3 低充電率復雜形狀異常體模型激電特性

為分析復雜地質結構對激電效應的影響，并評估算法在實際地質勘查中的適用性，基于蒙古昌圖錫力錳銀鉛鋅多金屬礦床 [33]的探測結果構建了一個三維模型。鉆探數據和激發極化反演結果顯示，該 礦 床 主 體 長 約 350m ，寬 約 170m ，平 均 厚 約50m ，傾角較小，圍巖因侵蝕作用呈現較弱的充電特性。據此，設計圖 11 所示階梯狀異常體模型，其厚度沿 x 方向由 60m 逐漸減小為 40m ，且每層階梯寬度為 50m 。發射裝置采用 210m×210m 的大定源回線，其中心位于坐標系原點。

圖10 三維極化/非極化模型電磁響應虛線對應負響應。

模型磁場正演結果見圖 12。由于異常體埋深和厚度均較小，地面接收信號的差異主要集中在0.1ms 到 0.4ms 。早期（圖 12a、圖 12b）的響應只能大致反映異常體的相對長寬特性，但由于厚度差異不明顯，難以有效區分。結合圖 13 所示的部分測點衰減曲線可見， 0.4ms 后異常體所在區域的信號逐漸擴散，差異逐漸減小。

為分析激電效應的影響，繪制 z=60m 所在平面的四個時刻 （0. 05，0.10，0.15，0.20ms） ）的電流密度分布（圖 14），可見四個時刻的電流密度分布均能夠清晰地反映異常體的形狀特征。在較早階段（0.1，0.2ms） ），電流方向與發射電流基本保持一致，顯示出明顯的渦旋特征，電流主要集中于充電率較高的異常體邊界，而中心區域渦旋效應較弱。隨時

間推移 （0.4ms） ），異常體區域內的電流強度逐漸弱于背景區域呈現出放電特征。至 1.00ms 時，放電過程產生的電流超過了傳導電流，形成逆向渦旋電流，其強度遠高于圍巖中的渦旋電流，但對整體響應的影響仍然較小。圍巖由于充電率較低，放電特性不明顯，未能在有效觀測時間觀察到放電過程特有的逆渦旋電流。

進一步分析充放電狀態。設置與圖 14 一致的模型參數，但目標體和圍巖的充電率均設置為零。分別對極化模型和非極化模型模擬計算 y=0 平面不同時刻的電流密度分布進行比較，結果見圖 15。與對照組（非極化模型）相比，充電狀態下異常體的電流密度明顯更高。觀察電流方向，可以發現在充電狀態下，極化模型與非極化模型的電流向下傳導特征明顯不同：極化模型的極化電流在早期出現水平方向的渦旋流動（圖 15），到 0.4ms 時刻，異常體頂部電流密度顯著下降至低于圍巖水平，預示著異常體即將進入放電階段；非極化模型則未出現這一特征。圖 15d 揭示了異常體的放電過程，放電產生的電流方向與傳導電流相反，圍巖周圍的電流均匯集于異常體區域。

圖 14 實際模型 z=60m 平面不同時刻的電流密度分布圖中白色方框為異常體輪廓，箭頭所指方向為電流方向。

4 結論

本文基于電場分量的時域有限體積法，采用隱式時間步長迭代方式，實現了含激電效應的瞬變電磁三維正演模擬。對于 Padé 級數法中中心頻率的選取，定量分析了一維空間的逼近誤差，并基于誤差信息實現了中心頻率的自適應選取，所選取的中心頻率可適用于頻率相關系數變化較大的區域。通過一維、三維模型的數值實驗驗證了本文算法的精度和有效性。對比分析了極化體模型和非極化體模型在地面的電磁響應特征，發現極化體模型具有更明顯的響應異常。最后，基于實際生產數據進行建模和模擬，結果表明地下電流主要集中于異常體邊界，且充電率高的模型表現出更為明顯的充、放電現象。

圖15 極化模型（上）與非極化模型（下）不同時刻 y=0 平面的電流密度分布

參 考 文 獻

［ 1］ LI H， XUE G Q， HE Y M. Decoupling induced po? larization effect from time domain electromagnetic data in a Bayesian framework[J]. Geophysics， 2019， 84（6）： A59?A63.

［2］ 李昊錦， 毛玉蓉， 王新宇， 等. 地—井瞬變電磁三維 各向異性響應特征研究[J]. 石油地球物理勘探， 2024， 59（5）： 1184?1196. LI Haojin， MAO Yurong， WANG Xinyu， et al. Sur? face?borehole transient electromagnetic responses in 3D anisotropic media [J]. Oil Geophysical Prospecting， 2024， 59（5）： 1184?1196.

［3］ BORGATTI H， MELO A T， KNAUER L G. Arti? sanal and small ? scale gold mining in the Southern Es? pinhaco Range， Brazil， using time?domain electromag? netic induction： prospecting， efficiency， and environ? mental aspects[J]. Brazilian Journal of Geology， 2023， 53（3）： 349?358.

［4］ 張瑩瑩 . 電性源瞬變電磁法綜述[J]. 物探與化探， 2021， 45（4）： 809?823. ZHANG Yingying. Review on the study of grounded source transient electromagnetic method[J]. Geo? physical and Geochemical Exploration， 2021， 45（4）： 809?823.

［5］ 馬振軍， 底青云， 薛國強，等 . 地—空瞬變電磁法電 阻率成像研究與應用[J]. 地球物理學報， 2021， 64（3）： 1090?1105. MA Zhenjun， DI Qingyun， XUE Guoqiang， et al. The research and application of resistivity imaging of semi?air? borne transient electromagnetic method[J]. Chinese Jour nal of Geophysics， 2021， 64（3）：1090?1105.

［6］ 殷長春， 繆佳佳， 劉云鶴，等 . 時間域航空電磁法激 電效應對電磁擴散的影響[J]. 地球物理學報， 2016， 59（12）： 4710?4719. YIN Changchun， MIAO Jiajia， LIU Yunhe， et al. The effect of induced polarization on time ? domain air? borne EM diffusion[J]. Chinese Journal of Geophysics， 2016， 59（12）： 4710?4719.

［ 7］ HOHMANN G W， NEWMAN G A. Transient elec? tromagnetic responses of surficial， polarizable patches [J]. Geophysics， 1990， 55（8）： 1098?1100.

［ 8］ WAIT J R， DEBRORUX P. Induced polarization in electromagnetic inductive schemes[J]. Geophysical Prospecting， 1984， 32（6）： 1147?1154.

［9］ 王少杰， 周磊， 謝興兵， 等 . 基于差分進化算法的瞬 變電磁一維反演[J]. 石油地球物理勘探， 2024， 59（2）： 343?351. WANG Shaojie， ZHOU Lei， XIE Xingbing， et al. Transient electromagnetic one ? dimensional inversion based on differential evolution algorithm[J]. Oil Geo? physical Prospecting， 2024， 59（2）： 343?351.

［10］ ZASLAVSKY M， DRUSKIN V， KNIZHNERMAN L. Solution of 3D time?domain electromagnetic problems using optimal subspace projection[J]. Geophysics， 2011， 76（6）： F339?F351.

［11］ 周建美， 景旭， 魯凱亮， 等. 基于重啟多項式Krylov 子空間模型降階的全波形瞬變電磁三維正演[J]. 地 球物理學報， 2023， 66（6）： 2646?2657. ZHOU Jianmei， JING Xu， LU Kailiang， et al. Three?dimensional full?time transient electromagnetic modeling using restarting polynomial Krylov sub? space model order reduction method[J]. Chinese Jour? nal of Geophysics， 2023， 66（6）： 2646?2657.

［12］ MARCHANT D， HBBER E， OLDENBURG D W. Three?dimensional modeling of IP effects in time?domain electromagnetic data[J]. Geophysics， 2014， 79（6）： E303? E314.

［13］ CAI H Z， HU X Y， XIONG B， et al. Finite ? element time?domain modeling of electromagnetic data in general dispersive medium using adaptive Padé series[J]. Com ? puters and Geosciences， 2017， 109： 194?205.

［14］ JI Y J， WU Y Q， GUAN S S， et al. 3D numerical modeling of induced?polarization electromagnetic re? sponse based on the finite?difference time?domain method [J]. Geophysics， 2018， 83（6）： E385?E398.

［15］ SMITH R S， WALKER P W， POLZER B D， et al. The time?domain electromagnetic response of polarizable bodies?an approximate convolution algorithm[J]. Geo? physical Prospecting， 1988， 36（7）： 772?785.

［16］ 王 洪 華 ， 王 敏 玲 ， 張 智 ，等 . 基 于 Padé 逼 近 的 Cole? Cole 頻散介質 GPR 有限元正演[J]. 地球物理學報， 2018， 61（10）： 4136?4147. WANG Honghua， WANG Minling， ZHANG Zhi， et al. Simulation of GPR in Cole?Cole dispersive me? dia by finite element method based on the Padé ap? proximations[J]. Chinese Journal of Geophysics， 2018， 61（10）： 4136?4147.

［17］ COMMER M， PETROV P V， NEWMAN G A. FDTD modelling of induced polarization phenomena in transient electromagnetics[J]. Geophysical Journal International， 2017， 209（1）： 387?405.

［18］ REKANOS I T， PAPADOPOULOS T G. An auxilia? ry differential equation method for FDTD modeling of wave propagation in Cole?Cole dispersive media[J]. IEEE Transactions on Antennas and Propagation， 2010， 58（11）： 3666?3674.

［19］ HBER E， ASCHER U M， OLDENBURG D W. In? version of 3D electromagnetic data in frequency and time domain using an inexact all ? at ? once approach[J]. Geophysics， 2004， 69（5）： 1216?1228.

［20］ PELTON W H. Mineral discrimination and removal of inductive coupling with multi?frequency IP?reply[J]. Geophysics， 1984， 49（9）： 1557?1558.

［21］ GE J C， EVERETT M E， WEISS C J. Fractional diffusion analysis of the electromagnetic field in frac? tured media， Part 2： 3D approach[J]. Geophysics， 2015， 80（3）： E175?E185.

［22］ WEEDON W H， RAPPAPORT C M. A general method for FDTD modeling of wave propagation in ar? bitrary frequency ? dispersive media[J]. IEEE Transac? tions on Propagation， 1997， 45（3）： 401?410.

［23］ OLDENBURG D W， HABER E， SHEKHTMAN R. Three dimensional inversion of multisource time domain electromagnetic data[J]. Geophysics， 2013， 78（1）： E47?E57.

［24］ HYMAN J M， SHAHKOV M. Mimetic discretiza? tions for Maxwell′s equations[J]. Journal of Computa? tional Physics， 1999， 151（2）： 881?909.

［25］ 彭榮華， 胡祥云， 韓波， 等. 基于擬態有限體積法的 頻率域可控源三維正演計算[J]. 地球物理學報， 2016， 59（10）： 3927?3939. PENG Ronghua， HU Xiangyun， HAN Bo， et al. 3D frequency?domain CSEM forward modeling based on the mimetic finite ? volume method[J]. Chinese Journal of Geophysics， 2016， 59（10）： 3927?3939.

［26］ YEE K. NNMERICAL solution of initial boundary value problems involving Maxwell’s equations in iso? tropic media[J]. IEEE Transactions on Antennas and Propagation， 1966， 14（3）： 302?307.

［27］ 劉亞軍 . 時間域電磁法三維各向異性正演及分析 [D]. 武漢： 中國地質大學（武漢）， 2020. LIU Yajun. Three?Dimensional Forward Modeling and Analysis of Time?Domain Electromagnetic Data for An Anisotropic Medium[D]. Wuhan： China University of Geosciences （Wuhan）， 2020.

［28］ AMESTOY P R， DUFF I S， L’Excellent J Y， et al. A fully asynchronous multifrontal solver using dis? tributed dynamic scheduling[J]. SIAM Journal on Ma? trix Analysis and Applications， 2001， 23（1）： 15?41.

［29］ NABIGHIAN M N. Electromagnetic Methods in Ap? plied Geophysics[M]. Society of Exploration Geophy? sics， 1988.

［30］ 李建慧， 劉樹才， 朱自強，等 . 矩形回線激發的電磁 場與磁場的對稱關系研究[J]. 中南大學學報（自然科 學版）， 2010， 41（2）： 638?642. LI Jianhui， LIU Shucai， ZHU Ziqiang， et al. Rela? tionship between electromagnetic field and magnetic field′s symmetric excited by rectangular loop[J]. Journal of Central South University （Science and Technology）， 2010， 41（2）： 638?642.

［31］ NEWMAN G A， HOHMANN G W， ANDERSON W L. Transient electromagnetic response of a 3?dimen? sional body in a layered earth[J]. Geophysics， 1986， 51（8）： 1608?1627.

［32］ SEIDEL M， TEZKAN B. 1D Cole?Cole inversion of TEM transients influenced by induced polarization[J]. Journal of Applied Geophysics， 2017， 138： 220?232.

［33］ 馬一行， 顏廷杰， 何鵬，等. 激發極化法在覆蓋區礦 產勘查中的應用——以內蒙古昌圖錫力錳銀鉛鋅多 金 屬 礦 勘 查 為 例 [J]. 物 探 與 化 探 ， 2019， 43（4）： 709?717. MA Yihang， YAN Tingjie， HE Peng， et al. Explora? tion in coverage area： A case study of the Changtuili Mn?Ag?Pb?Zn polymetalic ore deposit lnner Mongolial [J]. Geophysical and Geochemical Exploration， 2019， 43（4）： 709?717.

（本文編輯：劉海櫻）

作 者 簡 介

楊威 副研究員，1984 年生；2007 年獲成都理工大學信息工程專業學士學位，2011 年獲成都理工大學地球探測與信息技術專業碩士學位，2024 年獲成都理工大學地球探測與信息技術專業博士學位；現就職于中石化西北油田分公司勘探開發研究院，主要從事地震資料處理及地震和電磁方法應用研究。