一种适于存在极性反转的微震初至到时拾取方法

毕丽飞，曾志毅，张建中，黄忠来，芮拥军，刁 瑞

（1.中国石油化工股份有限公司胜利油田分公司科技处，山东东营257001；2.海洋科学与探测技术教育部重点实验室，中国海洋大学海洋地球科学学院，山东青岛266100；3.中国石油化工股份有限公司胜利油田分公司物探研究院，山东东营257001）

水力压裂微震监测技术是一项通过观测注水压裂过程中所诱发的微震事件来监测和评估裂缝发育情况的技术[1-4]。微震事件识别和到时拾取是微地震数据处理的关键环节，直接影响对裂隙成像的可靠性[5]。通常情况下，水力压裂产生的微震事件能量十分微弱，微震观测资料的信噪比较低。低信噪比资料的微震事件识别及其到时拾取是微震监测技术的重要研究课题。

微震监测记录的数据量巨大，人工识别并拾取微震事件初至到时的工作量极大，且容易引入人为误差。微震事件到时拾取是根据有效信号与背景噪声在能量、偏振、频谱和统计等方面的特性差异来实现的。常见的初至到时拾取方法基于单道记录的算法，其中以长短时窗能量均值比（STA/LTA）算法[6-9]及其改进方法为代表[10-12]。该类方法计算速度快，可以满足实时处理要求，但是对低信噪比微震资料的应用效果不佳。近年来，一些学者综合利用微震信号与环境噪声的多种特征差异，提高了识别和拾取效果。如：宋维琪等[13]利用小波分解与Akaike信息准则相结合，有效压制噪声引起的局部效应，提高微震事件初至拾取的准确性；吴治涛等[14]联合小波变换与偏振分析实现微震事件的全自动识别和拾取P 波初至；谭玉阳等[15]通过综合地震信号与环境噪声在能量、偏振和统计方面存在的差异，提出一种微震事件初至拾取算法，有效提高了抗噪能力和拾取精度。

GIBBONS等[16]的研究表明，具有相似震源位置及破裂机理的地震事件通常在记录上表现出相似的波形特征。一些学者发展了基于事件多道记录相似特征的初至拾取算法[17-19]，如：谭玉阳等[20]通过多道相似系数实现微震事件的识别，避免了由于某些道数据的信噪比过低而遗漏微震事件的情况；VANDECAR 等[21]利用多道互相关与最小二乘法相结合，确定远震震相的相对到时，提升了拾取效率和抗噪性；DE MEERSMAN 等[22]和AKRAM 等[23]分别基于互相关的迭代算法优化了微震事件的初至拾取结果；魏梦祎等[24]提出一种基于波形互相关的微震事件自动拾取方法，通过叠加获取高信噪比参考道，提高了初至拾取精度。在拾取微震事件的过程中，以上方法需要对各道微震事件记录进行叠加得到信噪比更高的叠加参考道，但微震事件信号的初至极性与震源机制相关，各道微震事件的初至极性会不一致，即初至极性具有正负之分[25-28]，致使通过波形直接叠加得到的叠加参考道的信噪比反而降低，影响了对微震事件的识别和到时拾取精度。

本文提出了一种适于初至极性反转的微震到时拾取方法。通过计算任意两道之间的延时互相关，获得各道微地震的相对到时；利用相邻道相乘后再叠加得到参考道，避免了道间波形相似系数（Semblance）和直接波形叠加得到参考道时无法有效克服极性反转影响的困难。本文首先介绍了微震事件识别和初至拾取的方法原理和技术流程，然后分别利用合成数据和实际资料对方法的性能进行了测试和分析。

1 方法原理

1.1 相对走时计算

取一时窗内的微震记录，利用互相关函数计算任意两道微震记录之间的互相关系数序列。时窗内应包含同一微震事件在各道记录的P 波或S波初至，时窗长度应大于微地震波的时长以及最短和最长偏移距检波器之间的到时差，可根据射孔信号时长及最短与最长偏移距检波器的到时差长度进行估计。互相关系数计算公式为：

式中：xi和xj分别为地震道i和j的微震记录；W为时窗长度；k为xi和xj的时间间隔；ci，j（k）为互相关系数。考虑到初至极性存在反转的情况，对互相关系数ci，j（k）取绝对值，其最大值cmaxij对应的k为xi和xj的到时差，记为Δti，j。两个地震道i和j包含的同一微震事件的到时差Δti，j可表示为：

式中：ti和tj分别为地震道i和j的微震初至到时。对所有地震道记录进行两两互相关处理，可得到两两记录中微震事件的到时差。以4个地震道（i＝1，2，3，4）为例，可得到6个到时差Δt1，2，Δt1，3，…，Δt3，4。为了避免方程组求解过程的不适定性，增加一个约束方程∑ti＝0，使相对到时的均值为0。这样有下列方程：

若参与计算的微震记录有N道，则需利用（1）式进行N（N－1）/2次互相关函数计算，得到N（N－1）/2个两道间的到时差数据。（3）式写成矩阵形式为：

式中：A为0，1和－1组成的稀疏系数矩阵；Δt为到时差向量；t为到时向量，也是待求向量。（4）式的最小二乘解为：

利用（5）式求出t后，再用t中的各个分量对各道记录进行时差校正。时差校正后，时窗内各道记录的微震初至波在时间上将趋于一致。

1.2 参考道获取

利用（5）式求出的是相对于时差校正后各道记录平均到时的相对到时。求微震事件在各道记录上的绝对到时，还需要求出各道的平均到时，即微震事件在叠加道或参考道上的初至到时。求参考道及其微震事件的绝对到时的常用做法是：先将时窗内时差校正后的各道记录波形直接进行叠加，形成高信噪比的参考道，再在参考道上拾取绝对到时。但是，当时窗内微震事件的初至极性存在反转时，对校正后的各道微震记录进行直接叠加，会使微震记录振幅相互抵消，反而降低了参考道记录的信噪比，从而影响微震事件绝对到时的拾取精度。有的方法通过对各道的初至极性进行识别或通过震源机制反演，将初至极性校正为相同极性后，再叠加得到参考道[23]，这样做将会增加计算量。为了解决这个问题，我们利用相邻道零延迟相乘再叠加的方法获取参考道（图1），计算式为：i＝1

式中：N是相邻地震道数的配对数；W是时窗长度；w是时窗内采样点序号；x′i（w）表示时差校正后的地震道i在采样点w处的幅值；X（w）为获得的参考道在采样点w处的幅值。图1a中，时差校正后微震初至被“拉平”，其中左边记录（黑色地震道）与右边记录（红色地震道）的初至极性相反。对时差校正后的微震记录进行相邻道的零延迟相乘，初至极性相同（同为正或负极性）的相邻道记录的零延迟相乘后的幅值均为正值；而初至极性相反（一正一负）的情况仅出现在不同初至极性区域的分界处，如图1b中第8道（红色）所示。幅值为负的道的数量很少，对整个相乘道的叠加值影响很小，而随机噪声的极性无一定规律，叠加后振幅减弱，从而得到信噪比更高的参考道，如图1c所示。

图1 相邻道零延迟相乘再叠加得到参考道示意a存在初至极性相反的地震道；b 相邻道零延迟相乘道；c参考道

1.3 绝对到时拾取

由相邻道零延迟相乘再叠加得到的参考道，虽然没有完整的微震信号波形，但是在微震信号初至时刻附近有较大的振幅或能量变化，因此，本文利用STA/LTA 算法判断参考道是否含有有效微震信号以及确定其初至到时。STA/LTA 算法的计算式为：

式中：X（w）为参考道记录；WL和WS分别为长、短时窗的长度；m为采样时刻；STA 和LTA 分别为短时窗和长时窗内的能量均值。在初至到时附近，STA 比LTA 变化更快，得到时窗内的STA/LTA会出现明显的极大值，一般近似地将STA/LTA 曲线的极大值点作为初至到时点。若时窗内的STA/LTA 的极大值超过所设阈值[9]，则认为该时窗内含微震事件，并且同时将STA/LTA 曲线的极大值点所对应的时刻取为该微震事件在参考道的初至到时，记为T0；否则认为该时窗内不含微震事件。本文以各个时窗的长短窗能量比值的均值为依据，将该均值的某个倍数（以下称为长短窗能量比值均值倍数）作为微震事件的识别阈值，这样，就可在每个时窗内自适应地设置阈值，从而避免固定阈值所带来的识别错误和遗漏。

获得参考道的初至到时T0，再结合（5）式中各道的相对到时ti，则各道微震事件的初至绝对到时Ti可由下式得到：

1.4 算法流程

本文的微震事件初至拾取算法流程如图2所示，关键步骤概括如下。

1）确定滑动时窗大小及滑动步长，设置用于判别微震时间的长短窗能量比值均值的倍数b。

2）提取时窗内的微震记录，计算该时窗内的道间互相关系数，确定微震事件的到时差Δti，j，并利用（5）式计算各道微震事件的相对到时ti。

3）利用得到的相对到时，对各道微震信号进行时差校正，对时差校正后的微震记录进行相邻道零延迟相乘再叠加处理，得到参考道。

4）利用STA/LTA 算法得到参考道的长短时窗能量比曲线，并计算能量比均值与设置的倍数b的乘积，即微震事件识别阈值R。判断该时窗能量比极大值F是否大于微震事件识别阈值R，若是，进行步骤5）；否则返回步骤2），进行下个时窗的处理。

5）拾取STA/LTA 能量比曲线的极大值对应时刻T0，再利用（8）式求取时窗内各道微震事件初至的绝对到时Ti。

图2 微震事件初至到时拾取算法流程

2 合成数据测试

首先利用合成数据来检验本文方法的应用效果。图3 为合成微震记录，时间长度为3s，采样率为1 ms。图3a为无噪声微震记录，含有3个有效微震事件，从第1个到第3个微震事件的振幅不断减小，其中第1个微震事件初至极性在第8～16道上发生反转。图3b为添加均值为0，标准差为0.4的高斯白噪的微震记录，其中第3个微震事件信号已经完全淹没在噪声中。分别采用常规方法和本文方法进行测试，选用的滑动时窗长度和滑动步长分别为120个采样点和15个采样点，取各个时窗参考道的长短窗能量比均值的3.5倍，作为微震事件的识别阈值。

图4a和图4b分别为采用各道记录波形线性叠加参考道和本文的相邻道零延迟相乘叠加参考道得到的STA/LTA 能量比曲线（蓝色曲线）及自适应阈值曲线（红色曲线）。可以看出，图4a所示的微地震记录直接叠加的STA/LTA 能量比曲线只有1个明显的峰值（红色箭头所指），指示了第2个微震事件，而对能量最强但存在初至极性反转的第1个微震事件和信噪比较低的第3个微震事件却没反映出来；图4b所示的曲线有3 个明显的峰值（红色箭头所指），3个微震事件都被识别出来。可见，本文方法不但适于初至极性存在反转的情况，对低信噪比资料也有很强的适应性。

图3 合成微震记录

图4 采用不同方法识别参考道微震事件的STA/LTA 能量比曲线（蓝色）及自适应阈值曲线（红色）

下面用本文方法对微震事件初至到时进行拾取，关键步骤的处理结果如图5 所示。图5a 为根据图3b中3个微震信号所在的时窗分别提取出来的经时差校正后的含微震时间记录，与理论微震事件的到时一致，说明微震事件相对到时计算的准确性。图5b是相邻道零延迟相乘结果，在一定程度上压制了不相干的背景噪声。图5c是对图5b各道进行叠加的结果，即参考道，具有明显的振幅极值，消除了初至极性反转对参考道的不利影响。图5d 为利用STA/LTA 算法得到的参考道的STA/LTA 能量比曲线及其极大值对应的参考道的初至到时T0（红色虚线所示）。

图5 用本文方法对微震事件初至到时拾取过程的相关结果（左、中、右分别对应图3中的第1、第2和第3个微震事件）

利用得到的各道相对到时ti和参考道初至到时T0并用（8）式计算出微震事件在各道的初至到时Ti。图6为采用本文方法拾取微震事件初至到时的误差（拾取到时与理论到时之差）。从图6可以看出，整体拾取误差较小，说明采用本文方法可以得到较为准确的微震事件初至的绝对到时。同时我们采用基于波形线性叠加参考道的识别方法拾取第2个微震事件初至到时，其相关结果和拾取误差如图7所示，与采用本文方法得到的第2个微震事件初至到时拾取结果仅相差1 ms（1个采样点长度），说明本文方法的拾取精度与常用拾取方法的拾取精度相当。但是，采用常用方法只能识别出不存在极性反转的第2个微震事件，未能识别出含有极性反转的第1个微震事件和信噪比很低的第3个微震事件，而采用本文方法可以识别出这3个微震事件。可见，采用本文方法不仅可以较为准确地识别和拾取低信噪比微震事件，且能较好地克服初至极性反转所带来的不利影响，具有更强的适用性。

图6 采用本文方法拾取微震事件初至到时的误差（拾取到时减去理论到时）

图7 采用基于波形线性叠加参考道的识别方法拾取第2个微震事件到时的相关结果及拾取误差

3 实际资料处理

图8 实际野外检波器布设位置

采用本文方法对某区实际地面微震监测资料进行处理和分析。野外布设42个检波器，位置如图8所示。图9为一段经过带通滤波后的地面微震监测记录，时间长度为4 s，时间采样间隔为1 ms。采用本文方法和常用方法处理该段实际资料时选取的参数为：滑动时窗为200个采样点，滑动步长为15个采样点，判别微震事件的阈值设置为该时窗内长短窗能量比均值的3.5倍。图10a和图10b分别为采用微地震记录波形线性叠加方法和本文方法得到的参考道STA/LTA 能量比曲线（蓝色）及自适应阈值曲线（红色），可以看出，微地震记录波形线性叠加参考道的STA/LTA 能量比曲线没有出现明显的极大值，而本文方法的STA/LTA 能量比曲线出现了2个超过识别阈值的明显的峰值（图10b中红色箭头所示），检测出该段微震记录中的2个微震事件。

图9 野外实际微震监测记录

利用本文方法对该段微震记录的2个微震事件的初至到时进行了拾取。图11为实际微震信号的拾取过程的相关结果。图11a为时差校正后的微震信号记录，可以看出，微震事件波形基本上被校正拉平，其中第1个微震事件中的第5～12道和第35～40道的微震初至极性发生了反转，波形线性叠加得到的参考道振幅会相互抵消，这是利用常规波形线性叠加参考道方法无法识别出微震事件的主要原因。第2个微震事件记录的信噪比较低，个别道初至同相轴未被完全拉平，也导致波形线性叠加得到的参考道振幅较小，漏识了该微震事件。本文方法利用相邻道零延迟相乘再叠加得到参考道，突出了被拉平的各道初至同相轴的作用，受部分未拉平道的影响相对较小。图11b为利用相邻道零延迟相乘的结果，有效压制了背景噪声，提高了信噪比。图11c为由相邻道零延迟相乘再叠加形成的参考道。图11d为利用STA/LTA 算法对参考道进行处理得到的STA/LTA 能量比曲线，以及由该曲线的极大值点得到参考道的初至到时T0（图中红色虚线所示）。最后把各道的相对到时ti和参考道的绝对到时T0相加，得到了微震事件在各道的初至绝对到时Ti，如图12所示（红点为各道微震事件的初至到时时刻）。

图10 采用不同方法得到的参考道STA/LTA 能量比曲线（蓝色）及自适应阈值曲线（红色）

图11 实际微震事件初至到时拾取过程的相关结果（左列对应915～1 115 ms时窗记录；右列对应2 565～2 765 ms时窗记录）

为进一步说明本文方法对实际资料中微震事件的识别率及初至拾取精度，以震源扫描定位算法[28]得到微震事件数目作为标准，分别采用本文方法和常用的基于波形线性叠加的微震事件识别方法对同一段压裂的微地震监测数据进行识别和拾取，结果如表1所示。从表1可以看出，采用本文方法的误拾率和漏拾率都小于采用波形线性叠加参考道方法的结果，尤其是漏拾率远小于后者，说明采用本文方法能识别出更多的有效事件，具有更高的微震事件识别能力。

图12 采用本文方法对某段监测记录的微震事件到时的拾取结果（红点为初至到时时刻）

由于采用扫描定位方法确定的初至到时受到速度模型精度的影响，因此本文在采用扫描定位法确定到时的基础上进行人工拾取，再以人工拾取结果为准确值，计算本文方法和波形线性叠加参考道方法拾取微震初至到时的误差。波形线性叠加参考道方法共识别出26个微震事件，这里分别计算了采用本文方法和波形线性叠加参考道法拾取这26个微震事件的到时误差，统计的拾取误差的频率如图13所示。从图13可以看出，采用波形线性叠加参考道方法的拾取误差主要分布在2～21 ms，而采用本文方法的拾取误差主要集中在2～12 ms。总体上，采用本文方法的拾取误差要低于采用波形线性叠加参考道方法的误差。

表1 采用不同方法识别微震事件的结果

图13 采用不同方法拾取初至到时的误差率

4 结论

常用的多道相似系数法由于在识别有效微震事件时受初至极性反转的影响，因而会导致微震事件的遗漏。本文根据微震事件的波形相似特征，并用相邻道零延迟相乘再叠加的方法取代常用的波形线性叠加方法，获取更高信噪比的参考道，提出了一种从多道微震监测记录中识别和拾取微震到时的方法。该方法不但能够压制不相干噪声，从低信噪比监测资料中获取有效微震事件信息，而且能够克服初至极性反转的不利影响。合成数据和实际资料应用结果表明，本文方法具有较强的微震事件识别能力和较高的微震到时的拾取精度，对低信噪比资料和存在极性反转资料有较强的适应性。

在判别参考道各时窗是否存在微震事件时，以各个时窗的长短窗能量比均值为依据，将该均值的某个倍数作为判别的阈值，这样自适应地设置阈值，避免了固定阈值带来的识别误差。

采用本文方法能较好地区别具有波形相似性的微震记录与不规则噪声干扰。但在实际资料中，除背景噪声外，还可能包括一些相干噪声，如反射波和折射波及续至波等，这些噪声会对互相关时差估计和参考道产生影响，导致对微震事件识别困难及对到时拾取产生误差，需要通过分析压裂区域微震事件的初至曲线特征来进一步研判和去除错误的拾取结果。