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

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

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

0 引言

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

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

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

1 理論方法

1.1 邊值問題

激發(fā)極化法中，地下點電流源產(chǎn)生的電位滿足微分方程［9］

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

其中σb滿足

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

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

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

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

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

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

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

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

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

因此

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

式（12）的通解為

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

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

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

1.2 有限單元法

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

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

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

1.3 關鍵環(huán)節(jié)分析

1.3.1 背景場求解思路

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

1.3.2 波數(shù)選取方法

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

1.3.3 迭代算子和波數(shù)域電場的求解

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

1.3.4 三維數(shù)值模擬算法流程

本次三維數(shù)值模擬采用以下流程：

（1）對求解域采用均勻網(wǎng)格剖分，并利用采樣定理得到波數(shù)域采樣點。

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

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

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

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

2 算例分析

2.1 算法正確性驗證

為了驗證本文算法的正確性，利用均勻半空間背景下的低阻、高極化率球體模型進行測試。模型如圖1所示。采用均勻網(wǎng)格剖分計算區(qū)域，節(jié)點總數(shù)為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 模型的電位解析解和數(shù)值解（左）及其相對誤差（右）

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

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

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

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

模型如圖3 所示，兩個異常體中心點在地面上的投影位于直角坐標系原點，鉆孔井口G坐標為（10 m，0，0），其他參數(shù)詳見圖中標注。分別利用五種裝置（二極裝置、三極裝置、偶極裝置、積分裝置、差分裝置）沿井孔測量。采用均勻網(wǎng)格剖分計算區(qū)域，節(jié)點總數(shù)為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 所示。可以看出，供電電極放置于不同位置時，五種裝置得到的極化體視電阻率和視極化率曲線形態(tài)大體相似，即礦體附近均出現(xiàn)異常。這種現(xiàn)象可依據(jù)“電偶極子”疊加原理進行分析，即礦體產(chǎn)生的異常場可近似等效為多個“電偶極子”組合疊加的結果［38］。此外，低阻異常體產(chǎn)生的異常遠大于高阻異常體產(chǎn)生的異常，說明低阻異常體對電流更敏感。另外，從圖中還看出，對于同一供電位置，差分裝置測得的異常幅值遠大于其他裝置，積分裝置和偶極裝置產(chǎn)生的異常幅值次之，三極裝置再次之，二極裝置測得的異常幅值最小，說明通過物理差分可有效提高分辨率。因三極裝置觀測的是電位差，所以對介質(zhì)的分辨率高于二極裝置；積分裝置和偶極裝置相當于在三極裝置的基礎上對供電源進行了一次差分，所以它們的分辨率也高于三極裝置。盡管積分裝置和偶極裝置的分辨率相同，但是積分發(fā)射裝置是由兩組偶極發(fā)射裝置同向供電組合而成，所以供電電流更大，信噪比更高；差分裝置則是在偶極裝置的基礎上對供電源再次進行差分計算，所以分辨率進一步提高。綜上所述，差分裝置勘探效果優(yōu)于其他裝置，可為野外勘探測量裝置的選擇提供參考。2.2.2 不同背景介質(zhì)下的高極化礦體模型

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

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

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

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

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

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

3 結論

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

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

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

附錄A 有限元單元分析

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

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

式中

其中l(wèi)為單元長度。

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

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

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

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

式中

式中

式中

式中

式中

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