高階自適應有限元三維直流電阻率正演方法及其在沁水盆地煤氣層壓裂監測中的應用

趙 寧 黃明衛 申亞行 陶德強 秦 策*

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

0 引言

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

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

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

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

1 正演理論

1.1 控制方程

假設在地面A點放置一個電流強度為I的點電源(圖1)，空間電位分布滿足

·(σU)=-IδA

(1)

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

(2)

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

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

圖1 點源示意圖

(3)

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

(4)

式中r是由點A指向邊界Γ的方向向量。

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

(5)

1.2 有限元離散

將區域Ω離散為一系列非結構化六面體單元，每個單元e上電位的近似解為

(6)

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

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

(7)

keue=fe

(8)

(9)

(10)

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

KU=F

(11)

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

1.3 形函數的生成

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

(12)

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

(13)

(14)

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

m=(p+1)d

(15)

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

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

(16)

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

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

1.4 后驗誤差估計方法

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

(17)

(18)

計算，其中l是單元中面數目。

得到所有單元的誤差估計后，對誤差前10%的單元進行加密得到新的網格，然后在新的網格上重新計算有限元解，直到迭代次數或自由度個數達到設定值，迭代結束。h型自適應有限元算法流程如圖4所示。

圖4 h型自適應有限元算法流程圖

1.5 懸掛點上的約束條件

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

圖5 八叉樹加密原理

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

圖6 懸掛點示意圖自由度4是懸掛點

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

(19)

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

(20)

當通過1.3節中的步驟求得自由度0、1和2上的形函數表達式后，可以得到

(21)

上式即為懸掛點4上的約束條件。

2 數值算例

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

2.1 正確性驗證

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

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

(22)

式中R為測點與點源的距離。

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

表1 第1次、第3次和第5次網格剖分方案下的自由度個數、測點視電阻率最大相對誤差及耗時統計

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

2.2 收斂速度

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

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

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

表2 不同p值時網格自由度個數和所有測點視電阻率平均相對誤差統計

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

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

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

2.3 異常體模型

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

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

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

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

然后，在地表平行于x軸布置41條測線，測線范圍為y=-200～200m，測線間距為10m，沿x方向測點的布設同圖11a。供電電極極距仍為400m。截取視深度z=200m的數據繪制視電阻率平面圖(圖11b)。可以看到，大致由x軸-100～100m、y軸-50～50m所圍區域是一個低阻區，其中心位置與模型異常體的中心位置吻合，覆蓋范圍與模型也基本一致。

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

2.4 合成數據模擬

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

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

圖12 煤層氣模型示意圖

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

(23)

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

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

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

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

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

3 結論

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

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

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

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

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