基于指紋和相似性閾值方法的水力壓裂過程微地震事件檢測研究

李秋雨,孟曉波,陳海潮,陳欣星,陳信宇,王麗瑋

(1.成都理工大學地球勘探與信息技術教育部重點實驗室,成都理工大學,四川成都610059;2.中國石油大學非常規油氣研究所,油氣資源與勘探國家重點實驗室,中國石油大學(北京),北京102249;3.中國石油國際勘探開發有限公司中油國際(蘇丹)6區項目公司,北京100034)

近年來,隨著非常規油氣資源勘探開發技術逐漸成熟,水力壓裂技術在非常規油氣領域中得到快速發展,水力壓裂通過井筒向目標儲層注入高壓液體,使地層巖石破裂同時釋放能量級別較低的微地震信號[1-3],這些微地震信號的位置對應壓裂過程中巖石破裂的位置。微地震監測通過檢測識別微地震信號,進一步反演得到微地震事件的位置分布,能夠對壓裂過程中產生的裂縫范圍和裂縫發育情況進行實時監測和反饋,有效分析壓裂作業效果,及時調整壓裂方案[4-5]。

微地震事件的檢測識別是微地震數據處理的關鍵環節之一。由于微地震信號能量較弱,受噪聲影響較大,微地震事件檢測前首先需進行去噪處理,去噪效果也會極大程度影響微地震事件的識別和定位效果[6-10]。通常來說,微地震臺站記錄到的信號是從震源通過各種路徑傳播到各個接收臺站的地震波,因此臺站記錄包含了地震波在震源處的信息,以及地震波傳播到臺站路徑受到的影響。同態反褶積去噪能夠在不改變事件頻率特性的情況下清除地震波從震源到接收臺站之間因傳播路徑所產生的干擾,主要包括反射系數,使用該方法能夠提高微地震事件的信噪比,便于后續檢測[11]。微地震檢測方法目前應用較廣泛的一種是基于能量的長短時窗比值(short-term average/long-term average,STA/LTA)方法,STA/LTA方法快速簡單,利用微地震信號和背景噪聲振幅(能量)的差異,計算滑動長短時窗的振幅(能量)比值,提取微地震初至。STA/LTA方法能夠有效地拾取高信噪比微地震事件,但是弱信號檢測精度較低[12-16]。另一種是基于波形相似性的模板匹配方法,該方法通過比較模板地震波形和實時波形之間的相似性,識別微地震事件[17-18]。模板匹配方法能夠檢測到非常低信噪比的微地震信號[19-20],但是該方法依賴于模板的選擇,在模板數量多時,計算其相關系數需要大量時間。YOON 等[21]提出了一種新的波形匹配的思路,將波形轉換為指紋,利用局部敏感哈希比較相似度,匹配相似指紋,并命名為指紋和相似性閾值(fingerprint and similarity thresholding,FAST)方法。FAST方法通過將連續數據分成N個重疊窗口,并將每個窗口對應波形制作為指紋,壓縮了大量信息,對比每一個窗口指紋的相似性,大大提高了計算效率,其檢測能力與模板匹配方法相當。但是在應用于采樣頻率較高的微地震數據時,由于重疊窗口過多,對比指紋相似度依然要比較N(N-1)次,體現不出快速計算的優勢。FAST方法是對單個臺站地震數據處理,需要利用多臺站震相關聯算法將多個臺站識別結果關聯起來,能夠有效剔除誤觸發事件[22-23]。

采用單一方法進行微地震事件檢測時各有適用條件和局限性。近年來,利用各種算法優勢的綜合分析方法備受關注。本文以FAST方法為基礎,提出了一種微地震事件檢測方法,旨在提高計算效率的同時,還能夠識別出低信噪比的微地震事件。本文方法首先利用同態反褶積去噪方法對微地震數據去噪,結合了STA/LTA、模板匹配和FAST方法的優勢,利用STA/LTA方法獲得模板事件,再利用FAST方法將模板事件和連續波形制作為指紋,將指紋進行匹配獲取微地震事件;進一步運用“帶噪聲的基于密度的聚類”(Density-Based Spatial Clustering of Application with Noise,DBSCAN)[24]對多個微地震臺站進行震相關聯,去除錯誤檢測到的事件。

1 方法原理

本文方法結合了同態反褶積去噪、STA/LTA、FAST和DBSCAN方法,其數據處理流程如圖1所示。

圖1 本文方法的數據處理流程

1.1 同態反褶積去噪

本文采取同態反褶積去噪方法[10]分離出微地震信號在傳播路徑中受到干擾而產生的噪聲。假定傳播通道作為一個線性時不變系統(LTI),用卷積對震源產生的信號和擾動進行建模,可以表示為:

r[n]=s[n]?hT[n]

(1)

式中:r[n]為臺站記錄的地震信號;s[n]為震源處的地震信號;hT[n]為反射系數,表示沿路徑產生的干擾;?表示卷積算子。同態反褶積,通過將卷積模型轉換為求和,使震源信號和反射系數分離,在理論上可以得到分開的震源信號和反射系數。

同態反褶積去噪方法首先將長度為n的離散時間信號r[n]分成m個窗口的離散數據,每個窗口的長度為N,利用公式(2)和公式(3)計算每一個窗口內時間信號的離散傅里葉變換(DFT)。

(2)

WN=e-j(2π/N)0≤k≤N-1

(3)

對R[k,m]取自然對數得到L[k,m]:

L[k,m]=ln(R[k,m])

(4)

(5)

上述過程是對(1)式分窗做傅里葉變換以及取對數,所以公式(1)變為:

(6)

(7)

S[k,m]=eL′[k,m]

(8)

(9)

將S[n,m]按列合并之后得到的數據即為去噪處理后的微地震數據。

1.2 STA/LTA獲取模板

STA/LTA方法是地震信號初至拾取的常規有效方法,能夠快速、準確地獲取高信噪比的微地震事件[11-15]。STA/LTA方法原理如圖2所示,對于采集到的微地震記錄,計算兩個連續滑動的長、短時間窗口之間平均能量之比。當微地震信號到達時,長短時窗能量比會有一個突變,對于大于預定的判斷閾值k的時間點,將該時刻標記為震相到時。

圖2 STA/LTA方法原理示意

對于時間點i,STA/LTA比值計算公式如下:

(10)

(11)

式中:S代表滑動短時窗的長度;L代表滑動長時窗的長度;xi代表時間點i處微地震數據的振幅;CF為特征函數,表示微震數據能量的變化情況;STAi和LTAi分別為微震信號在i時刻的短時窗和長時窗平均能量值。短時窗的窗長取值應與待測信號的周期有關,取值太短會誤觸發微地震事件,取值太長會漏觸發微地震事件;長時窗的窗長取值與信號的噪聲水平有關,取值范圍相對固定。

STA/LTA方法適用于高信噪比的微地震數據識別,對信噪比的微地震數據識別并不準確。因此將識別到的信噪比較高的信號當作模板事件,進一步使用FAST方法以其結果為指紋,識別出低信噪比的微地震信號。

1.3 FAST方法

YOON等[21]提出的FAST方法(以下簡稱Yoon方法)是將地震數據轉換為指紋,壓縮了信息,能夠快速地處理長時間、高采樣頻率的地震數據(分析一個星期的連續波形數據用時不到2h)[19]。同時,FAST方法還能夠檢測到STA/LTA方法檢測不到的低信噪比微地震事件。

以STA/LTA方法獲取的模板事件為例,采用FAST方法對地震信號數據進行如下處理。首先對模板事件的時間序列數據(圖3a)進行短時傅里葉變換,計算得到復數離散傅里葉變換的功率(振幅的平方),地震事件在頻譜圖中以瞬態高能的形式出現(圖3b)。對這個頻譜圖像進行二維哈爾小波變換,小波變換是多分辨率分析工具,地震信號在全分辨率小波系數中具有高能量(圖3c),有利于進一步識別圖像。小波變換對數據進行壓縮的同時,能夠最大程度地保留圖像信息。

圖3 FAST方法對微地震信號數據處理流程a 某微地震事件的時間序列數據; b 對a進行短時傅里葉變換之后的頻譜; c b的二維哈爾小波變換的小波變換結果; d 對c圖像進行壓縮的指紋圖

接下來選擇丟棄一部分小波系數對數據繼續壓縮。由于接收到的連續信號大部分都是噪聲,并且通常情況下地震信號的能量大于背景噪聲,因此對小波圖像中的哈爾小波系數矩陣進行z-score的標準化。

(12)

式中:Hij為小波變換后第i行j列的小波系數;ui為小波矩陣第i行的均值;σi為第i行的標準差。標準化的目的是為了更好地區別信號和背景噪聲。保留與標準化的哈爾小波系數平均值偏差最大的前50%的哈爾小波系數。將大于平均值的正向偏差的小波系數設置為+1,表示正數;將小于平均值的負向偏差設置為-1,表示負數;將舍棄的50%的哈爾小波系數設置為0。為了便于使用后續的局部敏感哈希算法有效搜索相似指紋并減少存儲所用的比特數,用兩個比特來表示標準哈爾系數的符號:-1變為01,0變為00,1變為10。因此每個指紋的比特是哈爾系數的2倍。處理后,哈爾小波變化結果轉換為只包含0和1且稀疏(大多為0)的“指紋”(圖3d),白色代表1,黑色代表0。再利用FAST方法,將STA/LTA方法獲取的高信噪比微地震事件波形制作為指紋,便于快速進行相似性搜索。

表1 最小哈希數組

與上述過程一致,依據模板事件以選取時窗長度0.7s為例,將臺站收集的一段10min的微地震數據分成若干個重疊的0.7s的時窗波形數據,時窗間隔為0.01s。這樣,10min的數據分成了59931段波形數據,將這些數據經過圖3所示的處理過程轉換為指紋,微地震事件指紋與模板指紋會更相似,比較這些指紋與模板指紋的相似度,實現對微地震事件的識別。

最小哈希方法和局部敏感哈希方法可以用來比較指紋之間的相似度。最小哈希方法是一種快速判斷兩個集合是否相似的方法。通過多個哈希函數hi(x)作用于同一個指紋,每個哈希函數將指紋映射為一個整數。如表1所示,事件X和事件Y的指紋通過選取6個哈希函數hi(x)映射之后,分別轉換為6個整數。

最小哈希的一個重要屬性就是兩個指紋映射為同一整數的概率等于指紋的杰卡德(Jaccard)相似度。杰卡德相似度定義為:

(13)

事件X和事件Y相似性越高,指紋被映射為同一整數的概率越高。最小哈希將兩個事件X,Y的二維指紋映射為了兩個一維向量,最小哈希函數越多,映射之后的一維向量越能夠保留X和Y之間的相似度,而當指紋總數足夠多時,映射為一維向量相似度比較的計算量也越大。

局部敏感哈希方法將指紋對應的一維數組分為多個條帶,任一條帶相同時,放入一個哈希桶,放入同一哈希桶的次數越多,指紋間相似的可能性就越大,從而減小計算量。如表2所示,將事件X和Y的6個哈希值分成3個條帶,每個條帶由2個哈希值組成,當且僅當兩個條帶中每個哈希值相等時,將兩個條帶放到同一哈希桶中。表2中h(X)和h(Y)的第1個條帶同為[45,23],因此它們放入同一哈希桶。對于第2個條帶,[14,11]和[21,11]是不完全相同的,因此將它們放入不同的哈希桶中,同理,h(X)和h(Y)的第3個條帶放入同一哈希桶。最終統計X和Y的指紋在同一桶內的概率作為兩個事件的相似度。例如表2中指紋X和指紋Y的最終相似度為2/3。概率越大,多個條帶中的向量相等的可能性就越大,兩個指紋相似度越高。

表2 相似度匹配

Yoon方法是通過比較采用FAST方法對全波形數據處理得到的所有指紋的兩兩相似性進行微地震事件的識別。一般認為,微地震事件會重復發生,對應指紋具有相似性,而噪聲是隨機的,不同時刻的噪聲指紋不存在相似性。事實上,井下接收到的連續信號大部分是噪聲,這樣的比較主要是噪聲和噪聲之間的比較,是沒有必要的。本文以STA/LTA方法獲取的微地震事件為模板并壓縮為指紋,然后與全波形數據處理得到的所有指紋比較相似度,從而減少比較次數,提高了計算效率。

1.4 DBSCAN方法

在水力壓裂過程中,井下微地震檢測臺站的位置距離裂縫破裂位置較近,通常會在一個短時間內記錄到同一微地震事件。比較模板指紋和每個記錄波形窗口的指紋相似度,相似度最高的波形窗口的初至到時和模板事件的初至到時對齊,從而檢測得到每個臺站的震相到時。將不同臺站的震相到時進行關聯,可以減少錯誤檢測。DBSCAN聚類方法是比較有代表性的基于密度的聚類算法,該方法可以在有噪聲的數據中發現任意形狀的聚類。DBSCAN聚類算法主要有兩個參數,即ε和minPts,其中,參數ε描述了核心點的鄰域大小,minPts描述了某一樣本距離為ε的鄰域中樣本個數的閾值。DBSCAN聚類方法原理如圖4所示,如果A的鄰域ε內至少有minPts個樣本點,那么A為核心對象。A中的不滿足鄰域ε內至少有minPts個樣本點的稱為非核心對象(圖4中A的鄰域內的黑色點)。如果核心對象在另一核心對象鄰域內,則兩個核心對象密度可達(圖4中通過黑色箭頭相連接的紅點)。所有密度可達的核心對象與其內部的非核心對象共同組成一個聚類簇,如圖4中通過聚類最終形成3個聚類簇A、B、C。在微地震震相關聯的過程中,點對應的為微地震信號在每個臺站的到時信息,任何兩個點之間的距離為不同臺站之間的到時差。ε定義為一個微地震事件在不同臺站到時差的最大值,考慮到臺站位置和壓裂區域的相對距離,minPts定義為檢測到同一微地震事件的最小臺站個數,如果臺站個數小于minPts,定義為噪聲。

圖4 DBSCAN聚類方法原理示意(圖中,minPts為4,紅點為核心對象(其鄰域至少有4個樣本點),黑點為非核心對象。如果核心對象在另一核心點鄰域內,則兩個核心對象密度可達,圖中由黑色箭頭連起來的核心對象組成了密度可達的A、B、C這3個聚類簇)

龍峰等[24]在利用DBSCAN的過程中指出在地面臺站數據對微地震事件定位的過程中,至少需要4個臺站才能完成定位。minPts越大,檢測到微地震事件的正確率越高,后續定位結果越可靠,但同時也會丟失一部分事件。本次研究基于深井中的10個臺站,同樣選擇minPts為4,在保證正確率的基礎上盡可能地保留更多微地震事件。通過計算得到微地震事件在各個臺站上的震相到時差在0.4s以內,此時設定ε為0.4s,即通過DBSCAN方法,在0.4s內有4個不同的臺站都檢測到了這一地震事件時,才判定這一事件為有效的微地震事件。

2 數值模擬

為了驗證本文方法檢測微地震事件的可行性,首先采用合成地震數據測試方法性能。本次合成數據基于四川盆地威遠頁巖氣開發水平井第19段壓裂井附近的10個井下臺站記錄得到的10個實際微地震事件(彼此互相關系數較低),該10個微地震事件位置和震源機制解如圖5所示,可以看出10個微地震事件的震源機制解相互不一致。調整10個微地震事件的振幅,并隨機插入對應臺站記錄的一段噪聲數據中的不同時間段合成171個微地震事件,不同的臺站保留事件的到時差,由同一個事件生成的微地震事件具有相同的位置和震源機制解。圖6a為合成的10個臺站中編號ZA20臺站的微地震數據,圖6b展示了ZA20臺站的微地震信噪比分布。

圖5 臺站與微地震事件位置(圖中,沙灘球為10個微地震事件位置和震源機制解,紅色五角星為10個臺站的位置)

圖6 ZA20臺站微地震數據(a)及其信噪比(b)

對10個臺站的合成數據進行測試以驗證本文方法的實用性。對每個臺站的微地震數據進行相同的數據處理,以ZA20臺站為例。首先利用同態反褶積去噪,并通過計算信噪比RSN評價同態反褶積去噪效果:

(14)

式中:PS和PN分別代表信號和噪聲的有效功率。圖7展示了同態反褶積去噪處理前、后ZA20臺站在第514.55s的效果。圖7a為去噪前的微地震信號,圖7b為同態反褶積去噪后的微地震信號。經試驗,去噪后微地震信號信噪比平均提升了6dB以上。

圖7 同態反褶積去噪前(a)、后(b)的信號

對去噪之后的數據,利用STA/LTA方法提取出強信號事件,便于后續制作為模板指紋利用FAST方法進行處理。其中STA/LTA方法的參數S和L分別設置為0.1s和0.4s,滑動步長為0.005s,觸發閾值λ設置為1.4。圖8為圖6a所示ZA20臺站波形計算得到的STA/LTA值,紅色虛線為設置的閾值,在該閾值條件下,可以檢測到90個微地震事件,降低閾值時,可以得到更多微地震信號,但同時增大了檢測錯誤率。分析圖8可以看到,STA/LTA值設定在1.3左右時,會使檢測數量突增,可能會帶來部分噪聲。我們將閾值設定為1.3時,檢測到146個事件,其中有12個是誤檢測,錯誤率為8%;另外有37個事件沒檢測到。為確保模板事件的可靠性,本文選取閾值為1.4時STA/LTA方法檢測到的信噪比大于2.0的事件作為模板,挑選得到80個事件。

圖8 STA/LTA拾取結果

采用FAST方法進一步進行地震信號相似性搜索,可以搜索到額外的低信噪比信號。FAST方法能夠充分地檢測到與模板事件相似的事件,如圖9所示,以STA/LTA獲取到的119.526～119.626s的波形為模板波形(第1道),搜索到了不同時間的其它相似波形,圖中展示了加上模板的其它6個波形,圖9中從上往下事件波形的信噪比依次為17.9,1.7,1.1,18.0,14.6,16.8,12.1dB。其中第2道、第3道地震波形信噪比小于2.0,計算的STA/LTA值分別為1.04和0.99,為STA/LTA方法未能檢測到的地震事件。

圖9 FAST方法波形匹配結果

利用FAST方法對每個臺站進行搜索后,共得到資料段內10個臺站的2126次觸發(其中可能包括部分虛假檢測)。最后利用DBSCAN方法對10個臺站得到的所有震相進行關聯,排除檢測到的噪聲信號。采用ε=0.4s,minPts=4進行DBSCAN聚類計算,其中171次觸發參與了聚類,171次地震事件全部檢測到,未參加的離群點為隨機噪聲。圖10為某一微震事件在第504～520s的聚類效果,可以看出,504～520s發生了一次聚類,在506s處發生了一次錯誤檢測,但是離群點并不構成聚類條件,沒有對檢測構成影響。圖11展示了514.50s左右的細節,10個臺站均有檢測并構成一個聚類簇。在FAST方法設置閾值較低的情況下,這16s內存在一次離群點,但噪聲的觸發都未能構成聚類簇,說明在低閾值的情況下,依然沒對檢測效果產生影響。經過本文方法處理后,171次地震事件均被完全檢測?？梢?本文方法對于微地震事件的檢測是可行的。

圖10 基于DBSCAN方法的事件聚類

圖11 基于DBSCAN的事件聚類細節(514.50s左右)

3 實際應用

實際微地震數據為四川盆地威遠頁巖氣開發水平井19段壓裂井微地震數據集。本次微地震監測20級井中檢波器布設在一直井段,深度范圍為2270～2405m,每個臺站在深度方向間隔15m,各站點分布位置如圖12所示。

圖12 地震臺站位置分布(圖中,藍色豎線表示不同的壓裂段的射孔事件,紅色豎線表示本文數據的第19級壓裂段的射孔事件位置,紫色倒立三角表示放置在井中的20個地震檢測臺站)

采用本文方法對壓裂段在2014年11月10號7∶50至11∶50由井下10個臺站采集的數據(數據以10min為間隔,劃分為21個時間段)進行處理并與STA/LTA方法和模板匹配方法進行比較。井下數據的采樣間隔為2000Hz。本方方法選取的參數如表3所示,模板匹配方法所用波形模板與本文方法一致。

表3 輸入參數

采用STA/LTA方法處理后總共識別出754個微地震事件,從中挑選信噪比大于2.0的作為模板事件用于本文方法和模板匹配方法。處理后,本文方法共計識別出2604個微地震事件,模板匹配方法共計識別3052個事件。圖13展示了壓裂起始3種方法對于壓裂段微地震數據檢測事件數量分布。采用本文方法和模板匹配方法檢測的微地震事件數量均遠遠超過STA/LTA方法的檢測數量。本文方法在多數時間段內與模板匹配方法的檢測結果一致,但在第2,3,5,7,17這5個時間段內,模板匹配檢測到的微地震數量比本文方法的微地震數量多。

圖13 STA/LTA方法、模板匹配方法和FAST方法實際微地震數據處理結果對比

為進一步分析本文方法檢測微地震事件數量少于模板匹配方法的原因,取FAST方法和模板匹配方法在壓裂時間段775.8～780.0s(即圖13中第2個時間段)內的檢測結果進行分析,結果如圖14所示。采用FAST方法和模板匹配方法均檢測到776.0s和777.8s的震相,信噪比較高,分別為13.8和18.9,模板匹配方法還在777.4,778.8,779.5s處檢測到3個事件,這3個事件信噪比較低,分別為1.32,0.78和1.18。分析可知,對于信噪比小于1.5的微地震事件,本文方法檢測能力低于模板匹配方法,可能是由于FAST方法在將波形信息壓縮為指紋加快運算效率的同時,損失了部分波形信息。但是本文方法對微地震數據的處理速度比模板匹配方法的處理速度快,在intel(R)core(TM)i5-10400cpu上模板匹配方法處理本次壓裂數據用時5.0h,本文方法處理只用了1.5h。在模板數量更多的情況下,本文方法的效率比模板匹配方法的效率會更高。

圖14 FAST方法和模板匹配方法在775.8～780.0s的檢測結果

將實際壓裂過程4.0h的數據分成2879981個重疊窗口,采用本文方法時,是對754個模板波形和2879981個窗口波形的指紋進行相似性比較,比較次數為754×2879981;采用Yoon方法時,是對2879981個波形的指紋進行兩兩比較,計算2879981×2879980次相似性。另外STA/LTA方法獲取模板的時間很短,通常來說,連續時間波形劃分得到的窗口個數比模板事件波形的個數多得多?？梢?本文方法優化了Yoon方法,提高了計算效率,對于微地震數據的實時監測處理具有意義。

4 結論

針對微地震信號弱、數據量大的特點,本文提出了一種基于FAST方法并結合同態反褶積去噪、STA/LTA以及DBSCAN方法的水力壓裂過程中微地震事件檢測方法。理論測試驗證了該方法的有效性:能夠有效地去除噪聲,利用STA/LTA方法快速得到模板事件,基于FAST方法比較模板事件指紋和記錄波形指紋的相似度,能夠快速、精確地獲取各個微地震事件對應的震相到時信息,可供后續地震定位、速度結構成像等分析使用。在實際數據應用中,相較于常規STA/LTA方法,本文方法檢測能力更高;另外與模板匹配檢測方法相比,本文方法在應用于長時間的微地震數據,尤其是模板事件較多的時候,表現出更高的檢測效率,可以用于壓裂過程中的實時監測。