王旭林, 桂志先, 王 鵬*, 易娟子
(1.長江大學非常規油氣湖北省協同創新中心, 油氣資源與勘探技術教育部重點實驗室,武漢 430100;2.中國石油天然氣股份有限公司西南油氣田分公司重慶氣礦,重慶 400021)
頁巖油氣作為一種新型能源,因其資源十分豐富,越來越被人們所重視。頁巖油氣儲層具有低孔隙度、低連通性、低滲透率的物性特性,通常需要采取增產措施來改善其巖石物理特性[1]。由于頁巖油氣藏的特殊性,在勘探開發中影響頁巖油氣井產能的因素有很多。研究表明,儲層壓裂改造體積是影響頁巖油氣藏生產效果的關鍵因素之一[2-3]。前人研究發現被壓裂改造油氣藏的生產力與被壓裂改造體積正相關,改造的體積越大,其相應的增產效果越明顯[4-5]。因此對壓裂改造體積(stimulated reservior volume,SRV)的準確估算已經成為頁巖壓裂研究領域的重難點問題。前人分別采用通道長度與通道寬度來表征裂縫擴展的長度和寬度,水平井長度表征SRV[6]。上述理論模型方法均采用立方體來擬合SRV,但包括較多的無效范圍,所以擬合效果較為粗略[7]。溫慶志等[8]結合壓裂縫網絡形態,通過研究均質地下儲層中主裂縫周邊的應力場發現,地應力場方向改變的區域形成縫網,縫網區域大致為橢球體,因此可以近似地認為壓裂形成的縫網為橢球體。同時,微震事件的空間展布特征也被用來估算SRV[9]。因此,可將壓裂形成的微震事件的空間展布近似處理為橢球體,來估算儲層壓裂改造體積。邵媛媛等[7]基于地震事件震源點將壓裂改造體積擬合成橢球體,建立了儲層改造體積最小橢球模型,通過對比其他方法得出最小體積橢球模型獲得了更為詳細的SRV幾何結構。但上述改進方法只能得出儲層改造體積的近似形狀與各種參數對改造體積估算的影響,尚未給出比較完整的橢球體積估算方法。在實際生產過程中,多段壓裂施工作業會出現各段壓裂裂縫網絡存在交錯疊合的情況,并且由于多種因素的影響,單段壓裂裂縫可能會呈現出較為復雜的網絡形態,即表現為微震事件的空間展布形態復雜。因此,使用單一的模型很難解決上述復雜的壓裂裂縫網絡的壓裂改造體積的估算。
本文針對實際壓裂改造過程中存在的多段壓裂改造體積存在重疊和單段壓裂形成的復雜裂縫網絡兩類情況進行分析,估算對應的壓裂改造體積。對多段壓裂體積重疊的情況,首先利用DBSCAN(density based spatial clustering of application with noise)算法對各段微震事件分別進行處理,消除異常的微震事件,再利用MVEE算法估算出各段微震事件對應的最小體積封閉橢球體,并獲得橢球體的相關參數(橢球體中心點、特征值和特征向量等),然后對壓裂作業區進行網格化,結合各個橢球體參數,使用仿射變換判斷各個網格點與每個橢球體的歸屬關系,進而統計出各段壓裂改造體積間的重疊部分;針對單段壓裂形成的復雜裂縫網絡形態的情況,采用DBSCAN算法對微震事件進行分組,將復雜的微震事件展布轉化為多個展布相對簡單的分組,以各個分組的微震事件為基礎,計算對應的最小體積橢球體,并消除各個橢球體間的重疊部分,進而估算出單段壓裂的整體改造體積。本文還利用兩組實際微震監測數據,對上述兩類情況進行分析與估算,以驗證本文方法的可行性與準確性。
最小體積封閉橢球算法(MVEE)[10]是通過尋找一個n維最小封閉橢球以包含給定的所有點集合,并獲取相應橢球的信息,如橢球軸半徑、橢球中心、橢球體積等。MVEE問題在數據挖掘、模式識別等領域應用十分廣泛。
任意一個橢球可表示為
E(a,A)={x∈R3|(x-a)TA(x-a)≤1}
(1)
式(1)中,a∈R3為該橢球體的中心點;A∈R3×3為對稱正定矩陣,該矩陣實現橢球體向單位球體的仿射變換,矩陣A的特征值(λ1、λ2、λ3)平方根倒數分別為橢球體半軸長,即實現仿射變換的縮放處理;特征向量實現仿射變換的旋轉處理。式(1)也可表示為
E=Mtranslation×Mrotation×Mscaling×C
(2)
式(2)中,Mtranslation、Mrotation和Mscaling分別為平移、旋轉和縮放三種仿射變換;C為單位球體;E表示最終的橢球體。Mrotation與矩陣A的特征向量矩陣互為逆矩陣;Mtranslation由中心點a控制。
以微震事件的展布特征為基礎獲得的最小體積封閉橢球體的體積可表示為

(3)
求解最小體積封閉橢球問題轉可化為最優化問題:
minA[-(ΔA)1/2], (x-a)TA(x-a)≤1
(4)
由于-(ΔA)1/2不是一個凸函數,所以式(4)的優化模型不是一個凸優化模型,不具備凸優化的性質。利用KY算法[11]對之進行凸優化處理,并最終實現對最小體積封閉橢球的求解,獲得橢球體的中心點a與矩陣A。
實際微震監測數據處理會受到各種因素的影響,例如速構建度模型的誤差、微震事件初至信息無法準確拾取等,這些因素都會影響到最終的微震事件定位精度。同時,由于頁巖儲層的非均質性、天然裂縫的發育程度等因素的影響,壓裂形成的裂縫結構更為復雜。因此,需要對微震事件進行適當的預處理,消除定位異常的影響,并對復雜壓裂裂縫的分組處理。
DBSCAN是一種基于密度的空間聚類分析算法,能將具有足夠密度的區域劃分為簇,在具有噪聲的空間數據中發現任意形狀的簇[12]。該算法被廣泛應用于各個領域的數據處理與分析,例如提高雷電定位算法的抗干擾能力[13]、數據挖掘[14]、分析電力工程數據的完整性[15]和巖石物理實驗數據的處理[16]。其中,Eps和MinPts是整個算法最為關鍵的參數,其含義分別為單簇內樣點間的最大距離(又稱聚類半徑)和單簇包含的最小樣點數。并由此引申出核心點、直接密度可達、密度可達和密度相連等相關概念。若某個樣點的Eps臨域內,臨近樣點數大于或等于Minpts,則該樣點就是核心點;若某樣點到達核心點的距離小于Eps,就說明從該樣點到核心點直接密度可達;對于一個樣點序列,樣點間依次滿足直接密度可達,即直接密度可達滿足傳遞性,則這些樣點是密度可達;若兩個樣點對于同一樣點密度可達,就認為這兩個樣點是密度相連[17]。
DBSCAN的處理流程可描述為:①在整個樣點數據集中隨機挑選出一個樣點,并確保該樣點為核心點;②以該樣點為起點,尋找到所有密度相連的數據點,并將尋找結果建立新簇;③重新掃描樣點數據集(不包括之前已有簇中的任何樣點)尋找新的核心點,再重復步驟②,直到沒有新的核心點為止[18]。原始樣點數據集中沒有被包含在任何簇中的樣點就構成異常點。
對于微震事件空間展布為簡單形態的情況,使用DBSCAN算法可以得到微震事件的核心區域,排除異常微震事件定位點對后續處理的影響;對于微震事件空間展布為復雜形態的情況,可結合其他信息,使用DBSCAN算法對微震事件進行分組處理,再基于這些分組數據使用MVEE算法,計算出各個分組對應的橢球體,從而估算出復雜裂縫網絡對應的壓裂改造體積。
對于多段壓裂改造,每段壓裂形成的裂縫網絡往往存在重疊,對應的各段壓裂形成的微震事件在空間展布也存在彼此重疊。利用橢球體模型對單段的壓裂改造范圍進行描述時,可以獲得該橢球體的中心點a和矩陣A;對于多段壓裂的處理,可以得到一系列的中心點(a1,a2,…,an)和矩陣(A1,A2,…,An)。由于各段壓裂形成的微震事件在空間中可能會彼此重疊,以微震事件空間展布特征為基礎求得的橢球體也會出現重疊。為消除橢球體重疊區域對多段壓裂整體壓裂改造體積估算的影響,處理流程如下。
(1)對矩陣(A1,A2,…,An)進行處理,獲得每個矩陣的特征值矩陣Mi_λ和特征向量矩陣Mi_v(i=1,2,…,n)。由上述分析可知,矩陣A的特征值的平方根倒數是對應橢球體的各個半軸長。將所有特征值轉化為半軸長,尋找到最大半軸長Lmax。
(2)以每個橢球體中心點為中心,以Lmax為半徑建立網絡,并令Lstep為單個網格大小。判斷各個網格點g與所有橢球體之間的歸屬關系:

(5)
若該網格點位于橢球體的表面,經過逆向處理后,正好位于單位球體的表面,即D到坐標原點的距離為1;若該網格點位于橢球體的內部,經過逆向處理后,也位于單位球體的內部,即D到坐標原點的距離小于1;若該網格點位于橢球體的外部,經過逆向處理后,也位于單位球體的外部,即D到坐標原點的距離大于1。將所有橢球體的相關參數均代入式(5)中,就可以判斷出該網格點與所有橢球體間的歸屬關系:該網格點可能只歸屬于當前橢球體,或歸屬于多個橢球體,或不歸屬于任何橢球體。
(3)當以第一個橢球體為被分析對象時,統計只歸屬于第一個橢球體的網格點;再對下一個橢球體進行分析,此前已被分析處理的橢球體不再參與計算,統計出只歸屬于當前橢球體的網格點數;依次對所有橢球體進行處理,對每步統計出的網格點數進行累加,并換算成體積。此時獲得數值就為消除重疊區域后的整體壓裂改造體積。
實例數據來自重慶涪陵地區的頁巖氣田,為某頁巖氣儲層的水平井壓裂改造的微震監測項目。由于頁巖氣大規模開發的需要,對該頁巖氣儲層實施了水平井多段水力壓裂處理。該水平井井軌跡為南北向,并采用臨井井中觀測方式對壓裂過程實施微震監測。提取三個緊鄰壓裂段的微震監測結果進行分析,各段微震事件數量分別為50、150和130,相應空間分布圖像如圖1所示。在圖1中,X軸正向指向正東,Y軸正向指向正北。由各段微震事件的分布特征可見,微震事件的分布整體呈現東西向分布,與壓裂井井軌跡走向(南北向)近似正交。各段的微震事件分布形態相對簡單,但在各段微震事件之間存在較大程度的重疊。
首先,對這些微震事件進行分段處理,用 DBSCAN算法尋找到各段微震事件分布的核心區域,排除掉異常點。DBSCAN算法中的Eps和MinPts兩個參數分別被設置為50和3,其物理意義是簇內兩個臨近微震事件間的最大距離是50 m和簇內的至少應包含3個微震事件。再由核心區域的微震事件結合MVEE算法,獲得最小體積封閉橢球體及其相關參數,例如橢球體的中心,各個半軸長和橢球體的旋轉矩陣等。各段處理結果如圖2所示,以第2段為例,該段對應的橢球體主軸方向與微震事件的分布特征相一致,橢球體積大小可近似描述微震事件所占據的范圍。將各段分析結果進行綜合分析(圖3),各個橢球體之間也存在明顯重疊現象,各段計算獲得的壓裂改造體積之和顯然不能作為整體的壓裂改造體積,需對重疊區域進行剔除。
分別以各個橢球體為中心,以橢球體的最大半軸為半徑建立網格,結合各個橢球體的參數,判斷各個網格點的歸屬情況。逐個處理完所有橢球體后,就可獲得重疊區域大小和整體壓裂改造體積大小。如表1所示,各段壓裂改造體積分別為255×104、406×104和317×104m3,直接對各段壓裂改造體積求和為978×104m3,但消除重疊區域208×104m3,整體的壓裂改造體積為770×104m3。若對多段壓裂形成微震事件直接進行整體處理,計算獲得的橢球體會包含更多的區域,造成壓裂改造體積估值過大,并且橢球體主軸方向可能與壓裂裂縫分布形態不一致。

表1 各段壓裂改造體積及整體壓裂改造體積
由于多種因素的影響,單段壓裂可能會形成較為復雜的裂縫網絡結構。與此相對應的微震事件的空間分布形態就無法由單一的橢球體來加以描述。進而應將復雜的分布形態分解為多個子區域進行處理。被分析的實測數據來另一頁巖氣儲層的壓裂改造項目,提取該項目的一個壓裂井段的微震監測結果進行分析。該微震監測結果包含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
基于微震事件的空間展布形態,結合DBSCAN算法,Eps和MinPts兩個參數也被設置為50和3,對該段微震事件進行分組處理,共劃分為四組。對各組微震事件分別采用MVEE算法估算出各自對應的最小體積封閉橢球體(圖5)。由圖5可知,各橢球體之間不存在重疊,整體壓裂改造體積可表示為各橢球體體積之和1 101×104m3(表2)。本文方法也可自動計算出該段整體壓裂體積為1 079×104m3,該結果與直接求和結果略有差異。產生該差異的主要原因是所選取的網格大小不同。對微震事件進行分組處理時,可結合其他信息來提高分組的合理性和準確性,例如微震事件的相對能量強弱、微震事件的發震時刻前后關系等。

圖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
頁巖儲層的壓裂改造體積估算對評價壓裂改造效果和預測壓裂后的產能具有重要的意義。結合前人成果,針對兩種微震事件分布形態下的壓裂改造體積估算展開研究,得到以下結論。
(1)對于單段微震事件分布形態單一的情況,使用橢球體可近似描述壓裂裂縫的縫網形態,并使用MVEE算法來計算出相應橢球體的各項參數,估算出相應的壓裂改造體積。考慮到實際監測數據會受到各種噪聲干擾,使用DBSCAN算法消除異常點的干擾。
(2)對于單段微震事件分布形態復雜的情況,首先結合其他信息對微震事件進行劃分,將復雜的分布特征拆分為多個分布特征相對簡單的分組;再對各分組采用MVEE算法估算出其壓裂改造體積,并消除各組對應橢球體間存在的重疊區域;最終獲得該壓裂段的整體壓裂改造體積。
(3)選擇合適處理參數有助于對壓裂改造體積進行準確估計,例如聚類半徑和簇內樣點數的選擇會影響聚類分析結果,進而影響最小體積封閉橢球的形態。因此,相關參數的選擇應根據被分析數據的特征進行優化選擇。