張東曉， 尉玉龍， 李三百， 陳云天， 徐加放

(1.南方科技大學(xué) 工程學(xué)院，廣東深圳 518055； 2.中國石油大學(xué)(華東)石油工程學(xué)院，山東青島 266580；3.鵬城實驗室，廣東深圳 518055)

非常規(guī)頁巖油氣資源的有效開發(fā)得益于水力壓裂技術(shù)的應(yīng)用，并且中國非常規(guī)頁油氣儲層非均質(zhì)性較強，儲層物性較差[1]，因此這類油氣資源的主要開采手段為水平井分段多簇射孔壓裂[2-7]。在進(jìn)行射孔和壓裂時存在壓裂段的選擇不合理和壓裂效率低下兩大問題[8]。Cipolla等[9]對100口水平井做了測試發(fā)現(xiàn)，由于壓裂段的地質(zhì)力學(xué)參數(shù)性質(zhì)的差異導(dǎo)致大約40%射孔簇?zé)o效，并且只有約20%的射孔簇與油氣產(chǎn)量有直接關(guān)系。巖石力學(xué)性質(zhì)的差異對裂縫的發(fā)育也有很大影響。Tang等[10]研究了非均質(zhì)性對水力壓裂效果的影響。模擬結(jié)果表明不同程度的非均質(zhì)性對裂縫的發(fā)育有很大影響。趙金洲等[11-12]通過對水平井的分段多簇壓裂研究發(fā)現(xiàn)合理的劃分壓裂段能增大縫網(wǎng)波及的區(qū)域，使儲層有較好的改造效果。吳百烈等[13]對致密氣水平井的分段多簇裂的簇間距和射孔位置進(jìn)行了優(yōu)化研究，結(jié)果表明合理的參數(shù)設(shè)計會降低射孔之間的應(yīng)力干擾。因此為了提高射孔和壓裂效率，需要合理劃分壓裂段和選擇射孔壓裂位置，降低壓裂段內(nèi)儲層性質(zhì)的差異。筆者運用k-means聚類算法對儲層特征進(jìn)行聚類，再根據(jù)聚類結(jié)果指導(dǎo)壓裂段的劃分和射孔壓裂位置的選擇。

1 試驗流程

提出一個優(yōu)化射孔壓裂段劃分和選擇的工作流程。首先通過數(shù)據(jù)驅(qū)動的方法，選取井段的滲流和地質(zhì)力學(xué)參數(shù)作為k-means聚類算法的輸入，輸出是沿井段分布的不同類別；其次訓(xùn)練好的k-means模型可直接預(yù)測其他井沿井段分布的不同聚類類別，避免在大量數(shù)據(jù)集上進(jìn)行模型的重復(fù)訓(xùn)練和尋找不同聚類類別的中心，減少耗時并提高計算效率；再通過FracTHM[14-15]進(jìn)行壓裂過程模擬評價不同類別的壓裂效果；最后根據(jù)不同類別沿井段的分布情況和不同類別的壓裂效果指導(dǎo)射孔壓裂段的劃分和位置的選擇。

2 試驗方法

2.1 k-means聚類

經(jīng)典的k-means聚類方法采用歐氏距離作為樣本間距離的度量方法。因為k-means聚類是一種基于歐氏距離的聚類算法，將所選數(shù)據(jù)集在空間上劃分成不同的類別，保證每一類內(nèi)部的樣本之間的距離最近，同時使不同聚類類別之間的距離盡可能遠(yuǎn)。從物理意義角度看，同一聚類類別井段的儲層滲流性質(zhì)和地質(zhì)力學(xué)性質(zhì)差異較小，不同聚類類別井段的儲層滲流性質(zhì)和地質(zhì)力學(xué)性質(zhì)差異較大。給定m×n的樣本X= {x1,x2,...,xn}，則樣本之間的距離dij為

(1)

式中，xki和xkj分別為第i個樣本和第j個樣本的第k個特征。

定義樣本到分類中心距離的和為最終損失函數(shù)，即

(2)

2.2 FracTHM模擬

FracTHM模擬器[14-15]是一個熱-水-力(THM)完全耦合的模擬頁巖氣開發(fā)過程中的水力壓裂的三維模型，采用混合有限體積法和有限元法的全新壓裂模擬方法，其主要由初始化解耦模型、全耦合模型和斷裂模型3部分組成。

3 試驗數(shù)據(jù)

試驗數(shù)據(jù)來自川南地區(qū)的9口水平井，為A1～A9。每口井的測井曲線包括自然伽馬、縱波時差和橫波時差測井曲線及儲層解釋參數(shù)黏土體積分?jǐn)?shù)。

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

由于測井儀器直接測得的數(shù)據(jù)中包含異常值，并且有些井段的部分?jǐn)?shù)據(jù)點缺失，因此在進(jìn)行試驗之前需要對數(shù)據(jù)進(jìn)行預(yù)處理，包括缺失值處理、異常值處理和數(shù)據(jù)的標(biāo)準(zhǔn)化處理。本文中直接對部分特征缺失的樣本進(jìn)行刪除處理；通過箱型圖進(jìn)行異常值檢測和處理，標(biāo)準(zhǔn)化處理公式為

(3)

3.2 參數(shù)選擇

輸入?yún)?shù):儲層特征為孔隙度、滲透率和黏土體積分?jǐn)?shù)；地質(zhì)力學(xué)特征為彈性模量和泊松比。以上輸入?yún)?shù)為儲層解釋參數(shù)，不能通過測井的方式直接獲取，需要基于每口井中直接測得的測井曲線根據(jù)解釋模型計算得到。本文中根據(jù)井A1～A9中的自然伽馬、縱波時差和橫波時差測井曲線和黏土體積分?jǐn)?shù)計算對應(yīng)孔隙度、滲透率、彈性模量和泊松比。

(1)有效孔隙度計算。首先根據(jù)自然伽馬測井曲線計算泥質(zhì)含量 (Vsh)：

Vsh=0.33(22IGR-1).

(4)

其中

式中，IGR為標(biāo)準(zhǔn)化的自然伽馬測井曲線值；G為實際自然伽馬測井曲線值；Gmin和Gmax分別為每口井的自然伽馬測井曲線的最小值和最大值。

根據(jù)聲波時差測井曲線計算有效孔隙度(φt)：

(5)

式中，Δtma和Δtf分別為巖石骨架和鉆井液聲波時差；Cp為地層壓實系數(shù)；Δt為目的層聲波時差測井值。

(2)滲透率計算。根據(jù)每口井的有效孔隙度和黏土體積分?jǐn)?shù)，滲透率由Kozeny-Carmen模型[16]計算得到。

(3)彈性比模量和泊松比計算。根據(jù)每口井中的縱波時差和橫波時差測井曲線值計算彈性模量和泊松比，計算方法同Chen等[17]所用方法。

4 k-means聚類試驗

首先在所選的訓(xùn)練集A1～A5井訓(xùn)練k-means模型得到不同的聚類類別，并結(jié)合物理意義定性分析聚類類別的壓裂效果；其次用訓(xùn)練好的k-means模型對測試集A6～A9進(jìn)行預(yù)測。

4.1 k-means模型訓(xùn)練

4.1.1 計算聚類趨勢

霍普金斯統(tǒng)計量為描述數(shù)據(jù)的聚類趨勢的量。如果數(shù)據(jù)集沒有明顯的聚類趨勢，霍普金斯統(tǒng)計量的值接近0.5；如果數(shù)據(jù)集聚類趨勢明顯，則其值接近1。計算5口井的霍普金斯統(tǒng)計量為0.94，說明數(shù)據(jù)集有明顯的聚類趨勢。

4.1.2確定聚類數(shù)目k值

采用“肘方法”確定k值，即計算數(shù)據(jù)集中所有樣本點到聚類中心的距離之和的平方(SSE)，根據(jù)SSE隨k的變化選取k值。圖1為“肘方法”確定k值示意圖。隨著k的增加，SSE減小，但當(dāng)k等于5時，SSE變化緩慢，如圖1中紅色圓圈內(nèi)的值為k值。因此k-means將把數(shù)據(jù)集分為5個聚類類別。

圖1 “肘方法”確定k值

4.1.3 聚類結(jié)果

圖2為測試集A1～A5井標(biāo)準(zhǔn)化后的各參數(shù)在不同聚類類別上的數(shù)據(jù)分布箱型圖，可見不同類別的參數(shù)差異較大。類別3的孔隙度、滲透率和黏土體積分?jǐn)?shù)較低，泊松比的值處于相對中間的位置，且彈性模量最大，所以類別3最容易壓裂。類別5的孔隙度和滲透率較高，但黏土體積分?jǐn)?shù)最大，彈性模量最小，泊松比的較高，因此不易壓裂，壓裂效果最差。對比類別1和2：類別1有較高的孔隙度和滲透率，但黏土體積分?jǐn)?shù)較高，并且彈性模量低，因此類別1的壓裂效果較差。對比類別4和類別2：孔隙度、滲透率和黏土體積分?jǐn)?shù)的分布基本一致，但類別4的彈性模量較低，且泊松比較高，因此類別4的壓裂效果較差。

圖2 A1～A5井標(biāo)準(zhǔn)化后的參數(shù)在不同聚類類別上的數(shù)據(jù)分布箱型圖

4.2 k-means模型預(yù)測

經(jīng)過訓(xùn)練的k-means模型在測試集A6～A9進(jìn)行預(yù)測。圖3為A6～A9井標(biāo)準(zhǔn)化后的參數(shù)在不同聚類類別上的數(shù)據(jù)分布箱型圖。對比圖3與圖2，各參數(shù)在不同聚類類別的分布趨勢一致，說明訓(xùn)練好的k-means模型在測試集上有很好的預(yù)測性和較強的泛化能力，可以直接預(yù)測同區(qū)塊內(nèi)的水平井的聚類類別，從而避免模型的重復(fù)訓(xùn)練和重新尋找每個聚類類別的中心，減少耗時。在具體工程實施時可參考已壓裂井段的工程參數(shù)對相同類別的井段進(jìn)行壓裂。

圖3 A6～A9井的標(biāo)準(zhǔn)化后的參數(shù)在不同聚類類別的數(shù)據(jù)分布箱形圖

5 FracTHM模擬試驗

通過數(shù)據(jù)驅(qū)動的聚類算法將井段分成了5類，然而這種基于地層性質(zhì)的分類結(jié)果并不一定反映壓裂的效果。借助FracTHM模擬試驗分析不同類別井段的壓裂效果，進(jìn)一步構(gòu)建基于壓裂效果的聚類模型。

5.1 不同井的相同聚類類別的壓裂過程模擬

本次試驗?zāi)M包括訓(xùn)練集和測試集在內(nèi)的9口井的5個聚類類別的壓裂過程。在該試驗中，孔隙率、滲透率、彈性模量和泊松比取每個類別的平均值。

表1為9口井的不同聚類類別在100 s時的壓裂面積統(tǒng)計。可見A1～A9的5個類別經(jīng)過100 s的壓裂過程，聚類類別3的壓裂面積最大，壓裂效果最好。聚類類別5的壓裂面積最小，壓裂效果最差，其次壓裂效果由好到差依次為聚類類別2、聚類類別4和聚類類別1，并且聚類類別1、2和4的壓裂效果差別較小。同一聚類類別在不同井中的壓裂面積差異較小，并且在訓(xùn)練集和測試集上，同一類別的壓裂面積差異較小，因此訓(xùn)練好的k-means模型在預(yù)測集上有很好的預(yù)測性，可用于預(yù)測同一區(qū)塊其他井的聚類類別。模擬結(jié)果同4.1節(jié)中的定性分析結(jié)果一致。

表1 9口井的不同聚類類別在100 s時的壓裂面積統(tǒng)計

在模擬試驗的基礎(chǔ)上，進(jìn)一步將k-means獲得的井段分類根據(jù)壓裂效果合并。通過FracTHM對A1～A9井的不同類別的壓裂過程進(jìn)行模擬，對比不同井的不同聚類類別的壓裂面積，類別1、2和4的壓裂效果差別較小，因此根據(jù)不同聚類類別的實際壓裂效果，數(shù)據(jù)驅(qū)動的k-means的聚類結(jié)果存在冗余，可以把類別1、2和4并為一類。因此一共有C1、C2和C3三個大類，其中C1包括聚類類別3；C2包括聚類類別1、2、4；C3包括聚類類別5。圖4為A6～A9井的3個類別沿水平井段的分布。

圖4 A6～A9測試井的3個類別沿水平井段的分布

5.2 壓裂分段和射孔位置的優(yōu)化

如圖5所示，聚類類別在井段中的分布主要有兩種情況：聚類類別有較長的連續(xù)分布長度；不同類別交叉出現(xiàn)且連續(xù)分布長度較短。首先模擬聚類之前井段的壓裂過程；其次根據(jù)不同類別的分布情況重新進(jìn)行壓裂段的劃分和射孔壓裂位置的選擇，射孔壓裂效率為

(6)

式中，S為成功壓裂的射孔數(shù)目；T為總的射孔數(shù)目。

5.2.1 聚類類別有較長的連續(xù)分布長度時壓裂過程模擬

A6井2 900～3 100 m水平井段的類別C2和C3具有較長的連續(xù)分布長度，對這一井段進(jìn)行壓裂過程模擬。

圖5為聚類之前A6井2 960～3 060 m的水平井段壓裂模擬示意圖。圖5(a)中黑色箭頭為射孔和壓裂的位置，編號1～5對應(yīng)的壓裂位置分別為2 915、2 955、2 995、3 035和3 075 m。圖5(b)中藍(lán)色區(qū)域為成功壓裂的裂縫形態(tài)，黑色標(biāo)記點為射孔點，但未成功壓裂。可見只有兩個位置成功壓裂。由式(6)計算射孔壓裂效率為40%，并且成功壓裂的裂縫形態(tài)差異較大。這是由于壓裂段內(nèi)的儲層性質(zhì)差異較大，只有部分射孔處的裂縫成功起裂并擴(kuò)展。

圖5 聚類之前A6井2 900～3 100 m的水平井段壓裂模擬示意圖

由于C2和C3有較長的連續(xù)分布長度，因此將井段劃分為兩段進(jìn)行壓裂，使連續(xù)分布長度較長的類別處在單獨的壓裂段。圖6為聚類之后A6井2 960～3 060 m水平井段的壓裂模擬結(jié)果示意圖。在第1段中主要為C2，射孔編號1～4對應(yīng)的壓裂位置分別為2 935、2 955、2 975、2 995 m。在第2段為C1，射孔編號1～3對應(yīng)的壓裂位置分別為3 045、3 065和3 085 m。圖6(b)中藍(lán)色區(qū)域為成功壓裂的裂縫形態(tài)，可見重新劃分的兩段壓裂段都成功壓裂，并且每段中的裂縫形態(tài)發(fā)育面積較大。由式(6)計算兩段的射孔壓裂效率都為100%。

圖6 聚類之后A6井2 960～3 060 m水平井段的壓裂模擬結(jié)果示意圖

通過對比A6井的2 900～3 100 m處水平段在聚類前后的壓裂效果，在同一壓裂條件下，當(dāng)不同類別在井段中具有較長的連續(xù)分布長度時，把連續(xù)分布長度較長的類別處在單獨的壓裂段的壓裂效果明顯優(yōu)于聚類之前的整段的壓裂效果，不但提高射孔壓裂效率，而且裂縫形態(tài)也有更大的面積。

5.2.2 不同類別交叉出現(xiàn)且連續(xù)分布長度較短時壓裂過程模擬

A9井在3 000～3 100 m處3個類別交叉分布，并且每一類在井段中具有較短的連續(xù)分布長度，本次試驗對這一井段進(jìn)行壓裂過程模擬。首先模擬A9在3 000～3 100 m水平井段在聚類之前的壓裂過程；其次根據(jù)A9在3 000～3 100 m水平井段聚類類別的分布情況重新選擇射孔壓裂位置，再進(jìn)行壓裂過程的模擬。

圖7為聚類之前A9井在3 000～3 100 m處的水平井段射孔壓裂位置示意圖。圖7(a)中射孔編號1～5對應(yīng)的壓裂位置分別為3 010、3 030、3 050、3 070和3 090 m；圖7(b)中黑色標(biāo)記點為未成功起裂的射孔點?？梢钥吹皆?50 s時，只有其中3個位置成功壓裂，由式(6)計算射孔壓裂效率為40%，而且成功壓裂的兩個位置的裂縫形態(tài)差異較大，這是由各射孔壓裂之間的相互干擾導(dǎo)致的。

圖7 聚類之前A9井在3 000～3 100 m水平井段的壓裂模擬結(jié)果示意圖

根據(jù)A9井在3 000～3 100 m水平井段的聚類類別分布，設(shè)置該水平井段的射孔壓裂位置。根據(jù)不同類別的壓裂效果，把射孔壓裂位置設(shè)置在壓裂效果好的類別中，或設(shè)置在單一類別中。因此方案1把射孔壓裂位置設(shè)置在壓裂效果好的C1。方案2把射孔壓裂位置設(shè)置在C2。圖8為不同方案的射孔壓裂模擬結(jié)果示意圖。方案1中射孔編號1～3對應(yīng)的壓裂位置分別為3 012、3 070和3 090 m；方案2中射孔編號1～4對應(yīng)的壓裂位置分別為3 019、3 040、3 061和3 095 m。重新選擇射孔壓裂位置后，每個方案中的裂縫形態(tài)發(fā)育也基本一致，由式(6)計算兩種方案的射孔壓裂效率都為100%。由于方案1只在壓裂效果好的C1中進(jìn)行壓裂，方案1的裂縫面積大于方案2的，因此方案1優(yōu)于方案2。

圖8 聚類之后A9井在3 000～3 100 m水平井段的壓裂模擬結(jié)果示意圖

通過對比A9在3 000～3 100 m處水平段在聚類前后的壓裂效果可知，當(dāng)不同類別在距離較短的井段中交叉分布時，在同一壓裂條件下，把射孔壓裂位置設(shè)置在壓裂效果較好的類別中，射孔壓裂效率得到了提高，并且裂縫形態(tài)面積增大。

6 全井段的多段射孔壓裂設(shè)計

對A8井的全井段進(jìn)行多段射孔壓裂的設(shè)計和模擬。對各分段進(jìn)行壓裂過程的模擬時，分別選取每一段中的100 m范圍內(nèi)的井段進(jìn)行壓裂模擬示范。表2為射孔壓裂的位置和對應(yīng)儲層參數(shù)。

表2 A8井射孔壓裂位置

圖9為A8井的壓裂段在250 s時的裂縫形態(tài)?？梢钥吹礁鞫紊淇仔识紴?00%，并且各段中的裂縫形態(tài)能均勻發(fā)育并有較大面積。因此根據(jù)聚類類別在水平井段中的分布優(yōu)化壓裂段劃分和射孔壓裂位置的選擇，保證了裂縫的成功擴(kuò)展。

圖9 A8井的壓裂段在250 s時的裂縫形態(tài)

7 結(jié)束語

基于歐式距離的k-means聚類算法根據(jù)儲層滲流參數(shù)和地質(zhì)力學(xué)參數(shù)對儲層進(jìn)行聚類，劃分的不同聚類類別能夠有效識別可壓裂區(qū)域和易壓裂區(qū)域。訓(xùn)練好的k-means聚類模型具有泛化能力，可直接預(yù)測其他水平井的聚類類別，避免在大量數(shù)據(jù)集上進(jìn)行模型的重復(fù)訓(xùn)練和尋找不同聚類類別的中心，減少耗時，提高計算效率。根據(jù)提出的工作流程可以優(yōu)化壓裂段的劃分和射孔壓裂位置的選擇，當(dāng)不同類別的連續(xù)分布長度較長時，同一壓裂段內(nèi)盡可能只包含同一類別；當(dāng)不同類別在距離較短的井段中交叉分布時，把射孔壓裂位置設(shè)置在壓裂效果較好的類別中，有助于提高射孔壓裂的效率，保證所有裂縫的擴(kuò)展，從而獲得更好的壓裂效果。