基于南海巨厚塑性泥巖地層特征的鉆速預(yù)測(cè)模型

曾小龍，李 謙,2，魏宏超，陳嘉豪，朱海燕

(1.成都理工大學(xué) 環(huán)境與土木工程學(xué)院，四川 成都 610059；2.成都理工大學(xué) 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610059；3.煤炭科學(xué)研究總院，北京 100013；4.中煤科工西安研究院(集團(tuán))有限公司，陜西 西安 710077；5.成都理工大學(xué) 能源學(xué)院，四川 成都 610059)

我國(guó)南海油氣資源豐富，但大多埋藏在高溫高壓且地質(zhì)條件復(fù)雜的地層深處，特別是擁有強(qiáng)塑性和幾乎沒有滲透能力的巨厚泥巖地層以及被稱為世界海上三大高溫高壓地區(qū)之一的鶯瓊盆地。因此，在這樣復(fù)雜的工況下鉆進(jìn)效率較低，甚至引發(fā)次生災(zāi)害[1-2]。若能實(shí)現(xiàn)提前預(yù)測(cè)鉆速，將會(huì)對(duì)指導(dǎo)現(xiàn)場(chǎng)鉆進(jìn)參數(shù)優(yōu)化、穩(wěn)定甚至提高鉆速等起到極大的幫助[3]。由于地層條件的不確定性以及鉆進(jìn)設(shè)備、鉆井液性能等與鉆速之間并沒有明確的影響關(guān)系，使得鉆速大小的變化是非線性的，導(dǎo)致對(duì)鉆速的提前預(yù)測(cè)變得困難且具有挑戰(zhàn)性[4-5]。

早在幾十年前，人們就對(duì)開始了對(duì)鉆速預(yù)測(cè)的研究，但這些基于解析方程的鉆速預(yù)測(cè)模型大多是在簡(jiǎn)化某些因素后建立的，導(dǎo)致在實(shí)際運(yùn)用中具有局限性[6]。近年來，隨著機(jī)器學(xué)習(xí)算法的快速發(fā)展，其強(qiáng)大的性能引起了研究人員的強(qiáng)烈興趣，在通過比較基于預(yù)測(cè)方程的傳統(tǒng)預(yù)測(cè)模型和機(jī)器學(xué)習(xí)預(yù)測(cè)模型的性能，發(fā)現(xiàn)機(jī)器學(xué)習(xí)預(yù)測(cè)模型表現(xiàn)優(yōu)于傳統(tǒng)模型并能獲得較高的預(yù)測(cè)精度[7]。現(xiàn)在一些研究人員已經(jīng)在使用支持向量機(jī)回歸、人工神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林等機(jī)器學(xué)習(xí)算法應(yīng)用于鉆速的預(yù)測(cè)[8]，并且這種方法在實(shí)際線下鉆井中能夠獲得較好的效果[9]，例如C.Hedge 等[10]運(yùn)用機(jī)器學(xué)習(xí)算法預(yù)測(cè)鉆速指導(dǎo)實(shí)際生產(chǎn)，提高了12.56%的鉆進(jìn)效率，節(jié)省約30.12 h 的鉆進(jìn)時(shí)間。

雖然機(jī)器學(xué)習(xí)預(yù)測(cè)模型一般情況下能夠獲得可觀的預(yù)測(cè)結(jié)果，但有時(shí)僅使用單一的算法來構(gòu)建模型，會(huì)限制其性能和進(jìn)步[11]。而粒子群優(yōu)化算法(PSO)、遺傳算法(GA)、貝葉斯優(yōu)化算法(BO)等超參數(shù)優(yōu)化算法的出現(xiàn)，能夠改善這樣的情況，并進(jìn)一步提升模型的性能。在實(shí)際使用中，加入了優(yōu)化算法的模型往往擁有比未使用優(yōu)化算法的模型更高的性能[12]，這也證明了使用優(yōu)化算法的必要性。并非只有通過引入優(yōu)化算法這一種方法來提升模型的性能，也能夠結(jié)合多種機(jī)器學(xué)習(xí)算法來取長(zhǎng)補(bǔ)短。通過對(duì)不同的鉆進(jìn)工況選擇適宜的機(jī)器學(xué)習(xí)算法來建立的預(yù)測(cè)模型，同樣能獲得比單一算法構(gòu)建的模型更高的性能[13]。

使用多種機(jī)器學(xué)習(xí)算法結(jié)合優(yōu)化算法共同預(yù)測(cè)鉆速相比僅使用單一機(jī)器學(xué)習(xí)算法更具有優(yōu)勢(shì)[14]，因此，本文提出了由多種機(jī)器學(xué)習(xí)算法結(jié)合而成的混合模型。其主要由數(shù)據(jù)預(yù)處理、地層類型聚類識(shí)別和鉆速預(yù)測(cè)3 個(gè)步驟組成。數(shù)據(jù)預(yù)處理模塊會(huì)對(duì)原始測(cè)量數(shù)據(jù)進(jìn)行離群值刪除、降噪以及標(biāo)準(zhǔn)化處理，以減小最終預(yù)測(cè)的誤差；地層類型聚類和識(shí)別模塊首先通過對(duì)5 種地層參數(shù)進(jìn)行因子分析達(dá)到降維的目的，并利用因子分析的結(jié)果進(jìn)行聚類，最后利用聚類結(jié)果訓(xùn)練KNN分類模型，達(dá)到識(shí)別地層類型的作用；鉆速預(yù)測(cè)是通過對(duì)聚類結(jié)果分別建立利用貝葉斯優(yōu)化的隨機(jī)森林(BORF)來預(yù)測(cè)鉆速，以此來實(shí)現(xiàn)模型的自動(dòng)化調(diào)參和分類預(yù)測(cè)的目的。

1 地層概況

本文所選取的是中國(guó)南海某鉆采區(qū)域，該地區(qū)典型特征是存在巨厚泥巖覆蓋層，泥質(zhì)較純、微含灰質(zhì)且滲透性極低。其自上而下的地層分布為：樂東組主要由大套泥巖為主。鶯歌海組的上層區(qū)域主要由大套泥巖為主，夾薄砂層；下層區(qū)域以泥質(zhì)粉砂巖為主，夾層狀泥巖與薄層粉砂質(zhì)泥巖。黃流組的上、中層區(qū)域?yàn)榇筇谆疑鄮r夾灰色粉砂巖和泥質(zhì)粉砂巖；下層區(qū)域?yàn)闇\灰色厚層狀細(xì)砂巖。可見，隨著深度的增加，鉆進(jìn)工況也隨之變得復(fù)雜。其巖性占比如圖1 所示。

圖1 地層中不同巖性占比Fig.1 Proportion of different lithology

此處的巨厚泥巖覆蓋層在高圍壓下不僅會(huì)呈現(xiàn)出典型塑性特征而且抗壓強(qiáng)度還會(huì)顯著增大。如圖2 所示，對(duì)井下鉆采的巖心分別在30、45、50、55、60、75 MPa的圍壓下進(jìn)行巖石力學(xué)性能的測(cè)定。此次使用MTS815.03 型電液伺服巖石試驗(yàn)系統(tǒng)，以0.2 MPa/s的加載速率施加軸向荷載直至試件破壞，記錄試驗(yàn)全過程的軸向荷載和變形值。試驗(yàn)最終處理結(jié)果如圖3所示。隨著圍壓的增大，巖心抗壓強(qiáng)度增長(zhǎng)至176.69 MPa；在圍壓超過45 MPa 后，巖心出現(xiàn)了明顯的塑性特征。同時(shí)，通過所采巖心區(qū)域某3 口鉆井的實(shí)際鉆速變化(表1)可以觀察到，實(shí)鉆過程中，在高靜液柱壓力的條件下鉆速顯著下降，巖石力學(xué)性能測(cè)定試驗(yàn)的結(jié)果與之相符。

表1 鉆速-靜液柱壓力對(duì)應(yīng)Table 1 ROP -hydrostatic pressure correspondence

圖2 不同圍壓試驗(yàn)后的巖心Fig.2 Cores after tests under different confining pressures

圖3 不同圍壓下巖心的應(yīng)力應(yīng)變曲線Fig.3 Stress-strain curves of core under different confining pressures

2 方法原理

2.1 數(shù)據(jù)預(yù)處理

如圖4 所示，數(shù)據(jù)預(yù)處理作為整個(gè)過程的最初一步，其結(jié)果直接關(guān)系到最終建立的模型上限，本次預(yù)處理包括3 步：離群值刪除、降噪以及標(biāo)準(zhǔn)化。

圖4 混合模型原理Fig.4 Hybrid model principle

離群值與正常數(shù)值之間有較大的差異，如果不處理，將會(huì)對(duì)后續(xù)的分析與預(yù)測(cè)產(chǎn)生極大的影響。本次離群值處理將使用馬氏距離，其表示各樣本間的協(xié)方差距離[15]：

式中：R為馬氏距離；Z為協(xié)方差矩陣的多變量矩陣；a為均值向量；Σ為協(xié)方差矩陣。

由于本次使用的數(shù)據(jù)量和馬氏距離數(shù)量級(jí)較大，使得一些靠近合理數(shù)值的點(diǎn)無法判斷。因?yàn)轳R氏距離服從卡方分布，故可通過計(jì)算檢驗(yàn)水準(zhǔn)α來尋找出離群值，設(shè)定α=0.005 為閾值，低于0.005 的數(shù)據(jù)將被判定為離群值。

傳感器測(cè)量的數(shù)據(jù)受到各方面因素的影響，導(dǎo)致收集的數(shù)據(jù)包含較多的噪聲，直接應(yīng)用會(huì)產(chǎn)生較大的誤差，需要對(duì)原始數(shù)據(jù)進(jìn)行降噪處理[16]。降噪將使用S-G 濾波器，其首先對(duì)固定長(zhǎng)度的滑動(dòng)窗口內(nèi)y個(gè)數(shù)據(jù)進(jìn)行k-1 次多項(xiàng)式擬合：

設(shè)定數(shù)據(jù)值向量為Y，系數(shù)矩陣為A，自變量矩陣為X，殘差為ε，可將擬合方程改寫為：

最后經(jīng)過最小二乘法處理后得到最終的降噪數(shù)據(jù)：

若將原始數(shù)據(jù)直接用于后續(xù)的機(jī)器學(xué)習(xí)建模，數(shù)值較大的數(shù)據(jù)在模型中作用會(huì)更為凸顯，小數(shù)值的數(shù)據(jù)影響甚微，不能夠完全展現(xiàn)各數(shù)據(jù)對(duì)鉆速的影響，因此需要對(duì)原始數(shù)據(jù)使用標(biāo)準(zhǔn)化。本次使用Min-Max標(biāo)準(zhǔn)化對(duì)數(shù)據(jù)進(jìn)行處理，即：

2.2 地層類型聚類和識(shí)別

為了達(dá)到地層類型聚類與識(shí)別目的，將采用因子分析、K-Means++、KNN 這3 種方法實(shí)現(xiàn)。

因子分析通過利用指標(biāo)間的相關(guān)性系數(shù)矩陣來對(duì)各個(gè)指標(biāo)的相關(guān)性大小來進(jìn)行分類，使得最后的結(jié)果中各組內(nèi)部指標(biāo)間的相關(guān)性高，不同組間的相關(guān)性低，并將標(biāo)示各組數(shù)據(jù)基本結(jié)構(gòu)的新指標(biāo)稱作公共因子。因子分析原理構(gòu)建的數(shù)學(xué)模型可表示為：

式中：xm為第m個(gè)樣本；fn為第n個(gè)公共因子；cmn fn為構(gòu)建樣本xm的第n個(gè)公共因子fn的系數(shù)。

由于原始K-Means 算法中初始聚類中心的選取是在樣本中選擇K個(gè)點(diǎn)，這往往會(huì)導(dǎo)致算法陷入局部最優(yōu)，從而導(dǎo)致K-Means 算法具有較大的局限性。為獲得算法的全局最優(yōu)，提出了K-Means++算法[17](圖5)。其對(duì)于K-Means 算法的改進(jìn)主要在對(duì)初始中心的選取上，選取過程如下。

圖5 K-Means++算法原理Fig.5 Principles of K-Means++algorithm

(1) 設(shè)定最終的聚類中心個(gè)數(shù)為K后，隨機(jī)在數(shù)據(jù)集S內(nèi)選擇一個(gè)樣本點(diǎn)x1作為第一個(gè)聚類中心C1。

(2) 計(jì)算其余點(diǎn)到所有聚類中心的最短距離(由于本文使用的數(shù)據(jù)進(jìn)行了標(biāo)準(zhǔn)化，因此選擇歐氏距離)，計(jì)算式如下：

式中：xn為第n個(gè)樣本點(diǎn)；C1,2,···,m為聚物類中心；D(xn)為第n個(gè)樣本點(diǎn)最近的聚類中心的歐氏距離；d(xn,Cm)為xn到Cm的歐氏距離。

(3) 計(jì)算這些點(diǎn)被選擇為下一個(gè)聚類中心的概率P(xn)，計(jì)算式如下：

(4) 選取概率最高的樣本點(diǎn)將成為新的聚類中心，并重復(fù)步驟(2)和步驟(3)，直至能夠同時(shí)選取K個(gè)聚類中心。

聚類是在并沒有真實(shí)分類標(biāo)簽下開展的，因此無法直接判斷聚類結(jié)果的好壞，所以需要使用指標(biāo)來評(píng)價(jià)聚類效果。本文使用輪廓系數(shù)s作為指標(biāo)，表示為：

式中：A為所選樣本點(diǎn)和同簇內(nèi)其他點(diǎn)間的平均距離；B為所選樣本點(diǎn)與某一最近簇中所有點(diǎn)的平均距離。

其通過將某一對(duì)象與自身所在簇的相似程度和其他簇的相似程度來進(jìn)行評(píng)價(jià)聚類效果。輪廓系數(shù)的取值范圍在[-1,1]內(nèi)，當(dāng)值越接近于1 時(shí)，所選樣本點(diǎn)和簇內(nèi)點(diǎn)越相似，與簇外點(diǎn)越不相似。

地層類型識(shí)別使用的KNN 算法，其整個(gè)過程是以不同特征值之間的距離作為分類依據(jù)，即所選取的一個(gè)樣本在特征空間中的M個(gè)最鄰近樣本中大部分屬于占比最多的一類，該樣本也屬于這個(gè)占比最多的類別。

2.3 鉆速預(yù)測(cè)

混合模型內(nèi)鉆速的預(yù)測(cè)選擇使用隨機(jī)森林(圖6)。隨機(jī)森林的核心在于決策樹和Bagging 算法。對(duì)于回歸問題來說，最終預(yù)測(cè)結(jié)果由各個(gè)決策樹投票取均值決定。隨機(jī)森林相比于決策樹、支持向量機(jī)等機(jī)器學(xué)習(xí)算法能夠減少離群值、異常值的影響，并且能夠有效地預(yù)防過擬合現(xiàn)象的出現(xiàn)，并且對(duì)于處理多變量和大數(shù)據(jù)量有較大優(yōu)勢(shì)[18]。

圖6 隨機(jī)森林結(jié)構(gòu)Fig.6 Structure of Random Forest

隨機(jī)森林由眾多決策樹結(jié)合而成，而樹的多少、結(jié)構(gòu)以及分枝都會(huì)影響隨機(jī)森林的預(yù)測(cè)結(jié)果，因此本文將使用樹數(shù)目、樹的最大深度和樹的特征數(shù)3 個(gè)超參數(shù)對(duì)模型進(jìn)行優(yōu)化。對(duì)這些超參數(shù)的調(diào)整將會(huì)直接影響模型的結(jié)構(gòu)和預(yù)測(cè)結(jié)果，因此，需要尋找到最優(yōu)的超參數(shù)組合。貝葉斯優(yōu)化在解決像尋找最優(yōu)超參數(shù)這樣的黑箱問題時(shí)，能夠在較少的評(píng)估過程下尋找到復(fù)雜未知函數(shù)的最優(yōu)解[19]。貝葉斯優(yōu)化流程如圖7 所示，整個(gè)過程是利用高斯過程生成先驗(yàn)分布來逼近未知目標(biāo)函數(shù)的后驗(yàn)分布，然后根據(jù)分布選擇下一個(gè)采樣的超參數(shù)組合，以此重復(fù)并不斷更新，直至尋找到最優(yōu)的結(jié)果。

為了保證模型的泛化能力，在模型建立及超參數(shù)優(yōu)化時(shí)采用多重交叉驗(yàn)證。采用五折交叉驗(yàn)證，并將5 次的平均驗(yàn)證結(jié)果作為模型的最終性能(圖8)。

圖8 N 折交叉驗(yàn)證Fig.8 Principle of N-fold cross-validation

2.4 模型評(píng)價(jià)指標(biāo)

由于涉及到分類模型和回歸模型的建立，因此，需要兩類不同評(píng)價(jià)指標(biāo)。地層類型識(shí)別中建立的KNN分類模型將使用準(zhǔn)確率(ACC)作為評(píng)價(jià)指標(biāo)。依據(jù)定義，ACC 越大，模型預(yù)測(cè)結(jié)果越準(zhǔn)確；ACC 越小，模型的預(yù)測(cè)性能越差，其公式為：

式中：yture為預(yù)測(cè)正確的總數(shù)；n為總樣本量。

鉆速預(yù)測(cè)模型將使用均方根誤差(ERMS)、R2和平均絕對(duì)誤差(EMA)作為評(píng)價(jià)指標(biāo)，各指標(biāo)數(shù)學(xué)表達(dá)式為：

式中：yi為 真實(shí)值；為預(yù)測(cè)值；為平均真實(shí)值。

依據(jù)定義，ERMS越小，預(yù)測(cè)值與原始數(shù)據(jù)間的差距越小，模型性能越好。R2范圍在0～1，R2越接近1，則預(yù)測(cè)效果越好。EMA代表預(yù)測(cè)值到真實(shí)值的平均距離，因此平均絕對(duì)誤差越小，預(yù)測(cè)效果越好。

3 結(jié)果與討論

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

本文使用的全部數(shù)據(jù)來自于中國(guó)南海某區(qū)域10 口鉆井共計(jì)21 913 組實(shí)測(cè)數(shù)據(jù)，其中包括5 大類44 種不同的參數(shù)。若將全部參數(shù)用來預(yù)測(cè)建模，必然會(huì)使模型復(fù)雜度過高導(dǎo)致預(yù)測(cè)反應(yīng)速度較慢。因此，除選擇了鉆速作為目標(biāo)參數(shù)以及常用的深度、鉆壓、轉(zhuǎn)速、泵量和扭矩這些參數(shù)外[20-22]，額外選擇了入口泥漿密度、出口泥漿密度和地層特征參數(shù)：地震波速、孔隙壓力、破裂壓力、上覆壓力、地層巖性；總共4 類14 種參數(shù)，其詳細(xì)內(nèi)容見表2。各參數(shù)在省略單位后數(shù)量級(jí)的對(duì)比如圖9 所示，深度、泵量等參數(shù)數(shù)量級(jí)能達(dá)到上千，而孔隙壓力、入出口泥漿密度等參數(shù)最高僅達(dá)到個(gè)位數(shù)，各參數(shù)之間數(shù)量級(jí)差距巨大；并且部分參數(shù)存在一定的離群值，因此想要使用和分析這些數(shù)據(jù)，需要提前對(duì)原始數(shù)據(jù)進(jìn)行預(yù)處理消除這些隱患。

表2 鉆進(jìn)工藝參數(shù)Table 2 Drilling parameters

圖9 各參數(shù)箱型圖Fig.9 Box plots of various parameters

3.2 數(shù)據(jù)集劃分

選取鉆速作為目標(biāo)參數(shù)，深度、工藝參數(shù)和地層特征共3 類13 種參數(shù)作為輸入?yún)?shù)。隨機(jī)選取所有數(shù)據(jù)的80%作為訓(xùn)練集，剩余20%數(shù)據(jù)作為測(cè)試集。訓(xùn)練集會(huì)直接參與各模型的訓(xùn)練與建立，并且在訓(xùn)練階段使用五折交叉驗(yàn)證保證建立模型的泛化能力；而測(cè)試集直接用于最終模型性能測(cè)試。

3.3 數(shù)據(jù)預(yù)處理

首先進(jìn)行原始數(shù)據(jù)中離群值的去除，通過計(jì)算檢驗(yàn)水準(zhǔn)α，設(shè)定α閾值為0.005，最終確定的離群值如圖10 所示，此次判定為離群值的點(diǎn)共1 740 個(gè)，約占樣本總量的7.94%。可以觀察到被判定為離群值的某些點(diǎn)和合理數(shù)值的之間有較大差距，若不剔除離群點(diǎn)，必然會(huì)在后續(xù)的標(biāo)準(zhǔn)化、分析和預(yù)測(cè)時(shí)產(chǎn)生較大的誤差。

圖10 不同鉆速的馬氏距離Fig.10 Mahalanobis distance of different ROP

鉆速、鉆壓、轉(zhuǎn)速、扭矩這些數(shù)據(jù)受到傳感器、鉆進(jìn)工況和地質(zhì)環(huán)境等眾多因素影響，導(dǎo)致收集到的數(shù)據(jù)內(nèi)包含較多的噪聲數(shù)據(jù)，因此需要進(jìn)行降噪處理。圖11 為真實(shí)數(shù)據(jù)(紅色)與降噪數(shù)據(jù)(藍(lán)色)間的對(duì)比，降噪數(shù)據(jù)趨勢(shì)與真實(shí)數(shù)據(jù)趨勢(shì)變化一致，說明達(dá)到了良好的降噪效果。

圖11 各參數(shù)降噪前后對(duì)比Fig.11 Comparison of parameters before and after noise reduction

完成以上操作后，對(duì)14 種數(shù)據(jù)進(jìn)行Min-Max 標(biāo)準(zhǔn)化處理。以鉆速為例，圖12 為鉆速標(biāo)準(zhǔn)化前后對(duì)比，所有數(shù)據(jù)在處理后被縮小到[0,1]的區(qū)間內(nèi)，指標(biāo)本身的量綱得以消除，并且在處理后并沒有改變曲線形狀。

圖12 鉆速標(biāo)準(zhǔn)化前后對(duì)比Fig.12 Comparison of ROP before and after standardization

3.4 地層類別聚類和識(shí)別

在聚類時(shí)，將使用訓(xùn)練集中的地震波速、孔隙壓力、破裂壓力、上覆壓力和巖性5 種地層參數(shù)作為依據(jù)。如果直接使用這5 類指標(biāo)進(jìn)行分析，會(huì)造成分析的維度較多而增加工作量。如圖13 所示，這些參數(shù)間均具有一定的相關(guān)性，甚至孔隙壓力、破裂壓力和上覆壓力間的相關(guān)性系數(shù)均在0.9 以上，這必然會(huì)導(dǎo)致部分信息產(chǎn)生重復(fù)，使得分析更加復(fù)雜。因此，使用因子分析來對(duì)這些原始數(shù)據(jù)進(jìn)行降維和消除各指標(biāo)間的相關(guān)性是必要的。

圖13 各地層參數(shù)間皮爾遜相關(guān)性系數(shù)熱力圖Fig.13 Heat map of Pearson correlation coefficients between formation characteristics.

首先對(duì)5 種地層參數(shù)進(jìn)行因子分析，并確定公共因子的數(shù)目。這些提取出來的公共因子均能夠在一定程度上反映原數(shù)據(jù)中所包含的信息，但其所能反映的信息量未知，所以就需要累計(jì)方差貢獻(xiàn)率這樣的指標(biāo)來度量。累計(jì)方差貢獻(xiàn)率越大，公共因子所能反映的信息越接近于原始數(shù)據(jù)。一般要求累計(jì)貢獻(xiàn)率需達(dá)到85%，為了能夠盡可能獲得最大信息，此次分析中將累計(jì)方差貢獻(xiàn)率閾值設(shè)為90%。在公共因子數(shù)目為3 時(shí)，累計(jì)方差貢獻(xiàn)率達(dá)到98.563%，公共因子數(shù)目因此確定為3(表3)。

表3 累計(jì)貢獻(xiàn)率Table 3 Cumulative contribution rate

由式(6)確定各個(gè)變量在每一個(gè)公共因子下的得分，其結(jié)果見表4，公共因子1 主要是由孔隙壓力、破裂壓力和上覆壓力共同決定，公共因子2 由巖性所決定，公共因子3 由地震波速所決定。

表4 各變量因子得分Table 4 Scores of variable factors

以各組數(shù)據(jù)在3 個(gè)公共因子下的得分作為三維坐標(biāo)，在聚類中心數(shù)目從2 到20 的條件下分別進(jìn)行K-Means++聚類，并記錄每次聚類的輪廓系數(shù)，其結(jié)果如圖14 與圖15 所示，當(dāng)聚類中心數(shù)目為2 時(shí)，獲得的輪廓系數(shù)最大。

圖14 各聚類數(shù)的輪廓系數(shù)Fig.14 Silhouette coefficients for different numbers of clustering centers

因此，最終聚類結(jié)果可分為2 種地層類型。各地層詳細(xì)數(shù)據(jù)(圖16)，兩種地層類型之間存在明顯的區(qū)別：地層類型1 的地震波速、孔隙壓力和破裂壓力的上、下四分位數(shù)均小于地層類型2，且地層類型1 的地震波速的上、下四分位數(shù)間的差值206.19 m/s 比地層類型2 的小359.03 m/s，差距較大，但地層類型1 的孔隙壓力和破裂壓力的上、下四分位數(shù)差值比地層類型2的高出13.20 MPa 和12.09 MPa；地層類型2 的上覆壓力的下四分位數(shù)為80.71 MPa，比地層類型1 的上四分位數(shù)78.54 MPa 更大，但上、下四分位數(shù)間的差值比地層類型1 的小11.93 MPa。從巖性來看，地層類型1 的巖性是以泥巖和粉砂質(zhì)泥巖為主，而地層類型2的巖性是以粉砂巖、細(xì)砂巖和中砂巖為主的地層類型。此聚類結(jié)果同現(xiàn)場(chǎng)地層數(shù)據(jù)吻合，將在高圍壓下具備高塑性的泥巖單獨(dú)成組，符合實(shí)際需求。

圖16 地層特征參數(shù)箱型圖對(duì)比Fig.16 Comparison of formation characteristic parameter box plots

為實(shí)現(xiàn)地層識(shí)別，需要使用聚類結(jié)果訓(xùn)練KNN分類模型。選取訓(xùn)練集中的5 種地層特征參數(shù)作為輸入?yún)?shù)，地層類型作為目標(biāo)參數(shù)，并在模型訓(xùn)練時(shí)使用十折交叉驗(yàn)證，以ACC 為評(píng)價(jià)指標(biāo)。最終，以KNN模型在十次交叉驗(yàn)證的ACC 的平均值為結(jié)果，其ACC可達(dá)到0.997。下文將以聚類結(jié)果分別對(duì)各地層類型建立對(duì)應(yīng)的預(yù)測(cè)模型。

3.5 鉆速預(yù)測(cè)模型

為評(píng)估模型性能，將混合模型同貝葉斯優(yōu)化的隨機(jī)森林(BO-RF)、隨機(jī)森林(RF)、決策樹(DT)、支持向量回歸(SVR)和BP 神經(jīng)網(wǎng)絡(luò)(BP-ANN)進(jìn)行比較。其中，混合模型能夠進(jìn)行分類預(yù)測(cè)，貝葉斯優(yōu)化的隨機(jī)森林不具備分類預(yù)測(cè)；貝葉斯優(yōu)化的隨機(jī)森林和隨機(jī)森林對(duì)比進(jìn)行貝葉斯優(yōu)化性能的測(cè)試；隨機(jī)森林同決策樹、支持向量回歸和BP 神經(jīng)網(wǎng)絡(luò)進(jìn)行模型的預(yù)測(cè)性能對(duì)比。設(shè)定貝葉斯優(yōu)化初始節(jié)點(diǎn)為30；迭代次數(shù)為100；樹數(shù)目范圍為[1,100]，樹的最大深度范圍為[1,40]，樹的特征數(shù)目范圍為[1,20]。除混合模型和BO-RF 使用貝葉斯優(yōu)化外，剩余模型均在訓(xùn)練時(shí)使用網(wǎng)格搜索調(diào)整超參數(shù)。

通過對(duì)比圖17 中6 個(gè)模型在測(cè)試集下的預(yù)測(cè)結(jié)果，混合模型的預(yù)測(cè)結(jié)果更加貼合和集中于1∶1 line，并且偏離1∶1 line 的點(diǎn)即預(yù)測(cè)誤差較大的點(diǎn)更少，整體預(yù)測(cè)效果更好；BO-RF 相比RF 預(yù)測(cè)結(jié)果更加貼合1∶1 line，表明了貝葉斯優(yōu)化具有良好的超參數(shù)優(yōu)化能力。圖17 為混合模型、BO-RF、RF、DT、SVR 與BP-ANN 在經(jīng)過訓(xùn)練集訓(xùn)練后，在測(cè)試集下的R2、ERMS和EMA，通過對(duì)比，無論是R2、ERMS還是EMA，混合模型均優(yōu)于其他模型。表5 為BO-RF 與混合模型通過貝葉斯優(yōu)化尋找到的最佳超參數(shù)，以及混合模型和BO-RF 在訓(xùn)練集下使用貝葉斯優(yōu)化尋找最佳超參數(shù)的耗時(shí)。混合模型內(nèi)的BO-RF_1 和BO-RF_2 模型樹數(shù)目、最大深度以及特征數(shù)均小于BO-RF，讓混合模型相比于BO-RF 具有更小的復(fù)雜度，使得混合模型的貝葉斯優(yōu)化耗時(shí)低于BO-RF。并且由于混合模型是在地層聚類的基礎(chǔ)上對(duì)聚類結(jié)果分別建立預(yù)測(cè)模型，使得模型的建立更加細(xì)化且有針對(duì)性。因此，混合模型能夠獲得較好的預(yù)測(cè)性能和更快的反應(yīng)速度。

表5 各模型使用貝葉斯優(yōu)化獲得的最優(yōu)超參數(shù)和耗時(shí)Table 5 Optimal hyperparameters and time consumption of each model using Bayesian optimization

圖17 不同模型的最佳預(yù)測(cè)結(jié)果Fig.17 Optimal prediction results of each model

4 結(jié)論

a.為建立南海巨厚塑性泥巖地層特征的鉆速預(yù)測(cè)模型，使用K-Means++算法并基于其中5 種地層參數(shù)對(duì)原始數(shù)據(jù)進(jìn)行聚類分析，并在此基礎(chǔ)上建立了地層類型分類和針對(duì)不同地層類型分別建模的混合模型。基于本次使用的數(shù)據(jù)集，本文所提出的預(yù)測(cè)模型相比于常見的各類預(yù)測(cè)模型，具有較好的預(yù)測(cè)性能。

b.基于K-means++算法，找到了原始數(shù)據(jù)的潛在分類，即將數(shù)據(jù)集中的數(shù)據(jù)劃分成了兩種類型，其中一類是以泥巖和粉砂質(zhì)泥巖為主的地層類型，另一類是以粉砂巖、細(xì)砂巖和中砂巖為主的地層類型。此聚類結(jié)果與現(xiàn)場(chǎng)地層數(shù)據(jù)吻合，符合實(shí)際需求。

c.使用的貝葉斯優(yōu)化算法的BO-RF 無論從R2、ERMS、EMA還是預(yù)測(cè)結(jié)果上，均要優(yōu)于RF，表明了基于本數(shù)據(jù)集建立的鉆速預(yù)測(cè)模型中，貝葉斯優(yōu)化算法具有較好的超參數(shù)優(yōu)化性能。

d.使用的K-means++聚類算法僅在本次使用的數(shù)據(jù)集中表現(xiàn)良好，在其他數(shù)據(jù)集中表現(xiàn)未知，后續(xù)會(huì)進(jìn)一步開展對(duì)地層類型聚類的研究，尋找和建立更加具有普適性的地層聚類算法。