各向異性激發極化法三維數值模擬

胡代明 劉 昊 吳小平

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

2）中國長沙 410019 水能資源利用關鍵技術湖南省重點實驗室

3）中國合肥 230026 中國科學技術大學地球和空間科學學院

0 引言

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

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

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

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

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

1 激發極化模型正演

1.1 基本方程和邊值問題

在地球導電介質中，電勢的基本方程和邊值問題需在地空表面滿足Neumann邊界條件，在無窮遠邊界滿足混合邊界條件。公式如下

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

1.2 空間離散化

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

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

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

1.3 時間域激發極化響應

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

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

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

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

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

1.4 頻率域極化響應

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

式中，ρj0為第j種介質的電阻率；mj為充電率；τj表示體積極化效應的時間常數，通常在10-2—102s之間變化，是描述巖石結構的敏感指標；cj為頻率相關系數，數值區間為0.5—1.0。

當介質表現為各向異性時，Cole-Cole模型參數也可以表示為張量，因此可用于計算各向異性介質在不同頻率下的復電阻率響應ρ*(ω)，公式如下

FDIP中的頻率效應Fe被定義（Yang et al，2008）為

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

2 算法驗證

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

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

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

3 模型算例

3.1 時間域激發極化模型

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

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

設模型尺寸為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。表示各向異性系數，ρL表示橫向或平行層理方向電阻率，ρT表示縱向或垂直層理方向電阻率，異常體參數見表1。

表1 異常體參數Table 1 Parameters of anomalies

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

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

圖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 頻率域激發極化模型

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

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

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

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

4 結論

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