趙 寧 黃明衛(wèi) 申亞行 陶德強(qiáng) 秦 策*

(①河南理工大學(xué)物理與電子信息學(xué)院，河南焦作 454150； ②河南省瓦斯地質(zhì)與瓦斯治理重點(diǎn)實驗室——省部共建國家重點(diǎn)實驗室培育基地，河南焦作 454150； ③東方地球物理公司綜合物化探處，河北涿州 072751)

0 引言

煤層氣是世界公認(rèn)的三大非常規(guī)天然氣之一，儲層的壓裂改造是提高采收率的主要手段[1]。其中，微地震監(jiān)測技術(shù)可以通過儲層壓裂改造引發(fā)的震源點(diǎn)進(jìn)行成像，從而為儲層壓裂改造提供依據(jù)。然而，煤層氣水力壓裂區(qū)圍巖或煤層破碎嚴(yán)重、煤層較軟，致使微地震信號微弱，資料解釋困難[2]。由于電磁信號對導(dǎo)電流體和裂縫中的壓裂劑較敏感，Beskardes等[3]采用三維直流電法對煤層氣壓裂效果進(jìn)行了深入研究，其中高精度三維電法數(shù)值模擬是關(guān)鍵。

目前，三維直流電正、反演計算中應(yīng)用最廣的數(shù)值模擬方法主要有積分方程法[4-6]、有限差分法[7-8]和有限元法[9-13]。這三種方法各有優(yōu)缺點(diǎn)：積分方程法計算速度快，但只能應(yīng)用結(jié)構(gòu)化網(wǎng)格處理異常體，對復(fù)雜地質(zhì)體的模擬具有局限性；有限差分法精度不高，且只能應(yīng)用結(jié)構(gòu)化網(wǎng)格；有限元法精度較高，且可以結(jié)合非結(jié)構(gòu)化網(wǎng)格[14-17]處理復(fù)雜模型，但相較于前兩種方法，所需的計算內(nèi)存較高。近些年來，隨著計算機(jī)硬件計算能力的提高和數(shù)值計算方法的進(jìn)步，有限元法得到廣泛應(yīng)用[18]。

在三維直流電數(shù)值模擬中，有限元法通?？赏ㄟ^兩種策略提高解的精度，即加密網(wǎng)格和應(yīng)用高階形函數(shù)。為快速得到最優(yōu)的網(wǎng)格分布，目前正演通常采用h型自適應(yīng)加密算法[19-22]對網(wǎng)格進(jìn)行加密。該算法首先固定網(wǎng)格單元上形函數(shù)的階數(shù)(通常較低)，然后利用后驗誤差估計估算每個單元的誤差，最后對誤差較大的單元進(jìn)行加密。但是，若單元上解的光滑度較高，對單元進(jìn)行加密仍無法快速提高解的精度。為解決這一問題，Grayver等[23]在求解三維地電模型時，對網(wǎng)格中的所有單元應(yīng)用高階形函數(shù)，提高了解的精度，并證明高階自適應(yīng)有限元法可以用最少的自由度得到最精確的有限元解。

本文基于高階自適應(yīng)有限元實現(xiàn)了直流電阻率模型的正演。首先闡述了有限元算法的實現(xiàn)步驟和關(guān)鍵技術(shù)，包括高階形函數(shù)生成方法、后驗誤差估計技術(shù)和懸掛點(diǎn)上的約束條件；然后通過數(shù)值算例驗證了算法的精確性，并且對比了形函數(shù)階數(shù)不同時(p=1、2、3)的模擬結(jié)果；最后，將此方法應(yīng)用于沁水盆地南部某煤氣層壓裂監(jiān)測區(qū)的三維模型，通過三維直流電阻率模擬，驗證了方法的有效性。

1 正演理論

1.1 控制方程

假設(shè)在地面A點(diǎn)放置一個電流強(qiáng)度為I的點(diǎn)電源(圖1)，空間電位分布滿足

·(σU)=-IδA

(1)

式中：σ表示電導(dǎo)率；δA表示A點(diǎn)的Dirac delta函數(shù)；U表示電位。在地表Γs處，電流沿地表流動，所以電流密度的法向分量為0，表示為

(2)

式中n是邊界上外法向向量。

假定區(qū)域Ω內(nèi)，電性不均勻性對無窮遠(yuǎn)邊界Γ(圖1)處的電位分布不產(chǎn)生影響，即Γ處的電位為點(diǎn)電源所產(chǎn)生的電位，表示為

圖1 點(diǎn)源示意圖

(3)

式中：c是比例系數(shù)；r是A點(diǎn)到邊界Γ的距離。將式(3)對n求導(dǎo)，消去常數(shù)c，可得

(4)

式中r是由點(diǎn)A指向邊界Γ的方向向量。

綜上所述，三維直流電阻率模型邊界值問題可由以下偏微分方程描述

(5)

1.2 有限元離散

將區(qū)域Ω離散為一系列非結(jié)構(gòu)化六面體單元，每個單元e上電位的近似解為

(6)

式中：m是單元上自由度個數(shù)；ui是單元中第i個自由度上的電位；ue是由m個自由度上未知電位組成的向量；φi是第i個自由度上的形函數(shù)；φ是由m個形函數(shù)組成的向量。

對式(5)應(yīng)用伽遼金有限元法[24]，可得單元e上的弱解形式為

(7)

keue=fe

(8)

(9)

(10)

對所有單元分別組裝矩陣，可以得到系統(tǒng)線性方程組

KU=F

(11)

式中：K是一個n×n階的稀疏對稱正定矩陣，n是網(wǎng)格中自由度總數(shù)；U是網(wǎng)格中所有自由度上未知電位組成的n維向量；F是一個n維向量，描述源的分布。

1.3 形函數(shù)的生成

任意空間中、任意階數(shù)的形函數(shù)可由一維多項式的張量積表示[25]

(12)

j0(i)=i/(p+1)j1(i)=i×mod(p+1)j2(i)=i/(p+1)2

(13)

(14)

單元上形函數(shù)的個數(shù)m可由空間維度d和形函數(shù)的階數(shù)p確定

m=(p+1)d

(15)

以二維形函數(shù)(p=2)在單元上的分布為例(圖2)， 這些自由度上形函數(shù)的數(shù)學(xué)表達(dá)式為

圖2 二維形函數(shù)(2階)在單元上的分布示意圖

(16)

根據(jù)式(16)，分別繪制自由度0和自由度7上形函數(shù)在單元上的分布形態(tài)(圖3)。圖中黑色區(qū)域表示z=0平面，所以形函數(shù)為負(fù)值的區(qū)域在圖3中顯示不完整。可以看出，形函數(shù)所在自由度上函數(shù)值為1，在其他自由度上函數(shù)值均為0。

圖3 自由度0(a)和自由度7(b)上形函數(shù)在單元上的分布形態(tài)

1.4 后驗誤差估計方法

后驗誤差估計可以無需人工干預(yù)，以最快的速度得到最優(yōu)的網(wǎng)格分布。目前主要有兩類后驗誤差估計方法，即基于梯度恢復(fù)[15]和基于殘差[26]的后驗誤差估計方法。本文采用前者將單元中每個面上有限元解的梯度跳躍量[27]作為單元上的基本誤差

(17)

(18)

計算，其中l(wèi)是單元中面數(shù)目。

得到所有單元的誤差估計后，對誤差前10%的單元進(jìn)行加密得到新的網(wǎng)格，然后在新的網(wǎng)格上重新計算有限元解，直到迭代次數(shù)或自由度個數(shù)達(dá)到設(shè)定值，迭代結(jié)束。h型自適應(yīng)有限元算法流程如圖4所示。

圖4 h型自適應(yīng)有限元算法流程圖

1.5 懸掛點(diǎn)上的約束條件

本文采用八叉樹加密方法[28](圖5)對網(wǎng)格進(jìn)行自適應(yīng)加密，所以當(dāng)一個單元一個面被加密時，在相鄰兩個單元的公共面上會出現(xiàn)懸掛點(diǎn)。為保證有限元的連續(xù)性，需要在這些懸掛點(diǎn)上添加一些約束條件。為方便起見，以二維網(wǎng)格上應(yīng)用2階形函數(shù)為例說明懸掛點(diǎn)上所施加的約束條件(圖6)。三維空間及其不同階數(shù)形函數(shù)的情況與此類似。

圖5 八叉樹加密原理

圖6中，自由度0、1和2屬于單元M1，自由度3、4和5屬于單元M2。為保證單元M1與單元M2的公共面上有限元的連續(xù)性，需滿足

圖6 懸掛點(diǎn)示意圖自由度4是懸掛點(diǎn)

u0φ0+u1φ1+u2φ2=u3φ3+u4φ4+u5φ5

(19)

其中自由度2與自由度3、自由度1與自由度5表示同一自由度，即

(20)

當(dāng)通過1.3節(jié)中的步驟求得自由度0、1和2上的形函數(shù)表達(dá)式后，可以得到

(21)

上式即為懸掛點(diǎn)4上的約束條件。

2 數(shù)值算例

本文在程序?qū)崿F(xiàn)中使用了開源有限元框架deal.II[29-30]。所使用的計算平臺配備了Intel Xeon E5 2680 V4 CPU，包含14個處理器核心，主頻為2.4GHz，內(nèi)存為128GB。

2.1 正確性驗證

采用均勻半空間模型驗證形函數(shù)階數(shù)為3時算法的正確性。模型的計算區(qū)域為200m×200m×200m，點(diǎn)源位于原點(diǎn)O(0,0,0)，沿x軸正方向以2m的間隔布置20個測點(diǎn)。

由于此模型較為簡單，電位解析解[31]計算公式為

(22)

式中R為測點(diǎn)與點(diǎn)源的距離。

通過有限元模擬各測點(diǎn)的電位，再計算得到視電阻率曲線(圖7a)；通過與解析解對比，得到所有測點(diǎn)的相對誤差(圖7b)。從圖7a可以看出：距離點(diǎn)源越遠(yuǎn)，有限元的解越精確；隨著迭代次數(shù)的增加，有限元解快速收斂于解析解。從圖7b可以看到，隨著自由度個數(shù)的增加，所有測點(diǎn)的相對誤差大幅度降低。表1給出了第1次、第3次以及第5次網(wǎng)格剖分方案下的自由度個數(shù)、計算時間和測點(diǎn)的最大相對誤差，最大相對誤差從124.8%降至0.3%，證明本文算法具有很高的精度。

表1 第1次、第3次和第5次網(wǎng)格剖分方案下的自由度個數(shù)、測點(diǎn)視電阻率最大相對誤差及耗時統(tǒng)計

圖7 第1次、第3次和第5次網(wǎng)格剖分方案下的視電阻率曲線(a)及相對誤差(b)

2.2 收斂速度

應(yīng)用均勻半空間模型分析形函數(shù)階數(shù)p=1、2、3時自適應(yīng)有限元解的收斂速度。計算區(qū)域為500m×500m×500m，背景電阻率為100Ω·m，初始剖分網(wǎng)格為8000個單元。從源點(diǎn)S(0,0,0)處向地下注入電流強(qiáng)度為1A的電流，沿x軸正方向100～400m范圍內(nèi)等間隔布置31個測點(diǎn)，間距為10m。對不同的網(wǎng)格剖分方案分別計算各測點(diǎn)的有限元解，然后參照解析解計算所有測點(diǎn)的相對誤差，最終用所有測點(diǎn)的視電阻率平均相對誤差展示有限元解的收斂速度，結(jié)果如圖8所示。表2列出了形函數(shù)階數(shù)p=1、2、3時最終網(wǎng)格上的自由度個數(shù)和所有測點(diǎn)的電阻率有限元解的平均相對誤差。

從圖8和表2可以看出，對于粗糙網(wǎng)格，形函數(shù)的階數(shù)越高，有限元解的精度越高；有限元解的相對

圖8 形函數(shù)p=1、2、3時有限元解的收斂速度對比

表2 不同p值時網(wǎng)格自由度個數(shù)和所有測點(diǎn)視電阻率平均相對誤差統(tǒng)計

誤差隨著迭代次數(shù)的增加逐漸下降，p=1時有限元解曲線下降最慢，且最終網(wǎng)格上自由度個數(shù)最多為8256929，而所有測點(diǎn)的有限元解的平均相對誤差僅降至1.5×10-4；p=2時有限元解的收斂速度有所提升，最終網(wǎng)格上的自由度個數(shù)為3277344，大約是p=1時最終網(wǎng)格上自由度個數(shù)的5/8，所有測點(diǎn)有限元解的平均相對誤差降至1.4×10-5，即有限元解的精確度相較于p=1時提高了約11倍；p=3時有限元解的收斂性能最佳，最終網(wǎng)格上自由度個數(shù)為1968695，所有測點(diǎn)的有限元解的平均相對誤差為8.8×10-6，相較于p=2時，自由度個數(shù)減少了1308649，精度提高了約1.5倍。綜上所述，p=3時有限元程序可以用最少的自由度個數(shù)得到最精確的解。

圖9展示了形函數(shù)階數(shù)p=1、2、3時最終網(wǎng)格剖分?jǐn)?shù)目和誤差分布。眾所周知，在點(diǎn)源附近電勢變化劇烈，導(dǎo)致解的光滑度較低，在距離點(diǎn)源較遠(yuǎn)的區(qū)域中解的光滑度較高。從圖中可以看到，在距離點(diǎn)源較遠(yuǎn)、解比較光滑的區(qū)域中，形函數(shù)的階數(shù)越高，最終網(wǎng)格上有限元解的誤差越??；形函數(shù)的階數(shù)p越高，最終整體誤差越小。在距點(diǎn)源非常近的區(qū)域，p=1、2時網(wǎng)格得到了充分加密，而p=3時網(wǎng)格未得到充分加密。由于點(diǎn)源附近網(wǎng)格加密程度不同，從圖9c(右)可以清楚看到，p=3時在距離點(diǎn)源非常近的區(qū)域，有限元解的誤差較高；而p=1、2時點(diǎn)源附近的誤差較p=3時明顯降低。

圖9 p=1(a)、2(b)、3(c)時最終網(wǎng)格剖分方案(左)及誤差分布(右)

2.3 異常體模型

對低阻異常體模型(圖10)進(jìn)行模擬，以驗證算法的有效性。異常體是一個電阻率為10Ω·m的低阻長方體，位于坐標(biāo)原點(diǎn)正下方。背景電阻率為100Ω·m。

圖10 低阻異常體模型示意圖

采用對稱四極裝置模擬計算沿x軸的視電阻率擬斷面。x方向的計算范圍為-400～400m，測點(diǎn)間距為20m；z方向為100～600m。供電電極的極距由200m逐漸增加到1200m，測量電極的極距為相應(yīng)供電電極距的1/3。

首先，計算得到沿x軸的視電阻率擬斷面(圖11a)，可以清楚地看到，在x軸-100～100m范圍內(nèi)有一個低阻異常區(qū)，其寬度與模型異常體基本一致。

然后，在地表平行于x軸布置41條測線，測線范圍為y=-200～200m，測線間距為10m，沿x方向測點(diǎn)的布設(shè)同圖11a。供電電極極距仍為400m。截取視深度z=200m的數(shù)據(jù)繪制視電阻率平面圖(圖11b)?？梢钥吹剑笾掠蓌軸-100～100m、y軸-50～50m所圍區(qū)域是一個低阻區(qū)，其中心位置與模型異常體的中心位置吻合，覆蓋范圍與模型也基本一致。

圖11 x軸視電阻率擬斷面(a)及視深度為200m(b)的視電阻率平面圖黑色虛線為異常體邊界

2.4 合成數(shù)據(jù)模擬

為驗證高階自適應(yīng)有限元算法對較復(fù)雜模型的模擬效果，本文基于沁水盆地南部某煤層氣壓裂監(jiān)測資料，建立三維電性模型驗證本文算法的有效性。

沁水盆地南部煤層埋藏淺、厚度大、煤儲層中煤層氣含量高，且具有相對較高的滲透率，壓裂后富水地層較上覆和下伏地層通常表現(xiàn)為低阻特征[32-33]。模型的計算區(qū)域為2000m×2000m×1000m。煤層氣異常體在x、y和z軸方向的范圍分別為-76～92m、-50～50m、124～190m，異常體示意圖如圖12所示。模型的背景電阻率為100Ω·m，壓裂前異常體電阻率設(shè)為10000Ω·m，壓裂后為1Ω·m。

圖12 煤層氣模型示意圖

采用對稱四極電阻率測量裝置。測線沿x方向布設(shè)，測線范圍y為-80～80m，測線間隔5m，測點(diǎn)范圍x為-100～100m，間隔為5m。供電電極的極距為300m。分別計算煤層壓裂前、后的視電阻率分布，并從計算結(jié)果中提取x方向-100～100m、y方向-80～80m、z=150m的視電阻率數(shù)據(jù)。通過定量計算得到壓裂后的視電阻率相對異常

(23)

式中B、C分別為壓裂前、后的視電阻率。圖13a、圖13b分別為第一次和壓裂結(jié)束后的視電阻率相對異常分布。

第一次壓裂后，x方向-74～-40m范圍內(nèi)充滿壓裂劑，該區(qū)域電阻率約1Ω·m(即壓裂劑的電阻率)，其他范圍內(nèi)異常體的電阻率值仍很高，約為10000Ω·m。計算得到對稱四級裝置條件下z=150m平面的視電阻率后，通過定量計算得到了視電阻率相對異常圖(圖13左)。因為壓裂范圍較小，所以圖中左側(cè)區(qū)域顯示出小范圍的低阻異常(圖13中藍(lán)色的低阻區(qū)域)。

壓裂結(jié)束后，異常區(qū)域內(nèi)全部充滿壓裂劑，這些區(qū)域視電阻率主要體現(xiàn)為壓裂液的低阻(1Ω·m)特征(圖13右)。圖中的藍(lán)色低異常區(qū)域與模型(圖12)中的煤層氣分布區(qū)域基本一致。

分析圖13可知，采用本文算法并結(jié)合對稱四級電阻率測量裝置，可以定性地反演異常體的水平分布，計算結(jié)果與模型相吻合。因此，該方法在煤層氣壓裂監(jiān)測領(lǐng)域具有重要應(yīng)用價值。

圖13 第一次壓裂(左)及壓裂完成后(右)z=150m平面視電阻率相對異常分布圖

3 結(jié)論

本文應(yīng)用高階自適應(yīng)有限元算法實現(xiàn)了直流電阻率模型正演，算例表明該算法具有較高的精度。

(1)形函數(shù)階數(shù)p=3時，數(shù)值模擬有限元解收斂最快，可以用少量的自由度個數(shù)得到精確的有限元解。

(2)在點(diǎn)源附近電勢劇烈變化的區(qū)域，可通過加密網(wǎng)格提高解的精度；在離點(diǎn)源較遠(yuǎn)處，解析解比較光滑，可通過提高單元上形函數(shù)的階數(shù)提高解的精度，效果顯著。

(3)合成數(shù)據(jù)模型的模擬結(jié)果證明三維直流電阻率法是煤層氣儲層壓裂監(jiān)測的一項重要補(bǔ)充技術(shù)。

感謝河南理工大學(xué)對本項目提供了高性能計算資源！