孫致學，姜寶勝，肖 康，李吉康

（1.中國石油大學（華東）石油工程學院，山東青島 266580；2.非常規(guī)油氣開發(fā)教育部重點實驗室中國石油大學（華東），山東青島 266580；3.中國石油勘探開發(fā)研究院非洲研究所，北京 100083）

近年來，全球已在30余個盆地發(fā)現了基巖油氣資源，基巖油氣藏成為重要的勘探開發(fā)陣地。相對于常規(guī)沉積巖油氣儲層，基巖儲層中天然裂縫類型、產狀、特征參數的精準評價更加重要。由于該類油氣藏基巖致密、孔滲性極低，天然裂縫系統不僅控制了有效儲層發(fā)育程度和油氣儲量規(guī)模，同時也是油氣開采過程中的重要運移通道。裂縫開度是評價基巖油氣藏儲層質量的重要參數。相對于裂縫密度，裂縫開度對儲層有效滲透率的貢獻更為顯著，也是影響其產能的主控因素之一［1-3］。目前儲層天然裂縫開度預測方法主要分為2 大類：一類是直接觀察法，包括巖心觀測、露頭識別、電鏡觀測等；另一類是間接觀察法，包括成像測井、數值模擬、經驗公式、動態(tài)資料分析等。其中露頭識別法是獲得裂縫開度最直接的途徑［4］，但地表風化作用使裂縫充填特征發(fā)生顯著變化，影響測量結果，同時露頭區(qū)可能遭受后期的改造或掩埋，使得典型裂縫發(fā)育的露頭不容易獲得。巖心中包含著最為直觀、詳實的裂縫信息，但取心資料往往少且不連續(xù)，同時機械應力對巖心的破壞影響天然裂縫開度的測量。隨著斷層掃描機、陰極射線發(fā)光、核磁共振、三維激光掃描技術的發(fā)展，裂縫開度表征朝更微觀、更立體和更精細的方向發(fā)展［5-6］。但由于儀器探測能力的限制，無法系統、大規(guī)模表征天然裂縫開度。成像測井具有高分辨率且連續(xù)測量的特點，能夠直觀地反映裂縫信息，但由于測量成本高，導致獲得的數據非常有限。通過室內數值模擬裝置進行裂縫開度模擬分析，測量相對誤差較小，但裝備適用范圍有限，實驗參數難以獲得，無法真實還原地層條件。進行應力場有限元數值模擬需考慮地質體的巖石物理特征，所需參數較多，難以準確獲取。該方法可以在一定程度上預測裂縫分布，但由于模型并不能反映實際地層情況，導致誤差較大［7-8］。依據滲流力學原理，利用泥漿漏失數據建立裂縫泥漿漏失數學模型，根據鉆井資料進行裂縫開度計算也是當下研究的熱點。但該方法受漏失數據的限制，適用范圍有限。對大多數油田而言，現有資料中除了少量取心資料外，其余幾乎是常規(guī)測井資料，因此如何利用常規(guī)測井信息建立裂縫的測井響應機理模型，進而計算天然裂縫開度，是不得不面對的實際問題。目前應用測井數據解釋天然裂縫仍停留在定性分析水平上［9-10］。主要是由于傳感器捕獲的實時測井數據具有高維、非線性和高噪性的特點，難以建立與裂縫開度之間的量化關系。機器學習對于解決非線性問題具有先天優(yōu)勢，而集成學習是機器學習領域的研究熱點［11］，相較于單一機器學習算法，集成學習算法具有更高的精度和更顯著的泛化性能。

圖1 中非乍得某基巖潛山巖心裂縫照片Fig.1 Photos of fractured core samples from bedrock buried hill in Chad，Central Africa

1 樣本集構建

研究區(qū)位于B 盆地中非乍得西南部、中非剪切帶中段北側，大量巖心和井壁取心分析資料揭示B盆地基巖潛山巖性分為變質巖和巖漿巖2 大類13個亞類30多種巖石類型，主要由花崗巖、正長巖、閃長巖和二長巖等巖漿巖及混合花崗巖和片麻巖類等正變質巖構成［12］。通過對該潛山15口取心井（含新完鉆井）巖心裂縫形態(tài)、規(guī)模及典型特征進行觀察及描述，發(fā)現研究區(qū)油氣運移通道包括構造裂縫、網狀縫、張剪縫及沿縫溶蝕孔洞，以張剪縫為主（圖1）。

1.1 學習樣本特征

準確且全面的裂縫數據樣本是實現模型訓練的基礎。所用數據集包括測井數據（表1）及相應裂縫開度值，共2 140 組，由基巖潛山油藏巖心描述、關鍵井成像測井、巖礦薄片鑒定等數據組成。由裂縫開度分布（圖2）可知，樣本集中裂縫開度最小值為0.011 mm，最大值為0.544 mm，平均值為0.183 mm，標準差為0.087 mm，裂縫開度主要集中在0.126～0.258 mm。

表1 測井數據統計結果Table1 Statistics of well logging parameters

圖2 裂縫開度分布Fig.2 Fracture aperture distribution

1.2 Z-score標準化處理

在進行裂縫預測前，需將學習樣本進行Z-score標準化處理，即將其轉換為均值為0、方差為1 的分布，其表達式為：

如果一個特征的方差比其余特征的方差大許多個數量級，那么該特征將會主導整個目標函數，使得模型不能從其余特征學習到數據的特征。相對于min-max 歸一化方法，該方法不僅能夠去除量綱，還能夠均衡考慮所有維度的變量。

1.3 K均值聚類算法進行樣本去噪

樣本數據由于測量儀器及人為因素的干擾，不可避免的引入噪聲。為此采取利用K均值聚類算法進行數據過濾的思路，以去除冗余，提高學習樣本質量。

K均值聚類算法是基于距離的聚類算法，將距離作為相似性的評價指標，即對象之間的距離越近，相似度越大。而異常點通常距離中心點較遠，檢測異常點，從而進行樣本過濾［13］。假設輸入的樣本向量集合為：

K均值聚類算法具體步驟包括：①從輸入的樣本向量集合中隨機選取1 個向量作為第1 個簇中心點，簇中心集合記為center。②對于滿足條件的任意向量，計算與最近簇中心的距離。③計算每個向量被選為簇中心的概率，其表達式為：

④P(xj)最大時對應的向量就是新的簇中心，若新的簇中心改變則重復步驟②—③直到目標函數收斂，聚類結束。

K的取值對聚類算法的效果具有極大影響。若K取值過小，將導致數據粗化，在剔除異常點的同時會誤判正常數據，造成有效樣本丟失；若K取值過大，將致使聚類結果無法有效收斂，計算時間過長，導致無法有效篩選異常數據。為此采用手肘法來確定K值。手肘法的核心指標是誤差平方和，其表達式為：

該方法的核心思想是隨著K值的增加，樣本數據劃分更加精細，各個簇的聚合度逐漸提高，誤差平方和逐漸減小。當K值小于真實聚類數時，K值增加會顯著增強各個簇的聚合度，誤差平方和下降幅度變大。當K值接近真實聚類數時，提高K值各個簇的聚合度變小，誤差平方和變化幅度驟減（圖3）。

由圖3 可知，當K取值為5，即當聚類數為5 簇時，K均值聚類算法性能最優(yōu)，過濾異常值能力最強，因此本文聚類數取5。同時計算距離時容易受較大數據的影響而忽略取值較小的數據，需在聚類前進行Z-score 標準化處理。然后通過K均值聚類算法對樣本數據進行去噪，找出異常點72組。將異常點剔除，其余2 068 組樣本數據用于后續(xù)算法的訓練與測試。

圖3 誤差平方和、計算時間與聚類數的關系Fig.3 Relationship among sum of squared errors，calculation time and cluster number

2 預測算法建立

作為機器學習的最新技術，集成學習在智能計算和機器學習領域引起了廣泛關注。集成學習不是一種特定的模型而是一種思想，通過結合較簡單的基礎模型來構建強化模型。本文引入集成學習技術，將2種不同的基礎模型結合起來，生成一個更好的模型來預測裂縫開度。

2.1 支持向量回歸算法

鑒于地質認識及資料豐度的不確定性，以及特征之間具有復雜的非線性關系，應用傳統回歸模型不能較好地進行裂縫開度預測。而支持向量回歸算法可通過核函數將樣本數據映射到高維空間，解決非線性問題，同時該算法具有良好的穩(wěn)定性和泛化能力［14］。支持向量回歸算法可形式化為：

引入松弛變量ξi和，可將（5）式寫為：

引入拉格朗日算子ui≥0,≥0,ai≥0,≥0，其拉格朗日函數表達式為：

其中：

根據wolf 對偶的定義，在KKT 條件下得到拉格朗日對偶形式為：

支持向量回歸算法函數表達式為：

對于非線性問題，可通過非線性變換轉化為某個高維空間中的線性問題，即用核函數k(x,xi)替換可以實現非線性函數擬合，能較好處理非線性以及高維數的問題，可表示為：

2.2 XGBoost回歸算法

XGBoost 是由一系列回歸樹組成的強大的預測模型。其核心思想是不斷添加回歸樹，通過生成新樹來擬合前一棵樹的殘差。當訓練完成得到N棵回歸樹時，將每棵樹對應的分數加起來就是該樣本的預測值［15］，其表達式為：

XGBoost目標函數為：

為避免算法擬合過程中的過擬合，算法不能同時訓練所有回歸樹，因此利用固定訓練好的回歸樹，依次添加一棵新樹來解決，假設步驟t的預測值用表示，（12）式可以寫為：

將其進行二階泰勒展開為：

則（15）式可以改寫為：

2.3 基于嶺回歸的集成學習算法

本文所提的裂縫開度預測集成算法以XGBoost回歸算法和支持向量回歸算法為基礎模型。每個基礎模型均接收輸入數據，并給出獨立的裂縫開度預測結果，這些預測結果均作為元特征，被饋送到元學習器中（本文的元學習器采用嶺回歸算法），并給出最終的裂縫開度預測結果（圖4）。

圖4 基于嶺回歸的集成學習算法Fig.4 Ensemble learning algorithm based on ridge regression

該算法為基礎學習器δg=1,2（g=1 為支持向量回歸，g=2 為XGBoost 回歸）對于H折交叉驗證中的每一個待預測的訓練樣本集合DH都有與之對應的訓練集預測結果集合ZHg。這樣的循環(huán)過程完畢后，對于每個基礎學習器而言，都有H對同質基礎學習器訓練集預測值，將其整合為D2g。再將所有基礎學習器D2g整合作為元特征定義為D2。將D2饋送到嶺回歸算法得到加權結果即為最終預測開度值?；趲X回歸的集成學習算法的最終表達式為：

在此基礎上，經K均值聚類算法去噪后，應用基于嶺回歸的集成學習算法進行裂縫開度預測，筆者將其定義為新型集成學習算法。

3 預測算法應用

3.1 模型參數優(yōu)化求解

機器學習算法參數的選擇直接決定了算法的性能。網格搜索法是當前應用最為廣泛的參數優(yōu)化算法。但該方法依靠窮舉所有參數進行優(yōu)化，計算成本過于龐大，同時對于連續(xù)數據需要等間取樣，不一定能取得全局最優(yōu)。故采用隨機搜索進行參數優(yōu)化，該方法主要原理是從指定的分布中采樣固定數量的參數設置。與網格搜索法相比，該方法在保障準確度的同時，顯著減少計算時間。

根據測試集上模型的均方根誤差值來判斷基礎模型最佳超參數。其中支持向量回歸算法的主要超參數為懲罰系數，XGBoost 回歸算法的主要超參數為最大深度，其超參數隨機優(yōu)化調參過程如圖5 所示。在搜索的過程中，超參數快速收斂，并找出最優(yōu)值。支持向量回歸算法的懲罰系數搜索范圍為0～20，最優(yōu)值為0.147，對應均方根誤差為0.113；XGBoost 回歸算法的最大深度搜索范圍為0～18，最優(yōu)值為13，對應均方根誤差為0.076。

3.2 模型應用評價分析

確定好模型參數之后，隨機選取80%經過Zscore 標準化處理后的樣本數據作為訓練集共1 712組，20%的樣本數據作為測試集共428 組來驗證模型效果。以均方根誤差（RMSE）和真實裂縫開度值與預測裂縫開度值間相關系數（R2）作為評價標準。將測試集分別代入訓練好的支持向量回歸算法、XGBoost 回歸算法及基于嶺回歸的集成學習算法中，計算測試集中真實裂縫開度與預測裂縫開度間相關系數（圖6）。

圖5 隨機搜索優(yōu)化調參過程Fig.5 Parameter adjustment optimization process based on random search

圖6 預測裂縫開度與真實裂縫開度交會圖Fig.6 Cross plot of measured and predicted apertures

由圖6 可以看出，3 種算法中，基于嶺回歸的集成學習算法的R2最高，達0.928。同時為探究K均值聚類降噪效果，將樣本數據饋送于基于嶺回歸的集成學習算法中進行訓練和測試，并與先前計算結果進行綜合對比（表2），發(fā)現4 組方法中K均值-基于嶺回歸的集成學習算法的RMSE 最小，R2最大。即該算法的預測裂縫開度值與真實裂縫開度值之間的偏差最小，支持向量回歸算法的RMSE 最大，R2最小。K均值聚類算法能夠對學習樣本進行有效降噪，去除冗余，提高了學習樣本的質量。

表2 各算法預測效果對比Table2 Comprehensive comparison results of prediction effects in various algorithms

圖7 樣本觀測值Fig.7 Sample observation values

為進一步綜合分析該算法的預測效果，將部分樣本真實值及各算法預測值進行可視化研究。從圖7 可以明顯看出，支持向量回歸算法和XGBoost回歸算法預測值整體在真實值上下波動，支持向量回歸算法總體變化平穩(wěn)，但對裂縫開度突變值檢測不明顯。XGBoost 回歸算法對數據敏感，部分數值波動較大。新型集成學習算法的計算結果緊密圍繞真實裂縫開度值波動，很好地結合了基礎算法的優(yōu)點，平衡了基礎算法的缺點，預測精度明顯提升。

4 結論

利用測井數據及其對應裂縫開度值，提出基于新型集成學習的基巖潛山油藏儲層裂縫開度預測算法。該算法先通過K均值將學習樣本進行聚類、降噪來提升學習樣本質量；以支持向量回歸算法和XGBoost 回歸算法為基礎模型，并利用隨機搜索進行基礎模型參數優(yōu)化。然后利用嶺回歸算法對優(yōu)化好的基礎模型進行集成組合。所提出的新型集成學習算法建立的裂縫開度預測模型彌補了單一回歸算法不穩(wěn)定的特點，提升了預測精度，能夠充分挖掘測井數據中蘊含的地質信息，為裂縫開度定量預測提供了新的思路。同時該方法可實現自動、快速優(yōu)化調參，具有廣泛的適用性。

符號解釋

x'——標準化處理后的樣本數據；x——樣本數據；μ——樣本均值；δ——樣本方差；K——聚類數，簇；X——樣本向量集合；j——樣本序號；n——樣本向量總數；center——簇中心集合；P(xj)——簇中心概率；D(xj)2——最近簇中心的距離；SSE——誤差平方和；i——樣本列號；p——Ci中的樣本點；Ci——第i個簇；mi——Ci中所有樣本的均值；w——法向量；b——位移項；C——懲罰系數；m——樣本數量；lε——不敏感損失函數；ε——軟間隔帶；f(x)——支持向量回歸算法函數；yi——開度實際值；ξi和松弛變量；和——在第i數據下不同的拉格朗日算子；α——拉格朗日算子αi的合集；拉格朗日算子的合集；——在第j數據下不同的拉格朗日算子；xi和xj——輸入的第i和第j個數據；k(x,xi)——核函數開度預測值；N——回歸樹的總數，個；fk——第k棵回歸樹算法；Lt——步驟t下的目標函數；t——步驟序號；——預測開度值與真實開度值的差；Ω(f)——懲罰項；γ——回歸樹分割的難度系數；T——回歸樹葉子節(jié)點個數，個；λ——L2正則系數；w′——回歸樹葉子節(jié)點權重——步驟t的預測值；Const——常數；gi——損失函數一階導數；hi——損失函數二階導數；Ij——第j個葉子的樣本集合；δg=1,2——基礎學習器；H——交叉驗證折數；DH——H折交叉驗證中的每一個待預測的訓練樣本集合；ZHg——與DH對應的訓練集預測結果集合；D2g——H對同質基礎學習器訓練集預測值；D2——元特征；Y——預測的最終開度值；W——元特征數據矩陣；β^——嶺回歸估計值；λL——嶺參數；I——單位矩陣。