基于非監(jiān)督學習的水力壓裂裂縫平面識別

馬孜卓 鄭憶康 薛清峰 翟鴻宇 雷興林

1)中國科學院地質(zhì)與地球物理研究所，中國科學院油氣資源研究院重點實驗室，北京 100029 2)中國科學院地球科學研究院，北京 100029 3)中國地震局地球物理研究所，北京 100081 4)日本產(chǎn)業(yè)技術綜合研究所，筑波 3058567

0 引言

近年來，隨著水力壓裂技術在頁巖氣開采中的廣泛應用，對于裂縫網(wǎng)的研究變得越來越重要。由于頁巖氣儲層具有低孔隙度和低滲透性的物性特征，因此只有極少數(shù)天然裂縫發(fā)育明顯的頁巖氣區(qū)塊才可以直接投入生產(chǎn)。大多數(shù)頁巖氣區(qū)塊需要通過水力壓裂施工達到預期產(chǎn)量(李新景等，2007；謝春輝等，2015；仝少凱等，2019)。通過裂縫網(wǎng)絡監(jiān)測技術，施工人員可以在施工的過程中針對壓裂效果進行有效評估，從而促進開采的進行。裂縫網(wǎng)絡監(jiān)測技術具有以下優(yōu)點：①有助于了解地下裂縫結構，確定裂縫大小和幾何特征，并能通過監(jiān)測實現(xiàn)壓裂過程新生裂縫與原生裂縫的有效區(qū)分；②有助于了解水力壓裂施工后的壓裂效果，確定裂縫是否與生產(chǎn)層相交并分析裂縫與天然裂縫是否相交；③可以采用裂縫優(yōu)化和生產(chǎn)經(jīng)濟評價來優(yōu)化壓裂設計，增加施工規(guī)模。因此，準確的裂縫監(jiān)測技術是促進頁巖氣開采的有效手段。

Ding等(2011)通過地質(zhì)方法、測井方法、地震方法、鉆孔方法和機械干擾應力測試，實現(xiàn)了對裂縫的有效識別。目前，最成熟的方法是利用測井數(shù)據(jù)來識別和評估裂縫(趙軍龍等，2012；賴錦等，2015；杜貴超等，2016；李麗慧等，2019；董春梅等，2020)。由于其有效性，該方法已被廣泛應用于石油和天然氣勘探(Yu et al，2017)。然而，由于井的覆蓋范圍有限，通過有限的井資料，無法獲得整個工作區(qū)完整的裂縫發(fā)育情況，因此基于測井數(shù)據(jù)的裂縫識別方法仍然具有局限性。隨著近年來微地震監(jiān)測技術的迅速發(fā)展，大規(guī)模高密度的微震監(jiān)測已成為水力壓裂施工的標準配置。如果可以根據(jù)微地震信號直接提取地下裂縫信息，不僅可以節(jié)省開發(fā)成本，而且能夠全面地了解裂縫的發(fā)展情況。微地震監(jiān)測技術能夠通過對水力壓裂過程中產(chǎn)生的微地震事件進行觀測和分析，從而監(jiān)測壓裂過程，其理論基于聲發(fā)射和地震學。施工人員在監(jiān)測區(qū)域布置多個探測器，在巖石破裂的過程中實時收集地震數(shù)據(jù)，這些監(jiān)測結果對于制定非常規(guī)油氣藏開發(fā)方案、提高產(chǎn)量和評價壓裂效果具有重要意義。傳統(tǒng)的微地震監(jiān)測通常側重于提取裂縫的破裂區(qū)信息(Eshiet et al，2013)，而利用微震數(shù)據(jù)直接提取裂縫信息的研究開展不多(Xue et al，2018)。

本文通過引入基于時空屬性的非監(jiān)督學習數(shù)據(jù)挖掘分析方法，分析水力壓裂過程中產(chǎn)生的微地震數(shù)據(jù)。結合微地震有效事件的時空屬性，進一步分析在水力壓裂中產(chǎn)生的裂縫幾何形態(tài)，確定裂縫的破裂平面信息，并采用水力壓裂巖石物理實驗手段和數(shù)據(jù)對該方法進行可行性驗證。相對于常規(guī)的地球物理觀測手段和反演方法，水力壓裂實驗提供了一個在巖樣尺度下進行裂縫平面分析研究的實驗平臺，其優(yōu)勢主要包括：溫度和應力加載條件可控，便于設計各種實驗，研究不同科學問題，可以相對真實地模擬地層中的水力壓裂過程(李霞穎等，2015；蘭恒星等，2017；申海萌等，2018)。

1 理論方法

非監(jiān)督學習算法OPTICS(Ordering Points To Identify the Clustering Structur)聚類算法使用三維空間屬性，是進行裂縫平面識別的核心算法，再通過引入時間緯度，可將算法從三維擴展至四維。本文提出的方法基于微地震事件位置數(shù)據(jù)，實現(xiàn)裂縫平面自動識別，具體由兩步實現(xiàn)：首先，獲得對應于單個斷裂面的微地震事件數(shù)據(jù)，采用基于時空屬性約束數(shù)據(jù)的聚類方法來確定這些數(shù)據(jù)；然后，使用最小二乘擬合方法來識別裂縫平面。

OPTICS聚類算法是基于密度的聚類算法(Ankerst et al，1999)，該算法的目標是將空間中的數(shù)據(jù)按照密度分布進行聚類，其思想與DBSCAN(Density-based Spatial Clustering of Applications with Noise)算法(Ester et al，1996)類似，但與DBSCAN算法思想不同的是，OPTICS聚類算法可以通過對密度的不同設置獲得不同的聚類結果，即通過OPTICS算法處理的聚類結果，原理上可以實現(xiàn)對任意密度聚類的提取，這是因為OPTICS算法的輸出基于樣本的有序隊列，從該隊列可以獲得任意密度的聚類結果。Xue等(2018)將DBSCAN算法應用在裂縫識別當中，取得了較好的裂縫識別效果。然而，在基于密度的DBSCAN算法中，半徑參數(shù)ε和密度閾值MinPts的給定對最終聚類結果影響較大。由此可知，通過設置全局參數(shù)進行整體密度聚類的算法其結構不甚合理。正確的算法邏輯應該是分析數(shù)據(jù)空間的稀疏性，針對不同密度簇選取不同的局部密度參數(shù)。該思想的簡單原理如圖1 所示，圖中A簇和B簇選擇的密度參數(shù)顯然與C簇不同；C簇中可以區(qū)分出3個不同的簇，分別為C1簇、C2簇和C3簇，且三簇的密度參數(shù)也不一樣，若強制設置一個密度參數(shù)必將導致其他簇的聚類出現(xiàn)偏差。OPTICS聚類算法是對DBSCAN算法的擴展，該算法對半徑參數(shù)ε不再敏感，在實際計算時，只要確定MinPts值，ε的輕微變化并不會影響聚類結果。OPTICS聚類算法并不直接產(chǎn)生類簇結果，而是通過生成一個增廣的簇排序進行聚類分析，例如，通過設置可達距離參數(shù)作為縱軸，設置樣本點輸出次序為橫軸，其坐標圖示例如圖2 所示。增廣的簇排序代表了各樣本點基于不同密度參數(shù)的聚類結構，根據(jù)該增廣的簇排序可以得到基于任何DBSCAN算法參數(shù)組合的聚類分析結果。

圖 1 密度分布不均勻的空間簇(據(jù)Ankerst等(1999)修改)

圖 2 基于OPTICS聚類分析算法計算得到的增廣簇排序X軸表示OPTICS算法處理對象的序列，Y軸代表可達距離；簇在坐標軸中表述為山谷，山谷越深，簇越緊密；不成簇的點代表噪聲，其不形成任何山谷(據(jù)Ankerst等(1999)修改)

OPTICS聚類算法包含2個核心概念： 核心距離(core-distance)與可達距離(reachability-distance)。

定義樣本集D，假設樣本x∈D，對于給定的ε和MinPts，使x成為核心點的最小鄰域半徑大小稱為x的核心距離(Ankerst et al，1999)。核心距離的數(shù)學表達式為

(1)

可達距離根據(jù)核心距離來定義。對于核心點x的鄰點y，若鄰點y到點x的距離大于點x的核心距離，則其可達距離為鄰點y到點x的實際距離；若鄰點y到點x的距離小于點x的核心距離，則其可達距離為點x的核心距離(Ankerst et al，1999)。可達距離的數(shù)學表達式為

(2)

其中，d(x，y)表示為鄰點y到點x的實際距離。

上述2個核心概念的示意圖見圖3，假設MinPts=4、半徑為ε，則x點的核心距離為d(3，x)，點1、點2和點3的核心距離均為d(3，x)，點4的可達距離為d(4，x)。

圖 3 2個核心概念示意圖(據(jù)Ankerst等(1999)修改)

將時間屬性引入聚類算法，基于時空的非監(jiān)督學習聚類分析方法是對基于空間的非監(jiān)督學習聚類分析的擴展。其使用密度作為對象之間相似性的度量，并將時間-空間聚類視為由低密度區(qū)域分割的一系列高密度連接區(qū)域(Birant et al，2007)。在DBSCAN算法的基礎上，Birant等(2007)發(fā)展了ST-DBSCAN算法，該算法考慮了時間維度的信息。ST-OPTICS算法的主要擴展是使用了時空鄰域方法，該方法類似于時空掃描統(tǒng)計方法(Tango，2010；Pereira et al，2015)。采用時空鄰域方法，ST-OPTICS實現(xiàn)了對時空實體密度的有效估計。引入時間屬性后的時空鄰域區(qū)域示意圖見圖4，時空實體STi的時空鄰域被定義為以半徑參數(shù)(ε)為直徑、2倍時間窗口(2ΔT)為高度的圓柱體。密度概念在時空聚類里指的是在圓柱時空相鄰區(qū)域中的實體STi數(shù)量。ST-OPTICS 算法中的核心距離和可達距離，與OPTICS算法中的定義相同(Ester et al，1996)。由于增加了時間維度，ST-OPTICS算法需要3個參數(shù)：半徑參數(shù)ε、時間窗口ΔT和密度閾值MinPts，其中前2個參數(shù)用來確定時空鄰域。ST-OPTICS算法繼承了OPTICS算法的優(yōu)良特性，例如適應任意形狀的時空簇、不需要數(shù)據(jù)分布預設以及抗噪聲的能力，這一系列優(yōu)良特性是選擇ST-OPTICS算法來進行時空聚類的原因。

圖 4 引入時間屬性后的時空鄰域區(qū)域示意圖(據(jù)Xue等(2018)修改)

假設輸入?yún)?shù)為半徑參數(shù)ε、時間窗口ΔT、密度閾值MinPts以及定位事件點集合C，則ST-OPTICS算法的計算步驟為：

(1)對定位事件點計算并增加核心距離與可達距離屬性，初始化對象結果序列；

(2)循環(huán)集合C中的每個對象點p；

(3)如果對象點p在ε和2ΔT鄰域范圍內(nèi)，其密度滿足參數(shù)閾值MinPts，則判定對象p為核心對象，計算對象點p鄰域內(nèi)其他點的核心距離與實際距離的數(shù)值，取2個數(shù)值中的較大值作為可達距離；

(4)將核心對象與可達距離加入對象結果集合中；

(5)導出對象結果集合，繪制以定位事件點屬性為橫軸、可達距離為縱軸的簇排序，根據(jù)簇大小、鄰域確定ST-OPTICS參數(shù)。

ST-OPTICS算法的工作流程如圖5 所示，其中C為定位事件點集合；Q為有序隊列，元素按照可達距離排序，可達距離最小的在隊首；E為結束隊列，即最后輸出結果的點集的有序隊列。

圖 5 ST-OPTICS聚類算法的詳細工作流程

得到結果隊列后，從結束隊列中按順序取出對象點。如果該對象點的可達距離小于等于給定的半徑參數(shù)ε，則該點屬于當前聚類，處理結束；反之，則跳轉至上文所述的計算步驟(2)。如果該點的核心距離大于給定的半徑參數(shù)ε，則推斷該點為噪聲，可以略過該點；如果該點的核心距離小于等于ε，則該點屬于新的聚類，跳至步驟(1)結束隊列遍歷結束，直至算法結束。最終得到聚類分析結果，即結果序列、每個節(jié)點的可達距離和核心距離。

將上述基于時空的非監(jiān)督學習聚類分析方法具體到裂縫平面識別問題，首先，將考察對象選擇為獲取得到的準確微地震事件信息；然后，利用微地震事件的位置信息以及發(fā)震時間信息進行微地震數(shù)據(jù)時空密度聚類。通過聚類分析方法得到的不同微地震事件聚類簇對應于不同裂縫平面上的微地震事件，進而采用最小二乘擬合方法，可以實現(xiàn)對裂縫平面的求取。

在特征值最小二乘平面擬合方法中，假設平面方程為

ax+by+cz-d=0

(3)

其中，a、b、c為平面方程的單位法向量，d為坐標原點到平面的距離(d≥0)。

假設對某一平面進行掃描得到了n個數(shù)據(jù)點{(xi，yi，zi)，i=1，2，…，n}，根據(jù)式(3)，則任意一數(shù)據(jù)點(xi，yi，zi)到該平面的距離為

di=|axi+byi+czi-d|

(4)

為了獲得最佳擬合平面，需要在滿足a2+b2+c2=1的約束條件下獲得最小殘差，即

e=∑i(axi+byi+czi)2→min

(5)

利用拉格朗日乘數(shù)法求解函數(shù)極值

f=∑i(axi+byi+czi-d)2-λ(a2+b2+c2-1)

(6)

其中，λ為拉格朗日系數(shù)。

分別對a、b、c、d求偏導，并將偏導數(shù)取零，得到矩陣A

(7)

其中，Δxi=xi-∑Δxi/n，Δyi=yi-∑Δyi/n，Δzi=zi-∑Δzi/n。

利用下式求解矩陣A的特征向量I

|A-λI|=0

(8)

矩陣A有3個實數(shù)特征值，最小特征值λmin對應的特征向量即為平面方程的參數(shù)a、b、c，利用重心點可求得參數(shù)d，由此得到擬合的平面方程表達式(官云蘭等，2008；王峰等，2011)。

將微地震事件點投影到上述擬合得到的平面上，根據(jù)投影所得的散點坐標位置確定投影區(qū)域，并計算投影面積，即可得到裂縫平面的大小(李紅梅，2015)。

水力壓裂裂縫平面識別方法可總結如下：①獲得高信噪比的微地震數(shù)據(jù)；②獲得高精度的微地震事件位置和發(fā)震時間；③使用所提出的聚類分析算法獲得分類結果；④使用最小二乘擬合等擬合方法獲得精確的斷裂平面。

2 室內(nèi)水力壓裂實驗驗證

在提出水力壓裂裂縫平面自動識別算法后，需要進行方法可行性的驗證實驗，以驗證方法的有效性。由于在水力壓裂施工中，微地震的產(chǎn)生和裂縫平面的識別結果均缺乏直接檢驗結果的手段，因此，可以選擇水力壓裂巖石物理實驗作為方法的驗證手段。通過水力壓裂巖石物理實驗，將CT掃描數(shù)據(jù)結果與裂縫平面識別結果進行比較。如果CT掃描結果與算法分析結果相似，則可以驗證所提出的水力壓裂裂縫平面自動識別算法的有效性。本次實驗的巖芯為直徑50mm、長度125mm的圓柱狀致密砂巖。

實驗前的CT掃描巖芯結果如圖6 所示。為了監(jiān)測水力壓裂過程中產(chǎn)生的聲發(fā)射信號，在巖芯的表面放置了24個壓電陶瓷傳感器(PZT)。巖芯采用硅膠涂層密封，將測試樣本與圍壓油分離。實驗后，將破裂的巖芯樣本進行CT掃描分析，進而得到實驗后的CT掃描巖芯結果，如圖7 所示，該CT掃描結果清晰地揭示了巖石內(nèi)部的斷裂破壞情況。

圖 6 實驗前CT掃描巖芯結果

圖 7 實驗后CT掃描巖芯結果

掃描分析后，利用記錄到的全波形聲發(fā)射數(shù)據(jù)來獲得聲發(fā)射事件位置和發(fā)震時間信息(圖8)，結合巖芯斷裂后的CT掃描結果，對本文提出的裂縫平面自動識別算法進行驗證。利用圖8 中所有聲發(fā)射事件位置進行平面擬合，得到一個擬合平面(圖9)。本文所提算法的計算結果如圖10 所示，可以看到通過聚類分析，得到了2個聚類結果，其中紅色點和綠色點所示的分類結果對應著巖芯中的2個裂縫，黑色點表示被剔除掉的噪聲點。通過提取圖10 顯示的主要信息，根據(jù)紅色點進行裂縫平面擬合，獲得斷裂平面1(圖11(a))；根據(jù)綠色點進行裂縫平面擬合，獲得斷裂平面2(圖11(b))。

圖 8 聲發(fā)射事件3D定位結果

圖 9 聲發(fā)射事件擬合平面結果

圖 10 聲發(fā)射事件的ST-OPTICS時空聚類結果

圖 11 裂縫平面1(a)和裂縫平面2(b)的識別結果

為了直觀地比較CT掃描結果與裂縫平面識別結果，將二者進行疊合顯示，結果如圖12 所示，其中藍色虛線代表直接利用聲發(fā)射事件定位結果進行平面擬合得到的裂縫平面位置，紅色和綠色實線代表通過聚類分析后得到的2個裂縫平面位置。由圖12 可見，CT掃描得到的裂縫與通過聚類分析計算得到的2個裂縫平面在對應的切片上吻合良好，由此檢驗了本文算法的有效性。

圖 12 疊合顯示CT掃描結果與裂縫平面識別結果

3 結論

本文提出了一種用于微地震定位數(shù)據(jù)的非監(jiān)督學習裂縫平面識別算法，其主要優(yōu)點是利用微地震定位數(shù)據(jù)自身的時間與空間分布屬性，屬于一種完全由數(shù)據(jù)驅動的方法。此外，該算法具有抗噪聲特性，適用于水力壓裂微地震監(jiān)測環(huán)境。通過水力壓裂巖石物理實驗數(shù)據(jù)的驗證，巖石中的主要裂縫可以自動識別，結果是可以接受的。采用該算法可以確定巖芯的整體裂縫斷裂趨勢，雖然已經(jīng)通過巖石物理實驗數(shù)據(jù)驗證，但該算法仍然有很多限制。該算法對巖芯內(nèi)部宏觀的斷裂有較好的識別效果，然而對于位于巖芯邊緣的微裂縫，由于聲發(fā)射事件數(shù)據(jù)有限，目前該算法還無法得到可靠的結果。后續(xù)研究需要挖掘聲發(fā)射數(shù)據(jù)的地球物理特性，如地震矩張量和地震震級等，并將這些屬性添加到聚類算法中，以期得到更加可靠的裂縫識別效果。