利用微地震事件重構三維縫網

劉 星 金 衍 林伯韜 向建華 鐘 華

(①中國石油大學(北京)石油工程學院，北京 102200；②中國石油西南油氣分公司，四川成都 610065)

0 引言

隨著水力壓裂技術在頁巖氣開發過程中的廣泛應用，微地震監測技術逐漸成為實時監測裂縫擴展和壓裂效果評價的關鍵技術[1-4]。但是，目前基于微地震的縫網表征僅能通過事件點推測縫網的宏觀參數(如縫網尺度和方位)，從而建立基于一定假設的二維裂縫模型，缺乏壓裂改造后三維縫網的重構建模方法。

在水力壓裂施工中，在儲層中原始裂縫和新生裂縫周圍產生不同程度的應力集中，致使總體的應變能增加，當外力進一步增加到臨界值時，原始裂縫的弱面中產生微觀的屈曲和變形，導致裂縫開始擴展，局部產生應力釋放和松弛，儲層中的能量以彈性波的形式釋放、傳播，在地層中產生微地震信號[5-6]。微地震信號以空間事件點的形式被解釋、記錄，而事件點的空間分布特征和裂縫面存在一定的幾何對應關系，這為壓裂體積縫網的精細反演和建模提供了數學基礎。

在基于微地震事件點的裂縫定量化表征方面，Cai等[7]首次應用微地震監測得到巖石破壞時的事件點密度，并提出了事件點密度和巖石損傷評價指標的關系模型。Sherilyn等[8]通過微地震監測計算的裂縫位置、尺寸、形狀等參數校核油藏數值模擬的離散裂縫模型，以提高歷史擬合精度。Maxwell等[9]根據微地震事件反映的裂縫位置信息，結合不同裂縫類型的變形機理，模擬了不同的二維平面幾何模型的復雜縫網擴展和變形規律。趙爭光等[10]通過分析微地震監測過程中事件點的分布和破裂能量，計算了水力裂縫的二維動態延伸和分布范圍。Yu等[11]首先對微地震數據進行矩張量分析得到了裂縫的產狀參數，進一步結合霍夫變換取得裂縫的形狀參數，建立了較為可靠的復雜離散裂縫網絡模型。楊瑞召等[12]基于微地震事件點的信號強度數據建立了基于能量層析成像的裂縫反演模型，該模型能得到二維水力裂縫帶分布。張云銀等[13]利用Delaunay三角剖分算法計算了微地震事件點的空間分布體積，提出了儲層壓裂改造體積估算方法，可以得到改造縫網體的空間體積。

然而，由于頁巖壓裂改造后體積縫網的幾何形態和空間分布的復雜性，并且監測結果易受環境噪聲影響，因此微地震監測事件點中存在噪點[14]。已有的研究大多從定性角度利用微地震事件點分析體積縫網的二維建模，或者利用微地震事件點提取個別裂縫參數以校核基于一定假設的隨機離散裂縫網絡，但缺乏基于微地震事件的穩健、直接的三維縫網重構方法。本文基于隨機模擬一致性(Random Sample Consensus，簡稱RANSAC)開發了一種穩健的三維縫網重構方法(RFM3D)。首先通過對比室內真三軸水力壓裂實驗結果建立了隨機多邊形的單裂縫幾何模型，并通過alpha-shape方法開發了基于微地震事件的單裂縫形狀識別算法；為了消除事件點中噪點對裂縫重構的影響，采用RANSAC算法提取、計算裂縫面產狀，在得到單裂縫的全部信息后通過單裂縫疊加得到RFM3D縫網。為了驗證算法的穩健性，利用蒙特卡羅隨機模擬方法生成已知縫網，通過對縫網離散生成模擬微地震事件點，在模擬事件點中加入一定比例的噪點后進行RFM3D，以驗證算法的穩健性。模擬結果發現，RFM3D算法可在一定程度克服噪點影響，重構縫網與原始縫網符合度較高，具有較高的穩健性。

1 RFM3D原理

1.1 單裂縫幾何模型及識別算法

構建離散裂縫網絡模型的基礎是建立合理的單裂縫幾何模型，目前常用的單裂縫幾何模型主要有Bacher圓盤模型[15]、改進Bacher圓盤模型[16]、Possion圓盤模型[17]、隨機多邊形模型[18]等。由于Bacher圓盤模型、Possion圓盤模型與水力壓裂中的雙翼縫模型較一致，且數值模擬較簡便，因此這兩種模型得到廣泛應用。通過大型真三軸室內壓裂實驗發現，水力壓裂縫網中裂縫的幾何形狀符合隨機多邊形模型(圖1)，并且隨機多邊形模型具有表征準確、模型簡單、空間變換方便等優點[19]，故本文采用隨機多邊形模型進行縫網重構。

圖1 水力壓裂物模實驗及縫網描繪

圖2 隨機多邊形單裂縫幾何模型

在選擇了隨機多邊形模型(圖2)之后，基于裂縫產狀三要素(走向、傾向、傾角)構建單裂縫幾何模型。假設裂縫的走向為α1、傾角為α2、傾向為α3，裂縫面在三個坐標軸上的截距分別為a、b、c，通過幾何關系計算各個產狀參數。由截距關系得到裂縫的幾何方程為

(1)

裂縫走向為

(2)

其中α1∈[0,π]。同理得到裂縫的傾角為

(3)

(4)

其中α3∈[0,2π]。

在得到裂縫的產狀參數之后，根據微地震事件點結合alpha-shape方法開發了三維隨機多邊形裂縫識別算法[20](圖3)，alpha-shape方法是一種Delaunay三角剖分算法。對于一個二維平面上的點集P來說，alpha-shape通過遍歷P中的任意兩點p和q，當且僅當存在一個以p和q為弦、內部沒有其他點的alpha-shape圓盤時，認定p和q為P的二維形狀的頂點；在找到所有的頂點后，依次連接所有的頂點形成一個多邊形。

圖3 alpha-shape形狀識別算法示意圖

1.2 基于RANSAC的裂縫產狀識別方法

在水力壓裂施工中，除了周圍的噪聲源產生的噪點之外，在水力壓裂時也會有一部分天然裂縫在誘導應力場的作用下被“激發”，壓裂結束后這些處于“激發”狀態的裂縫也隨即閉合，該部分響應產生的事件點往往距井筒較遠，這些額外的事件點和環境噪聲形成的事件點形成了微地震監測信號中的“噪點”。對微地震事件點進行RFM3D的首要問題是如何克服噪點的影響、準確地識別裂縫的形狀和產狀，因此需要選擇穩健的回歸方法識別、提取裂縫信息。

RANSAC由Fisvhler 等[21]首次提出，主要為了解決被部分噪點污染的數據集的最小二乘估計失真問題。該算法基于隨機抽樣改進傳統的最小二乘法，大大增強了對有效點集的識別，在噪點比例超過10%時仍然能識別有效點集，是一種非常穩健的擬合算法，因此在圖像處理領域得到廣泛應用。對于單裂縫產狀的識別，RANSAC首先將所處理的點集分為噪點和有效點，通過隨機抽樣挑選一個隨機樣本；然后使用人為設定的閾值將樣本的噪點率控制在合理范圍內；最后將選擇的樣本代入給定的裂縫產狀幾何方程(式(1))，通過回歸得到幾何模型的對應產狀參數(圖4)。

1.3 RFM3D

基于微地震事件得到的位置信息，綜合裂縫面產狀識別和裂縫形狀識別開發了RFM3D(圖5)。

結合圖4，假設某個含有一定比例噪點的微地震事件點集合為S，其目標縫網中含有m條裂縫，裂縫幾何模型為f(x;α)，包含產狀和形狀參數向量α1,…,αm，RFM3D迭代計算步驟如下。

圖4 RANSAC算法識別單裂縫產狀示意圖

圖5 基于微地震數據的RFM3D流程

(1)從微地震事件集合S中隨機抽取一個大小為n的樣本，并將樣本進行裂縫產狀模型最小二乘回歸得到模型產狀參數αtest。

(2)設定距離閾值t，計算剩余所有微地震事件點到該裂縫面的距離，由統計得到在閾值范圍內的點集Stest?S，當Stest的大小滿足要求時，認為選定的點集符合裂縫的擬合條件。

(3)重復步驟(1)，再次隨機選擇一個大小為n的樣本，得到閾值范圍內的標識點集Sin和產狀參數標識向量αc，如果|Sin|<|Stest|，令Sin=Stest，則αc=αtest。

(4)設定點集閾值為T，重復步驟(1)～(3)，直到|Sin|≥T為止，經統計得到總迭代次數N，進而得到首個單裂縫產狀參數向量α11=αc，并將Sin代入alpha-shape形狀識別算法，得到裂縫的形狀參數向量α12，通過組合得到首個裂縫幾何模型的全部產狀和形狀參數向量α1。

(5)從原始事件點集S中去除Sin得到剩余微地震事件點集Sres，對Sres重復步驟(1)～(4)得到第二個裂縫模型的全部參數向量α2。

(6)重復步驟(1)～步驟(5)m次，直止得到所有裂縫的幾何參數向量α1,…,αm，算法運行結束。

在隨機抽樣過程中，如果p表示在N次抽樣中至少有一次抽到完全沒有噪點的樣本的概率，那么可以得到迭代次數N

(5)

式中：p為成功抽到完全為有效點樣本的概率；ω為有效點占所有事件點的比例；n為每次隨機抽樣的樣本大小。

式(5)表明，隨著ω的減小、n的增大，迭代次數隨之增加。但是在實際微地震監測中一般很難得到確切的噪點比例(噪點個數與有效點的比值)，因此也無法得到準確的迭代次數N，可以采取估計的方式近似用每次得到的標識點集的大小和事件點集大小的比值作為每一次迭代過程的ω，根據每次的計算結果調整總的循環迭代的上限N，可以在大規模縫網重構中有效減少計算時間，提高計算效率。

2 模型穩健性定量分析

2.1 模型穩健性模擬分析

為了驗證RFM3D算法的穩健性，采用蒙特卡羅方法生成裂縫數量為5、10、15、20的四種模擬縫網，采用隨機離散的方法將模擬縫網離散為模擬微地震事件點集，加入一定比例的、由均勻隨機過程生成的噪點，通過對模擬離散微地震事件點進行RFM3D得到新縫網，通過對比原始模擬縫網和重構縫網分析算法的穩健性。

圖6為不同裂縫數目m的縫網模擬結果。由圖可見，RFM3D能很好地克服噪點影響，重構縫網和模擬縫網的相似度較高，算法具有較好的穩健性。

圖7為裂縫數目為10的不同噪點比例的RFM3D結果。由圖可見，在保持縫網不變的條件下，對原始縫網離散化處理后(圖7a)，分別添加5%(圖7b)、10%(圖7c)、15%(圖7d)、20%(圖7e)、25%(圖7f)的噪點比例進行RFM3D。模擬結果表明，隨著噪點比例增加，RFM3D縫網(圖7b～圖7f)與原始縫網(圖7a)相似性降低，重構效果逐漸變差，并且當噪點比例大于10%時(圖7d～圖7f)，RFM3D縫網的不確定性明顯增加，當噪點比例不大于10%時(圖7a～圖7c)，RFM3D效果較好，表明RFM3D算法可以克服約10%的噪點影響，該結果在后續的模擬中會得到進一步證實。

圖6 不同裂縫數目m的縫網模擬結果(噪點比例均為10%)

圖7 裂縫數目為10的不同噪點比例的RFM3D結果

2.2 模型魯棒性分析方法

為了定量分析RFM3D算法對噪點的穩健性(圖6a)，首先使用蒙特卡羅方法在1m×1m×1m的立方體中隨機生成模擬縫網，結合測繪學空間形狀相似性理論[22]，提出了基于單裂縫的平均距離的相似性指標ADI(Average Distance Index)。

如圖8所示，以多邊形單裂縫ABCD和單裂縫A′B′C′D′為例，采用“五點法”計算ADI，其方法是尋找兩個縫網中的裂縫對應關系，然后通過計算五點對應距離的平均值作為兩個裂縫的ADI(圖6)，最后求取模擬縫網和重構縫網所有裂縫的ADI均值作為縫網ADI，其計算步驟如下：

(1)假設模擬縫網和重構縫網的個數均為m，首先確定兩種縫網的單裂縫對應關系，形成m個裂縫對。

(2)分別計算第i個裂縫對的頂點平均距離Lfi和裂縫中心點的距離Dfi。

圖8 多邊形裂縫相似性度量示意圖

(3)重復步驟(1)～(2)遍歷所有裂縫對，得到重構縫網和模擬縫網的平均單裂縫ADI

(6)

由式(6)可知，在一定的噪點比例下，ADI的范圍為0～1，其中縫網ADI越大，重構縫網和模擬縫網間的差異就越大，RFM3D的重構效果越差；反之，縫網ADI越小，重構縫網和模擬縫網越接近，RFM3D的效果越好，其中ADI=0表示重構縫網和原始縫網完全重合。根據ADI的不同取值范圍對應的噪點比例定義為第一、第二和第三臨界點，并且以此劃分了四個評價區間，據此給出了四種重構效果(較好、一般、較差、最差)(表1)。通過設定不同的噪點比例和裂縫數目進行模擬，以縫網ADI所在區間定量評價算法的重構效果。

表1 ADI評價區間劃分結果

2.3 縫網相似性模擬分析

在定量評價RFM3D算法的穩健性過程中，其核心是通過蒙特卡羅生成模擬縫網樣本，計算不同噪點、不同裂縫數目的ADI。圖9為RFM3D算法噪點比例—ADI曲線。由圖可見：①隨著噪點比例增加，ADI基本上落在0～0.5范圍，雖然由于隨機抽樣的原因致使計算結果略有波動，但ADI總體呈非線性增長，表明RFM3D效果越來越差，噪點對重構結果的影響也越來越大，重構縫網的不確定性增加。②即使模擬縫網中裂縫數目不同，ADI值在噪點比例大于200%時達到最大并保持穩定，表明重構縫網和原始縫網的相似性較差，說明當原始微地震事件點集中時若噪點比例高于200%，有效點集完全被噪點淹沒，重構縫網也完全失真，不具有參考價值，即將200%作為RFM3D算法的噪點比例上限。③對噪點比例—ADI曲線進行非線性回歸發現，ADI隨噪點比例滿足logistic增長模式(縫網ADI先增加后減小并趨于零，中間存在增速極值點，并且極值點位于第二臨界點附近)。因此對于一個特定微地震事件點的降噪處理來說，需要將噪點比例至少降低到第二臨界點以下，重構縫網形態才會逐漸趨于穩定，若將噪點比例降低到第一臨界點以下才能得到較為可靠的重構結果。④噪點比例—ADI曲線隨著裂縫數目的增加逐漸變得陡峭，ADI增長速度明顯加快，并很快趨于穩定，說明在同等條件下裂縫數目越多，噪點對算法的影響越大，需要將噪點比例降至更低才能得到可靠的重構結果。

圖10為m—噪點比例曲線。由圖可見：①第一臨界點隨著裂縫數目m增加呈下降趨勢，在m=1以后趨于穩定，并且噪點比例穩定在約10%，即使裂縫數目進一步增加，第一臨界點也幾乎保持不變。由于當噪點比例小于第一臨界點時RFM3D結果是穩定、可靠的，因此RFM3D算法的噪點比例下限為10%，說明在一般條件下對不同裂縫數目的縫網重構來說，RFM3D算法至少可克服10%的噪點影響，這和前文的不同噪點比例的RFM3D結果一致(圖7)。②隨著裂縫數目增加，三個臨界點的值也逐漸減低，表明裂縫數目越多，RFM3D重構效果越差，當m=9時，第三臨界點降低到50%以下，說明在使用RFM3D算法重構含有較多裂縫的目標縫網時需將噪點比例降至更低，否則無法得到可靠的結果。

圖9 RFM3D算法噪點比例—ADI曲線

設定噪點比例范圍為0～200%(步長為2%)，縫網中裂縫數目分別為3、5、7、9、10(步長為2)，編寫MATLAB程序對每種縫網進行100次RFM3D，分別計算重構縫網和模擬縫網的ADI

圖10 m—噪點比例曲線

2.4 模型穩健性定量分析

綜上所述：RFM3D算法是一種較穩健的縫網重構方法，且ADI隨噪點比例呈logistic曲線增長；在一般條件下，RFM3D算法的噪點比例下、上限分別為10%、200%；隨著裂縫數目增加，RFM3D算法受噪點的影響越來越嚴重。因此，在使用RFM3D算法重構大規模的縫網時必需對微地震事件點進行降噪處理，使噪點比例控制在10%以下，才能得到穩定、可靠的縫網重構結果。

3 工程應用分析

3.1 工程地質背景

YY1井位于川中隆起區的川西南低陡褶帶，鉆遇志留系龍馬溪組頁巖儲層，平均孔隙度為5.28%，基質滲透率均值分布范圍為(4～5)×10-5mD，含氣孔隙度平均值為2%，含氣飽和度平均值約為50%，總體上儲層物性較好。井下成像測井和巖心資料表明，YY1井中的天然裂縫主要為構造裂縫以及超壓填充裂縫，其中構造裂縫形成的張開縫多為高角度垂直縫，超壓裂縫多為微細裂縫，空間呈網狀分布，延伸距離較小，大多數被礦物充填(圖11a)；該區最大水平主應力方位約為135°(圖11b)，YY1井所在地區天然裂縫總體發育方向(圖11c)與該區最大主應力方向夾角較大，在水力壓裂過程中地應力使高角度天然裂縫閉合，因此不利于天然裂縫的張開[23]。由YY1井三個壓裂段微地震監測結果(圖12 )可見：①由于地應力的屏蔽作用，未出現遠井天然裂縫的激發響應產生的微地震事件點，均勻分布的事件點有利于RFM3D，在一定程度上降低了的影響，事件點沿井分布較為集中，未出現高導流的長段水力裂縫，并且在3個施工段分布相對均勻；④對于每一級壓裂來說，平均縫長和平均縫寬約為200m。為了評估YY1井壓裂效果，對YY1井的微地震事件點進行RFM3D，并結合地應力和天然裂縫分析結果對RFM3D結果進行對比、驗證。

圖11 YY1井地應力方位和裂縫方位[24]

圖12 YY1井三個壓裂段微地震監測結果

由孤立事件點引起的誤差；②YY1井微地震監測為3級接收，現場實時處理并定位的有效微地震事件點共1800個；③從微地震事件震級及分布特征看，事件在壓裂施工全程均有發生，受地應力屏蔽作用

3.2 RFM3D結果分析

由微地震監測結果(圖12)和RFM3D結果(圖13)可知，水力壓裂形成了一定規模的體積縫網；由裂縫與微地震事件對比圖(圖13b)可見，重構的三維縫網和原始微地震事件貼合度較高，每個壓裂段重構的裂縫延伸范圍僅限于該壓裂段對應的微地震信號的展布范圍。由于水平最大地應力方向和天然裂縫的夾角較大，因此天然裂縫對水力裂縫的延伸具屏蔽作用，造成水力裂縫的延伸范圍僅限于井筒周圍，未出現遠井區域的孤立微地震事件點，同時也未出現個別的高導流水力裂縫，整體上微地震事件和壓裂裂縫沿井分布較為均勻。

圖14為YY1井三維縫網產狀分布。由圖可見：

圖13 YY1井RFM3D結果(a)及裂縫與微地震事件對比圖(b)

圖14 YY1井三維縫網產狀分布

水力裂縫的優勢走向范圍為120°～150°，平均值為135°，垂直井筒并與水平最大地應力方向一致，RFM3D結果和由地應力分析得到的水力裂縫擴展方向一致；裂縫的傾向集中在240°，與最大主應力方向垂直；裂縫的優勢傾角范圍為60°～90°，表明裂縫經壓裂改造后主要為高角度垂直縫，與成像測井分析結果一致，說明井筒周圍的天然裂縫對水力裂縫的起裂和擴展具重要影響。總體來看，重構后的縫網和水力壓裂實驗中的裂縫擴展規律一致[25]，驗證了RFM3D算法的可行性。

4 結束語

本文基于隨機模擬一致性和alpha-shape方法開發了頁巖壓裂體積縫網重構算法(RFM3D)，通過蒙特卡羅縫網模擬定量分析了算法對噪點的穩健性，并用于川中地區YY1井的微地震縫網重構，得到以下認識：

(1)通過室內壓裂實驗發現，可用隨機多邊形模型描述實際裂縫形狀。

(2)RFM3D算法易于匹配復雜的裂縫幾何模型，在一般條件下該算法至少能克服10%的噪點干擾，可較準確地重構壓裂縫網的幾何形態，算法具有較好的穩健性。

(3)隨著噪點比例增大，重構相似性指標(ADI)隨噪點比例滿足logistic增長模式；隨著縫網中裂縫數目的增多，ADI臨界點也隨之降低。因此，在重構大規模體積縫網時噪點比例應嚴格控制在10%以下才能得到穩定、可靠的結果。

對于頁巖壓裂體積縫網重構算法來說，裂縫數目和噪點比例是影響縫網重構精度的主要因素。本文從微地震點和裂縫面的幾何關系出發，提出了一種較為可行的重構方法，該研究也可用于基于縫網的頁巖產能預測和儲層改造體積的計算。