各向異性激發(fā)極化法三維數(shù)值模擬

胡代明 劉 昊 吳小平

1）中國長(zhǎng)沙 410019 中國電建集團(tuán)中南勘測(cè)設(shè)計(jì)研究院有限公司

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

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

0 引言

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

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

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

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

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

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

1.1 基本方程和邊值問題

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

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

1.2 空間離散化

利用Galerkin方法對(duì)基本方程和邊值問題進(jìn)行簡(jiǎn)化（Ciarlet，2002），公式如下

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 算法驗(yàn)證

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

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

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

3 模型算例

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

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

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

設(shè)模型尺寸為5000 m×5000 m×5000 m，立方體異常位于x—y平面原點(diǎn)下方，埋深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

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

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

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

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

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

圖5 頻率域激發(fā)極化觀測(cè)模型及其有限元網(wǎng)格(a)測(cè)量線布置；(b)平坦地形非結(jié)構(gòu)化網(wǎng)格；(c)地臺(tái)地形非結(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 不同地形下地表觀測(cè)的頻率效應(yīng)等值線(a)平坦地形；(b)地臺(tái)地形；(c)地塹地形；(d)異常正上方測(cè)線的頻率效應(yīng)Fig.6 Frequency effect of surface observation under different terrains

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

4 結(jié)論

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