基于指纹和相似性阈值方法的水力压裂过程微地震事件检测研究

李秋雨,孟晓波,陈海潮,陈欣星,陈信宇,王丽玮

(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方法,本文方法检测能力更高;另外与模板匹配检测方法相比,本文方法在应用于长时间的微地震数据,尤其是模板事件较多的时候,表现出更高的检测效率,可以用于压裂过程中的实时监测。