孫予舒,黃 蕓,梁 婷,季漢成,向鵬飛,徐新蓉
(1.中國石油大學(xué)(北京)地球科學(xué)學(xué)院,北京 102249;2.中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室,北京 102249;3.中國石油華北油田分公司勘探開發(fā)研究院,河北任丘 062550)
巖性識別是儲層評價中重要的研究內(nèi)容之一。目前識別巖性最可靠的方法為鉆井取心,然而,受取心成本較高以及取心率不高等問題的限制,常常不能做到全井段取心。因此,在實際勘探過程中,利用測井資料間接對目的層進行巖性解釋成為了研究儲層的重要手段。
目前利用測井?dāng)?shù)據(jù)進行巖性解釋的方法主要分為傳統(tǒng)的測井巖性解釋方法,以及基于機器學(xué)習(xí)的智能化方法[1]。前者包括:基于測井曲線響應(yīng)特征的定性解釋方法、基于測井響應(yīng)方程的定量解釋方法、圖版法;后者則包括:支持向量機、神經(jīng)網(wǎng)絡(luò)和分類決策樹等方法。對于碳酸鹽巖來講,由于其礦物成分較為單一,測井巖電響應(yīng)特征不明顯[2]。同時,由于各測井曲線間存在著大量的信息冗余,導(dǎo)致其相關(guān)性較高,因此僅運用測井曲線資料難以全面而準(zhǔn)確地識別巖性。因此,傳統(tǒng)的測井解釋方法在判斷復(fù)雜碳酸鹽巖巖性上有很大的局限性,大部分方法僅能識別出灰?guī)r和白云巖或3~4 種巖性[3-4]。隨著機器學(xué)習(xí)技術(shù)的不斷發(fā)展,采用機器智能化學(xué)習(xí)的手段識別復(fù)雜碳酸鹽巖巖性成為了研究的熱點,不少學(xué)者已經(jīng)采用向量機[5-6],神經(jīng)網(wǎng)絡(luò)[7-9],模糊識別[10-11]和傳統(tǒng)決策樹方法[12-14]識別復(fù)雜碳酸鹽巖巖性或碳酸鹽巖成巖相,然而這些方法普遍采用單一學(xué)習(xí)器進行學(xué)習(xí),不能對錯誤樣本進行二次學(xué)習(xí),對復(fù)雜碳酸鹽巖的巖性識別具有一定的局限性。近年來,Chen 等[15]在梯度提升決策樹算法(GBDT)的基礎(chǔ)上進行了改進,提出了一種設(shè)計高效、靈活并且可移植性強的最優(yōu)分布式?jīng)Q策梯度提升庫XGBoost。其原理是通過弱分類器的迭代計算從而實現(xiàn)準(zhǔn)確的分類效果,具有高效性和擴展性強的特征[16]。其優(yōu)勢主要在于在代價函數(shù)中引入了正則化項,控制了模型的復(fù)雜度,降低了模型的方差,防止模型過擬合。同時該算法會在完成一次迭代之后,將葉子節(jié)點的權(quán)重與縮減系數(shù)相乘,來進一步縮減某一個決策樹對整體的影響,從而使得后面有更多的學(xué)習(xí)余地,使其判斷的準(zhǔn)確性較同類方法得到有效提升。同時,該方法之前也鮮用于測井解釋。鑒于此,筆者以測井、錄井資料為基礎(chǔ),采用XGBoost 系統(tǒng)對廊固凹陷北部奧陶系建立碳酸鹽巖識別模型,并以最終分類結(jié)果的準(zhǔn)確率作為評價標(biāo)準(zhǔn),驗證XGBoost 算法應(yīng)用于測井巖性識別的可行性,以期為復(fù)雜碳酸鹽巖巖性的測井識別提供新思路。
GBDT 算法由決策樹和梯度提升兩部分組成,是Friedman[17]提出的一種Boosting 算法。該算法通過讓每一輪迭代得到的損失函數(shù)沿著梯度方向下降來構(gòu)造一個弱分類器函數(shù),然后把多個弱分類器的結(jié)果以一定權(quán)重組合形成強分類器作為最終的預(yù)測輸出[18],其訓(xùn)練過程如圖1 所示。
圖1 GBDT 的訓(xùn)練原理(據(jù)文獻[17]修改)Fig.1 GBDT training principle
XGBoost 是GBDT 算法的優(yōu)化。其特點在于模型能自動利用CPU 進行多線程并行計算,提高運算速度,并且對損失函數(shù)進行泰勒公式二階展開使得預(yù)測精度更高,在損失函數(shù)后面增加正則項,可以約束損失函數(shù)的下降和模型整體的復(fù)雜度[19]。XGBoost 的目標(biāo)函數(shù)T為
式中:l為損失函數(shù);k為分類回歸樹的個數(shù),個;i為模型的預(yù)測值;yi為樣本xi的分類標(biāo)簽;Ω為正則懲罰項函數(shù);fk為第k個樹的模型。
式中:J為每個分類回歸樹葉子節(jié)點的數(shù)量;ω為該個樹葉子節(jié)點的權(quán)重之和;γ和λ 為懲罰系數(shù),為常數(shù),在具體應(yīng)用中可以調(diào)節(jié)。
對于第t輪迭代,模型的目標(biāo)函數(shù)為
式中:ft(xi)表示第t個分類回歸樹;c為常數(shù)項。
對式(3)進行二階泰勒展開
式中:gi表示的一階導(dǎo)數(shù);hi表示的二階導(dǎo)數(shù)。
化簡后,最終目標(biāo)函數(shù)的形式為
本次研究數(shù)據(jù)來源于廊固凹陷北部奧陶系碳酸鹽巖地層,通過層序和巖相古地理分析,整體為潮坪相沉積;研究區(qū)地層經(jīng)歷了加里東運動、海西—印支運動、燕山運動和喜山運動等多次構(gòu)造作用,風(fēng)化剝蝕作用強烈。因此導(dǎo)致該區(qū)地層經(jīng)歷的成巖作用較為復(fù)雜,巖性變化較為頻繁,形成的巖石類型較為多樣。根據(jù)研究區(qū)取心井段的巖心觀察和薄片鏡下鑒定結(jié)果,識別出該研究區(qū)地層主要發(fā)育灰?guī)r、泥質(zhì)灰?guī)r、白云質(zhì)灰?guī)r、灰質(zhì)白云巖、泥質(zhì)白云巖、白云巖等6 種碳酸鹽巖類型(圖2)。其中,含油氣儲集層主要發(fā)育在白云巖和白云質(zhì)灰?guī)r段。所以,對于研究區(qū)地層巖性有效的識別是儲層評價的重要步驟,對后續(xù)進行油氣儲層預(yù)測具有重要意義。
圖2 研究區(qū)碳酸鹽巖巖性鏡下識別(a)白云質(zhì)灰?guī)r,安探2 x 井,5 186.0 m,下馬家溝組上段,單偏光;(b)泥晶灰?guī)r,安探101 x 井,5 815.0 m,上馬家溝組上段,單偏光;(c)泥質(zhì)灰?guī)r,安探4 x 井,峰峰組上段,4 805.3 m,單偏光;(d)灰質(zhì)白云巖,務(wù)古1 井,4 865.0 m,峰峰組下段,單偏光;(e)泥質(zhì)白云巖,安探3 井,4 806.2 m,峰峰組下段,單偏光;(f)粉晶白云巖,安探5 x 井,5 965.0 m,亮甲山組,單偏光。其中藍(lán)色部分為鑄體,紅色為茜素紅染色的方解石Fig.2 Microscopic identification of carbonate lithology in the study area
本文采用研究區(qū)內(nèi)目的層7 口井共2 101 個數(shù)據(jù)點作為巖性識別模型的數(shù)據(jù)集,其中將89 個取心樣品和537 個隨井XRD 數(shù)據(jù)作為本次模型的訓(xùn)練集,訓(xùn)練集中各個巖性的樣本比例如表1 所列,其余1 475 個隨井XRD 數(shù)據(jù)作為測試集來驗證XG‐Boost 算法的巖性識別效果。
在建立巖性識別模型時,選取的測井曲線類型會直接影響預(yù)測精度,故須要進行測井屬性優(yōu)選。本文調(diào)研了前人在碳酸鹽巖儲層相關(guān)研究中常用的測井曲線[20-24],選取了對碳酸鹽巖巖性比較敏感的聲波時差(AC)、自然伽馬(GR)、光電吸收截面指數(shù)(PE)、補償中子(CNL)、補償密度(DEN)、深側(cè)向電阻率(RLLD)等6 種測井曲線作為模型樣本的屬性值(表2)。結(jié)合研究區(qū)不同類型碳酸鹽巖典型巖性測井響應(yīng)特征(圖3),挑選研究區(qū)內(nèi)每種巖性中最符合巖性描述的深度段具有的測井響應(yīng)特征制成圖版,使其具有一定的代表性。從圖3 可以看出,不同巖性的測井響應(yīng)特征不同,泥質(zhì)灰?guī)r和泥質(zhì)白云巖的GR值與其他4 種巖性相差較大,其典型測井曲線值為80~130 API,灰?guī)r、灰質(zhì)白云巖、白云質(zhì)灰?guī)r和白云巖的GR值均較小,為0~50 API;AC,DEN,PE這3 條測井曲線在不同巖性中響應(yīng)特征的區(qū)分度較差,形態(tài)上也沒有太大波動,具體來講,研究區(qū)灰?guī)r的DEN值稍低,PE值稍高,白云巖的PE值稍低;從CNL上看,泥質(zhì)白云巖和灰質(zhì)白云
巖的值稍高,為8%~10%,其次是白云巖,為5%~8%,灰?guī)r、泥質(zhì)灰?guī)r和白云質(zhì)灰?guī)r的值十分接近,均為0~5%;RLLD曲線特征主要呈現(xiàn)為泥質(zhì)灰?guī)r和泥質(zhì)白云巖的相對低值和白云質(zhì)灰?guī)r、灰?guī)r和白云質(zhì)灰?guī)r的相對高值,白云巖的RLLD值位于這兩類巖性之間。
表1 訓(xùn)練集不同巖性樣本比例Table 1 Proportion of different lithological samples in the training set%
表2 研究區(qū)不同巖性原始測井響應(yīng)特征統(tǒng)計Table 2 Original logging response characteristics of different lithologies in the study area
圖3 研究區(qū)不同類型碳酸鹽巖典型測井響應(yīng)曲線特征Fig.3 Characteristics of typical logging response curves of different types of carbonate rocks in the study area
和砂巖的巖性識別不同,碳酸鹽巖通常由于礦物成分相似,非均質(zhì)性較強,使用單一曲線往往很難識別其中的過渡類型。雖然泥質(zhì)灰?guī)r和泥質(zhì)白云巖的典型測井特征與其他類型碳酸鹽巖有所區(qū)別,但從整體來看,研究區(qū)不同類型碳酸鹽巖的測井曲線區(qū)別不是十分明顯,取值區(qū)間較為重疊,說明單一使用某一測井曲線數(shù)據(jù)無法達到識別碳酸鹽巖巖性的目的。為了進一步探究測井曲線的相關(guān)性,通過交會分析各類測井響應(yīng)參數(shù)對于儲層的敏感性,并用置信橢圓來加以表征(圖4)。橢圓越窄,說明2 種測井曲線在該巖性中的相關(guān)性較好,對于巖性區(qū)分有效性較強。在同一個二維測井交會圖中,不同巖性的橢圓離散,說明該二維交會圖對巖性的區(qū)分度較好。研究樣本結(jié)果表明,在研究區(qū)內(nèi)的碳酸鹽巖樣品中,CNL和PE,CNL和RLLD以及GR和DEN的置信橢圓較窄(長軸與短軸之比大于2),總體上看,不同巖性之間的置信橢圓分散程度一般。因此,二維測井曲線交會圖版不能將復(fù)雜碳酸鹽巖巖性很好地區(qū)分開,須要采用機器學(xué)習(xí)的方法進一步識別研究區(qū)碳酸鹽巖巖性。
圖4 二維測井響應(yīng)參數(shù)置信橢圓交會圖(90%置信區(qū)間)圖中置信橢圓表示90%的點所落到的區(qū)域,橢圓越窄,說明2 種測井曲線相關(guān)性越好,橢圓越分散,說明測井響應(yīng)參數(shù)對巖性的判別越敏感Fig.4 Two-dimensional logging response parameter confidence ellipse intersection graph
由于XGBoost 在進行建模時不能識別字符串類型的文本,為了方便編程處理以及避免簡單用數(shù)值替換后數(shù)值的大小對模型產(chǎn)生影響,本文對巖性標(biāo)簽采用了One-HotEncoding 進行了編碼,通過One-Hot 編碼將巖性的n個取值轉(zhuǎn)換為n個二元特征,以稀疏矩陣的形式儲存在模型的標(biāo)簽列中。在進行算法識別巖性時,不同類型的測井曲線具有不同的量綱和數(shù)量級,其差異性會對模型的識別精度產(chǎn)生影響,故對6 種測井曲線進行了標(biāo)準(zhǔn)化處理,由于不同井的測井曲線屬性值大小不一,也不存在統(tǒng)一的上下限標(biāo)準(zhǔn),因此本文采取標(biāo)準(zhǔn)分?jǐn)?shù)的方法進行標(biāo)準(zhǔn)化,其計算式為
式中:z為標(biāo)準(zhǔn)化處理后的數(shù)據(jù);x為樣本數(shù)據(jù);μ為樣本數(shù)據(jù)均值;σ為樣本數(shù)據(jù)的標(biāo)準(zhǔn)差。
2.4.1 模型參數(shù)設(shè)置
通過XGBoost 算法建立巖性識別模型,其中對模型性能影響較大的參數(shù)主要為:迭代次數(shù),即生成決策樹的個數(shù)t,若迭代次數(shù)太大則容易使得模型過擬合,從而導(dǎo)致模型的泛化能力下降;學(xué)習(xí)率eta,用于控制每一次迭代的步長,提高模型的穩(wěn)定性,學(xué)習(xí)率太高會降低模型識別的正確率,太低則會影響模型的運算速度;最大樹深度Dmax與子節(jié)點中最小樣本權(quán)重和Wmin,如果一個葉子節(jié)點的樣本權(quán)重和小于設(shè)置的Wmin,則此次葉子節(jié)點的拆分過程結(jié)束。該參數(shù)用于控制模型的復(fù)雜度,樹太淺會降低模型的準(zhǔn)確率,樹太深則會使得模型過擬合,降低模型的泛化能力;葉子節(jié)點所需的最小損失函數(shù)的下降值γ,γ值越大,算法越保守,計算時間越長;以及用于增加模型隨機性的參數(shù)隨機采樣訓(xùn)練樣本的比例S,若將S設(shè)置為0.5 表示XGBoost 將隨機從整個樣本集合中抽取50%的子樣本建立樹模型,通過調(diào)整S的大小,可以提高模型的穩(wěn)定性,增加最終識別的正確率。
由于本次模型所用數(shù)據(jù)規(guī)模較小,故采用k-折交叉驗證的方法對模型的參數(shù)進行調(diào)優(yōu)。k-折交叉驗證的過程是將數(shù)據(jù)集分成k份,輪流將其中k-1份作為訓(xùn)練數(shù)據(jù),剩余的1 份作為測試數(shù)據(jù)進行試驗。每次試驗都會得出相應(yīng)的測試正確率,k次測試正確率的平均值作為最終的測試正確率[25],在本次調(diào)優(yōu)中根據(jù)樣本量較小的特點將k的值取5。在模型參數(shù)為默認(rèn)值的情況下(eta=0.3,Dmax=6,Wmin=1,S=1,γ=1)對迭代次數(shù)進行調(diào)優(yōu)(圖5),從圖中可以看出,隨著迭代次數(shù)的增加,模型測試集的對數(shù)似然損失值逐漸下降,迭代到420 次左右逐漸穩(wěn)定在0.27 左右;模型的錯誤率隨著迭代次數(shù)的增加先下降后上升,迭代到660 次左右后錯誤率取得最小值,為0.112,以選取最小誤差和迭代精度不再出現(xiàn)明顯變化為原則,最終選定巖性識別模型的迭代次數(shù)為660 次。
在確定迭代次數(shù)的基礎(chǔ)上,對模型的其他參數(shù)進行調(diào)優(yōu),其中對eta,S,γ的調(diào)優(yōu)步長設(shè)置為0.1,Dmax,Wmin的調(diào)優(yōu)步長設(shè)置為1,并按照eta,Dmax,Wmin,S和γ的順序依次進行調(diào)節(jié),最終調(diào)節(jié)完成后模型的參數(shù)為eta=0.1,Dmax,=10,Wmin=6,S=0.5,γ=0.1。
圖5 模型的對數(shù)似然損失值(a)和錯誤率(b)隨迭代次數(shù)的變化Fig.5 Logloss value(a)and error rate(b)of the model vary with the number of iterations
2.4.2 實驗結(jié)果及方法對比
通過交叉檢驗對模型參數(shù)調(diào)節(jié)后,進行模型的訓(xùn)練,并采用測試集樣本對建立的巖性識別模型預(yù)測效果進行驗證,結(jié)果如表3 所列。
表3 巖性識別模型預(yù)測結(jié)果Table 3 Prediction results of lithology identification model
從表3 可以看出,模型對灰?guī)r和白云巖的識別效果均較好,識別準(zhǔn)確率分別達到了91.76%和90%,其次是泥質(zhì)灰?guī)r、白云質(zhì)灰?guī)r和泥質(zhì)白云巖,這3 種巖性的識別率均達到了85%以上,灰質(zhì)白云巖的識別準(zhǔn)確率較低,僅為81.25%。對灰質(zhì)白云巖識別率較低的問題,本文采取了多種方式調(diào)節(jié)參數(shù)均不能提高其準(zhǔn)確率,這可能是因為所選取的測井曲線對灰質(zhì)白云巖的敏感程度較低。若想進一步優(yōu)化,可以采取構(gòu)造測井曲線的特征參數(shù)等特征工程進一步提高巖性識別的準(zhǔn)確率[26-27]。
同時將采用XGBoost 算法的巖性識別結(jié)果與采用分類決策樹、支持向量機方法所得到的巖性識別結(jié)果進行比較(表4)。其中,分類決策樹采用C4.5 算法,表4 中的運行時間為訓(xùn)練和測試模型的總時間,準(zhǔn)確率是測試的結(jié)果。為進一步證明XG‐Boost 算法對碳酸鹽巖巖性識別的優(yōu)越性,本文在測試集中選取了XRD 測試數(shù)據(jù)密度較大的安探4 X 井進行了單井預(yù)測驗證(圖6)。與XRD 測試數(shù)據(jù)的對比結(jié)果顯示,XGBoost 方法對于碳酸鹽巖巖性識別較為準(zhǔn)確,尤其是在巖性快速變化及曲線變化特征不明顯的井段,XGBoost 相較于SVM 和決策樹方法可以更準(zhǔn)確地作出響應(yīng),識別率更高。
表4 不同模型的性能比較Table 4 Comparison of performance of different models
圖6 安探4 x 井單井巖性識別效果圖Fig.6 Lithology identification results of well Antan 4x
對于廊固凹陷北部奧陶系碳酸鹽巖地層,XG‐boost 識別模型對于巖性識別的準(zhǔn)確率為三者最高,達到了88.18%(表4)。其次是分類決策樹,巖性識別的準(zhǔn)確率為76.64%,支持向量機的性能較弱,識別的準(zhǔn)確率為73.59%(表4)。從時間上來看分類決策樹和支持向量機的耗時相近,分別為2.27 s 和2.45 s,XGBoost 由于采用多個弱分類器組合的方式進行機器學(xué)習(xí),其運算相對于其他2 種算法耗時稍長,為3.52 s,由于其采用并行運算,在進行較大模型性訓(xùn)練和測試時并不會大幅度增加運算時間。
綜合來看,采用XGBoost 算法構(gòu)建的復(fù)雜碳酸鹽巖巖性識別模型,通過采用對多個弱學(xué)習(xí)器的線性組合以及在損失函數(shù)后面增加正則項的方式,在稍微增加了運算時長的情況下提升了其巖性識別的準(zhǔn)確率,為復(fù)雜碳酸鹽巖巖性的測井識別方面提供了新的思路。
(1)XGBoost 算法的優(yōu)勢主要在于在代價函數(shù)中引入正則化項,控制了模型的復(fù)雜度,降低了模型的方差,防止模型過擬合。同時在完成一次迭代之后,將葉子節(jié)點的權(quán)重,乘上縮減系數(shù),來縮減一顆樹的影響,從而使得后面有更多的學(xué)習(xí)的余地,可以很好地解決常規(guī)測井巖性識別方法中過擬合以及準(zhǔn)確性不高的問題,使其在利用測井信息在巖性識別的準(zhǔn)確性上較同類方法得到有效提升。
(2)模型對灰?guī)r和白云巖的識別效果均較好,識別準(zhǔn)確率分別達到了91.76%和90%,其次是泥質(zhì)灰?guī)r、白云質(zhì)灰?guī)r和泥質(zhì)白云巖,這3 種巖性的識別率均達到了85%以上,灰質(zhì)白云巖的識別率較低,僅為81.08%,且相較于傳統(tǒng)的機器學(xué)習(xí)模型,其準(zhǔn)確率提升了10%以上。驗證了XGBoost 算法應(yīng)用于測井巖性識別方面的可行性。
(3)數(shù)據(jù)的選取和預(yù)處理、模型的參數(shù)優(yōu)化對于模型的準(zhǔn)確率有著十分重要的影響。通過測井相關(guān)專業(yè)知識、置信橢圓的分析,確定了6 種與巖性較為相關(guān)的測井曲線作為模型的變量數(shù)據(jù),并通過交叉驗證和對模型參數(shù)的調(diào)節(jié)確定了初步最優(yōu)化模型,但參數(shù)調(diào)優(yōu)進行模型優(yōu)化存在優(yōu)化上限,若想進一步優(yōu)化巖性識別模型,須采用特征工程、模型組合等方法。
(4)基于XGBoost 算法的巖性識別模型對復(fù)雜碳酸鹽巖巖性的識別的準(zhǔn)確率優(yōu)于分類決策樹、向量機這2 種傳統(tǒng)的機器識別方法,且由于XGBoost 算法采用多線程和分布式計算的方法,使得訓(xùn)練時間大大縮短,因此該方法可以應(yīng)用于數(shù)據(jù)較為龐大的模型。