基于三維數值模擬的激電法差分裝置響應特征分析

凌嘉宣，陳志芳，鄧維，戴世坤，周印明，陳輕蕊

（1. 桂林航天工業學院計算機科學與工程學院，廣西桂林 541004；2. 浙江省錢塘江管理局勘測設計院，浙江杭州 310016；3. 中南大學地球科學與信息物理學院，湖南長沙 410083；4. 湖南科技大學計算機科學與工程學院，湖南湘潭 411201；5. 國防科技大學氣象海洋學院，湖南長沙 410073）

0 引言

激發極化法（簡稱激電法）作為電法勘探重要的分支方法，通常情況下比電阻率法和大地電磁法更適用于目標體與圍巖電阻率差異不大的浸染型金屬礦床的勘探［1］。根據數據采集方式，激電法測量方式包括地—地、地—井、井—地、井—井［2］，其中后三種方法是利用已有鉆孔探測井旁盲礦或礦體空間分布，確定礦體埋深、方位，追蹤并圈定礦化帶。這些方法勘探成本較高。地—地測量是在地面發射和接收，無需鉆孔，成本較低，但是不適用于深部礦產的勘測［3］。不同測量方法各有優缺點，因此開展激發極化法探測應根據礦區（工區）實際情況和地質任務選擇合適的測量方法。

正演是反演的基礎，為了能更好地進行精細反演，需開展激電法三維數值模擬方法研究。激電法數值模擬與電阻率法數值模擬具有相通性，即激電法正演需要進行兩次直流電阻率法的正演計算［4］。目前三維數值模擬方法主要包括積分方程法［5-8］、有限單元法［9-15］及有限差分法［15-18］，其他數值模擬方法均由這三種方法衍生得到。針對不同測量方法的激電三維數值模擬方法主要包括：Hohmann［19］和樸華榮等［20］基于積分方程法研究了地—地測量的三維極化體激發極化效應；張輝等［21］利用體積分方程法對同時存在激電效應和電磁效應的半空間三維體進行研究，并分析了激電效應和電磁效應并存對正演結果的影響；Tsourlos 等［22］和黃俊革［23］利用有限單元法開展了地—地三維激電數值模擬，驗證了采用偏導數矩陣與視極化率的關系式計算極化率的可行性；呂玉增［24］和王智［25］利用基于四面體網格交叉剖分的有限單元法開展了地—井、井—地三維激電數值模擬，系統分析了井旁不同形態礦體的激電異常特征；周峰等［26］利用有限差分法開展了井—井激電三維正演，并利用二次異常電位差曲線反演得到礦體的激發極化特征；李長偉［27］利用放射狀三棱柱單元結構網格的有限單元法實現了井—井激電法三維數值模擬；高文龍等［28］利用變步長有限差分法實現了井—井激電法三維數值模擬，并分析了侵入帶、井徑、不同電極系裝置及不同井液電阻率對激電效應結果的影響；李靜和等［29］為了高效、高精度探測重金屬污染的激發極化空間分布，提出了接觸式供電—地面觀測和接觸式供電—直接觀測兩種觀測系統，進行了接觸式激發極化法滲漏型目標探測的應用研究，解決了復雜施工環境條件的限制問題，同時提高了觀測信號的強度和觀測精度。

以上傳統的激電法三維數值模擬大多基于空間域，在計算復雜模型時存在計算量大、存儲需求高、效率較低等問題。本文基于二維傅里葉變換，將激電法三維數值模擬問題轉化為波數域一維數值模擬問題進行求解，可以避免大型稀疏線性方程組的求解，降低計算量和存儲需求，明顯提高計算效率。 此外，傳統激電測深裝置主要為二極裝置、三極裝置、偶極裝置等，這些裝置對地下介質的分辨率有限。本文通過改善數據采集裝置提高對地下電性異常體的分辨率，即引入差分裝置［30］和積分裝置。差分裝置指的是以兩組供電電流強度相等、方向相反的偶極發射裝置并列或共線作為供電裝置，利用偶極電極接收信號；積分裝置以兩組供電電流大小相等、方向相同的偶極發射裝置并列或共線作為供電裝置，也是利用偶極電極接收信號。本文基于井—井測量方式，對比傳統測深裝置、差分裝置及積分裝置對井旁極化礦體的分辨率。最后，對不同背景介質下二極裝置和差分裝置得到的視電阻率和視極化率響應特征進行分析，比較各自的優勢。

1 理論方法

1.1 邊值問題

激發極化法中，地下點電流源產生的電位滿足微分方程［9］

式中：σ為介質電導率；U為總電位；I為電流強度；δ為狄拉克函數；ω為點電流源A到計算區域的立體角，當點電流源在地面時，ω=2π，當點電流源在地下時，ω=4π。總電位U可分解為背景電位Ub和異常電位Ua，即U=Ub+Ua；總電導率σ可分解為背景電導率σb和異常電導率σa，即σ=σb+σa。因此，式（1）可改寫為

其中σb滿足

將式（3）代入式（2），可得

因電位U與電場E滿足E=-?U，散射電流密度J=σbE，故式（4）改寫為

對于水平層狀介質模型，將式（5）沿三個方向展開，可得

式中J=［Jx，Jy，Jz］，Jx，Jy，Jz為空間域三個方向上的散射電流。

對式（6）沿水平方向進行二維傅里葉變換，可得

從式（7）可看出，方程右端Jx、Jy、Jz為未知，該式為的非線性方程，無法直接求出。因此，本文提出以下求解思路：①首先令總電場等于背景場Eb（n表示迭代次數，此時n=1），求出散射電流J；②對散射電流J進行二維正傅里葉變換，得到波數域散射電流的三分量并代入式（7），求得波數域異常電位；③利用波數域異常電位與異常電場的關系式求出波數域異常電場三分量，并對和進行二維反傅里葉變換，得到空間域異常電位Ua和異常電場分量Eax、Eay、Eaz，再與背景場相加可得總電位和總電場；④計算殘差，若ε小于期望值，則輸出結果，否則，根據迭代算子修改得到新的總電場，并利用總電場計算散射電流，重復步驟②～步驟④直至滿足迭代結束條件。

根據前述求解激電法三維數值模擬問題的思路，需給出方程組的邊界條件、求解方程組的方法、背景場、波數域電場與電位的關系式、波數選取方法、迭代算子等。

因式（7）沿x和y方向為波數域，由傅里葉變換特性可知，x和y方向的邊界條件自動滿足，所以只需給出z方向的邊界條件。假設z軸向下為正方向，地面為上邊界zmin，此處電流密度法向為零，滿足

選取異常體下方某位置為下邊界zmax，假設供電電極位于該區域之外，總電位U和背景電位Ub均滿足拉普拉斯方程

因此

對應的波數域異常電位滿足微分方程

式（12）的通解為

式中C和D為常數。因無窮遠處（z→+∞）異常電位Ua等于零，有

對式（14）求導，下邊界zmax上的異常電位Ua滿足

式（7）、式（8）和式（15）即組成三維激電法的邊值問題。

1.2 有限單元法

對于三維激電法的邊值問題，可采用基于二次插值的一維有限單元法求解。利用伽遼金法［8］得到邊值問題等價有限單元方程

式中：M為垂向剖分單元總數；Ni為第e個單元的第i個節點的二次插值形函數，這里i=1，2，3。式中各項的表達式見附錄A。式（16）可表示為矩陣形式的線性方程組

從求解三維激電異常電位邊值問題的思路可看出，本文利用二維傅里葉變換對復雜數值模擬問題進行降維求解，即將三維數值模擬問題轉換成為不同波數域中的一維數值模擬問題進行求解，可大幅度降低計算量及對存儲的需求，提高計算效率。

1.3 關鍵環節分析

1.3.1 背景場求解思路

本文算法基于異常場，因此需先確定背景場。在給出邊值問題時，分別假設背景介質為均勻半空間、層狀模型、電阻率隨深度連續變化等不同情況。電阻率隨深度連續變化的特殊情況就是層狀介質，即垂向節點電阻率值連續變化。當背景為均勻半空間時，源在任意位置所產生的背景場的求解方法詳見文獻［1］。當背景為層狀介質時，本文推導了點電流源在任意位置產生的背景場，詳見附錄B。

1.3.2 波數選取方法

本文采用基于標準FFT法的二維傅里葉變換，波數的選取基于采樣定理，具體方法詳見文獻[32-33]。需要說明的是，二維反傅里葉變換可以采用基于標準FFT法或高斯-FFT法［33-34］。

1.3.3 迭代算子和波數域電場的求解

針對求解電位的邊值問題，首先根據背景場求出電位和電場的近似解，然后利用算子進行迭代求解，最終獲得高精度數值解。本文采用的迭代算子［35-36］為

1.3.4 三維數值模擬算法流程

本次三維數值模擬采用以下流程：

（1）對求解域采用均勻網格剖分，并利用采樣定理得到波數域采樣點。

（2）計算背景電場，并求出相應散射電流，利用二維傅里葉變換得到波數域散射電流。

（3）令邊值問題式（7）等式右端的波數散射電流等于前一個步驟得到的值，求解方程組得到波數域異常電位，進而得到波數域異常電場。

（4）利用二維反傅里葉變換得到空間域異常電位和異常電場，加上背景場得到總場。

（5）將更新后的總電場與先前總電場進行對比，判斷是否滿足殘差要求。若滿足，則輸出總電位和電場；否則，利用迭代算子修改總電場，并根據修改的總電場計算波數域散射電流，重復步驟（3）～步驟（5），直至滿足殘差要求。

2 算例分析

2.1 算法正確性驗證

為了驗證本文算法的正確性，利用均勻半空間背景下的低阻、高極化率球體模型進行測試。模型如圖1所示。采用均勻網格剖分計算區域，節點總數為101×101×81，水平方向剖分間隔為1 m，垂直方向剖分間隔為0.5 m。算法利用Gauss-FFT法進行二維傅里葉反變換。采用二極裝置進行測量，電流強度為1 A，供電電極置于點A（-10 m，0，0）。測點沿地面x方向布設。

圖1 低阻、高極化球體模型示意圖

采用解析解［1］作為理論值，兩種算法的總電位和相對誤差如圖2 所示。可以看出，兩種算法得到的總電位曲線吻合較好，最大相對誤差為0.54%，驗證了本文算法的正確性。

圖2 圖1 模型的電位解析解和數值解（左）及其相對誤差（右）

2.2 差分裝置激電異常響應特征

為了研究均勻半空間背景下激電差分裝置對高極化率異常體的異常響應特征及不同背景條件下差分裝置的異常響應特征，分別討論以下兩種情況。

2.2.1 高阻/低阻、高極化礦體模型

對于井—井勘探，為了探測井旁盲礦，為了節約成本，通常在單孔井中進行測量，即供電電極和測量電極均放置在同一口井中［2］。在金屬礦產勘探中，黃鐵礦、黃銅礦、方鉛礦等礦產多表現為低阻、高極化特性，鉻鐵礦、黑鎢礦及侵染狀硫化礦等通常表現為高阻、高極化特性［37］。為了研究高阻或低阻的高極化礦體在差分裝置下的激發極化響應特征，不失一般性，可利用簡化的金屬礦模型——兩個大小相同的棱柱異常體模型進行測試。

模型如圖3 所示，兩個異常體中心點在地面上的投影位于直角坐標系原點，鉆孔井口G坐標為（10 m，0，0），其他參數詳見圖中標注。分別利用五種裝置（二極裝置、三極裝置、偶極裝置、積分裝置、差分裝置）沿井孔測量。采用均勻網格剖分計算區域，節點總數為61×51×81，水平方向剖分間隔為2 m，垂直方向剖分間隔為1 m。另外，為了分析電極位于不同位置時盲礦的激電異常特征，分別將供電電極對A-B置于井口（中間點坐標：（10.0 m，0，1.5 m））、井中（中間點坐標（10 m，0，50 m））和井底（中間點坐標（10 m，0，90 m）），使用極距為1 m 的接收電極對M-N沿鉆孔自上而下布設。

圖3 兩個棱柱異常體模型示意圖

供電電極位于井口、井中和井底時得到的視電阻率和視極化率曲線如圖4 所示。可以看出，供電電極放置于不同位置時，五種裝置得到的極化體視電阻率和視極化率曲線形態大體相似，即礦體附近均出現異常。這種現象可依據“電偶極子”疊加原理進行分析，即礦體產生的異常場可近似等效為多個“電偶極子”組合疊加的結果［38］。此外，低阻異常體產生的異常遠大于高阻異常體產生的異常，說明低阻異常體對電流更敏感。另外，從圖中還看出，對于同一供電位置，差分裝置測得的異常幅值遠大于其他裝置，積分裝置和偶極裝置產生的異常幅值次之，三極裝置再次之，二極裝置測得的異常幅值最小，說明通過物理差分可有效提高分辨率。因三極裝置觀測的是電位差，所以對介質的分辨率高于二極裝置；積分裝置和偶極裝置相當于在三極裝置的基礎上對供電源進行了一次差分，所以它們的分辨率也高于三極裝置。盡管積分裝置和偶極裝置的分辨率相同，但是積分發射裝置是由兩組偶極發射裝置同向供電組合而成，所以供電電流更大，信噪比更高；差分裝置則是在偶極裝置的基礎上對供電源再次進行差分計算，所以分辨率進一步提高。綜上所述，差分裝置勘探效果優于其他裝置，可為野外勘探測量裝置的選擇提供參考。2.2.2 不同背景介質下的高極化礦體模型

圖4 供電電極AB 位于不同位置時的電阻率（上）和極化率（下）曲線

為了進一步研究不同背景下礦體的激發極化特征，建立圖5所示低阻、高極化異常體模型。異常體中心位置在地面的投影位于坐標原點。為了便于對比，對此模型分別利用二極裝置和差分裝置進行測量。點源（兩組偶極源中心位置）位于（-20 m，0，0），接收電極M、N極距為1 m，沿鉆孔向下布設。假設背景介質的極化率為0，電阻率分為三種情況：（a）均勻半空間，電阻率為1000 Ω ?m；（b）三層水平介質，自上而下各層介質電阻率為100、200、100 Ω·m，第一、第二層的厚度分別為10、20 m；（c）連續介質，電阻率隨深度變化滿足關系式計算區域大小為120 m×100 m×50 m，剖分網格總數均為61×51×51，水平方向剖分間隔均為2 m，垂直方向剖分間隔為1 m。三種不同背景下得到的電阻率和極化率見圖6。

圖5 低阻、高極化異常體模型示意圖

圖6 不同背景介質下的電阻率（上）和極化率（下）曲線

對比不同背景介質下兩種裝置的視電阻率曲線，可見差分裝置能更好地反映背景介質和異常體的垂向電性特征和激發極化特征。

需要說明的是，不同背景介質下同一種裝置得到的視極化率曲線形態相似，說明無激電效應的背景介質對極化率影響較小。在實際生產中，對復雜背景介質下極化礦體的勘查可綜合考慮電阻率和極化率兩個參數，以便更好地分析礦體及圍巖的空間分布。

3 結論

（1）將空間域與波數域相結合能有效求解激電法三維數值模擬問題，為探索激電法數據的高效、高精度的數值模擬提供一條新思路。

（2）與其他裝置相比，差分裝置對異常體的視電阻率/視極化率異常響應更明顯，分辨率更高。但因該裝置是兩組發射電流大小相同、方向相反的偶極發射裝置組合而成，相對于其他裝置，向地下供電能量更低，需增大供電電流，所以在實際應用中，需要綜合考慮地質任務及經濟成本，選擇合適的測量裝置及參數。

（3）視電阻率能有效反映目標體及圍巖的電性特征，而視極化率能反映目標體的激電效應特征，綜合分析這兩個參數可更精確地解釋異常體的空間分布及地電異常特征。

附錄A 有限元單元分析

z軸方向可利用節點二次插值函數表示單元e的電位、電流密度及電導率，其形函數表達式為

式中j、p、m表示單元內節點編號。因此，式(16)的積分項為

式中

其中l為單元長度。

式(16)中的最后一項表示z=zmax處節點的值，此時

附錄B 任意位置點源層狀介質電位表達式推導

令笛卡爾坐標系原點位于地面，點電流源置于層狀介質中任意位置，則源所在地層內測點和非源所在地層內測點的電位和電場分別為

式（B-1）和式（B-2）中：j表示源所在地層序號；i表示非源層序號；系數A和B的表達式根據點電流源在不同位置分為五種情況。假設地下包含n層介質，則點電流源位于第1 層、第2 層、第3～第n-2 層、第n-1層、第n層五種情況下A和B的取值分別見式（B-3）～式（B-7）。

式中

式中

式中

式中

式中

上面公式適用于地下介質至層數n>5 的情況。當n<5時，可將某層分解為電導率相等的多層，即可采用以上公式。若需計算多源下任意位置的電位和電場，可分別計算每個點電流源產生的場值，進行矢量疊加即可。求解點源在層狀介質任意位置產生的電位和電場可采用濾波法。傳統的濾波法需要逐個測點遞推計算層系數，對于層狀和節點較多的情況，計算速度慢。因此，本文采用優化的濾波法求解，即先計算電位和電場表達式中與z軸深度方向的相關系數，然后計算所有波數下與z軸深度方向無關的層系數和某一平面節點的電場／電位系數并進行存儲，后續計算任意深度上測點的場值時，只需調用存儲數據，將其與z軸相關系數相乘并進行累加即可，如此可快速計算背景場值。模型測試結果表明，采用此方法比優化前計算速度提高約20倍。