孫致學(xué)，姜寶勝，肖 康，李吉康

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

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

圖1 中非乍得某基巖潛山巖心裂縫照片F(xiàn)ig.1 Photos of fractured core samples from bedrock buried hill in Chad，Central Africa

1 樣本集構(gòu)建

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

1.1 學(xué)習(xí)樣本特征

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

表1 測井?dāng)?shù)據(jù)統(tǒng)計(jì)結(jié)果Table1 Statistics of well logging parameters

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

1.2 Z-score標(biāo)準(zhǔn)化處理

在進(jìn)行裂縫預(yù)測前，需將學(xué)習(xí)樣本進(jìn)行Z-score標(biāo)準(zhǔn)化處理，即將其轉(zhuǎn)換為均值為0、方差為1 的分布，其表達(dá)式為：

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

1.3 K均值聚類算法進(jìn)行樣本去噪

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

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

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

④P(xj)最大時(shí)對應(yīng)的向量就是新的簇中心，若新的簇中心改變則重復(fù)步驟②—③直到目標(biāo)函數(shù)收斂，聚類結(jié)束。

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

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

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

圖3 誤差平方和、計(jì)算時(shí)間與聚類數(shù)的關(guān)系Fig.3 Relationship among sum of squared errors，calculation time and cluster number

2 預(yù)測算法建立

作為機(jī)器學(xué)習(xí)的最新技術(shù)，集成學(xué)習(xí)在智能計(jì)算和機(jī)器學(xué)習(xí)領(lǐng)域引起了廣泛關(guān)注。集成學(xué)習(xí)不是一種特定的模型而是一種思想，通過結(jié)合較簡單的基礎(chǔ)模型來構(gòu)建強(qiáng)化模型。本文引入集成學(xué)習(xí)技術(shù)，將2種不同的基礎(chǔ)模型結(jié)合起來，生成一個更好的模型來預(yù)測裂縫開度。

2.1 支持向量回歸算法

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

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

引入拉格朗日算子ui≥0,≥0,ai≥0,≥0，其拉格朗日函數(shù)表達(dá)式為：

其中：

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

支持向量回歸算法函數(shù)表達(dá)式為：

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

2.2 XGBoost回歸算法

XGBoost 是由一系列回歸樹組成的強(qiáng)大的預(yù)測模型。其核心思想是不斷添加回歸樹，通過生成新樹來擬合前一棵樹的殘差。當(dāng)訓(xùn)練完成得到N棵回歸樹時(shí)，將每棵樹對應(yīng)的分?jǐn)?shù)加起來就是該樣本的預(yù)測值［15］，其表達(dá)式為：

XGBoost目標(biāo)函數(shù)為：

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

將其進(jìn)行二階泰勒展開為：

則（15）式可以改寫為：

2.3 基于嶺回歸的集成學(xué)習(xí)算法

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

圖4 基于嶺回歸的集成學(xué)習(xí)算法Fig.4 Ensemble learning algorithm based on ridge regression

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

在此基礎(chǔ)上，經(jīng)K均值聚類算法去噪后，應(yīng)用基于嶺回歸的集成學(xué)習(xí)算法進(jìn)行裂縫開度預(yù)測，筆者將其定義為新型集成學(xué)習(xí)算法。

3 預(yù)測算法應(yīng)用

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

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

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

3.2 模型應(yīng)用評價(jià)分析

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

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

圖6 預(yù)測裂縫開度與真實(shí)裂縫開度交會圖Fig.6 Cross plot of measured and predicted apertures

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

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

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

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

4 結(jié)論

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

符號解釋

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