馬先林 ，周德勝，蔡文斌，李憲文，何明舫

1.西安石油大學石油工程學院，陜西 西安 710065 2.西部低滲--特低滲油藏開發(fā)與治理教育部工程研究中心，陜西 西安 710065 3.中國石油長慶油田公司油氣工藝研究院，陜西 西安 710018

引言

傳統(tǒng)的多段壓裂水平井產(chǎn)能預測模型，大多采用機理驅(qū)動的方式建立，建模方法包括解析、半解析及數(shù)值模擬方法等。解析法基于均質(zhì)儲層、單相滲流及簡單平板狀雙翼對稱裂縫等假設，以油氣滲流力學為基礎，應用壓力疊加原理和位勢能理論等得出產(chǎn)能計算公式，計算效率高[1-5]；而對于非均質(zhì)較強、多相滲流、復雜裂縫網(wǎng)絡的儲層，應用油氣藏數(shù)值模擬法進行產(chǎn)能預測，在精細油氣藏描述的基礎上，應用連續(xù)介質(zhì)或離散模型描述人工裂縫和天然裂縫，通過模擬油氣在儲層中滲流過程，計算分段壓裂水平井的產(chǎn)能[6-9]。

傳統(tǒng)的產(chǎn)能預測方法對儲層特征和復雜縫網(wǎng)幾何形態(tài)表征不但存在較強的不確定性，而且還引入了一些理想化假設和計算上的簡化，使得產(chǎn)能預測結(jié)果可能偏離實際。另外，在滲流方程求解時，大多數(shù)采用網(wǎng)格加密及迭代計算等手段，導致建模周期長，計算效率低下。事實上，分段壓裂水平井的生產(chǎn)規(guī)律受到地質(zhì)、鉆完井、壓裂施工及生產(chǎn)制度等諸多因素的影響，因此，過多的理想化假設以及求解方法的復雜性限制了機理驅(qū)動產(chǎn)能預測方法在現(xiàn)場的應用與推廣。

隨著人工智能的興起，特別是機器學習技術為多級壓裂水平井產(chǎn)能預測提供了新的途徑。與傳統(tǒng)的機理驅(qū)動方法相比，該方法通過對油氣大數(shù)據(jù)的挖掘，建立壓裂井產(chǎn)能與儲層、鉆完井和壓裂施工等參數(shù)之間的相關模型，能夠迅速地對壓裂效果進行評價，從而進行壓裂優(yōu)化設計?；跀?shù)據(jù)驅(qū)動的建模方法已經(jīng)在致密油氣藏開發(fā)中得到了初步的應用[10-14]，體現(xiàn)了在產(chǎn)能預測方面獨特的優(yōu)勢。但所建立的預測模型大多是“黑盒子”模型，即輸入一組相關參數(shù)值能夠獲得產(chǎn)能的預測值，卻不能解釋模型內(nèi)部的預測機制，從而降低了模型的可信度與實用性。為此，許多學者對提高機器學習模型可解釋性進行了研究，并提出了一些解釋方法[15]。機器學習產(chǎn)能建模是純數(shù)據(jù)驅(qū)動的方法，預測精度依賴于樣本數(shù)據(jù)的個數(shù)和質(zhì)量，小樣本會導致預測失真和泛化能力降低。

本文首先介紹致密氣藏水平井產(chǎn)能預測機器學習建模流程和機器學習原理，其次討論對預測模型進行解釋的SHAP（SHapley Additive exPlanations）方法原理，并以蘇里格氣田東區(qū)分段壓裂水平井為例，對方法的有效性和實用性進行驗證。

1 水平井產(chǎn)能預測機器學習建模流程

利用機器學習建立致密氣藏分段壓裂水平井產(chǎn)能預測模型，包括6 個步驟，如圖1 所示。

圖1 水平井產(chǎn)能機器學習建模流程圖Fig.1 Workflow for modelling horizontal well productivity using machine learning

（1）原始數(shù)據(jù)收集。數(shù)據(jù)集包括主要影響參數(shù)及其壓裂后水平井的產(chǎn)能評價指標，其中，影響參數(shù)包括地質(zhì)、工程等因素，產(chǎn)能可以是無阻流量或產(chǎn)氣量等。

（2）數(shù)據(jù)預處理。先進行數(shù)據(jù)清洗、缺失數(shù)據(jù)插補、數(shù)據(jù)降維及數(shù)據(jù)轉(zhuǎn)換等[16]，再將預處理后的數(shù)據(jù)集劃分為訓練集和測試集，劃分比例一般取70%～30%或80%～20%。

（3）機器學習建模。應用訓練集數(shù)據(jù)建立產(chǎn)能模型，主要尋找機器學習模型中最優(yōu)的超參數(shù)取值，常用的優(yōu)化方法包括網(wǎng)格搜索、多折交叉驗證及自動學習等。

（4）產(chǎn)能預測模型評價。利用測試集數(shù)據(jù)評估產(chǎn)能預測模型精度，常用精度評價指標包括決定系數(shù)（R2）、平均絕對誤差（MAE）及均方誤差（MSE）。根據(jù)評價指標，選擇在測試集上具有最高預測性能的機器學習算法所建立的產(chǎn)能預測模型。

（5）模型解釋?；诮⒌淖顑?yōu)產(chǎn)能預測模型，利用SHAP 方法對產(chǎn)能預測進行全局和局部解釋。

（6）模型應用。應用建立的產(chǎn)能預測模型，評價壓裂效果，優(yōu)化新井的壓裂參數(shù)設計。

2 機器學習原理

2.1 機器學習方法

機器學習是一門多領域交叉技術，其本質(zhì)是一種特殊的算法。通過分析大數(shù)據(jù)，發(fā)現(xiàn)數(shù)據(jù)內(nèi)部潛在的模式并應用這些模式進行預測[17]。假設一個訓練數(shù)據(jù)集D包含n個學習樣本，每個樣本有m個影響參數(shù)和一個輸出參數(shù)，即

選擇相應的機器學習算法在此數(shù)據(jù)集建立一個預測模型?y=f(x,θ)，優(yōu)化算法的超參數(shù)θ，使得以下誤差函數(shù)值最小

式中：f--某一機器學習算法建立的模型；θ--該算法的超參數(shù)；L(yi,f(xi,θ))--誤差函數(shù)。

典型的機器學習算法包括人工神經(jīng)網(wǎng)絡、支持向量機、決策樹及隨機森林及梯度提升樹等。

2.2 機器學習模型的可解釋性

機器學習模型可解釋性為用戶提供了理解機器學習模型一個接口，既是機器學習模型的代理，又是一種解釋模型的方法[18-19]。模型可解釋性分為兩類：事前可解釋性和事后可解釋性。事前可解釋性是指在建模時，采用可解釋性好的模型或設計有解釋能力的模型，使模型自身具有解釋能力。事后可解釋性是對已建立的機器學習模型進行解釋，獨立于機器學習建模過程，靈活性較強。根據(jù)可解釋的范圍，事后可解釋性又分為全局和局部兩種，全局可解釋用于理解模型內(nèi)部的工作原理，對模型整體能力的解釋[20]；局部可解釋用于理解機器學習模型針對單個樣本的預測過程和依據(jù)，對樣本的預測結(jié)果進行解釋[21]。

本文是對建好的水平井產(chǎn)能模型預測結(jié)果進行解釋，屬于事后可解釋性，應用SHAP 可解釋技術。SHAP 法通過計算每一個輸入變量對預測的貢獻值進行模型解釋，是一種加性解釋方法。比如，使用SHAP 法解釋樣本x*的機器模型預測值?y=f(x*)時，預測值f(x*)可以分解成

式中：φ0--預測模型f(x) 在數(shù)據(jù)集上的平均預測值；--第j個輸入變量對樣本x*預測的貢獻值，即SHAP 值；M--輸入變量的數(shù)量。

對于3 個輸入變量的預測問題，式（3）的示意流程見圖2，紅色表示SHAP 值是正的，藍色則表示SHAP 值是負的。圖2 表明，SHAP 值表示某一個樣本的預測值各個輸入變量的貢獻大小[22]，φ1表示變量1 用于預測后引起當前預測值的變化，是正影響，導致預測值增加，同理，φ2為變量2 加入后當前預測值的增加量，但是，當變量3 加入預測后，導致當前預測值減少，減少量為φ3。

圖2 3 個輸入變量SHAP 方法Fig.2 SHAP method with three input variables

3 應用實例

3.1 數(shù)據(jù)收集

收集的原始數(shù)據(jù)包括蘇里格氣田東區(qū)598 口分段壓裂水平井，產(chǎn)能影響參數(shù)包括7 個地質(zhì)和完井參數(shù)以及15 個壓裂施工輸入?yún)?shù)，輸出參數(shù)為無阻流量，具體參數(shù)描述見表1。

表1 產(chǎn)能影響參數(shù)Tab.1 Influential parameters of well productivity

由于大部分水平井沒有進行測井，不能分析孔隙度、滲透率、楊氏模量及泊松比等參數(shù)對產(chǎn)能的影響。

3.2 數(shù)據(jù)預處理

3.2.1 缺失值處理

在數(shù)據(jù)集中，儲層水平段長度的缺失值多達25.00%，破裂壓力有13.00% 的數(shù)值缺失，遠大于5.00% 安全最大閾值插補的要求，因此，這些變量在建模時不予以考慮。其他變量缺失值的分布如圖3所示。由圖3a 可以看出，大部分變量的缺失值低于2.00%。

圖3b 顯示了變量缺失值的分布，其中，最上面的紅色豎條表示入井總液量有0.25%數(shù)值缺少，第二行的紅豎條表示液氮量缺失0.50%，第三行的兩個紅豎條表示液氮量和施工最高壓力同時缺失達1.26%，第四行的紅豎條表示垂深缺失2.02%，最下行表示95.97%數(shù)據(jù)沒有缺失。由于數(shù)據(jù)缺失較少，故決定刪除數(shù)據(jù)表中有缺失值的那些井，剩余水平井數(shù)量為574 口。

圖3 變量缺失值分析Fig.3 Missing data analysis

3.2.2 異常值處理

原始數(shù)據(jù)中若存在不合理的值，會影響模型預測精度。采用基于箱線圖的異常值檢測方法，利用數(shù)據(jù)中的上、下四分位數(shù)（Q3、Q1）及四分位距（IQR），四分位距定義為上四分位數(shù)與下四分位數(shù)的差值，即IQR=Q3-Q1，變量的值大于Q3+1.5IQR和小于Q1-1.5IQR為異常值。在剔除數(shù)據(jù)集中的異常值后，用于機器學習建模的水平井為532 口。

3.2.3 描述性數(shù)據(jù)分析

描述性數(shù)據(jù)分析采用圖表等方式對數(shù)據(jù)進行統(tǒng)計性描述，是建模的重要內(nèi)容。

圖4 顯示了2012--2020 年平均無阻流量和主要壓裂施工參數(shù)的變化趨勢，隨著平均最大排量、加砂量和壓裂段數(shù)的增加，水平井分段壓裂后的增產(chǎn)效果明顯。

圖4 平均無阻流量和壓裂施工參數(shù)的變化趨勢Fig.4 Trend analysis of average absolute open flow potential and hydraulic fracturing treatment

圖5 給出了主要壓裂參數(shù)與無阻流量的箱線圖，從圖中可以看出，無阻流量與水平段長度、壓裂段數(shù)、加砂量、排量和入井總液量均存在正相關關系，采用水力橋塞分段壓裂工藝后顯著地提高了水平井的產(chǎn)能。

圖5 壓裂參數(shù)影響無阻流量的箱線圖Fig.5 Boxplots for effect of hydraulic fracturing treatment parameters on AFOP

3.2.4 相關性分析

若同時使用相關性較強的輸入?yún)?shù)，建模不僅會增加模型訓練的時間，而且影響模型可解釋性。

圖6 給出了變量間的Pearson 相關系數(shù)，與無阻流量正相關最強的參數(shù)是加砂量（R=0.50），其次是有效儲層段長度、壓裂段數(shù)及入井總液量等參數(shù)。完井井深與水平段長具有極強的線性正相關性（R=0.92），表明完井較深的水平井具有較長的水平段長度；水平段長同有效儲層水平段長、壓裂段數(shù)、入井總液量也具有較強的線性正相關；砂量與入井總液量也具有很強的線性相關性（R=0.91）。

圖6 Pearson 相關系數(shù)矩陣Fig.6 Pearson correlation matrix

為減少變量間的線性相關性，建立以下輸入?yún)?shù)：（1）段間距=水平段長度/壓裂段數(shù)，不再使用壓裂段數(shù)；（2）加砂強度=砂量/水平段長度，代替加砂量；（3）用液強度=入井總液量/水平段長度，代替入井總液量。

采用逐步向后回歸，對輸入?yún)?shù)進行了選擇，確定建模的輸入?yún)?shù)為水平井水平段長、有效儲層水平段長、完井井深、段間距、加砂強度、用液強度、改造工藝及最大排量。

3.3 產(chǎn)能預測模型的構(gòu)建

經(jīng)過預處理后的數(shù)據(jù)集按80%～20% 的比例隨機劃分成訓練集與測試集，采用人工神經(jīng)網(wǎng)絡（ANN）、支撐向量機（SVM）、隨機森林（RF）和梯度提升樹（GBDT）算法在訓練集學習，構(gòu)建無阻流量預測模型。

3.3.1 模型超參數(shù)優(yōu)化

采用網(wǎng)格搜索和十折交叉驗證方法，在訓練集上對4 種機器學習算法的超參數(shù)進行優(yōu)化。以GBDT 算法為例，GBDT 中的主要超參數(shù)有決策樹的棵數(shù)、決策樹的深度和決策樹的深度。首先確定這些超參數(shù)的范圍，再應用十折交叉驗證方法進行驗證，超參數(shù)決策樹的棵數(shù)、決策樹的深度及決策樹的深度取值分別為43、8 和2。GBDT 模型預測結(jié)果如圖7 所示。

由圖7 可以產(chǎn)出，無阻流量小于100×104m3/d時，模型預測值比較接近實際值。然而，當無阻流量大于100×104m3/d 時，模型預測精度不理想，主要原因是水力橋塞分段壓裂水平井只占數(shù)據(jù)集的19%，而在這些井中無阻流量大于100×104m3/d 的達57%，增加水力橋塞分段壓裂水平井數(shù)據(jù)有助于改進模型預測精度。

圖7 GBDT 模型預測結(jié)果Fig.7 Prediction results of GBDT model

3.3.2 模型評價

由于目前沒有一個機器算法能夠絕對優(yōu)于其他算法，主要通過評價機器學習模型在測試數(shù)據(jù)集的預測性能，選擇最適合的算法，為此，應用建立的機器學習模型對測試集中的106 口水平井的無阻流量進行了預測，表2 是這些模型的預測性能對比。

平均絕對誤差和均方誤差越小，模型預測值與實際值越接近；決定系數(shù)越接近1，模型預測值與實際值越接近。由表2 可以看出，GBDT 算法在測試集上的預測結(jié)果優(yōu)于其他幾種學習算法，故對用GBDT 建立的預測模型進行解釋。

表2 機器學習模型性能對比Tab.2 Comparison of performances of machine learning model

3.4 產(chǎn)能預測模型解釋

3.4.1 模型全局解釋

圖8 給出了GBDT 模型的解釋，壓裂施工參數(shù)中，改造工藝和最大排量對無阻流量的影響最顯著，2018 年以前，蘇里格氣田采用裸眼封隔器和水力噴砂分段壓裂改造工藝，排量大部分集中在2～4 m3/min，施工排量低，而且施工段數(shù)受限。2018年以后，改用水力橋塞分段壓裂改造方式，排量提高到6～10 m3/min，施工段數(shù)也有所增加，無阻流量大于100×104m3/d 的水平井明顯多于其他兩種改造工藝。有效儲層段長度和完井井深對無阻流量也有較大的影響，井越深，儲層壓力越大。

圖8 影響因素重要性全局分析Fig.8 Global interpretation of input variables

圖8b 展示了每個影響因素SHAP 值的分布，圖中，每個數(shù)據(jù)點代表一口壓裂水平井，顏色表示變量的值，從藍色到紅色表示變量數(shù)值由低到高的變化。

正（負）的SHAP 值表示影響參數(shù)與無阻流量呈正（負）相關，例如，增大排量，排量的SHAP 值會變大，導致無阻流量增大；而無阻流量則隨著段間距增加而減少。

3.4.2 模型局部解釋

圖9 給出了一口水平井的局部解釋。該井的無阻流量是所有輸入?yún)?shù)貢獻的總和，基準值為GBDT 模型預測的無阻流量平均值，等于50.88×104m3/d。紅色導致該平均值增加，而藍色表示使其下降。有效儲層段長、完井井深和用液強度的SHAP 值分別為9.891，5.294 和1.578；最大排量、改造工藝和加砂強度的SHAP 值分別是-6.898，-5.582 和-1.745，段間距和水平段長的影響也是負的，它們的SHAP 值太小，在圖中無法顯示。最終，該壓裂水平井的無阻流量預測值為53.03×104m3/d，接近實際值59.93×104m3/d。

圖9 某水平井無阻流量預測結(jié)果解釋Fig.9 Local interpretation of a fractured horizontal well

3.4.3 變量相關解釋

SHAP 方法還能揭示輸入變量之間相關性對無阻流量的影響。

如圖10a 所示，當采用水力橋塞分段壓裂工藝時，改造工藝的SHAP 值迅速增加，說明水力橋塞分段壓裂工藝與無阻流量呈正相關；而且隨著壓裂水平段增長，進一步增加水力橋塞分段壓裂工藝的SHAP 值，因此，采用長壓裂水平段的橋塞分段壓裂能夠增加無阻流量。

圖10 變量相關性解釋Fig.10 Interpretation of variable correlation

圖10b 表明，當段間距小于100 m 時，段間距的SHAP 值迅速增加，導致無阻流量增大，而且提高加砂強度有助于段間距的SHAP 值增大，進一步增加無阻流量。但是，隨著段間距的增大，卻引起其SHAP 值減少，因此，段間距與無阻流量呈負相關。

4 結(jié)論

（1）提出了一種基于機器學習算法進行致密氣藏分段壓裂水平井產(chǎn)能預測以及模型解釋的方法，該方法是對現(xiàn)有機理驅(qū)動產(chǎn)能預測方法的有效補充，具有多類型數(shù)據(jù)綜合和較強的預測能力，能夠有效地提高分段壓裂水平井產(chǎn)能預測效率和精度，實用性強。

（2）學習模型全局可解釋性有助于理解梯度提升樹算法產(chǎn)能預測機制，明確影響蘇里格氣田東區(qū)分段壓裂水平井無阻流量主要壓裂施工參數(shù)是排量和改造工藝。

（3）SHAP 值不僅有助于分析單一壓裂水平井產(chǎn)能的影響因素，還能理解輸入變量之間的交互作用對產(chǎn)能的影響。

（4）首次應用可解釋機器學習算法對蘇里格氣田東區(qū)分段壓裂水平井產(chǎn)能評價，可為該區(qū)的壓裂優(yōu)化設計和開發(fā)方案制定提供依據(jù)。