基于MCMC 的頁巖氣藏壓裂縫網不確定性反演法

賈品，封皓然，王蘇冉，王志凱，牛烺昱，程林松

1 中國石油大學(北京)油氣資源與工程全國重點實驗室，北京 102249

2 中國石油大學(北京)石油工程學院，北京 102249

3 中海油研究總院有限責任公司，北京 100028

0 引言

我國頁巖氣資源分布廣泛且儲量豐富，可采資源量達15×1012～30×1012m3，近些年先后設立了涪陵、長寧—威遠和昭通等國家級頁巖氣示范區，3500 m以上淺儲層地質壓裂工程一體化技術基本成型，實現了我國淺層頁巖氣的工業化生產；3500 m以下深儲層頁巖氣也得到有效的嘗試性開發，發展潛力巨大。不論對于淺層還是深層頁巖氣，由于頁巖基質以微納米級孔隙為主，儲層滲透率極低、物性更差，能夠形成三維立體復雜裂縫網絡的水平井體積壓裂技術成為了頁巖氣經濟有效開發的主要手段。其中，準確反演壓裂縫網參數是壓裂效果評價、壓裂優化設計及產能設計的關鍵問題，也一直是頁巖氣藏壓裂及開發工程方面的熱點問題。

頁巖氣藏壓裂縫網的反演通常是基于一個適合于目標區的滲流數學模型，通過不斷調整縫網參數對實際生產動態數據或其相關變形量進行擬合從而實現縫網關鍵參數的推算。目前常見的線性流分析法、圖版擬合法及生產歷史擬合法均屬于此類方法的范疇。線性流分析法通過繪制擬壓力規整化產量和物質平衡擬時間的雙對數曲線，找出曲線上的線性流段，擬合線性流段的斜率和截距獲得縫網參數。國內外大多數研 究 都 采 用Wattenbarger[1]、Bello和Wattenbarger[2]、Al-Ahmadi和Wattenbarger[3]、Brown[4]、Stalgorova和Mattar[5]等提出的線性流模型。Nobakht和Matter[6]、Clarkson等[7]、Qanbari和Clarkson[8]、Xu等[9]通 過 確定的原始地層壓力，儲層厚度，壓裂段間距(通常認為段間全改造)，基質及裂縫的壓縮系數，基質及裂縫的孔隙度，及氣體吸附解析參數等參數對壓裂縫長度、導流能力及改造區大小、滲透率進行反演。圖版擬合法是基于縫網流動模型繪制裂縫關鍵參數影響下的產量敏感性圖版，運用理論圖版擬合實際生產，從而反演縫網參數。由于簡化的物理模型，線性流分析法、圖版擬合法能夠反演的參數有限。Moghadam[10]等和Abdulal等[11]應用圖版擬合法反演了頁巖氣井儲層參數。生產歷史擬合法多用于復雜縫網形態的近似表征，Li和Zhu[12]研究了水平井壓裂過程的溫度不穩態響應，Zhang和Zhu[13]進一步提出了裂縫長度和導流能力反演方法。結合微地震和巖心數據，Cipolla等[14]、Weng等[15]提出了一種利用微地震事件位置校準地質力學模型并生成縫網的方法，Gamboa等[16]和Sun等[17]進一步通過考慮頁巖氣水排采動態建立了一種壓裂縫網展布形態解釋方法。基于生產歷史擬合的縫網參數反演存在歷史擬合數據量大、耗時耗力及適用性差等不足。

確定性反演需要相對準確的輸入參數，才能得到合理的反演值。由于滲流數學模型本身包含了眾多因素，例如儲層參數、流體性質及非線性滲流機理，這都給頁巖氣壓裂縫網反演帶來了極大的不確定性。因此，頁巖氣藏壓裂縫網反演應該不是一個確定性過程，而是由于基礎參數的不確定性引起的一個不確定過程。上述的3 類方法均是確定性反演方法，即給定了基礎參數后，便反演出了一套縫網參數，人們還無法明確縫網反演參數的多組合性及其概率分布，難以客觀并全面評估頁巖氣壓裂改造效果。近幾年，基于馬爾可夫—蒙特卡洛方法(MCMC)的不確定反演方法已經廣泛應用于巖石物理領域[18-19]，但在頁巖氣藏領域還僅停留在基于MCMC的不確定性產能評價，針對壓裂縫網參數不確定性的反演研究鮮少，亟需提出一種適用于頁巖氣藏壓裂縫網參數的不確定性反演方法。

本文針對頁巖氣藏壓裂縫網參數在歷史擬合過程中的不確定性，基于馬爾可夫—蒙特卡洛方法(MCMC)的原理，提出一種適用于頁巖氣藏的壓裂縫網反演方法。通過隨機生成反演所需的基礎參數樣本集，構建縫網參數解釋模型，聯合線性流分析法和圖版擬合法，獲取縫網反演參數的多組合性及其概率分布，評價由于基礎參數的不確定性對壓裂縫網參數解釋結果的影響程度，對頁巖氣藏壓裂效果評價及合理開發對策制定提供較全面且實際的指導。

1 壓裂縫網確定性反演

如上所述，頁巖氣藏壓裂縫網的不確定性反演要基于壓裂縫網的確定性反演，即通過隨機抽樣方法產生若干基礎參數的樣本集，針對每個樣本集進行對應的壓裂縫網確定性反演，然后收集若干基礎參數樣本下的反演結果，從而得到壓裂縫網參數的組合及概率分布。以下首先給出頁巖氣藏壓裂縫網參數的確定性反演方法，然后闡述壓裂縫網反演所需基礎參數的不確定性，從而建立不確定性反演方法，最后進行實例分析。

1.1 頁巖氣藏壓裂水平井物理模型

如圖1 所示，頁巖氣藏一口水平井經體積壓裂后，形成了復雜的裂縫網絡系統，即壓裂改造體(SRV)。整個儲層可等效為由人工裂縫、壓裂改造區(內區)和未改造區(外區)三個區域組成的系統。人工裂縫考慮為等長雙翼縫；內區考慮壓裂改造形成的復雜縫網系統，將其處理成高滲透區域；外區由于幾乎沒有受到壓裂改造，將其處理成低滲透區域。其他基本假設為：

圖1 頁巖氣藏壓裂縫網物理模型示意圖Fig.1 Schematic diagram of physical model of fracture network in shale gas reservoir

(1)儲層為水平等厚，水平井位于儲層中心，裂縫貫穿整個儲層；

(2)考慮氣體及巖石的壓縮性，不考慮重力作用；

(3)各區域內氣體流動都視為線性流，即外區系統流體線性流入內區系統，再由內區系統線性流入人工裂縫，最后由人工裂縫線性流入水平井筒(如圖1 所示)。

1.2 壓裂縫網確定性反演方法

頁巖氣主要以游離態、吸附態和溶解態3 種形式賦存于頁巖儲層中，其中游離態和吸附態占比最高。Al-Ahmadi和Wattenbarger[20]針對考慮頁巖氣吸附解析及非線性流動機理的三線性流模型做了大量工作，本文在此基礎上構建基于產量遞減分析的壓裂縫網參數反演方法。對于圖1 所示的物理模型，其在定壓生產條件下拉普拉斯空間的無因次產量解為[21]：

式中，f(s)定義為：

其中，qD為無因次產量，yFD為無因次裂縫半長，yeD為無因次外區寬度，λIF為內區與人工裂縫竄流系數比，λIo為內區與外區竄流系數比，ηD為無因次導壓系數，ω為彈性儲容系數，s為拉普拉斯變量。以上所有變量均為無因次，其具體定義可參考筆者文章[22]。

大量礦場頁巖氣壓裂井開發動態表明生產早期頁巖氣流動通常會出現較為明顯的線性流特征，生產中后期流態以邊界控制流為主[23]，直線分析方法是分析線性流及邊界控制流的重要技術手段。直線分析法基于線性流或邊界控制流階段壓力及產量與時間的局部解，通過繪制特定變量的坐標圖，利用斜率或截距推算壓裂縫網參數，如表1 所示。本文首先采用直線分析方法定量分析頁巖氣產量動態數據，利用線性流獲取人工裂縫半長與改造區滲透率關系利用邊界控制流獲得改造體積大小，再以這些值為初值或約束條件利用圖版擬合法反演得到人工裂縫半長、人工裂縫導流能力、改造區滲透率及改造區寬度參數。

表1 直線分析法變量及反演參數表Table 1 Table of variables and inversion parameters of linear analysis method

實際礦場生產中，頁巖氣壓裂井往往是以變產量或變井底流壓生產，生產動態分析需引入物質平衡擬時間和擬壓力規整化產量分析變工作制度下生產動態。其中，物質平衡擬時間項為：

其中，Wg是累積產氣量；Vpm是改造區系統的總孔隙體積，m3；Smgi為初始含氣飽和度；Bgi是原始地層壓力下的氣體體積系數，Bg是平均壓力下的氣體體積系數，是壓力的函數;pi為原始地層壓力，MPa。平均壓力的求解可以使用牛頓迭代法，是改造區系統的總孔隙體積可以表示為：

其中h為儲層厚度，m；xf為人工裂縫半長，m；這里人工假設半縫長初值用于計算物質平衡擬時間，與解釋出的結果對比迭代。Ls為改造區寬度，m；通過統計水平井長度和壓裂段數可折算得到。

擬壓力歸整化產量項為：

式中，ψi為原始地層壓力下對應的擬壓力，MPa2/mPa·s；ψw為井底流壓對應的擬壓力，MPa2/mPa·s。以一口體積壓裂水平井數值模擬模型驗證上述線性流分析法和圖版擬合法聯合解釋壓裂縫網參數的流程，模型示意如圖2 所示，水平段長度為1500 m，壓裂25 段，段間完全改造。人工裂縫滲透率為5000 mD，人工裂縫半長為102.5 m，改造區滲透率為5×10-4mD，改造區寬度為30 m，壓裂液僅賦存于人工主縫中，基質內不含水。

圖2 數值模型示意圖Fig.2 Schematic diagram of numerical model

基于線性流分析法和圖版法聯合解釋壓裂縫網參數，具體解釋流程為：(1)通過實際產氣量及井底流壓動態數據計算物質平衡擬時間和歸一化擬壓力，繪制兩者的雙對圖版曲線并識別流動階段，如圖3 所示。圖3 中可以明顯看出線性流，邊界控制流兩個流動階段，在邊界流之前，由于數據點時間間隔過短，沒有明顯的雙線性流階段，通過流態的劃分可以進一步驗證三線性流的流動階段劃分合理；(2)提取雙對數圖版線性流階段曲線特征，擬合直線段斜率mL(圖4)，利用表1 中的關系式確定人工裂縫半長與改造區滲透率關系，進一步提取雙對數圖版中的邊界控制流階段曲線特征，擬合直線與x軸的截距(圖5)，利用表1中的關系式確定改造體大?。?3)以線性流擬合關系式及邊界控制流得到的改造體大小為約束，不斷調整人工裂縫半長、人工裂縫滲透率、改造區滲透率及改造區寬度參數，擬合全生產歷史動態下的物質平衡擬時間和規整化產量圖版，反演縫網參數，如圖6 所示。反演參數見表2，可以看出壓裂縫網的4 個關鍵參數反演值與數值模型值幾乎一致，吻合度均在94.4%以內，驗證了確定性反演方法的準確性和可靠性。

表2 數值模型參數反演參數表Table 2 Numerical model parameter inversion parameter table

圖3 生產不同階段流態劃分Fig.3 Flow regime division in different production stages

圖4 線性流直線分析圖Fig.4 Linear flow line analysis diagram

圖5 邊界控制流直線分析圖Fig.5 Boundary -dominated flow analysis diagram

圖6 產量遞減圖版擬合結果Fig.6 Production decline type curve fitting

2 基礎參數不確定性分析

壓裂縫網反演過程中所需基礎參數的不確定性影響著線性流分析和圖版擬合的參數解釋結果。頁巖氣藏壓裂縫網參數反演的不確定性主要來源于以下兩方面：一方面由于輸入參數多且適用范圍有限，如天然巖心及測井資料獲取的儲層參數本質上僅代表了近井筒附近儲層物性，而無法準確反映儲層整體性質；另一方面復雜縫網多參數反演易出現局部最優解，圖版擬合的多解性強。以下給出了頁巖氣藏壓裂縫網反演中的不確定性參數分析。

頁巖儲層參數不確定性：頁巖儲層參數主要包括儲層有效厚度、原始地層壓力、原始地層溫度、含氣飽和度、基質孔隙度。頁巖縱向不同埋深及平面井間參數存在差異，基于目標區儲層物性評價資料可確定該類參數的分布及范圍。

頁巖吸附解析及流動參數不確定性：頁巖不同井區的干酪根及初始含氣量存在差異，且含氣量隨著儲層埋深增加而增加。另外，由于不同區域儲層頁巖基質孔隙直徑變化大，氣體流動參數在不同位置處也不同。因此，頁巖吸附解析及流動參數不確定性分別體現在蘭氏體積、蘭氏壓力和擴散系數中。

頁巖壓裂縫網基礎參數不確定性：頁巖儲層經過大型水力壓裂后，人工裂縫半長、人工裂縫滲透率、改造區滲透率及改造區寬度參數需要反演，但裂縫壓縮系數、改造區孔隙度、改造區壓縮系數需要作為反演的基礎參數提前提供，該類參數具有不確定性。通過統計目標區塊的各井的原始地層壓力，儲層有效厚度，朗格繆爾體積，朗格繆爾壓力，水平井長度，基質孔隙度，基質壓縮系數，裂縫壓縮系數，選擇合適的概率密度函數(表3)進行擬合，得到針對該區塊的不同儲層參數的概率密度函數。

表3 常用概率密度函數Table 3 Common probability density function

3 基于MCMC 的壓裂縫網不確定反演

本文基于馬爾可夫-蒙特卡洛方法(MCMC)對頁巖氣藏水平井壓裂縫網參數進行不確定性反演?；陧搸r氣藏壓裂水平井基礎參數，采用MCMC抽樣方法獲取儲層參數、吸附解析及流動參數、壓裂縫網基礎參數的概率分布，結合線性流分析法和圖版擬合法識別典型流動階段，擬合生產動態，獲取各反演參數的概率分布。MCMC抽樣方法的基本思想是采用Metropolis抽樣方法，構造一個遍歷的馬爾可夫鏈，從任意概率向量出發，通過多次迭代計算，獲得不變分布，使其成為所需要的抽樣分布。抽樣次數越大，抽樣分布越趨于一個穩定分布。MCMC抽樣基于貝葉斯理論，其表達式為：

① 生成基礎參數：設置儲層、流動及裂縫基礎參數的初始樣本θs，通常取各參數分布的中值；

② 基礎參數分布抽樣：根據區塊基礎參數分布概率(統計獲得)，產生新樣本θp；

③ 計算新樣本θp接收概率p：接受概率p的計算方法：

④ 判斷是否更新樣本：產生一個(0，1)隨機數r，若r＜p，接受新樣本θs=θp，反之θs=θs-1；

⑤ 判斷是否結束：若達到最大迭代步數，則抽樣結果統計，否則返回步驟②。具體的不確定性反演方法流程如圖7 所示。

圖7 頁巖氣藏壓裂縫網參數不確定性反演流程圖Fig.7 Flow chart of parameter uncertainty inversion of pressure fracture network in shale gas reservoir

4 實例應用

4.1 基礎參數不確定性分析

以我國西南地區某頁巖氣藏一口體積壓裂水平井為例進行分析，該氣藏埋深約3500 m，儲層有效厚度約30 m，儲層中深壓力為67.3 MPa，地層溫度375 K，初始含水飽和度約為0.45。巖心分析結果表明，頁巖基質孔隙度約為0.04，滲透率跨度較大(10-5mD～10-3mD)。該水平井分為26 段壓裂，每段壓裂3 簇，壓裂總段長1430 m，入地總液量53831 m3，入地總砂量2442.4 t。該頁巖氣藏壓裂水平井的產氣動態及井底流壓如圖8 所示。根據該地區臨井基礎資料分析，獲取了儲層參數、流動參數及壓裂縫網基礎參數的范圍及概率分布，如表4 所示。

表4 各相關參數取值Table 4 Values of relevant parameters

圖8 生產動態數據圖Fig.8 Production dynamic data fitting

4.2 壓裂縫網參數解釋

基于表4 中基礎參數的概率分布，利用MCMC隨機生成至少1000 個基礎參數樣本集合，再利用不同的基礎參數樣本結合等式(3)和(4)計算物質平衡擬時間及歸一化擬壓力，繪制出1000 個頁巖氣流動特征雙對數曲線。為了便于展示，圖9(a)給出了其中基于6 個基礎參數樣本的曲線，可以看出生產早期氣體流動表現出斜率為1/2 的線性流特征，表明改造區得到有效動用，改造區內氣體線性流入人工裂縫中；隨后出現斜率為1 的邊界控制流，改造區壓力不斷衰竭。如圖9(b)和9(c)所示通過擬合線性流階段和邊界控制流階段的直線段，分別可獲得人工裂縫半長與改造區滲透率關系以及改造體大小，進一步通過擬合物質平衡擬時間和歸一化產量圖版，如圖9(d)所示，進而反演獲得壓裂縫網在每個基礎參數樣本下的參數。

圖9 實際井壓裂縫網參數不確定性反演結果Fig.9 Uncertainty inversion results of actual well pressure fracture network parameters.(a) Flow regime division in different production stages (b) Linear flow line analysis diagram (c) Boundary-dominated flow analysis diagram (d) Production decline type curve fitting

4.3 壓裂縫網參數不確定性分析

基于圖9 所展示的壓裂縫網反演流程，獲得了人工裂縫半長、人工裂縫導流能力、改造區滲透率及改造區寬度參數組合下的不確定分布，其概率分布以及累積概率分布如圖10～13 所示。同時也給出了各反演參數的擬合值及P10、P50、P90 下的反演結果，見表5。

表5 W1 井裂縫參數不確定性反演小結Table 5 Summary of uncertainty inversion of fracture parameters in well W1

圖10 人工裂縫長度分布圖Fig.10 Distribution of hydraulic fracture length.(a) Uncertainty inversion probability distribution diagram (b) Cumulative probability distribution

圖11 人工裂縫導流能力不確定性分析結果Fig.11 Uncertainty analysis results of artificial fracture conductivity.(a) Uncertainty inversion probability distribution diagram (b) Cumulative probability distribution

圖12 改造區滲透率不確定性分析結果Fig.12 Analysis results of permeability uncertainty in the reconstruction area.(a) Uncertainty inversion probability distribution diagram (b) Cumulative probability distribution

圖13 改造區寬度不確定性分析結果Fig.13 Uncertainty analysis results of width of reconstruction area.(a) Uncertainty inversion probability distribution diagram(b) Cumulative probability distribution

從人工裂縫半長不確定性反演分析圖可以看出，反演得到的人工裂縫半長基本符合正態分布且跨度較大，同時累計頻率分布P10、P50、P90 對應的參數數值大小較為分散，表明其不確定性較強。

從人工裂縫導流能力不確定性反演分析圖可以看出，反演得到的人工裂縫導流能力基本符合正態分布且跨度較大，同時累計頻率分布P10、P50、P90 對應的參數數值大小較為分散，表明其不確定性較強。

從改造區滲透率不確定性反演分析圖中可以看出，反演得到的改造區滲透率也基本符合正態分布且跨度較大但整體范圍較小，同時累計頻率分布P10、P50、P90 對應的參數數值大小較為分散，表明其不確定性一般。

由改造區寬度不確定性反演分析圖中可以看出，反演得到的改造區寬度的概率分布較為集中且主要集中在45 m左右，且累計頻率分布P10、P50、P90 對應的參數數值大小也較為相近，表明其不確定性較弱。

綜合分析基于隨機生成的基礎參數樣本反演得到的縫網參數不確定性可得：該地區頁巖氣藏體積壓裂后人工裂縫半長和人工裂縫導流能力的不確定性較強，改造區滲透率不確定性一般，改造區寬度不確定性較弱。通過對壓裂裂縫參數不確定性的分析，即可為實際區塊中壓裂裂縫參數多解性提供參考。

5 結論

本文基于頁巖氣單井產能模型，首先構建了壓裂縫網參數確定性反演方法，然后結合馬爾可夫—蒙特卡洛方法(MCMC)形成了壓裂縫網參數不確定性反演方法，最后對我國西南地區一口頁巖氣井壓裂縫網的4 個關鍵參數進行了不確定性反演，得到了不同基礎參數下的壓裂裂縫反演參數的概率分布及累積概率分布，得到結論如下：

(1)構建了頁巖氣壓裂井確定參數下的儲層參數反演方法，結合直線分析法和圖版擬合法對人工裂縫半長、人工裂縫滲透率、改造區滲透率和改造區寬度進行確定性解釋，并結合數值模型論證了方法的可靠性；

(2)頁巖儲層縫網參數反演參數多、不確定性大，基于馬爾可夫-蒙特卡洛方法(MCMC)的頁巖氣藏壓裂縫網不確定性反演方法主要包括蒙特卡洛抽樣，典型流動階段分析及生產動態擬合等步驟，能快速獲取人工裂縫長度、導流能力、改造區滲透率及改造區裂縫寬度等關鍵縫網參數及其概率分布；

(3)針對我國西南地區某頁巖氣區塊一口體積壓裂水平井生產動態，對其進行不確定性反演，明確了壓裂縫網四個關鍵參數的不確定性排序，人工裂縫半長的擬合值為75 m，累計概率分布P10 為40.5 m，P50為69.3 m，P90 為104.3 m。人工裂縫導流能力的擬合值為125 mD，累計概率分布P10 為84.1 mD，P50 為166.9 mD，P90 為270.3 mD。改造區滲透率的擬合值為5.0×10-4mD，累計概率分布P10 為4.4×10-4mD，P50 為5.8×10-4mD，P90 為7.0×10-4mD。改造區寬度的擬合值為45.3 m，累計概率分布P10 為40.2 m，P50 為45.3 m，P90 為47 m。