一种微震震相到时自动拾取方法

，，，

(1.山东科技大学 计算机科学与工程学院，山东 青岛 266590；2.山东科技大学 山东省智慧矿山信息技术省级重点实验室，山东 青岛 266590)

微震的初至到时拾取对实现海量微地震数据的自动处理有重要意义[1]。由于微震监测环境复杂，通常采集到的微震信号具有信噪比较低、非平稳、非线性等特征。受煤矿井下机械振动的影响，微震信号常常夹杂大量的振动噪声，这些噪声对微震震相的初至到时拾取有很大干扰。因此，在机械振动噪声干扰下研究微震震相初至的快速准确拾取具有重要意义。

目前震相初至到时的自动拾取方法主要包括：长短时均值比方法(STA/LTA)、AIC法、PAI-S/K法、Hausdoff分形法等。其中，Stevenson等[2-4]提出了基于时间域的长短时均值比法(STA/LTA)； Akaike等[5-9]依据地震波到达前后数据统计的差别，提出了AIC准则；Saragiotis等[10]基于地震波形的偏斜度和峰度提出了PAI-S/K方法；Chang等[11]提出了基于Hausdorff分形维识别震相的方法；Coppens[12]提出了在不同时窗内进行能量比较的方法；王继等[4]应用单台Akaike信息准则(AIC)和多台最小二乘互相关方法，发展了震相自动精确检测技术，实现了震相初至的自动拾取； 马强等[13]综合STA/LTA方法及AIC 准则，提出了多步骤P波自动拾取方法；刘劲松等[14]通过分析STA/LTA，AIC，PAI-S/K等几种方法的原理及特点，提出了移动时窗峰度的快速算法和改进的峰度拾取初至算法。STA/LTA方法在拾取震相到时过程中受信噪比的影响较大，不易拾取低信噪比信号的准确到时；而AIC算法拾取到时没有选取时窗大小的规则，不易确定合适的窗长；分形维识别法、PAI-S/K方法、能量比法等都没有对低信噪比信号进行过多的研究。

由于煤矿井下环境复杂，采集到的信号还包含有大量机械振动噪声，而上述拾取方法虽具有一定的抗噪能力，但对于低信噪比微震信号不能很好的拾取有效到时。目前工程上对于低信噪比微震信号的震相拾取主要是在专业软件的支持下进行人工拾取，拾取工作效率低，且震相初至的判识精度取决于操作人员的经验。为了解决低信噪比对到时拾取精度影响的问题，提出了一种基于低信噪比微震信号震相初至自动拾取方法。该方法首先应用小波阈值降噪法对采集到的微震信号进行降噪预处理，消除部分外部噪声对到时拾取精度的干扰；然后用STA/LTA获取信号的粗略到时，以该到时为基础选取合适的时间窗口，在该事件窗口内使用AIC方法计算震相的精确到时，从而实现低信噪比微震信号的震相初至的自动拾取。

1 震相拾取原理

本文的震相拾取方法是基于小波阈值降噪法的多方法联合拾取，涉及到的震相识别方法有STA/LTA法及AIC方法。下面介绍两种识别震相的基本理论基础。

1.1 长短时平均(STA/LTA)方法

STA/LTA 方法是由 Stevenson 提出[2-4]，并应用于地震初至到时的判别，后来Allen 等对该方法进行了改进并用于检测地震事件和震相识别[15-17]。

STA/LTA是通过STA(short-term average)和LTA(long-term average)的比值，反映信号的幅度、频率等特征的变化。其工作原理图如图1所示，当地震信号到达时，由于短时窗STA内幅值平均值变化比长时窗LTA内变化剧烈，此时STA与LTA的比值会有一个突变，当其比值大于某一设定阈值λ时，则判定为有效信号[18]，且此剧变点被判定为震相初至点。

STA/LTA拾取到时，是通过计算长短时窗特征函数的均值得到的STA、LTA的值，并通过STA和LTA的比值得到阈值。计算公式为:

(1)

(2)

(3)

其中，公式(1)～(3)中的n指微震信号的采样时刻，sw为短时窗的长度，lw是长时窗的长度，α为设定的触发阈值,CF(i)是在时刻i关于微震信号的特征函数值，表示微震数据的振幅、能量或其变化。

图1 STA/LTA原理Fig.1 STA/LTA Principle

1.2 池赤准则(AIC)法

AIC方法的基本原理是求取地震信号AIC函数的局部最小值，Sleeman[19]提出了AR-AIC 准则，根据自回归过程将地震波形数据分成2个局部统计时段(图2)，AR-AIC函数表示为:

图2 AR-AIC模型Fig.2 AR-AIC Model

(4)

其中：k为两个局部统计时段分界点;p为AR过程阶数;l为地震波形数据长度;σ1,max、σ2,max分别为2个局部统计时段的拟合误差;C为一个常数。为了求出震相初至，必须求出该函数中AR模型的阶数和系数，该方法计算复杂度较高，不利于震相初至的实时拾取。

不同于AR-AIC模型，Maeda提出直接由地震波形数据计算AIC函数[20]，求取AIC函数的局部最小值(图3)，该值对应的位置即为震相初至,AIC函数表示为

AIC(n)=n·log{var(x[1,n])}+(L-n-1)·log{var(x[n+1,L])}

。

(5)

其中n的范围为数据窗口内所有的地震采样点,x(i)(i=1,2,…L)为地震波形离散数据。

图3 使用AIC法拾取震相到时Fig.3 Onset time picking of seismic phase with AIC

通过以上2种AIC方法识别信号的对比，可见后者不需要计算AR模型的阶数，即可直接求取AIC值。后者能满足震相初至拾取的实时性较高的要求，但是该算法需要在震相初至的附近寻找一个合适的时窗来计算AIC值，这是因为不同的时窗可能使AIC函数局部最小值出现的位置不同。为同一个微震波形数据在不同时窗内的波形，由于时窗设置不合理导致震相初至拾取错误。因此，如何选择合理的时窗是AIC法拾取准确震相初至的关键问题之一。

由于采集的微震信号包含有大量机械振动噪声，而机械振动噪声对微震震相到时的拾取是一项极大的挑战，因此在进行到时拾取之前要进行机械振动噪声的滤除。

1.3 小波阈值降噪法

为压制机械噪声对微震震相到时拾取的影响，使用小波阈值降噪法对低信噪比微震信号进行降噪预处理。

小波阈值降噪[21]是在小波变换的基础上发展起来的一种降噪方法。小波变换对于信号具有良好的去相关性，在实际环境下，有意义的信号一般是比较平稳的或者低频信号，但是含噪信号则为高频。基于含有噪声的信号在小波域的分布特性，可以通过将时域信号转换到小波域内，通过设定适当的阈值滤除噪声，从而达到对含有噪声的信号进行降噪处理的目的。

小波阈值降噪法可以采用硬阈值或软阈值的处理方法对小波系数做阈值处理[22]。硬阈值是将指定阈值和信号的绝对值对比，把绝对值小于或等于阈值的置为零，对于绝对值大于阈值的信号则不变。软阈值是将指定阈值和信号的绝对值对比，将小于或等于阈值的置为零，大于阈值的信号，设置为本身与阈值相减，这样数据点就会向零靠拢。硬阈值、软阈值的公式定义分别为公式(6)和公式(7)所示：

(6)

(7)

(8)

sgn(ωi,j)为符号函数，即：

(9)

2 微震震相初至拾取过程

2.1 算法流程

该微震震相初至拾取过程由三部分组成，分别为降噪、大致的到时拾取、精确拾取。

进行第一部分降噪处理时，采用了小波阈值降噪法，对要进行到时拾取的信号进行预处理，消除大部分噪声。第二部分进行大致到时拾取，用的是STA/LTA方法。由于STA/LTA方法拾取到时的局限性，单纯的用此方法得不到准确的结果，将得到的到时作为下一步处理的一个中间值。第三部分为到时的精确拾取，在前两部分的基础上，用AIC方法进行到时拾取，AIC的窗长根据第二部分得到的中间值确定。

微震信号到时拾取流程如图4。

图4 微震信号到时拾取流程Fig.4 Onset time picking process of microseismic signal

2.2 微震震相初至拾取算法描述

算法具体步骤如下：

Step 1: 获取初始微震信号x(t)；

Step 2: 用小波阈值对信号x(t)进行降噪预处理，得到降噪后的信号x′(t)；

Step 4: 利用设置好的STA/LTA算法对x′(t)进行到时拾取，得到一个大致到时t1；

(10)

当公式(10)的结果首次大于δ时，获取t1的值；

Step 5： 在STA/LTA算法拾取到时的基础上，确定AIC算法的时间窗长，以t1为对称点，向前向后分别取500个采样点；

L=[t1-500,t1+499]，

(11)

Step 6: 利用AIC算法对Step5中时窗内的信号进行微震信号的拾取，得到的到时为t

AIC(x′(t))=t1*log{var(x′[1,t1])}+(L0-t1-1)*log{var(x′[t1+1,Ln])}。

(12)

3 实验

实验数据来源于微震监测系统，从监测数据中随机选取一组含噪微震信号进行到时拾取实验。

信号采集频率设置为1 000 Hz，选取的采样点数为4 000，该组信号的人工拾取到时点设置在1 229 ms处。原信号的波形图与时频谱如图5。

图5 原信号波形及时频谱Fig.5 Waveform and frequency spectrum of the original signals

图6为用STA/LTA以及AIC方法拾取的到时结果图。用STA/LTA进行到时拾取时，选取的短时窗为150 ms，长时窗为900 ms，设置的阈值为4，拾取的到时点为1 425 ms。如图6所示，AIC方法理论拾取到时为1 231 ms，但是由于时窗选取太长，结果中的最低点可能不止一个。而且由于噪音的影响，STA/LTA、AIC方法拾取波形时，在拾取到时点附近分辨率不高，不能确定到时点。

图6 不同方法对原信号进行震相拾取Fig.6 Seismic phase picking of the original signals by different methods

针对上述出现的拾取到时不准确以及分辨率不高的情况，提出了基于小波阈值降噪的低信噪比微震震相初至拾取方法。首先用小波阈值法对这组含噪的信号进行降噪预处理，采用无偏风险估计原则(rigrsure)[23]，设置小波分解层数为4层[24]，选取的小波基函数为Symlets小波系的sym8[25]。得到降噪后的信号波形图(如图7)。

图7 降噪波形及频谱图Fig.7 De-noising waveform and frequency spectrum

用STA/LTA法对降噪后的信号进行到时拾取，短时窗选取为150 ms，长时窗为900 ms，阈值设置为4，得到一个大致的到时，如图8(b)所示。图8(c)为用STA/LTA方法确定的时窗，该时窗是(b)中拾取到的到时(1 325±500) ms所得，即825～1 825 ms之间。

图8 STA/LTA方法拾取降噪波形并确定时窗Fig.8 STA/LTA method of de-noising waveform picking and determination of the time window

图9 本文方法拾取的到时Fig.9 Seismic phase picking of signals by the model

图9为本文模型拾取的到时，即将图8(c)得到的时窗作为AIC方法的窗长进行到时拾取，得到最低点即到时点为1 231 ms。

表1为该实验的性能分析，由本文方法拾取到时较传统的STA/LTA方法运行速度平均降低了34.28%；相对于AIC方法运行平均速度提升了6.65%。

对50组微震信号进行实验分析，选取效果较好的实验结果进行不同方法拾取到时准确率比对。以人工拾取到时为准，误差±10 ms作为本文所有拾取到时的准确拾取，统计对比结果如表2所示。通过对比分析可见，本文方法拾取准确率比传统的STA/LTA方法平均提高9.78个百分点，比传统的AIC方法平均提高10.08个百分点。

表1 一组信号不同方法拾取到时性能对比Tab.1 Performance comparison of different onset time picking methods for a set of signals

表2 不同方法拾取到时准确率对比Tab. 2 Accuracy comparison of different onset time picking methods

4 结论与未来工作

通过对比STA/LTA、AIC以及本文方法在低信噪比环境下拾取微震震相拾取的准确率，得出如下结论：STA/LTA算法简单，拾取速度快，但在噪声干扰下震相到时拾取误差较大；AIC方法具有较强的抗噪能力，且能准确拾取到时，但在时窗过长情况下算法运行时间过长，不利于震相到时的实时拾取。本文方法使用STA/LTA法确定震相的大致到时，并以此为基础为AIC算法确定一个尽可能小的时窗，大大缩短了AIC算法迭代的次数，不仅具有较强的抗噪能力，而且震相到时的拾取准确率也较高。

使用小波阈值去噪能够降低机械振动噪声对微震信号的干扰，提高了震相到时拾取的准确率，但由于微震信号采集环境复杂，尚未对其他类型的噪声进行深入研究，但初步实验表明本文方法能够对随机非平稳噪声进行有效压制，下一步我们将进一步改进小波阈值降噪算法，使其适应更多类型的噪声压制，进而提高震相初至拾取的精度及实时性。

[1]张文东,马天辉,唐春安,等.锦屏二级水电站引水隧洞岩爆特征及微震监测规律研究[J].岩石力学与工程学报,2014,33(2):339-348.

ZHANG Wendong,MA Tianhui,TANG Chun’an,et al.Research on characteristics of rockburst and rules of microseismic monitoring at diversion tunnels in Jinping II hydropower station[J].Chinese Journal of Rock Mechanics and Engineering,2014,33(2):339-348.

[2]田优平.近震P波震相自动识别方法研究[D].北京：中国地震局地球物理研究所,2015.

[3]ANDREW H RANKIN.作为赋存于花岗岩中钨锡矿化勘探指南的流体包裹体异常-前景展望[J].岩石学报,2007,23(1):3-14.

RANKIN H A.Fluid inclusion anomalies as exploration guides for granite hosted Sn-W mineralization:Prospects for the future[J].Acta Petrologica Sinica,2007，23(1):3-14.

[4]王继,陈九辉,刘启元,等.流动地震台阵观测初至震相的自动检测[J].地震学报,2006,28(1):42-51.

WANG Ji,CHEN Jiuhui,LIU Qiyuan,et al.Automatic onset phase picking for portable seismic array observation[J].Acta Seismologica Sinica (English Edition),2006，28(1):44-53,115.

[5]AKAIKE H.Information theory and an extension of the maximum likelihood principle[C]//2nd International Symposium on Information Theory.Tsahkadsor,1971:267-281.

[6]贾瑞生,谭云亮,孙红梅,等.低信噪比微震P波震相初至自动拾取方法[J].煤炭学报,2015,40(8):1845-1852.

JIA Ruisheng,TAN Yunliang,SUN Hongmei,et al.Method of automatic detection on micro-seismetic P-arrival time under low signal-to-noise ratio[J].Journal of China Coal Society,2015,40(8):1845-1852.

[7]LEONARD M,KENNETT B L N.Multi-component autoregressive techniques for analysis of seismograms[J].Physics of the Earth Planetary Interiors,1999,113(1-4):247-264.

[8]SLEEMAN R,ORILD V E.Robust automatic P-phase picking:An online implementation in the analysis of broadband seismogram recordings[J].Physics of the Earth planetary Interiors,1999,113(1/2/3/4):265-275.

[9]MAEDA N.A method for reading and checking phase times in autoprocessing system of seismic wave data[J].Journal of Seismological Society of Japan,1985,38(3):365-397.

[10]SARAGIOTIS,C.HADJILEONTIADIS L,PANAS S M.PAI-S/K:A robust automatic seismic P phase arrival identification scheme[J].IEEE Transactions on Geosciences and Remote Sensing,2002,40(6):1395.

[11]常旭,刘伊克.地震记录的广义分维及其应用[J].地球物理学报,2002,45(6):839-846.

CHANG Xu,LIU Yike.The generalized fractal dimension of seismic records and its applications[J].Chinese Journal of Geophysics,2002,45(6):839-846.

[12]COPPENS F.First arrival picking on common offset trace collections for automatic estimation of static correction[J].Geophysical Prospecting,1985,33(8):1212-1231.

[13]马强,金星,李山有,等.用于地震预警的P波震相到时自动拾取 [J].地球物理学报,2013,56(7):2313-2321.

MA Qiang,JIN Xing,LI Shanyou,et al.Automatic P-arrival detection for earthquake early warning[J].Chinese Journal of Geophysics,2013,56(7):2313-2321.

[14]刘劲松,王赟,姚振兴.微地震信号到时自动拾取方法[J].地球物理学报,2013,56(5):1660-1666.

LIU Jinsong,WANG Yun,YAO Zhenxing.On micro-seismic first arrival identification:A case study[J].Chinese Journal of Geophysics,2013,56(5):1660-1666.

[15]刘晗,张建中.微震信号自动检测的STA/LTA算法及其改进分析[J].地球物理学进展,2014,29(4):1708-1714.

LIU Han,ZHANG Jianzhong.STA/LTA algorithm analysis and improvement of microseismic signal automatic detection[J].Progress in Geophysics,2014,29(4):1708-1714.

[16]吴治涛,李仕雄.STA/LTA算法拾取微地震事件P波到时对比研究[J].地球物理学进展,2010,25(5):1577-1582.

WU Zhitao,LI Shixiong.Comparison of STA/LTA P-pickers for micro seismic monitoring [J].Progress in Geophysics,2010,25(5):1577-1582.

[17]白添羊,吴顺川,王进进,等.岩石破裂声发射压缩波到时拾取方法及其优化改进研究[J].岩石力学与工程学报,2016,35(9):1754-1766.

BAI Tianyang,WU Shunchuan,WANG Jinjin,et al.P-onset picking methods of acoustic emission compression waves and optimized improvement[J].Chinese Journal of Rock Mechanics and Engineering,2016,35(9):1754-1766.

[18]姜福兴,尹永明,朱权洁,等.单事件多通道微震波形的特征提取与联合识别研究[J].煤炭学报,2014,39(2):229-237.

JIANG Fuxing,YIN Yongming,ZHU Quanjie,et al.Feature extraction and classification of mining microseismic waveforms via multi-channels analysis[J].Journal of China Coal Society,2014,39(2):229-237.

[19]SLEEMAN R,ORILD V E.Robust automatic P-phase picking:An online implementation in the analysis of broadband seismogram recordings[J].Physics of the Earth Planetary Interiors,1999,113(1/2/3/4):265-275.

[20]MAEDA N.A method for reading and checking phase times in autoprocessing system of seismic wave data[J].Journal of Seismological Society of Japan,1985,38(3):365-379.

[21]DONOHO D L.Denoising by soft thresholding[J].IEEE Transactions on Information Theory,1995,41(3):613-627.

[22]王宏强,尚春阳,高瑞鹏,等.基于小波系数变换的小波阈值去噪算法改进[J].振动与冲击,2011(10):165-168.

WANG Hongqiang,SHANG Chunyang,GAO Ruipeng,et al.An improvement of wavelet shrinkage denoising via wavelet coefficient transformation[J].Journal of Vibration and Shock,2011(10):165-168.

[23]王维,张英堂,任国全.小波阈值降噪算法中最优分解层数的自适应确定及仿真[J].仪器仪表学报,2009,30(3):526-530.

WANG Wei,ZHANG Yingtang,REN Guoquan.Adaptiveselection and simulation of optimal decomposition level in threshold de-noising algorithm based on wavelet transform[J].Chinese Journal of Scientific Instrument,2009,30(3):526-530.

[24]臧玉萍,张德江,王维正.小波分层阈值降噪法及其在发动机振动信号分析中的应用[J].振动与冲击,2009,28(8):57-60.

ZANG Yuping,ZHANG Dejiang,WANG Weizheng.Wavelet hierarchical threshold de-noising method and its application in vibration signal analysis of engine[J].Journal of Vibration and Shock,2009,28(8):57-60.

[25]钟建军,宋健,由长喜,等.基于信噪比评价的阈值优选小波去噪法[J].清华大学学报(自然科学版),2014,54(2):259-263.

ZHONG Jianjun,SONG Jian,YOU Changxi,et al.Wavelet de-noising method with threshold selection rules based on SNR evaluations[J].Journal of Tsinghua University (Science and Technology),2014,54(2):259-263.