胡代明 劉 昊 吳小平

1）中國長沙 410019 中國電建集團中南勘測設(shè)計研究院有限公司

2）中國長沙 410019 水能資源利用關(guān)鍵技術(shù)湖南省重點實驗室

3）中國合肥 230026 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院

0 引言

激發(fā)極化法是電法勘探的經(jīng)典方法之一，可分為時間域激電法（TDIP）和頻率域激電法（FDIP），在金屬礦產(chǎn)、非金屬礦產(chǎn)、油氣等資源、能源勘探中獲得廣泛而有效的應(yīng)用（傅良魁，1983；李金銘，2005；李帝銓等，2007；鄭冰，2015；向葵等，2016；趙榮春等，2016；孫仁斌等，2017；趙后越，2017；景強，2018）。近年來，激發(fā)極化法還不斷拓展到水文地質(zhì)、環(huán)境地質(zhì)監(jiān)測與應(yīng)用等領(lǐng)域（Revil et al，2010；Okay et al，2013；Schwartz et al，2014；Boaga，2017）。在金屬礦勘探中，激發(fā)極化法因其獨特的激發(fā)極化特性，比直流電阻率方法更有效，成為金屬礦勘探中不可或缺的方法。由某些礦物引起的激電異常通過電阻率難以分辨，而激發(fā)極化法與電阻率的結(jié)構(gòu)特點無關(guān)，利用該方法則可觀測到較強的充電率特征。在實際礦產(chǎn)資源勘探中起伏地形難以避免，與其他電磁勘探方法相比，利用激發(fā)極化方法觀測到的異常理論上基本不受地形影響，這是其眾所周知的另一大優(yōu)勢，當(dāng)然這只是在地下電阻率、極化率各向同性情況下導(dǎo)出的結(jié)論（傅良魁，1983；李金銘，2005）。

激發(fā)極化法與電磁法有著密切的聯(lián)系，其數(shù)值模擬主要基于“等效電阻率法”的時間域激發(fā)極化進行（Pelton et al，1978）?；诜€(wěn)定電流場方程，根據(jù)極化率和電阻率已知的公式可求得等效電阻率，利用等效電阻率去替換初始電阻率，可求得極化總場，一次場可由解析公式求得，然后求得時間域極化響應(yīng)。頻率域激發(fā)極化的現(xiàn)有模型如Cole-Cole模型和Dias模型（Dias，1968）。Cole-Cole模型和Dias模型描述了復(fù)電阻率，該復(fù)電阻率取決于頻率，通過引入其他參數(shù)，如時間常數(shù)、頻率相關(guān)系數(shù)以及充電率，可以得出有關(guān)極化材料詳細(xì)特征的結(jié)論，可由含電容器的等效電路表示。

在現(xiàn)有理論中，只基于各向同性介質(zhì)對復(fù)雜地形影響下的激發(fā)極化法進行了研究。然而，地下介質(zhì)的各向異性是客觀存在的，一些常見巖石具有明顯的各向異性。Linde等（2004）在石灰?guī)r地區(qū)測量的各向異性系數(shù)高達3.7和4.5。Schmutz等（2000）發(fā)現(xiàn)，直流電法和瞬變電磁法在資料反演解釋中必須引入較大的各向異性系數(shù)。Herwanger等（2004）指出，如果采用各向同性井間成像，將產(chǎn)生無法解釋的層狀條帶圖案。Li等（2013）對任意各向異性介質(zhì)中的可控源電磁法響應(yīng)進行了自適應(yīng)二維有限元分析，結(jié)果表明，如果忽略各向異性，資料解釋會導(dǎo)致儲層位置的錯誤估計。因此，對于激發(fā)極化法而言，各向異性的研究具有重要意義。

目前所知，僅Liu等（2017）基于各向異性介質(zhì)進行了激發(fā)極化法的數(shù)值模擬，結(jié)果表明，在不考慮電磁效應(yīng)的情況下，實際解釋中必須考慮各向異性。然而，其只計算了簡單各向異性介質(zhì)的響應(yīng)，未考慮地形的影響，難以滿足激發(fā)極化法勘探的需要。目前，各向異性和地形對激發(fā)極化資料解釋的綜合影響尚未見報道。本研究旨在發(fā)展一種適用于帶地形、非結(jié)構(gòu)有限元三維各向異性激發(fā)極化的數(shù)值模擬方法。

本文主要介紹了激發(fā)極化的基本方程、邊值問題、TDIP和FDIP響應(yīng)，討論了各向異性和地形條件下的激發(fā)極化響應(yīng)，并進行規(guī)律總結(jié)。

1 激發(fā)極化模型正演

1.1 基本方程和邊值問題

在地球?qū)щ娊橘|(zhì)中，電勢的基本方程和邊值問題需在地空表面滿足Neumann邊界條件，在無窮遠(yuǎn)邊界滿足混合邊界條件。公式如下

式中，σ表示電導(dǎo)率張量，δ為狄拉克三角函數(shù)，?為拉普拉斯算子，r是觀測位置，I為電流強度，Γ0是地空邊界，ΓR是無窮遠(yuǎn)邊界，r0為電流源位置，n表示法向分量。

1.2 空間離散化

利用Galerkin方法對基本方程和邊值問題進行簡化（Ciarlet，2002），公式如下

式中，Ω為目標(biāo)區(qū)域。通過模型離散化，電勢可以用插值函數(shù)u=N近似表示，公式如下

式中，Ue為單元節(jié)點電勢向量，N為形函數(shù)，Ne為單元總數(shù)目U分別為體單元矩陣和面單元矩陣。將所有局部矩陣組合成一個全局系統(tǒng)方程，公式如下

1.3 時間域激發(fā)極化響應(yīng)

極化率是時間域激發(fā)極化法的核心參數(shù)，極化率的定義（Soueid Ahmed et al，2018）為

式中，η∈[0,1]；σ0為直流電法中介質(zhì)的原始電導(dǎo)率，與初始電阻率ρ0成反比，即ρ0=1/σ0。因此，極化率可表示為

由于文中研究的是三維各向異性介質(zhì)，故等效電阻率ρ∞(r)、初始電阻率ρ0(r)和極化率η(r)均為三階張量，則

時域激發(fā)極化模型需要計算2次，利用等效電阻率法可實現(xiàn)時域激發(fā)極化數(shù)值模擬（Seigel，1959）。視極化率ηs(r)可表示為二次電勢Us(r)與總電勢U(r)之比，公式如下

因此，視極化率可以通過一次電勢和總電勢來求解。

1.4 頻率域極化響應(yīng)

在FDIP中，當(dāng)介質(zhì)為極化時，可以用Cole-Cole模型在多頻率下的復(fù)電阻率ρ*(ω)來描述介質(zhì)的各個部分（Pelton et al，1978），公式如下

式中，ρj0為第j種介質(zhì)的電阻率；mj為充電率；τj表示體積極化效應(yīng)的時間常數(shù)，通常在10-2—102s之間變化，是描述巖石結(jié)構(gòu)的敏感指標(biāo)；cj為頻率相關(guān)系數(shù)，數(shù)值區(qū)間為0.5—1.0。

當(dāng)介質(zhì)表現(xiàn)為各向異性時，Cole-Cole模型參數(shù)也可以表示為張量，因此可用于計算各向異性介質(zhì)在不同頻率下的復(fù)電阻率響應(yīng)ρ*(ω)，公式如下

FDIP中的頻率效應(yīng)Fe被定義（Yang et al，2008）為

式中，ρL為低頻下確定的視電阻率，ρH為高頻下確定的視電阻率。

2 算法驗證

為了驗證算法的正確性，使用一個兩層各向異性模型，參數(shù)見圖1所示。該兩層水平層狀各向異性（VTI）模型具有一個等效的各向同性層狀模型，其具有解析解（Pek?en et al，2014）。在圖1所示等效的各向同性層狀模型中，第一層厚度λh，各向同性介質(zhì)平均電阻率為100 Ω·m；第二層為均勻半空間，各向同性介質(zhì)平均電阻率為10 Ω·m。使用對稱四極裝置進行測量，測線最大距離為180 m，電位電極在地表上間隔2 m分布。將兩層各向異性模型的數(shù)值模擬結(jié)果與兩層各向同性模型的解析解結(jié)果進行對比，最大相對誤差小于0.5%（圖2），從而驗證了算法的準(zhǔn)確性。

圖1 一個兩層各向異性模型Fig.1 A two-layer anisotropic model

圖2 視電阻率的數(shù)值解和解析解以及相對誤差對比結(jié)果Fig.2 Comparison of apparent resistivity and relative error between numerical solutions and analytical solutions

3 模型算例

3.1 時間域激發(fā)極化模型

圖3所示為在平坦地形模型下方有一個深度為150 m的正方體異常，研究各向同性（λ=1）、水平層狀各向異性VTI（λ=2,4）和垂直層狀各向異性TTI（λ=2,4）條件下，該異常體的時間域激發(fā)極化響應(yīng)結(jié)果。

圖3 各向異性的立方體異常模型示意Fig.3 Schematic diagram of anisotropic cube anomaly model

設(shè)模型尺寸為5000 m×5000 m×5000 m，立方體異常位于x—y平面原點下方，埋深150 m，大小為100 m×100 m×100 m。該模型的背景圍巖（圖3中A）電阻率為100 Ω·m，各向同性極化率為η1=η2=η3=0.01 mV/V。表示各向異性系數(shù)，ρL表示橫向或平行層理方向電阻率，ρT表示縱向或垂直層理方向電阻率，異常體參數(shù)見表1。

表1 異常體參數(shù)Table 1 Parameters of anomalies

在地表采用對稱四極裝置測視極化率，電位電極間距5 m，電流源電極最大距離為1000 m，最小距離為10 m，電流強度為1 A。以各向同性異常為例，設(shè)初始電阻率為10 Ω·m，極化率為0.2 mV/V，計算得到異常體等效電阻率為12.5 Ω·m，即可用等效電阻率來計算三維有限元建模中的總電勢U(r)。

圖4(a)—(c)顯示了由對稱四極裝置獲得的視極化率剖面圖，在平坦地形中，當(dāng)異常的電阻率和極化率分別表現(xiàn)為各向同性和各向異性時，在時間域激發(fā)極化響應(yīng)中的視極化率剖面表現(xiàn)不同。當(dāng)異常為各向同性時（λ=1），觀測的異常體視極化率約1%，與背景圍巖相近。而當(dāng)VTI異常為各向異性時（λ=2），觀測的最大視極化率約為1.6%；當(dāng)VTI異常表現(xiàn)出較強的各向異性時（λ=4），觀測的最大視極化率大于2.2%。與各向同性結(jié)果[圖4(a)]相比，各向異性結(jié)果[圖4(c)]具有更清晰的極值中心，異常位置定位更準(zhǔn)確。由圖4(d)—(e)可見，將η2和η3的值進行交換后，對于任意各向異性TTI介質(zhì)，模擬結(jié)果均無法定位異常位置。由圖4(f)可見，在各向異性VTI（λ=4）條件下，視極化率對異常體的反應(yīng)最強；TTI介質(zhì)的視極化率在背景極化率值（1%）附近波動，說明在實際激發(fā)極化勘探中，在地表觀測不到實際存在于地下的TTI各向異性異常；同時，VTI介質(zhì)的各向異性越強，異常定位越容易，而TTI介質(zhì)在平坦地形下也很難被檢測到。

圖4 視極化率斷面(a)各向同性異常（λ=1）；(b)各向異性VTI（λ=2）；(c)各向異性VTI（λ=4）；(d)各向異性TTI（λ=2）；(e)各向異性TTI（λ=4）；(f)異常體正上方的測線視極化率；虛線表示異常的真實位置Fig.4 Cross-section of apparent polarizability

3.2 頻率域激發(fā)極化模型

為了探討地形對頻率域激發(fā)極化模型的影響，設(shè)計如圖5的示例模型。圖5中平坦地形、臺地地形、地塹地形的球體異常模型被用來研究FDIP響應(yīng)結(jié)果。模型尺寸為1000 m×1000 m×500 m。在地表采用中間梯度裝置觀測，電位電極分布范圍為200 m×480 m，測點間隔10 m，電流電極間距800 m。球體異常半徑10 m，位于x—y平面中心，異常頂部距離地表為30 m。背景圍巖參數(shù)如下：電阻率ρx=ρy=ρz=100 Ω·m，充電率為mx=my=mz=0.05，時間常數(shù)為τx=τy=τz=1 s，頻率相關(guān)系數(shù)為c=0.25；球體異常參數(shù)如下：電阻率ρx=ρy=ρz=10 Ω·m，充電率為mx=my=mz=0.05，時間常數(shù)為τx=τy=τz=1 s，頻率相關(guān)系數(shù)為c=0.25，2個計算頻率分別為2 Hz和2/13 Hz。平坦地形、臺地地形和地塹地形的頻率效應(yīng)等值線見圖6。

圖5 頻率域激發(fā)極化觀測模型及其有限元網(wǎng)格(a)測量線布置；(b)平坦地形非結(jié)構(gòu)化網(wǎng)格；(c)地臺地形非結(jié)構(gòu)化網(wǎng)格；(d)地塹非結(jié)構(gòu)化網(wǎng)格Fig.5 Frequency domain induced polarization observation model and its finite element mesh

圖6 不同地形下地表觀測的頻率效應(yīng)等值線(a)平坦地形；(b)地臺地形；(c)地塹地形；(d)異常正上方測線的頻率效應(yīng)Fig.6 Frequency effect of surface observation under different terrains

由圖6可見，可在地面觀測到球體異常的頻率效應(yīng)響應(yīng)，以及異常體正上方地表中心位置出現(xiàn)橢球形的頻率效應(yīng)Fe最大區(qū)域，該區(qū)域與圖中黑色虛線表示的異常的真實位置較為一致，根據(jù)該頻率效應(yīng)值可以判斷異常體為高極化異常；臺地地形的最大頻率效應(yīng)最弱（紅線），地塹地形的最大頻率效應(yīng)最強，平坦地形最大頻率效應(yīng)則介于二者之間[圖6(d)]。因此，臺地地形對真實異常的響應(yīng)存在減弱效果，地塹地形對真實異常的響應(yīng)存在增強效果，但3種地形對異常真實位置的反映均較為準(zhǔn)確。

4 結(jié)論

在金屬礦產(chǎn)勘查過程中，地下極化介質(zhì)分布具有任意性和各向異性特征。采用非結(jié)構(gòu)有限元法建立了考慮各向異性的三維激電模型。在此基礎(chǔ)上，首次研究了各向異性對激電介質(zhì)的綜合影響。在TDIP結(jié)果中，真實介質(zhì)中各向異性的變化對觀測到的地表激電數(shù)據(jù)有較大影響。VTI異常在平坦地形下容易被準(zhǔn)確定位，各向異性越大，信號越強，定位越準(zhǔn)確。而TTI介質(zhì)即使在平坦地形下也很難被探測到。在實際勘探中，在激發(fā)極化數(shù)據(jù)反演工作中應(yīng)考慮介質(zhì)的各向異性特征。在FDIP結(jié)果中，臺地地形對異常結(jié)果響應(yīng)存在減弱作用，地塹地形對異常響應(yīng)則存在增強效果。