基于半監督高斯混合模型與梯度提升樹的砂巖儲層相控孔隙度預測

魏國華，韓宏偉，劉浩杰，李明軒，袁三一*

(1.中國石化勝利油田分公司物探研究院，山東東營 257000； 2.中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249)

0 引言

孔隙度(φ)作為儲層預測和油藏描述的重要物性參數之一，可以為含油氣性預測、儲層品質評估和儲量計算等提供較為可靠的參考依據[1-3]?？紫抖仁艿讲煌刭|因素的影響，包括構造位置、埋藏深度、成巖程度、沉積環境和巖性變化等。通過鉆井取心并進行巖石物理分析獲取孔隙度是最為直接、準確和有效的方法[4-5]，但由于取樣和測試成本太高而不適合于大規模應用。在利用測井資料預測孔隙度方面，通常采用基于多種巖石物理模型的方法[6-7]。Luo等[8]在Xu-White模型的基礎上提出了改進的巖石物理模型以預測儲層的孔隙度。Li等[9]基于Gassmann方程和多孔巖石物理模型建立了新的三維巖石物理模板以定量表征孔隙度。Wang等[10]基于軟孔隙度模型與Gassmann方程實現了儲層孔隙度的預測。除了采用巖石物理模型外，還存在基于經驗關系[11]、統計關系以及引入貝葉斯理論的孔隙度預測模型[12]，如通過數據驅動的方式，利用大量測井資料建立測井曲線與孔隙度之間統計性、經驗性的巖石物理關系，以實現孔隙度預測[13]；在巖石物理理論的基礎上，根據貝葉斯理論實現孔隙度的概率估計[14-15]。

隨著人工智能和大數據挖掘等技術廣泛應用于石油行業，人工智能算法開始應用于測井孔隙度的預測。該類方法選擇測井數據作為輸入，以專家精細解釋的孔隙度曲線作為標簽，利用人工智能數據驅動類算法建立多種測井屬性與孔隙度之間的非線性關系，在此基礎上再進行孔隙度預測。支持向量機[16-17]、隨機森林(Random Forest，RF)[18-19]、極度梯度提升樹(Extreme Gradient Boosting, XGBoost)[20]、深度神經網絡[21-22]、卷積神經網絡[23]以及長短期記憶網絡[24-26]等數據驅動類方法能建立高維度的多尺度特征與孔隙度之間的復雜聯系，但是不同巖性的測井曲線與孔隙度之間的非線性關系相差較大，僅利用數據驅動的人工智能算法可能對不同巖性的非線性關系擬合效果較差。同時，在人工解釋砂巖儲層孔隙度時，專家通常僅解釋儲層段的孔隙度而無法獲取全井段的孔隙度曲線。此外，高昂的鉆井成本導致智能孔隙度預測方法存在標簽數量不足和測井樣本空間分布不均衡等問題，極大地制約了智能化孔隙度預測精度的提高以及其進一步的推廣。因此，需要在孔隙度預測時考慮巖相控制，通過巖相劃分砂巖與泥巖，再分別進行砂巖段和泥巖段的孔隙度預測。為此，本文提出了一種基于半監督高斯混合模型與梯度提升樹(Gradient Boosting Decision Tree，GBDT)的相控孔隙度預測方法。首先利用少量具巖相標簽的測井數據確定高斯混合模型的初始聚類中心及對應的巖相類別；其次利用大量無標簽測井數據優化高斯混合模型，實現砂巖與泥巖的正確劃分；再次基于地質認識將泥巖孔隙度解釋為固定的極小值，從而后續只開展砂巖孔隙度預測；然后將巖石物理方法導出的孔隙度先驗信息和測井敏感屬性作為梯度提升樹算法的多源輸入信息，通過學習統計性巖石物理關系建立砂巖孔隙度的計算模型；最終根據巖相結果將砂巖段與泥巖段的孔隙度進行組合得到相控孔隙度。實際工區的測井資料測試表明該方法應用效果較好。

1 方法原理

1.1 測井相控孔隙度智能預測流程

如圖1所示，本文提出的基于半監督高斯混合模型與梯度提升樹的相控孔隙度預測方法的主要步驟包括：測井數據預處理、基于半監督高斯混合模型的巖相分類、基于梯度提升樹與巖石物理模型約束的相控孔隙度預測。具體流程為：①對測井數據進行巖性敏感屬性與孔隙度敏感參數的篩選，再對測井巖性敏感屬性進行高通濾波處理，消除低頻趨勢的同時保留高頻細節信息；②利用有巖相標簽的敏感屬性建立巖相分類的初始模型，并結合無巖相標簽的敏感屬性構建半監督高斯混合模型，進行巖相(砂巖與泥巖)的劃分；③根據上一步的巖相結果將測井曲線分成砂巖段和泥巖段兩部分，再分別開展孔隙度預測。

本文的目標儲層為常規的砂巖儲層，因此可將泥巖段孔隙度統一設置為一個極小的常數。針對砂巖段孔隙度預測，首先根據前期優選的孔隙度敏感參數，建立測井敏感參數與孔隙度之間的巖石物理模型；再基于巖石物理模型粗糙估算孔隙度；然后將估算的粗糙孔隙度與其他孔隙度敏感參數(如聲波時差AC、井徑CAL等)作為集成學習算法的多源輸入信息，訓練梯度提升樹模型，預測砂巖段的孔隙度；最終將砂巖孔隙度與泥巖孔隙度進行組合即可得到相控孔隙度。

1.2 基于高通濾波的測井數據預處理

為降低相控孔隙度預測的難度，本文將相控孔隙度預測分為巖相分類解釋和孔隙度預測兩個任務。專家在對測井曲線進行巖相解釋時，通常只關注測井曲線的某一部分并利用該部分測井曲線的差異判斷巖相。這一過程相當于在測井曲線中設置了一個時窗，通過時窗內的曲線變化判斷局部巖性。將該時窗沿著測井曲線進行滑動，可以實現整條測井曲線的巖性解釋。但是，一些常規的機器學習方法(如聚類算法、決策樹算法和支持向量機算法等)無法根據局部的測井曲線差異判斷巖性，而是需要考慮整條測井曲線的差異。盡管圖2中的自然伽馬(GR)和自然電位(SP)參數與巖性具有較好的對應關系，但是將整條GR和SP曲線交會時，卻無法區分巖性(圖3a)。

多種濾波方法已廣泛應用于測井曲線的處理與解釋[27-29]。如圖2所示，GR和SP隨著深度的增加存在低頻變化趨勢，這可能緣于地層的沉積環境的變化[30]。為了使機器學習算法在訓練過程中更好地實現巖性的分類，本文使用高通濾波對低頻趨勢進行處理。由圖2可見，高通濾波處理后的測井曲線消除了低頻趨勢，保留了高頻的局部信息。通過交會圖(圖3)分析，高通濾波處理后砂巖、泥巖的重合部分減少(圖3b)，說明處理后的砂巖和泥巖的GR、SP測井響應差異更大。因此，后續采用處理后的GR和SP參數作為巖相劃分的敏感屬性。

圖2 高通濾波處理前(左)、后(右)的測井曲線對比

圖3 高通濾波處理前(a)、后(b)GR與SP交會分析

1.3 基于半監督高斯混合模型的巖相分類

聚類方法已廣泛應用于測井數據巖性解釋。本文使用高斯混合模型對高斯濾波處理后的GR和SP曲線進行巖性分類。初始模型的選取是無監督聚類算法的關鍵，其質量的好壞不僅影響效率，而且影響精度。因此，本文提出了一種基于半監督高斯混合模型的巖相分類方法。即在傳統高斯混合模型的基礎上，增加具有巖性解釋的測井曲線，計算帶有巖性標簽的初始模型，實現一種半監督的聚類算法。在理論上，通過多個高斯分布線性組合可以擬合出地層巖性的分布[31]。假設測井參數樣本x中每個簇樣本的特征服從多元高斯分布，則x的概率密度函數為

(1)

式中：p(·)表示概率函數；K為樣本中簇的總數；ck為第k個高斯模型(或簇)對應的系數；N(·)表示高斯分布；μk為第k個簇的樣本分布的均值；Σk表示第k個簇的樣本分布的協方差；Dk表示數據的維度；上標T表示轉置。均值μk描述的是不同巖性對應測井參數的均值，而協方差Σk描述的是每一簇樣本的分布形態，反映測井參數的分布情況。

根據式(1)可計算出測井參數每個樣本的概率密度函數，通過最大化所有樣本的概率密度函數的乘積(最大似然估計法)可求得模型的最佳參數。從本質上講，最大化所有樣本的概率密度函數的乘積等同于最大化所有樣本對數概率密度函數的和，即最大化

lgL1(μk,Σk)

(2)

可獲得高斯混合模型中的最佳均值和方差。式中：L1是高斯混合模型的目標函數；J表示樣本總個數。式中有眾多相加項不能通過偏導直接獲得其最大值[32]，因此在高斯混合模型的參數求解中廣泛選用最大期望(EM)算法，通過迭代的思想求解高斯混合模型參數。具體步驟如下。

(1)基于具有專家巖性解釋的測井敏感特征結果初始化K個高斯分布的均值μk與協方差Σk，并隨機初始化混合系數ck。

(2)遍歷所有測井曲線x的所有樣本點，計算第i個樣本點xi(i=1，2，…，J)屬于第k個高斯分布的概率

γi,k=p(xi|zi=k)

(3)

式中：zi表示xi所屬的類；d為xi的維度。

(3)按照下式更新μk′和Σk′，即不斷優化不同巖相的聚類中心位置

(4)

(5)

(4)重復步驟(2)和步驟(3)，直到高斯混合模型參數收斂為止，從而找到不同巖相對應的聚類中心，完成對測井數據的巖相分類。

1.4 基于梯度提升樹的砂巖孔隙度預測

在利用半監督高斯混合模型劃分砂巖和泥巖的基礎上，設置泥巖孔隙度為固定常數0.01，后續采用梯度提升樹算法預測砂巖孔隙度。梯度提升樹算法和隨機森林算法都是基于決策樹的集成算法。在隨機森林算法的基礎上，梯度提升樹算法“提升”了其內部使用的多個決策樹之間的相關性，并使用梯度算法快速求解目標函數。

(6)

式中：FM(xW)為預測的孔隙度,xW為輸入的砂巖段測井數據；M為決策樹總個數；T(xW;θm)表示第m個決策樹，其中θm為第m個決策樹的參數。模型采用向前分步算法，具體步驟如下。

(1)首先確定初始提升樹

F0(xW)=0

(7)

(2)采用向前分步算法，得到第m步的模型，即

Fm(xW)=Fm-1(xW)+T(xW;θm)

(8)

式中：Fm-1(xW)、Fm(xW)分別表示第m-1步、第m步建立的梯度提升樹模型。

(3)利用經驗風險最小化確定下一棵決策樹的參數

(9)

L2[y,Fm(xW)]=[y-Fm(xW)]2

(10)

結合式(8)，可得

L2[y,Fm(xW)]=[y-Fm-1(xW)-T(xW;θm)]2

(11)

由于決策回歸器模型的自身特點，即使數據中的輸入與輸出之間存在著復雜關系，多個回歸器的線性組合也可以很好地擬合孔隙度敏感參數與孔隙度之間的非線性關系。由待測試數據輸入建立的梯度提升樹孔隙度預測模型后，即可獲得砂巖段的孔隙度結果。最終得到的相控孔隙度為砂巖孔隙度與泥巖孔隙度的組合。

2 實例應用

將本文相控孔隙度預測方法應用于中國東部D油田的18口井(分別命名為W1～W18井)。測井段深度為1500～2100 m，巖性主要為砂巖和泥巖。測井曲線主要包括AC、GR、CAL和SP。18口井已做了巖性解釋和孔隙度解釋。選擇其中的W10井進行展示(4條常規測井曲線和專家解釋的孔隙度曲線、地層巖性)，如圖4所示。選擇W1井的測井巖性敏感屬性和標簽計算半監督聚類算法的初始模型。選擇16口井(W1～W16井)的測井數據訓練相控孔隙度預測模型(即訓練集)，其余2口井(W17、W18井)測試相控孔隙度預測模型的性能(即測試集)。

圖4 W10井的測井曲線及巖性解釋結果

2.1 數據預處理與敏感屬性篩選

在人工解釋時，專家通常僅解釋儲層段(即砂巖)的孔隙度，因此需先劃分巖性，再根據巖性的識別結果分段預測孔隙度。如圖5a所示，對16口井的測井曲線進行交會分析。根據不同測井曲線對砂巖與泥巖的區分程度，本文選擇了GR、SP作為測井巖性敏感屬性。

根據訓練集中的巖性解釋結果提取儲層段(即砂巖段)的測井曲線。圖5b顯示了16口井儲層段測井曲線的分布與相關性?；谇€相關性結果，選擇砂巖段的AC和SP作為砂巖孔隙度敏感屬性，它們與孔隙度的相關系數分別為0.95和0.42。

在訓練和測試過程中，對測井曲線進行標準化，其定義為

(12)

式中：X、X*分別為標準化前、后的測井曲線；Xmean為測井曲線的平均值；Xstd為測井曲線的標準差。本文將測井敏感屬性、砂巖孔隙度敏感屬性以及砂巖孔隙度進行歸一化處理。

2.2 基于半監督高斯混合模型的巖性識別

如圖5a所示，雖然巖性敏感屬性可以較好地區分砂巖與泥巖，但是兩者之間存在重疊部分，需要對歸一化后的巖性敏感屬性進行高通濾波處理，從而提高敏感屬性的巖性區分效果。

圖5 W1～W16井的測井曲線分布與相關性分析

選擇W1井的巖性敏感屬性與巖性解釋結果計算得到帶有標簽的初始模型，將初始模型與其余測井的巖性敏感屬性(W2～W16井)輸入半監督高斯混合模型中進行迭代，并使用W17、W18井的巖性敏感屬性(濾波處理后)測試模型性能。同時，使用W17井未經濾波處理的巖性敏感屬性進行半監督高斯混合模型的巖性識別。然后，將其與常規的機器學習模型(K均值算法、支持向量機、隨機森林和K近鄰算法)進行對比 (圖6)。

預測結果的準確率如表1所示。巖性敏感屬性通過高通濾波處理后，巖性識別的準確率有了明顯提升，解決了1500～1800 m巖性識別效果較差的問題(圖6)。同時，本文采用的半監督高斯混合模型對巖性的識別準確率均達到了94%以上，高于常規的機器學習模型。

圖6 W17井(a)、W18井(b)濾波處理前(左)、后(右)的不同方法巖性預測結果對比

表1 W17、W18井不同方法預測結果的準確率對比

針對砂巖儲層，將訓練集16口井的砂巖段聲波時差與專家解釋的真實孔隙度進行交會分析(圖7，黑線表示最小二乘法擬合得到的孔隙度計算模型)，可見兩者之間具有較好的線性關系。通過最小二乘法建立的孔隙度預測模型為

圖7 16口井的砂巖段孔隙度與聲波時差的交會分析

φ=AC×0.16-29

(13)

根據該模型，利用訓練集16口井砂巖段的AC可以初步估算砂巖孔隙度。將最小二乘法估算的孔隙度與砂巖孔隙度敏感屬性作為多源信息輸入梯度提升樹，并以砂巖段的真實孔隙度作為標簽，建立砂巖孔隙度精確預測模型。

應用相關系數作為預測結果的評價標準，即

(14)

本文方法預測的W17和W18井的砂巖孔隙度與實際結果(圖8)的相關系數均為0.97；而最小二乘法使用式(13)在2口井上估計的砂巖孔隙度與實際結果的相關系數分別為0.95、0.94。由于本文所使用的測井數據中聲波時差與砂巖孔隙度的相關系數高達0.95，因此使用最小二乘法和梯度提升樹算法預測砂巖孔隙度精度相當。本文將最小二乘法估計的孔隙度、聲波時差和井徑三種參數(即多源信息)共同作為梯度提升樹算法的輸入，因此與最小二乘法相比，梯度提升樹算法提高了砂巖段孔隙度預測準確率。

圖8 W17井(a)、W18井(b)本文方法砂巖孔隙度預測結果與實際結果對比

將砂巖的孔隙度預測結果與泥巖的孔隙度進行組合，可得W17、W18井整段的相控孔隙度預測結果(圖9)。W17、W18井本文方法預測的孔隙度與實際結果較吻合，兩者相關系數分別為0.77、0.84。雖然砂巖段的孔隙度預測結果精度較高，但是由于巖相控制過程中砂、泥巖劃分中存在誤差，因此整段相控孔隙度的預測結果與實際結果的相關系數略有下降。

在不考慮巖性控制的情況下，利用訓練數據直接使用梯度提升樹算法進行模型訓練，同樣利用W17、W18井2口井進行測試(圖9)，預測結果與實際結果的相關系數分別為0.80、0.30。由于W18井的預測結果在1500～1700 m的范圍內無法對巖性進行有效劃分，導致孔隙度預測結果更偏向于砂巖，因此導致整段預測結果與實際結果的相關系數大大降低。因此，在不考慮巖性控制的情況下，直接利用梯度提升樹模型進行擬合會導致預測精度與泛化能力下降。這可能需要深度學習構建更加深層、更加復雜的結構改善這一情況。同時，該結果說明本文提出的基于半監督高斯混合模型與梯度提升樹的相控孔隙度預測方法具有更高的準確率和更強的泛化能力。

圖9 兩口井不同方法的預測孔隙度與真實孔隙度的對比

3 結論

本文提出了一種基于半監督高斯混合模型與梯度提升樹的測井相控孔隙度預測方法，得出以下結論。

(1)測井曲線的低頻趨勢干擾地層巖相的正確劃分，經過高通濾波處理消除低頻趨勢后，不同巖相的測井響應差異更為明顯。因此，高通濾波預處理可為半監督高斯混合模型提供更好的分類特征，提升巖相分類的準確率。

(2)相比于無監督巖相分類算法和有監督巖相分類算法，半監督高斯混合模型能最大化利用有標簽數據和無標簽數據，引入待測試數據，優化簇類中心位置，從而獲得更好的巖相識別效果。

(3)不同巖相孔隙度的差異將增加智能孔隙度預測模型的建模難度。在考慮相控的情況下，孔隙度預測可簡化為巖性分類和砂巖孔隙度預測兩個子問題，緩解了孔隙度與測井敏感屬性之間非線性關系的復雜度。在巖相識別準確的前提下，相控方法預測的孔隙度與真實結果具有更高的吻合度。