基于最小體積封閉橢球算法的多段水力壓裂改造體積估算

王旭林， 桂志先， 王 鵬*， 易娟子

(1.長江大學(xué)非常規(guī)油氣湖北省協(xié)同創(chuàng)新中心, 油氣資源與勘探技術(shù)教育部重點實驗室,武漢 430100；2.中國石油天然氣股份有限公司西南油氣田分公司重慶氣礦,重慶 400021)

頁巖油氣作為一種新型能源,因其資源十分豐富,越來越被人們所重視。頁巖油氣儲層具有低孔隙度、低連通性、低滲透率的物性特性,通常需要采取增產(chǎn)措施來改善其巖石物理特性[1]。由于頁巖油氣藏的特殊性,在勘探開發(fā)中影響頁巖油氣井產(chǎn)能的因素有很多。研究表明,儲層壓裂改造體積是影響頁巖油氣藏生產(chǎn)效果的關(guān)鍵因素之一[2-3]。前人研究發(fā)現(xiàn)被壓裂改造油氣藏的生產(chǎn)力與被壓裂改造體積正相關(guān),改造的體積越大,其相應(yīng)的增產(chǎn)效果越明顯[4-5]。因此對壓裂改造體積(stimulated reservior volume,SRV)的準確估算已經(jīng)成為頁巖壓裂研究領(lǐng)域的重難點問題。前人分別采用通道長度與通道寬度來表征裂縫擴展的長度和寬度,水平井長度表征SRV[6]。上述理論模型方法均采用立方體來擬合SRV,但包括較多的無效范圍,所以擬合效果較為粗略[7]。溫慶志等[8]結(jié)合壓裂縫網(wǎng)絡(luò)形態(tài),通過研究均質(zhì)地下儲層中主裂縫周邊的應(yīng)力場發(fā)現(xiàn),地應(yīng)力場方向改變的區(qū)域形成縫網(wǎng),縫網(wǎng)區(qū)域大致為橢球體,因此可以近似地認為壓裂形成的縫網(wǎng)為橢球體。同時,微震事件的空間展布特征也被用來估算SRV[9]。因此,可將壓裂形成的微震事件的空間展布近似處理為橢球體,來估算儲層壓裂改造體積。邵媛媛等[7]基于地震事件震源點將壓裂改造體積擬合成橢球體,建立了儲層改造體積最小橢球模型,通過對比其他方法得出最小體積橢球模型獲得了更為詳細的SRV幾何結(jié)構(gòu)。但上述改進方法只能得出儲層改造體積的近似形狀與各種參數(shù)對改造體積估算的影響,尚未給出比較完整的橢球體積估算方法。在實際生產(chǎn)過程中,多段壓裂施工作業(yè)會出現(xiàn)各段壓裂裂縫網(wǎng)絡(luò)存在交錯疊合的情況,并且由于多種因素的影響,單段壓裂裂縫可能會呈現(xiàn)出較為復(fù)雜的網(wǎng)絡(luò)形態(tài),即表現(xiàn)為微震事件的空間展布形態(tài)復(fù)雜。因此,使用單一的模型很難解決上述復(fù)雜的壓裂裂縫網(wǎng)絡(luò)的壓裂改造體積的估算。

本文針對實際壓裂改造過程中存在的多段壓裂改造體積存在重疊和單段壓裂形成的復(fù)雜裂縫網(wǎng)絡(luò)兩類情況進行分析,估算對應(yīng)的壓裂改造體積。對多段壓裂體積重疊的情況,首先利用DBSCAN(density based spatial clustering of application with noise)算法對各段微震事件分別進行處理,消除異常的微震事件,再利用MVEE算法估算出各段微震事件對應(yīng)的最小體積封閉橢球體,并獲得橢球體的相關(guān)參數(shù)(橢球體中心點、特征值和特征向量等),然后對壓裂作業(yè)區(qū)進行網(wǎng)格化,結(jié)合各個橢球體參數(shù),使用仿射變換判斷各個網(wǎng)格點與每個橢球體的歸屬關(guān)系,進而統(tǒng)計出各段壓裂改造體積間的重疊部分；針對單段壓裂形成的復(fù)雜裂縫網(wǎng)絡(luò)形態(tài)的情況,采用DBSCAN算法對微震事件進行分組,將復(fù)雜的微震事件展布轉(zhuǎn)化為多個展布相對簡單的分組,以各個分組的微震事件為基礎(chǔ),計算對應(yīng)的最小體積橢球體,并消除各個橢球體間的重疊部分,進而估算出單段壓裂的整體改造體積。本文還利用兩組實際微震監(jiān)測數(shù)據(jù),對上述兩類情況進行分析與估算,以驗證本文方法的可行性與準確性。

1 方法原理

1.1 最小體積封閉橢球

最小體積封閉橢球算法(MVEE)[10]是通過尋找一個n維最小封閉橢球以包含給定的所有點集合,并獲取相應(yīng)橢球的信息,如橢球軸半徑、橢球中心、橢球體積等。MVEE問題在數(shù)據(jù)挖掘、模式識別等領(lǐng)域應(yīng)用十分廣泛。

任意一個橢球可表示為

E(a,A)={x∈R3|(x-a)TA(x-a)≤1}

(1)

式(1)中，a∈R3為該橢球體的中心點；A∈R3×3為對稱正定矩陣,該矩陣實現(xiàn)橢球體向單位球體的仿射變換,矩陣A的特征值(λ1、λ2、λ3)平方根倒數(shù)分別為橢球體半軸長,即實現(xiàn)仿射變換的縮放處理；特征向量實現(xiàn)仿射變換的旋轉(zhuǎn)處理。式(1)也可表示為

E=Mtranslation×Mrotation×Mscaling×C

(2)

式(2)中，Mtranslation、Mrotation和Mscaling分別為平移、旋轉(zhuǎn)和縮放三種仿射變換；C為單位球體；E表示最終的橢球體。Mrotation與矩陣A的特征向量矩陣互為逆矩陣；Mtranslation由中心點a控制。

以微震事件的展布特征為基礎(chǔ)獲得的最小體積封閉橢球體的體積可表示為

(3)

求解最小體積封閉橢球問題轉(zhuǎn)可化為最優(yōu)化問題:

minA[-(ΔA)1/2]， (x-a)TA(x-a)≤1

(4)

由于-(ΔA)1/2不是一個凸函數(shù),所以式(4)的優(yōu)化模型不是一個凸優(yōu)化模型,不具備凸優(yōu)化的性質(zhì)。利用KY算法[11]對之進行凸優(yōu)化處理,并最終實現(xiàn)對最小體積封閉橢球的求解,獲得橢球體的中心點a與矩陣A。

1.2 DBSCAN算法

實際微震監(jiān)測數(shù)據(jù)處理會受到各種因素的影響,例如速構(gòu)建度模型的誤差、微震事件初至信息無法準確拾取等，這些因素都會影響到最終的微震事件定位精度。同時,由于頁巖儲層的非均質(zhì)性、天然裂縫的發(fā)育程度等因素的影響,壓裂形成的裂縫結(jié)構(gòu)更為復(fù)雜。因此,需要對微震事件進行適當?shù)念A(yù)處理,消除定位異常的影響,并對復(fù)雜壓裂裂縫的分組處理。

DBSCAN是一種基于密度的空間聚類分析算法,能將具有足夠密度的區(qū)域劃分為簇,在具有噪聲的空間數(shù)據(jù)中發(fā)現(xiàn)任意形狀的簇[12]。該算法被廣泛應(yīng)用于各個領(lǐng)域的數(shù)據(jù)處理與分析,例如提高雷電定位算法的抗干擾能力[13]、數(shù)據(jù)挖掘[14]、分析電力工程數(shù)據(jù)的完整性[15]和巖石物理實驗數(shù)據(jù)的處理[16]。其中，Eps和MinPts是整個算法最為關(guān)鍵的參數(shù),其含義分別為單簇內(nèi)樣點間的最大距離(又稱聚類半徑)和單簇包含的最小樣點數(shù)。并由此引申出核心點、直接密度可達、密度可達和密度相連等相關(guān)概念。若某個樣點的Eps臨域內(nèi),臨近樣點數(shù)大于或等于Minpts,則該樣點就是核心點；若某樣點到達核心點的距離小于Eps,就說明從該樣點到核心點直接密度可達；對于一個樣點序列,樣點間依次滿足直接密度可達,即直接密度可達滿足傳遞性,則這些樣點是密度可達；若兩個樣點對于同一樣點密度可達,就認為這兩個樣點是密度相連[17]。

DBSCAN的處理流程可描述為:①在整個樣點數(shù)據(jù)集中隨機挑選出一個樣點,并確保該樣點為核心點；②以該樣點為起點,尋找到所有密度相連的數(shù)據(jù)點,并將尋找結(jié)果建立新簇；③重新掃描樣點數(shù)據(jù)集(不包括之前已有簇中的任何樣點)尋找新的核心點,再重復(fù)步驟②,直到?jīng)]有新的核心點為止[18]。原始樣點數(shù)據(jù)集中沒有被包含在任何簇中的樣點就構(gòu)成異常點。

對于微震事件空間展布為簡單形態(tài)的情況,使用DBSCAN算法可以得到微震事件的核心區(qū)域,排除異常微震事件定位點對后續(xù)處理的影響；對于微震事件空間展布為復(fù)雜形態(tài)的情況,可結(jié)合其他信息,使用DBSCAN算法對微震事件進行分組處理,再基于這些分組數(shù)據(jù)使用MVEE算法,計算出各個分組對應(yīng)的橢球體,從而估算出復(fù)雜裂縫網(wǎng)絡(luò)對應(yīng)的壓裂改造體積。

1.3 重疊區(qū)域的處理

對于多段壓裂改造,每段壓裂形成的裂縫網(wǎng)絡(luò)往往存在重疊,對應(yīng)的各段壓裂形成的微震事件在空間展布也存在彼此重疊。利用橢球體模型對單段的壓裂改造范圍進行描述時,可以獲得該橢球體的中心點a和矩陣A；對于多段壓裂的處理,可以得到一系列的中心點(a1,a2,…，an)和矩陣(A1,A2,…，An)。由于各段壓裂形成的微震事件在空間中可能會彼此重疊,以微震事件空間展布特征為基礎(chǔ)求得的橢球體也會出現(xiàn)重疊。為消除橢球體重疊區(qū)域?qū)Χ喽螇毫颜w壓裂改造體積估算的影響,處理流程如下。

(1)對矩陣(A1,A2,…，An)進行處理,獲得每個矩陣的特征值矩陣Mi_λ和特征向量矩陣Mi_v(i=1,2,…，n)。由上述分析可知,矩陣A的特征值的平方根倒數(shù)是對應(yīng)橢球體的各個半軸長。將所有特征值轉(zhuǎn)化為半軸長,尋找到最大半軸長Lmax。

(2)以每個橢球體中心點為中心,以Lmax為半徑建立網(wǎng)絡(luò),并令Lstep為單個網(wǎng)格大小。判斷各個網(wǎng)格點g與所有橢球體之間的歸屬關(guān)系：

(5)

若該網(wǎng)格點位于橢球體的表面,經(jīng)過逆向處理后,正好位于單位球體的表面,即D到坐標原點的距離為1；若該網(wǎng)格點位于橢球體的內(nèi)部,經(jīng)過逆向處理后,也位于單位球體的內(nèi)部,即D到坐標原點的距離小于1；若該網(wǎng)格點位于橢球體的外部,經(jīng)過逆向處理后,也位于單位球體的外部,即D到坐標原點的距離大于1。將所有橢球體的相關(guān)參數(shù)均代入式(5)中,就可以判斷出該網(wǎng)格點與所有橢球體間的歸屬關(guān)系:該網(wǎng)格點可能只歸屬于當前橢球體,或歸屬于多個橢球體,或不歸屬于任何橢球體。

(3)當以第一個橢球體為被分析對象時,統(tǒng)計只歸屬于第一個橢球體的網(wǎng)格點；再對下一個橢球體進行分析,此前已被分析處理的橢球體不再參與計算,統(tǒng)計出只歸屬于當前橢球體的網(wǎng)格點數(shù)；依次對所有橢球體進行處理,對每步統(tǒng)計出的網(wǎng)格點數(shù)進行累加,并換算成體積。此時獲得數(shù)值就為消除重疊區(qū)域后的整體壓裂改造體積。

2 實例分析

2.1 多段簡單微震事件展布形態(tài)

實例數(shù)據(jù)來自重慶涪陵地區(qū)的頁巖氣田,為某頁巖氣儲層的水平井壓裂改造的微震監(jiān)測項目。由于頁巖氣大規(guī)模開發(fā)的需要,對該頁巖氣儲層實施了水平井多段水力壓裂處理。該水平井井軌跡為南北向,并采用臨井井中觀測方式對壓裂過程實施微震監(jiān)測。提取三個緊鄰壓裂段的微震監(jiān)測結(jié)果進行分析,各段微震事件數(shù)量分別為50、150和130,相應(yīng)空間分布圖像如圖1所示。在圖1中,X軸正向指向正東,Y軸正向指向正北。由各段微震事件的分布特征可見,微震事件的分布整體呈現(xiàn)東西向分布,與壓裂井井軌跡走向(南北向)近似正交。各段的微震事件分布形態(tài)相對簡單,但在各段微震事件之間存在較大程度的重疊。

首先,對這些微震事件進行分段處理,用 DBSCAN算法尋找到各段微震事件分布的核心區(qū)域,排除掉異常點。DBSCAN算法中的Eps和MinPts兩個參數(shù)分別被設(shè)置為50和3,其物理意義是簇內(nèi)兩個臨近微震事件間的最大距離是50 m和簇內(nèi)的至少應(yīng)包含3個微震事件。再由核心區(qū)域的微震事件結(jié)合MVEE算法,獲得最小體積封閉橢球體及其相關(guān)參數(shù),例如橢球體的中心,各個半軸長和橢球體的旋轉(zhuǎn)矩陣等。各段處理結(jié)果如圖2所示,以第2段為例,該段對應(yīng)的橢球體主軸方向與微震事件的分布特征相一致,橢球體積大小可近似描述微震事件所占據(jù)的范圍。將各段分析結(jié)果進行綜合分析(圖3),各個橢球體之間也存在明顯重疊現(xiàn)象,各段計算獲得的壓裂改造體積之和顯然不能作為整體的壓裂改造體積,需對重疊區(qū)域進行剔除。

分別以各個橢球體為中心,以橢球體的最大半軸為半徑建立網(wǎng)格,結(jié)合各個橢球體的參數(shù),判斷各個網(wǎng)格點的歸屬情況。逐個處理完所有橢球體后,就可獲得重疊區(qū)域大小和整體壓裂改造體積大小。如表1所示,各段壓裂改造體積分別為255×104、406×104和317×104m3,直接對各段壓裂改造體積求和為978×104m3,但消除重疊區(qū)域208×104m3,整體的壓裂改造體積為770×104m3。若對多段壓裂形成微震事件直接進行整體處理,計算獲得的橢球體會包含更多的區(qū)域,造成壓裂改造體積估值過大,并且橢球體主軸方向可能與壓裂裂縫分布形態(tài)不一致。

表1 各段壓裂改造體積及整體壓裂改造體積

2.2 單段復(fù)雜微震事件展布形態(tài)

由于多種因素的影響,單段壓裂可能會形成較為復(fù)雜的裂縫網(wǎng)絡(luò)結(jié)構(gòu)。與此相對應(yīng)的微震事件的空間分布形態(tài)就無法由單一的橢球體來加以描述。進而應(yīng)將復(fù)雜的分布形態(tài)分解為多個子區(qū)域進行處理。被分析的實測數(shù)據(jù)來另一頁巖氣儲層的壓裂改造項目,提取該項目的一個壓裂井段的微震監(jiān)測結(jié)果進行分析。該微震監(jiān)測結(jié)果包含297個微震事件,其空間分布如圖4所示。

圖2 各段壓裂的微震事件及其壓裂有效體積Fig.2 Microseismic event distributions and stimulated reservoir volume for each hydraulic fracturing stage

圖3 三段壓裂微震事件及各段壓裂有效體積Fig.3 Microseismic event distributions and stimulated reservoir volume of three stages

圖4 單段壓裂微震事件空間分布Fig.4 Microseismic event distributions in the single stage

基于微震事件的空間展布形態(tài),結(jié)合DBSCAN算法,Eps和MinPts兩個參數(shù)也被設(shè)置為50和3,對該段微震事件進行分組處理,共劃分為四組。對各組微震事件分別采用MVEE算法估算出各自對應(yīng)的最小體積封閉橢球體(圖5)。由圖5可知,各橢球體之間不存在重疊,整體壓裂改造體積可表示為各橢球體體積之和1 101×104m3(表2)。本文方法也可自動計算出該段整體壓裂體積為1 079×104m3,該結(jié)果與直接求和結(jié)果略有差異。產(chǎn)生該差異的主要原因是所選取的網(wǎng)格大小不同。對微震事件進行分組處理時,可結(jié)合其他信息來提高分組的合理性和準確性,例如微震事件的相對能量強弱、微震事件的發(fā)震時刻前后關(guān)系等。

圖5 單段壓裂微震事件及各組有效壓裂體積空間分布Fig.5 Microseismic event distributions and stimulated reservoir volume for each group in the single stage

第1組/m3第2組/m3第3組/m3第4組/m3整體m3593×10498×104363×10447×1041 079×104

3 結(jié)論

頁巖儲層的壓裂改造體積估算對評價壓裂改造效果和預(yù)測壓裂后的產(chǎn)能具有重要的意義。結(jié)合前人成果,針對兩種微震事件分布形態(tài)下的壓裂改造體積估算展開研究,得到以下結(jié)論。

(1)對于單段微震事件分布形態(tài)單一的情況,使用橢球體可近似描述壓裂裂縫的縫網(wǎng)形態(tài),并使用MVEE算法來計算出相應(yīng)橢球體的各項參數(shù),估算出相應(yīng)的壓裂改造體積。考慮到實際監(jiān)測數(shù)據(jù)會受到各種噪聲干擾,使用DBSCAN算法消除異常點的干擾。

(2)對于單段微震事件分布形態(tài)復(fù)雜的情況,首先結(jié)合其他信息對微震事件進行劃分,將復(fù)雜的分布特征拆分為多個分布特征相對簡單的分組；再對各分組采用MVEE算法估算出其壓裂改造體積,并消除各組對應(yīng)橢球體間存在的重疊區(qū)域；最終獲得該壓裂段的整體壓裂改造體積。

(3)選擇合適處理參數(shù)有助于對壓裂改造體積進行準確估計,例如聚類半徑和簇內(nèi)樣點數(shù)的選擇會影響聚類分析結(jié)果,進而影響最小體積封閉橢球的形態(tài)。因此,相關(guān)參數(shù)的選擇應(yīng)根據(jù)被分析數(shù)據(jù)的特征進行優(yōu)化選擇。