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

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

0 引言

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

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

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

1 理論方法

1.1 邊值問題

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

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

其中σb滿足

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

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

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

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

對(duì)式（6）沿水平方向進(jìn)行二維傅里葉變換，可得

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

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

因式（7）沿x和y方向?yàn)椴〝?shù)域，由傅里葉變換特性可知，x和y方向的邊界條件自動(dòng)滿足，所以只需給出z方向的邊界條件。假設(shè)z軸向下為正方向，地面為上邊界zmin，此處電流密度法向?yàn)榱悖瑵M足

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

因此

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

式（12）的通解為

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

對(duì)式（14）求導(dǎo)，下邊界zmax上的異常電位Ua滿足

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

1.2 有限單元法

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

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

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

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

1.3.1 背景場(chǎng)求解思路

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

1.3.2 波數(shù)選取方法

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

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

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

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

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

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

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

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

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

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

2 算例分析

2.1 算法正確性驗(yàn)證

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

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

采用解析解［1］作為理論值，兩種算法的總電位和相對(duì)誤差如圖2 所示?？梢钥闯觯瑑煞N算法得到的總電位曲線吻合較好，最大相對(duì)誤差為0.54%，驗(yàn)證了本文算法的正確性。

圖2 圖1 模型的電位解析解和數(shù)值解（左）及其相對(duì)誤差（右）

2.2 差分裝置激電異常響應(yīng)特征

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

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

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

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

圖3 兩個(gè)棱柱異常體模型示意圖

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

附錄A 有限元單元分析

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

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

式中

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

式(16)中的最后一項(xiàng)表示z=zmax處節(jié)點(diǎn)的值，此時(shí)

附錄B 任意位置點(diǎn)源層狀介質(zhì)電位表達(dá)式推導(dǎo)

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

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

式中

式中

式中

式中

式中

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