基于機器學習方法的海洋天然氣水合物水平井降壓開采模擬—優化耦合模型

辛 欣 王海彬 羅建男 于 涵 袁益龍 夏盈莉 朱慧星 陳 強

1.吉林大學地下水資源與環境教育部重點實驗室 2.自然資源部天然氣水合物重點實驗室·中國地質調查局青島海洋地質研究所 3.青島海洋科學與技術試點國家實驗室海洋礦產資源評價與探測技術功能實驗室

0 引言

天然氣水合物（以下簡稱水合物）以填充、膠結等多種形式賦存于沉積物中，通過降壓開采打破水合物形成的壓力平衡條件，可以使固態水合物分解為氣態甲烷和液態水[1]，從而使含水合物沉積層的抗剪強度和承載能力降低[2]，同時，降壓開采使含水合物沉積層壓力降低，有效應力增加?；谏鲜鲈颍祲洪_采易導致含水合物沉積層變形，從而產生海底面沉降、變形等地質災害[3-4]。因此，在水合物的開采過程中，如何保證在獲得較高產氣量的同時，將海底面沉降量控制在一定的安全范圍內，成為實現水合物安全、高效規?；_采的前提條件。

目前國內外學者多運用數值模擬的方法來評價不同開采方案下含水合物沉積層的產氣能力及海底面沉降量。Rutqvist等[5]在水合物數值模擬程序TOUGH+HYDRATE上鏈接巖土力學軟件FLAC3D，針對加拿大Mallik凍土區和美國阿拉斯加凍土區模擬預測了降壓開采水合物過程中的地層力學響應，結果表明影響含水合物沉積層剪切破壞程度的主控因素是降壓幅度。Kimoto等[6]利用化學—熱學—力學耦合模型，模擬開采水合物引起的沉積層力學變形過程，結果表明地層變形是基于水合物相態變化引起的松散沉積物應力強度降低。Gupta等[7]將含水合物沉積層力學場變化考慮為線彈性變形，基于有限元法建立了模擬水合物開采的流動、傳熱、力學多場耦合模型，并運用多個解析解和實驗數據對該模型進行了驗證。Jin等[8-9]利用自主開發的HYDRATEBiot程序模擬預測了南海神狐海域水合物在降壓開采過程中的產氣能力及海底面沉降量，結果表明在降壓開采初期，海底面發生快速沉降，而在開采后期，沉降速率變慢且趨于平緩。

上述研究均采用數值模擬手段預測了既定開采方案下含水合物沉積層的力學響應特征及其對甲烷累計產氣量的影響，但無法解決考慮含水合物沉積層力學變形影響下的水合物開采方案優選問題。而將數值模擬技術與運籌學中的優化技術相結合，構建模擬—優化耦合模型，既能預測水合物開采引起的多相流體滲流過程和溫度場、壓力場、應力場的變化趨勢，又能夠在滿足各種環境、經濟、技術等要求的前提下獲得最優的開采方案。因此，將數值模擬模型（以下簡稱模擬模型）與目標優化模型（以下簡稱優化模型）耦合的技術（Simulation and Optimization Technology）[10]為解決以水合物安全、高效開采過程為目標的開采方案優選問題提供了科學途徑。而目前關于應用模擬—優化耦合技術進行水合物安全、高效開采方案優選的報道還鮮見。

為了找到甲烷累計產氣量最優值與地層穩定性的關系，基于機器學習方法形成模擬—優化耦合技術，構建起水合物降壓開采傳熱—流動—力學數值模擬模型、可以替代數值模擬模型的機器學習模型和以甲烷累計產氣量最優為目標的混合整數非線性規劃優化模型；在此基礎上，選取南海神狐海域厚層Ⅱ類水合物藏W11站位為研究對象，獲得了海底面沉降量約束下的水合物儲層甲烷累計產氣量以及對應的最優開采方案參數，以期為水合物安全、高效開采方案的制訂提供支撐。

1 基于機器學習方法的模擬—優化耦合技術原理

模擬—優化耦合技術包括模擬模型和優化模型兩個部分，其中模擬模型是核心，能夠描述真實地質體輸入、輸出響應關系。建立水合物開采模擬模型包含以下3個步驟：①建立描述水合物儲層空間結構的概念模型；②建立水合物傳熱—流動—力學耦合數學模型；③運用數值模擬程序HYDRATEBiot求解。優化模型是基于運籌學理論的最優方案求解模型，能夠突破既定開采方案的限制并且求得全局最優解[10-11]。優化模型的構建方法包括線性規劃、非線性規劃、整數規劃等多種方法[10-11]。水合物開采引起的多相流體滲流及地層力學特征的響應過程具有典型的非線性特征，因此運用混合整數非線性規劃方法來構建優化模型，然后采用遺傳算法進行求解。

將模擬模型作為一個目標函數或約束條件編入優化模型中，在利用優化模型求解全局最優解時需要反復多次調用模擬模型，模擬模型的非線性特征越強，調用計算的負荷越大，嚴重時甚至出現無法運算的情況[12]，從而限制了模擬—優化耦合技術的應用。近年來，機器學習的理念逐漸滲透進入越來越多的研究領域，通過學習模擬模型計算得到的輸入—輸出特征響應關系，建立起模擬模型的替代模型，而該替代模型具備與模擬模型相同的計算能力，此即為采用機器學習方法的主要思路。建立了替代模型以后，優化模型對模擬模型的調用即可以更改為對替代模型的調用，由于替代模型結構簡單、求解容易，從而解決了計算負荷大、求解難度高等問題。替代模型的本質是一種機器學習方法，常用的建模方法包含多項式回歸法、人工神經網絡法及克里格法等[13]。

綜上，模擬—優化耦合技術的運用不但能夠考慮復雜的水合物分解行為（包含多相流體滲流行為、地層變形行為等），對于給定的約束條件（如海底面沉降量、開采時間、甲烷累計產氣量等），還能夠反演求得既定目標下的最優開采方案。因此，在考慮影響地層力學穩定性的限制條件下，采用模擬—優化耦合技術可以計算得到最優累計產氣量，主要步驟分為以下3步：①建立水合物開采模擬模型并求解；②對水合物開采模擬模型進行機器學習，即構建替代模型；③建立優化模型并求解。

2 水合物開采地層力學響應模擬模型

2.1 概念模型與網格剖分

南海北海神狐海域是我國海域水合物勘探開發的重點靶區[14]。該區域在新生代時期地層沉積速率較大，最大沉積厚度達到11 km，油氣資源豐富[15]。多期次構造運動使該區域具有泥底辟活動帶，在其上覆地層中形成高角度斷裂和垂向裂隙，為深部流體提供了流動通道，進而形成了特有的水合物成藏系統[16-18]。由GMGS3鉆探航次調查結果可知，W11站位海底面水深為1 312 m，水合物儲層最大厚度達80 m，平均水合物飽和度（SH）為30%，甲烷氣體體積分數超過99%。水合物儲層孔隙度為50%，滲透率為10 mD，地溫梯度為54.9 ℃/km。水合物儲層下伏地層無明顯的游離氣層特征，為Ⅱ類水合物藏[19]。

為求得W11站位水合物開采過程中的累計產氣量及地層力學響應特征，針對該站位建立水平井降壓開采概念模型。含水合物沉積層為水平延伸地層，由于傳熱—流動—力學耦合數學模型的求解對計算機內存消耗大，且在降壓開采過程中沿整個水平井筒的甲烷氣體壓力損失較小。因此考慮沿水平井段的網格尺寸為1 m。設置水合物儲層厚度為80 m，其上覆弱透水地層厚度為200 m，其下伏飽水地層厚度為20 m（圖1），模型頂板即為海底面。模擬區域在x方向上長度為500 m，網格尺寸隨著與生產井距離的增加而增加，分別為 1 m、2 m、3 m、4 m 和 5 m ；z方向上高度為300 m，針對水合物層劃分的網格尺寸為2 m，下伏飽水地層的網格尺寸也為2 m，上覆地層臨近水合物層的20 m范圍內網格尺寸為2 m，其余區域則為4 m。

圖1 W11站位天然氣水合物水平井開采傳熱—流動—力學耦合概念模型示意圖

2.2 初始條件與邊界條件

初始溫度場由海底面溫度（4.82 ℃）按地溫梯度（5.49 ℃/100 m）逐漸增加[19]，水合物儲層底部初始溫度為15.373 ℃、初始壓力為16.195 MPa，此時水合物儲層中固態水合物和液態水處于相平衡狀態，其中水合物飽和度為30%、含水飽和度為70%。上覆、下伏地層初始含水飽和度為100%，鹽度為3%。飽和水及不含水的含水合物沉積層導熱系數分別為3.1 W/（m·K）、1.0 W/（m·K）。

由于缺乏研究區相關站位的原位地應力數據，模型中初始有效地應力場由儲層骨架顆粒密度（2 600 kg/m3）和深度計算獲取，同時假設在x、y、z三個方向上初始地應力相等[2]。

將模型頂底板（其中模型頂板即為海底面）設定為可發生流體運移和熱量交換的定壓定溫邊界，其溫度和壓力值均為未開采前的初始值?？紤]模型的對稱性，側向邊界設定為隔水隔熱邊界。水平井為模型內邊界，內邊界壓力值即為水平井開采壓力值。

含水合物沉積層的相對滲透率采用Stone模型，如式（1）、（2）所示；毛細管壓力的計算采用van Genuchten模型，如式（3）所示。地質力學特征參數如表1所示[20-21]。

液相相對滲透率（KrA）計算式為：

氣相相對滲透率（KrG）計算式為：

式中SA表示液相飽和度；SirA表示殘余水飽和度，取0.60；nA表示液相衰減指數，取4.5；SG表示氣相飽和度；SirG表示殘余氣飽和度，取0.02；nG表示氣相衰減指數，取3.5。

毛細管壓力（pc）計算式為：

式中pc表示毛細管壓力，最大毛細管壓力為2.0×106Pa；p0表示毛細管進氣壓力，取 1.0×104Pa；m表示van Genuchten指數，取值為0.45。

表1 含水合物沉積層地質力學特征參數表

2.3 模擬模型的樣本空間

確定模擬模型樣本空間的目的是為以建立替代模型為目的的機器學習過程提供訓練樣本，樣本空間的覆蓋性越強、樣本數量越多，訓練效果則越佳，從而使替代模型具有與模擬模型相似的輸入輸出運算關系。

為確定不同水合物開采方案下甲烷累計產氣量與地層力學穩定性的平衡關系，以甲烷累計產氣量最佳為優化目標，以海底面沉降量小于某限值為約束條件，以降壓幅度、開采時間、布井位置（水平井段在垂向上的位置）、水平井段長度為自變量，構建樣本空間。依據水合物試采工程經驗來確定自變量取值（表2）[22]，得到樣本空間樣本總數為750組（自變量樣本數量的積）。

2.4 模擬模型求解

水合物傳熱—流動—力學耦合數學模型采用自主開發的數值模擬程序HYDRATEBiot求解，該程序是基于Biot固結理論建立的TOUGH+HYDRATE軟件框架的擴展程序。HYDRATEBiot既可以求解水合物分解引起的多相多組分系統中復雜的相態變化、流體輸運、熱量傳輸等過程，還可以求解考慮了孔隙水壓力和巖土體變形相互作用的含水合物沉積層力學特征響應過程[23]。此次采用該程序求解W11站位含水合物沉積層水平井降壓開采時的甲烷累計產氣量及地層力學響應特征。

表2 樣本空間自變量參數取值范圍統計表

2.5 模擬結果分析

2.5.1 降壓幅度和布井位置對單位長度水平段甲烷累計產氣量的影響

依據前述樣本空間自變量的參數取值范圍，設定開采時間為180天，考慮降壓幅度和布井位置變化，研究單位長度水平段甲烷累計產氣量和地層力學響應特征。由模擬計算結果可知，水平井降壓幅度為10 MPa，連續開采180天，不同布井位置下單位長度水平段甲烷累計產氣量分別為6.75 m3、7.82 m3和6.28 m3。在設定的開采時間（180天）內單位長度水平段甲烷累計產氣量隨降壓幅度升高而升高，且水平段布置在儲層中下部可以獲得相對較高的累計產氣量（圖2），這受控于下伏地層良好的熱動力學條件與水合物儲層的自封閉效應。一方面，對于厚層Ⅱ類水合物藏而言，水平井段在水合物儲層垂向上的位置越靠近水合物層底界（BSR），其熱動力學條件更有利于水合物發生分解；但另一方面，對于Ⅱ類水合物藏并非將水平井布置在鄰近BSR位置處最佳，在開采過程中還需考慮下伏地層涌水的影響。隨著開采進行，水平井下方水合物儲層中水合物逐漸分解，水平井與下伏透水地層取得良好的水力連通后將引起下伏地層水大量涌入水平井，從而導致水合物儲層內部的壓力降受到抑制，使水合物分解速率降低，不利于水合物的高效開采。因此，在以較大降壓幅度進行水合物開采時，水平段部署在儲層中下部為宜。

圖2 不同布井位置、降壓幅度下單位長度水平井段甲烷累計產氣量柱狀圖

2.5.2 海底面沉降量

水平井在Ⅱ類水合物藏中降壓開采30天時海底面沉降量與開采180天的總沉降量的占比介于15%～20%；開采60天時，不同開采制度下的海底面沉降量與總沉降量的占比介于40%～50%（圖3）。由此可見，水平井降壓開采Ⅱ類水合物藏時地層沉降主要發生在開采早期階段，地層的沉降主要是由于孔隙壓力降低導致儲層骨架壓縮，這是深部地層內部發生的一種體積應變。因此在降壓開采早期階段應采取必要的措施來防止井筒周圍地層產生強烈的變形作用進而導致嚴重的工程問題出現。

海底面沉降量隨降壓幅度增大而增大，以生產井位于水合物儲層中部為例，當降壓幅度依次為4 MPa、6 MPa、8 MPa、10 MPa、12 MPa，相應最大海底面沉降量分別為 0.035 m、0.060 m、0.083 m、0.101 m、0.118 m（圖3）。假定海底面沉降0.1 m為開采水合物保持地層長期穩定的臨界條件，那么針對屬于Ⅱ類水合物藏的W11站位，降壓幅度超過10 MPa的強降壓方案不能保證水合物儲層的長期、穩定開采。進行開采方案設計時，應結合考慮采取適當的措施來避免較大的海底面沉降量產生，同時需做好海底變形的多場監測。

圖3 不同開采時間下海底面沉降量占比及不同降壓幅度下最大海底面沉降量統計圖

3 替代模型的建立與驗證

3.1 替代模型的建立

采用徑向基函數人工神經網絡（Radial Basis Function Artificial Neural Network,簡寫為 RBFANN）方法建立模擬模型的替代模型。徑向基函數是對中心點徑向對稱衰減或增強的非線性高斯函數。RBFANN是一種前饋式神經網絡，包括輸入層、隱含層和輸出層，從輸入層到隱含層為非線性變換，從隱含層到輸出層為線性變換[24]?？傮w來說，RBFANN模型是一個黑箱模型，即不需要考慮復雜的物理過程，僅基于輸入輸出樣本數據集來構建內部映射過程，因此其特點是運算速度快，對于高度非線性問題具有較強擬合能力[24]。由此，RBFANN方法可以用于建立描述水合物開采的多相滲流非線性數學模型的替代模型。

利用前述750組樣本的輸入輸出數據作為替代模型的數據庫，隨機抽取600組為訓練樣本、150組為檢驗樣本，通過MATLAB基礎平臺編寫徑向基函數調用程序，對徑向基函數人工神經網絡進行訓練。在訓練過程中，將降壓幅度、開采時間、水平井段長度、布井位置作為替代模型的輸入變量，再分別將甲烷累計產氣量和海底面沉降量作為替代模型的輸出變量，進而建立起兩個替代模型，其中替代模型2的輸入變量不包括水平井段長度（表3），兩個模型訓練目標絕對誤差均為0.001。

表3 替代模型參數統計表 單位：個

3.2 替代模型的驗證

將基于RBFANN方法建立的替代模型的計算結果與150組檢驗樣本進行對比，如圖4、5所示，替代模型的計算精度較高。

運用R2（確定性系數）、MRE（平均相對誤差）兩個指標對替代模型的計算精度進行檢驗，其計算式如式（4）、（5）所示。R2數值越接近于1，表示替代模型的計算精度越高；MRE數值越接近于0，表示替代模型的計算精度越高。由檢驗結果可知，兩個替代模型的計算精度均較高（表4）。

圖4 替代模型與基于檢驗樣本的模擬模型計算結果對比圖（甲烷累計產氣量為輸出變量）

圖5 替代模型與基于檢驗樣本的模擬模型計算結果對比圖（海底面沉降量為輸出變量）

式中n表示樣本個數；i表示樣本序號；yi表示模擬模型響應值；表示替代模型響應值；表示模擬模型響應值的平均值。

表4 替代模型計算精度統計表

4 混合整數非線性規劃優化模型的建立與求解

4.1 優化模型的建立

4.1.1 目標函數和決策變量

設置甲烷累計產氣量為目標函數，降壓幅度、開采時間、水平井段長度、布井位置為決策變量，則目標函數式為：

式中Q表示甲烷累計產氣量，m3；L表示水平井段長度，m；Z表示與布井位置相關的整數變量（取值0或1）；p表示降壓幅度，MPa；t表示開采時間，d。

4.1.2 約束條件

約束條件包括每個決策變量的約束條件和附加約束條件，共考慮5個約束條件。

將海底面沉降量作為約束條件，其條件式為：

式中g表示海底面沉降量，m；M表示最大允許沉降量，m。

將水平井段長度作為約束條件，其條件式為：

式中L0表示水平井段長度最大值，取值為500 m。

將降壓幅度作為約束條件，其條件式為：

式中p0表示降壓幅度的最大值，取值為12 MPa。

將開采時間作為約束條件，其條件式為：

式中t0表示開采時間最大值，取值為180 d。

將布井位置作為約束條件，共考慮3種布井位置，Z1代表水平井段位于儲層中部，Z2代表水平井段位于儲層中下部，Z3代表水平井段位于儲層下部，上述3個變量取值為1，表示井存在，若取值為0，則表示井不存在，滿足的條件式為：

將式（6）～（11）聯立，則建立起優化模型，該優化模型具有兩個特點：①包含目標函數、決策變量、約束條件的兩個替代模型；②約束條件中含有整數變量。

4.2 優化模型的求解

基于MATLAB平臺編寫遺傳算法求解程序，求解前述混合整數變量的非線性規劃優化模型。求解步驟為：①選用二進制形式進行編碼；②設定遺傳算法參數（表5）；③通過選擇—交叉—變異—再選擇過程的反復迭代，使個體適應度不斷提高，最終得到滿足所有約束條件的最優解。在編碼、選擇、交叉、變異各過程中選用了不同的策略，如采用二進制形式進行編碼，選擇過程采用隨機遍歷抽樣，交叉方法采用單點交叉，采用二進制形式進行變異。基于上述策略，通過提高初始種群數和遺傳代數來保證遺傳算法的全局尋優能力。通過求解優化模型，即可獲得在海底面沉降量約束下的水合物儲層甲烷氣累計產氣量，以及對應的最優開采方案參數（表6）。

表5 遺傳算法參數統計表

4.3 優化模型計算結果與分析

如圖6所示，隨著最大允許沉降量數值增大，甲烷累計產氣量也增大，二者滿足正相關關系；當最大允許沉降量介于0.04～0.05 m時，是開采方案優化結果的轉折點。當最大允許沉降量較小時（小于0.045 m），影響甲烷累計產氣量和海底面沉降量的決策變量是開采時間，海底面沉降量隨開采時間增長而增大；當最大允許沉降量較大時（大于0.045 m），影響甲烷累計產氣量和海底面沉降量的決策變量是降壓幅度，海底面沉降量隨降壓幅度增大而增大。由此可知，因水合物開采引起海底面沉降主要發生在開采初期，如期望獲得較高甲烷氣產量同時控制最大海底面沉降量在0.05 m以下，在開采初期應該減小降壓幅度。

表6 優化模型計算甲烷累計產氣量與最優開采方案參數統計表

圖6 甲烷累計產氣量與海底面沉降量平衡關系曲線圖

以甲烷累計產氣量最高為目標，并且以最大允許沉降量為約束條件，優選水平井段部署在水合物儲層下部，這與模擬模型的計算結果（水平井段部署在水合物儲層中下部）不一致。原因在于模擬模型是在既定方案下，以甲烷累計產氣量最高為單一目標的求解結果，沒有考慮海底面沉降量的約束。另外，與模型設置的開采時間也有關，模擬模型和優化模型的水合物開采時間均為180天，如果生產時間再延長，布井位置的優選結果還會發生變化。

將表6列出的開采方案參數代入到前述模擬模型中進行求解，對比結果可知，優化模型的求解結果與模擬模型的計算結果相對誤差均低于±5%（表7），表明優化模型的計算結果可靠。

5 結論

1）模擬—優化耦合技術的關鍵是機器學習方法的運用，基于徑向基函數人工神經網絡方法而建立的替代模型計算精度較高，可以替代模擬模型來確定輸入輸出變量的關系，從而擺脫既定方案的限制，找到全局最優解。

2）模擬—優化耦合技術可以解決受含水合物沉積層力學響應特征影響的水合物開采方案優選問題。根據試采工程安全要求改變海底面沉降量最大允許值，可以計算得到對應的甲烷累計產氣量，以及降壓幅度、開采時間、井位布置、水平井段長度等最優開采方案參數。

表7 模擬模型和優化模型計算結果對比表

3）隨著最大允許沉降量增大，甲烷累計產氣量增大，二者滿足正相關關系；海底面沉降量隨開采時間增長而增大，也隨降壓幅度增大而增大；水合物開采引起海底面沉降主要發生在開采初期，為了獲得較高甲烷累計產氣量及較小海底面沉降量最大允許值，在開采初期必須減小降壓幅度。