水力压裂微地震监测技术是一项通过观测注水压裂过程所诱发的微弱地震事件来监测裂缝发育状况的技术，微地震事件初至拾取是该技术的一项关键处理步骤.由于目前常用的震源定位算法大都需要微地震事件的P波、S波到时信息，而人工拾取微地震事件初至又是一项十分费时费力的工作，因此，研究一种能够自动并准确提取微地震事件初至到时的方法是非常必要的.

微地震事件初至拾取是以地震信号和环境噪声的差异识别为基础(如振幅、频率成分、偏振、统计特征以及波形相关性等).目前，在微地震数据处理中常用的初至拾取算法主要借鉴了天然地震研究中的一些方法，例如STA/LTA(Ambuter and Solomon，1974; Stevenson，1976; Allen，1978，1982; McEvilly and Majer，1982; Baer and Kradolfer，1987; Earle and Shearer，1994; Withers et al.，1998; Chen，2005; 吴治涛和李仕雄，2010)、偏振分析(Defl and re and Dubesset，1992; Moriya and Niitsuma，1996; Anant and Dowla，1997; Soma et al.，2004; Moriya，2008，2009; 吴治涛等，2012)及AR-AIC方法(Takanami and Kitagawa，1988，1991，1993; Sleeman and van Eck，1999; Leonard and Kennett，1999; Leonard，2000; Zhang et al.，2003; 宋维琪和吕世超，2010)等.STA/LTA方法根据环境噪声振幅小、频带宽，而地震信号振幅大、频带窄的特点，通过计算一长一短两个滑动时窗内地震记录特征函数的平均值之比作为拾取地震波初至到时的依据；偏振分析方法利用了地震信号的振动轨迹近似线性，而环境噪声的振动轨迹无显著方向性这一特征差异，通过计算滑动时窗内局部信号的偏振度识别并提取地震波初至；AR-AIC方法是假设初至到时前后的两段信号可分别采用不同的AR模型表示，利用Akaike信息准则描述这两段信号与AR模型的匹配程度，认为AIC值最小点为地震波初至点.实践表明，实际微地震监测资料中存在大量信噪比较低的微地震事件，利用上述方法对这些事件进行初至拾取无法取得令人满意的效果.

研究表明，任何一种只依赖地震信号与环境噪声某一类特征差异的初至拾取方法都很难取得比较理想的效果(刘希强等，2009)，因此，可考虑通过综合地震信号与环境噪声的多种特征差异设计更加准确、可靠的算法.许多学者已做过相关的研究工作，并提出了不同的方法.Bai和Kennett(2000，2001)详细研究了地震信号与环境噪声在能量、偏振、瞬时频率以及自回归模型方面存在的特征差异；王继等(2006)提出了一种基于STA/LTA和AIC的两步法震相到时自动检测方法；刘希强等(2009)通过综合地震信号与背景噪声在能量、非高斯性、非线性以及偏振特征方面存在差异，提出了一种信噪综合差异特征量方法，简称EFGLP方法；赵大鹏等(2012)提出了一种以地震记录的峰度作为特征函数计算AIC值的Kurtosis-AIC方法；张唤兰等(2013)将基于能量比和AIC的两步法初至到时检测方法应用于微地震数据处理；吕世超等(2013)提出了一种利用偏振约束的能量比微地震初至拾取方法.

本文首先简要叙述了STA/LTA、偏振分析及AR-AIC方法的基本原理和实现步骤，然后在借鉴上述三种方法原理的基础上，设计了一种针对低信噪比微地震事件的初至拾取方法，称为SLPEA算法.该方法的基本思路是根据数据的信噪比构造不同检测函数来提取地震波初至到时.为了检验本文方法的可行性和可靠性，笔者分别对合成数据和实际资料进行了处理，并把处理结果与上述三种方法及手工拾取的结果进行了对比分析.

STA/LTA方法是目前常用的一种简单快速的初至拾取算法，它利用了地震信号与环境噪声的振幅特征差异来拾取地震波初至.STA/LTA比值函数的计算公式为

其中，N₁、N₂为长、短滑动时窗的长度，CF_i为特征函数.不同学者提出了不同形式的特征函数，目的是为了突出地震波到达时信号的变化规律.本文选用如下形式的特征函数计算STA/LTA比值函数：

其中X_i，Y_i，Z_i分别表示地震记录的三个分量.当地震波到达时，短时窗平均值(STA)比长时窗平均值(LTA)变化的更快，因此，二者的比值会出现一个明显的突跳.当R值超过了预先设定的触发阈值时，相应时刻即判定为地震波的初至到时.

地震数据通常采用三分量检波器进行采集，且地震信号与环境噪声的质点振动轨迹存在明显差别，这就为在三分量记录上拾取地震波初至提供了可能.利用偏振分析方法提取地震波初至到时的基本思路为：首先，利用一个滑动时窗内的局部三分量记录构造如下协方差矩阵：

其中E表示求信号序列的期望.求解该协方差矩阵的三个特征值 利用这三个特征值可计算得到记录的偏振度曲线：

P的范围为0~1，其中P=1表示信号为线性偏振(如地震信号)，P=0表示信号为圆偏振(如随机噪声).偏振分析方法选取P极大值点前曲线斜率最大的点作为地震波的初至点.

AR-AIC方法根据地震波初至前后波形数据的统计特征差异提出.该方法假设包含地震波初至的一段记录按初至点位置可划分为信号部分和噪声部分，从这两部分中各选取一小段记录，并利用AR模型将其表示为如下形式：

其中，X¹_i和X²_i分别为信号和噪声序列，M为AR模型阶数，a¹_j和a²_j为AR模型系数，ε¹_i和ε²_i为白噪声序列.(5)式和(6)式分别表示记录的信号模型和噪声模型.AIC值描述了记录中其余部分与这两个模型的匹配程度，其定义式为

其中，k为滑动分割点，它将原记录划分为信号部分和噪声部分；L为(对数形式)最大似然函数，其表达式为

AIC值最小点(或L最大点)表示由该点划分的两部分记录与信号、噪声模型的拟合度最佳，因此，该点被认为是两段稳定时间序列的分界点，即地震波的初至点.

上述三种方法均为利用地震信号与环境噪声的某一类特征差异进行初至拾取，实践表明在低信噪比条件下很难取得比较理想的效果.Bai和Kennett(2000，2001)认为通过综合地震信号与环境噪声的多种特征差异进行初至拾取将有助于降低噪声的影响并提高结果的准确性.根据这一思想，笔者在借鉴上述三种方法原理的基础上，设计了一种新的微地震初至拾取SLPEA算法.该算法的具体实现步骤如下：

(1)对于识别出的有效微地震事件，首先选取长度分别为N₁和N₂(N₁＞N₂)的两个滑动时窗，利用(1)式计算各道记录的STA/LTA比值函数R_i.然后，选取长度为N₃的一个滑动时窗，利用(4)式计算各道记录的偏振度函数P_i.

(2)对R_i和P_i进行归一化处理后将二者相乘，可得震相检测函数K_i.R_i和P_i的归一化过程可利用下式表示(以R_i为例)：

(3)选取长度为N₄的一个滑动时窗，计算如下边缘检测函数D_i：

D_i反映了K_i变化的快慢程度，它在K_i变化最快的点(即斜率最大的点)处存在一个局部极大值，因此，D_i的极值点位置比K_i更加接近地震波初至点.

(4)选取D_i最大的两个局部极值点，它们所对应的时刻可视为P波和S波初至到时的初始估计值，分别记为T′_P和T′_S.重复采用上述步骤对每道记录进行处理后可得到微地震事件P波和S波的初至走时曲线.对于信噪比较低的微地震事件，其初至走时曲线可能存在奇异值，因此，需要对它进行平滑处理，平滑后的初至到时分别记为T″_P和T″_S.计算T′_P与T″_P、T′_S与T″_S之间的绝对误差，记为ΔT_P、ΔT_S.

(5)给定误差阈值C₁，并判断每道记录的ΔT_P和ΔT_S是否超过该阈值.如果是，则在该道记录的边缘检测函数D_i中找到与T″_P(或T″_S)最接近的极值点位置，并对T′_P(或T′_S)进行更新；否则，保持T′_P(或T′_S)不变.

(6)对于每一道记录，分别以校正后的T′_P和T′_S为中心，向前(长度为N₅)、向后(长度为N₆)延拓一个短时窗，并利用下式计算短时窗内记录的信噪比S：

其中A_noise和A_signal+noise分别表示T′_P(或T′_S)前、后长度为N₆的两段记录的均方根振幅.

(7)截取短时窗内的部分地震记录X′_k，Y′_k，Z′_k及其对应的边缘检测函数D′_k，D′_k最大值点所对应的时刻记为T_D.利用(8)式分别计算X′_k，Y′_k，Z′_k的最大似然函数L^x_{k，L^yk，L^zk}，且

(8)给定信噪比阈值C₂，并判断S是否大于该阈值.如果是，则以T_L作为最终的初至拾取结果；否则，执行下一步.

(9)给定误差阈值C₃，并判断T_D与T_L之间的绝对误差是否小于该阈值.如果是，则以T_D与T_L的均值作为最终的初至拾取结果；否则，选取T_F作为最终结果.

笔者首先利用模型数据检验本文方法的可行性.图 1a、2a、3a为三组合成数据，它们均由一个按指数规律衰减的正弦波(模拟地震子波)与不同幅度的随机噪声叠加而成，其信噪比分别为12.4 dB、3.8、0.8 dB，时间采样间隔为0.0005 s.图 1b、2b、3b为不含噪声的合成记录，其中红色竖线表示地震波的真实初至到时0.2 s.图 1(c—h)、2(c—h)、3(c—h)分别为拾取初至到时的过程中计算得到的各类检测函数，所采用参数的取值见表 1.图中函数最大值点的位置用红色星号标出，其对应的时刻见表 2.

表 1(Table 1)

表 1 本文方法各参数的含义及取值 Table 1 Definitions and values of the parameters in the proposed method 参数 取值 含义 模型数据 实际资料 N1 100 150 用于计算长时平均值LTA的滑动时窗长度 N2 20 30 用于计算短时平均值STA的滑动时窗长度 N3 40 30 用于计算偏振度P的滑动时窗长度 N4 60 80 用于计算边缘检测因子D的滑动时窗长度 N5 无 120 以T′P、T′S为中心向前延拓的时窗长度 N6 无 40 以T′P、T′S为中心向后延拓的时窗长度 C1 无 0.01 用于校正初始到时的误差阈值 C2 10 20 信噪比阈值 C3 0.01 0.005 用于判断最终到时的误差阈值

表 1 本文方法各参数的含义及取值 Table 1Definitions and values of the parameters in the proposed method

表 2(Table 2)

表 2 图1—3中各类检测函数最大值点对应的时刻 Table 2 Time positions of the maximum points of the detector functions shown in Fig.1—3 组号 函数最大值点对应的时刻/s R P K D L Ave.* F 1 0.2095 0.2185 0.2090 0.2020 0.2000 0.2010 0.2005 2 0.2100 0.2190 0.2090 0.2040 0.1975 0.2010 0.2010 3 0.2120 0.2200 0.2190 0.2070 0.3235 0.2655 0.2025 注：Ave.表示函数D和L最大值点对应时刻的平均值.

表 2 图1—3中各类检测函数最大值点对应的时刻 Table 2Time positions of the maximum points of the detector functions shown in Fig.1—3

图 1

Fig. 1

图 1 (a) 第一组合成数据； (b) 不含随机噪声的合成数据； (c) STA/LTA比值函数R； (d) 偏振度函数P； (e) 震相检测函数K； (f) 边缘检测函数D； (g) 最大似然函数L； (h) 初至检测函数F Fig. 1 (a) The first synthetic dataset； (b) Synthetic dataset without random noise； (c) STA/LTA ratio function R； (d) Polarization attribute function P； (e) Phase detector function K； (f) Edge detector function D； (g) Maximized likelihood function L； (h) Arrival detector function F

对于第一组合成数据，通过对比图 1(c—h)可以看出，R、P、K三者最大值所对应的时刻均明显滞后于地震波的初至到时，而D、L、F最大值点的位置则十分接近地震波的真实初至点.由表 2中的数据进一步可知，L最大值点对应的时刻恰好为地震波的真实初至到时0.2 s，而D和F最大值点对应的时刻则与真实到时分别相差了0.002 s和0.0005 s(即4个和1个采样点).该结果表明对于信噪比较高的地震记录(例如信噪比大于10 dB)，根据L得到的初至到时比采用其他几类检测函数得到的结果更加准确.

对于第二组合成数据，由图 2(c—h)以及表 2中数据可知各类检测函数最大值点的位置与第一组合成数据相差不大，主要区别在于L最大值点的位置偏离了地震波的真实初至点，二者相差了0.0025 s(即5个采样点).进一步分析表 2中的数据可知，D和L最大值点对应时刻的平均值与F最大值点对应的时刻相同，且最接近地震波的真实初至到时.该结果表明当D和L最大值点的位置接近时(例如二者的距离小于20个采样点)，可通过计算二者的平均值作为地震波初至到时的估计值，从而无需计算F.

图 2

Fig. 2

图 2 (a) 第二组合成数据； (b) 不含随机噪声的合成数据； (c) STA/LTA比值函数R； (d) 偏振度函数P； (e) 震相检测函数K； (f) 边缘检测函数D； (g) 最大似然函数L； (h) 初至检测函数F Fig. 2 (a) The second synthetic dataset； (b) Synthetic dataset without random noise； (c) STA/LTA ratio function R； (d) Polarization attribute function P； (e) Phase detector function K； (f) Edge detector function D； (g) Maximized likelihood function L； (h) Arrival detector function F

对于第三组合成数据，通过图 3g可以看出，L最大值点的位置与地震波真实初至点的位置相差较大，表明L受环境噪声的影响较大，其初至拾取结果的准确度随着信噪比的降低而明显下降.由表 2中的数据可知，F最大值点的位置最接近地震波的初至点，表明F具有较强的抗噪性，在低信噪比条件下仍然能够较准确地描述地震波初至点的位置.

图 3

Fig. 3

图 3 (a) 第三组合成数据； (b) 不含随机噪声的合成数据； (c) STA/LTA比值函数R； (d) 偏振度函数P； (e) 震相检测函数K； (f) 边缘检测函数D； (g) 最大似然函数L； (h) 初至检测函数F Fig. 3 (a) The third synthetic dataset； (b) Synthetic dataset without random noise； (c) STA/LTA ratio function R； (d) Polarization attribute function P； (e) Phase detector function K； (f) Edge detector function D； (g) Maximized likelihood function L；(h) Arrival detector function F

为了深入研究本文方法的准确性并设置合理的参数取值，笔者对1000组具有不同信噪比的合成数据进行了处理.图 4为数据信噪比与不同检测函数拾取误差(即函数最大值点对应的时刻与真实初至到时之差)的关系图.由该图可知，当数据的信噪比大于22 dB时，除了极个别的点外，L和F的拾取误差均为0，表明利用二者均能够得到准确的初至到时；对于信噪比在10~22 dB范围内的数据，F的拾取误差为1个采样点，而L的拾取误差除少数点外仍然为0；当信噪比降低到10 dB以下时，L的拾取误差显著增大，而F的拾取误差仍然保持在一个可以接受的范围内，其最大值仅为5个采样点.通过上述分析可以得出，对于信噪比大于10 dB的合成数据，根据L得到的初至到时更加准确；而对于信噪比小于10 dB的合成数据，根据F得到的初至到时误差较小.该结论证明了笔者将信噪比阈值C₂设为10 dB(见表 1)是合理的.此外，通过图 4还可以看出根据L得到的初至到时一般早于地震波的真实到时(其拾取误差为正值)，而根据D得到的初至到时则晚于真实到时(其拾取误差为负值)，表明在一定的误差范围内二者的平均值更加接近真实的初至到时.

图 4

Fig. 4

图 4 合成数据信噪比与初至拾取误差关系散点图 Fig. 4 Scatter diagram of SNRs versus picking errors for synthetic datasets

本文采用的实际资料为某油田对一口水平井进行的11段水力压裂施工的微地震监测数据.在采集该资料时所使用的观测系统为一个15级的井下检波器串，其级间距为10 m，时间采样间隔为0.0005 s.实际监测时，该观测系统被布设在压裂井附近的一口直井中进行数据采集.

笔者从该资料中选取了4个具有不同信噪比的有效微地震事件，分别利用本文提出的SLPEA算法、STA/LTA、偏振分析以及AR-AIC方法对其进行了初至拾取，结果分别如图 5—8所示.图中的红色和绿色竖线分别表示P波和S波的初至到时.通过对比图 5—8可知，本文方法的应用效果最好，能够比较准确地拾取出4个微地震事件的P波和S波初至到时；STA/LTA方法与偏振分析方法的应用效果相似，其中前三个微地震事件的初至拾取结果较准确，而最后一个事件中某几道记录的P波初至到时则存在较大误差；AR-AIC方法的应用效果最差，其中每道记录的P波、S波初至到时均存在不同程度的误差，造成这一结果的原因可能是在进行初至拾取之前，笔者对实际资料进行了带通滤波处理，从而对记录的波形信息产生一定程度的破坏.

图 5

Fig. 5

图 5 本文方法初至拾取结果 (a) 事件1； (b) 事件2； (c) 事件3； (d) 事件4； 其中红色和绿色短线分别表示拾取的P波和S波初至. Fig. 5 Arrival picking results of the proposed method (a) Event 1； (b) Event 2； (c) Event 3； (d) Event 4； where the red and green vertical bars denote the picked P- and S-wave arrival times.

图 6

Fig. 6

图 6 STA/LTA方法初至拾取结果 (图注说明同图 5) Fig. 6 Arrival picking results of the STA/LTA method

图 7

Fig. 7

图 7 偏振分析方法初至拾取结果 (图注说明同图 5) Fig. 7 Arrival picking results of the polarization method

图 8

Fig. 8

图 8 AR-AIC方法初至拾取结果 (图注说明同图 5) Fig. 8 Arrival picking results of the AR-AIC method

为了定量对比几种方法的准确性，笔者对4个微地震事件进行了手工初至拾取(结果如图 9所示)，并以手工拾取结果作为参照，与上述四种方法的处理结果进行对比分析.图 10为图 5—8所示的初至拾取结果与手工拾取结果的误差分布直方图，误差的统计特征参数见表 3.由表 3中的数据可知，本文方法处理结果的绝对误差平均值仅为1.33×10^-3s，方差为3.21×10^-6s²，二者均小于其他三组结果的同类参数；误差在±0.005 s范围内的到时个数占总数的95.8%，在4组结果中为最大，证明了本文方法的应用效果优于其他三类方法.图 11为4组结果的绝对误差与记录信噪比的关系图，其中斜线表示对每组结果的误差进行线性拟合后得到的趋势线.通过该图可以看出，随着记录信噪比增大，4种方法的拾取误差均逐渐减小，但在信噪比相同的条件下，本文方法的拾取误差最小.

图 9

Fig. 9

图 9 手工初至拾取结果 (图注说明同图 5) Fig. 9 Manual arrival picking results

图 10

Fig. 10

图 10 初至拾取误差分布直方图 (a) 本文方法； (b) STA/LTA方法； (c) 偏振分析方法； (d) AR-AIC方法. Fig. 10 Histograms of the picking errors (a) The proposed method； (b) STA/LTA method； (c) Polarization method； (d) AR-AIC method.

图 11

Fig. 11

图 11 实际数据信噪比与初至拾取误差关系散点图 Fig. 11 Scatter diagram of SNRs versus picking errors for real datasets

表 3(Table 3)

表 3 初至拾取误差统计特征参数 Table 3 Statistical parameters of the picking errors 名称 绝对误差平均值/s 绝对误差方差/s 2 绝对误差满足以下范围的到时个数占总数的百分比/% ≤0.001 s ≤0.002 s ≤0.003 s ≤0.004 s ≤0.005 s 本文方法 1.33×10 -3 3.21×10 -6 60.0 80.8 90.0 94.2 95.8 STA/LTA方法 3.01×10 -3 5.25×10 -5 32.5 64.2 81.7 89.2 92.5 偏振分析方法 8.12×10 -3 2.82×10 -4 3.33 15.0 36.7 59.2 72.5 AR-AIC方法 8.93×10 -3 6.29×10 -5 23.3 31.7 35.0 38.3 40.0

表 3 初至拾取误差统计特征参数 Table 3Statistical parameters of the picking errors

本文提出了一种针对低信噪比微地震事件的初至拾取方法——SLPEA算法.该方法通过综合地震信号与环境噪声在振幅、偏振以及统计特征等方面存在的差异，根据数据的信噪比构造不同的检测函数来提取地震波初至到时，克服了传统方法抗噪性弱的缺点，提高了结果的准确性.实际资料的处理结果表明，在相同信噪比条件下本文方法的应用效果优于传统方法，得到的初至到时更加接近手工拾取的结果.本文方法的主要缺点是计算量较大、耗时较长，因此，在应用于实时处理方面还有待进一步研究.