一种基于dwt_sta/lta的含噪信号p波初至峰度拾取方法
【技术领域】
[0001] 本发明属于信号处理技术领域，尤其是涉及一种矿山微震信号P波初至时刻联合 拾取方法。
【背景技术】
[0002] P波初至拾取是矿山微震识别、震源定位及震源机制解释中最为基础和重要的一 步。目前主要采用人工拾取P波到时，但其易受个人经验、情绪等因素影响，且拾取数据量 大时非常耗时。鉴于上述原因，国内外提出了一系列P波初至自动拾取方法，但到目前为止 还没有任何一种算法能够在不同的震源环境、震中距范围、噪声背景下进行微震初至的一 致性检测，更不用说对后续震相初至的有效检测和拾取。
[0003] 常见的P波拾取方法包括长短时窗平均值比法、分形维数法、自回归模型、人工神 经网络法、高阶统计量法等。长短时窗平均值比法利用短时窗平均值与长时窗平均值之比 来反映信号幅值的变化，并取比值大于某一设定阈值时为P波到时。长短时窗平均值比法 具有算法简单、计算速度快、平均值特征函数多样等特点，但其拾取低信噪比信号困难，且 阀值较大时，则可能拾取不到P波初至或增大拾取误差；阀值较小时，则可能过早的被噪音 触发。分形维数法认为噪音与信号叠加时分形维数发生变化，并以此作为P波拾取的依据。 该算法抗噪性能好，拾取精度较高，但其对插值的准确性、时窗和步长的依赖性很高，稍有 不慎就会严重影响拾取结果，且其计算速度较慢。自回归模型将微震信号分为两个局部统 计时段，并取自相关最小值点作为P波到时。该算法拾取精度较高，但其对低信噪比和尾部 震荡信号拾取不稳定，且其采用最小值点作为P波到时，因此对纯噪音信号仍会拾取到时 (HaijiangZhang等，2003)。人工神经网络法采用峰值振幅、时窗内均方根振幅比、峰值与 其前后峰值的包络斜率及噪声与信号的比值等作为神经网络输入因素，人工神经网络法是 一种较为综合的P波初至拾取方法，对微震波形适应性强，但其参数计算工作量大，学习速 度慢，实现复杂。高阶统计量法（PAI-S/K法）由Saragiotis等（2002,2004)提出，取峰度 (或偏度）最大值点作为P波初至时刻，该算法借助于信号只包含噪音时，峰度值和偏度值 趋于0,而当P波到达时峰度和偏度值开始增加，并当滑动窗口内包含部分微震信号时达到 最大值。该算法拾取精度高，但未考虑尖刺和尾部震荡信号对拾取的影响，由此可能产生较 大的拾取误差，且其采用最大值拾取P波到时，对纯噪音信号仍会拾取到时。
[0004] 可见现有的矿山微震信号P波初至拾取方法存在很大的局限，需要研究一种拾取 精度高、稳定性好的自动拾取方法。
【发明内容】
[0005] 本发明所要解决的技术问题是提供一种改进的含噪信号P波初至峰度拾取方法， 解决峰度法（PAI-K法）拾取低信噪比、刺突、尾部震荡及纯噪音信号不稳定的技术问题，该 矿山微震信号P波初至拾取方法适用性强、准确性高。
[0006] 针对矿山微震信号P波初至PAI-K法拾取精度较高，但对低信噪比、刺突、尾部震 荡及纯噪音信号拾取不稳定的技术问题，本发明提出了W-S/L-K拾取法，该方法借助于小 波分解高频细节信号能够较好的保留P波初至信息以及STA/LTA法能很好地克服刺突和尾 部震荡的影响并确定P波大致到时，进而采用权重分析排除上述干扰，极大地降低了错误 拾取率，增强了P波拾取的准确性。
[0007] -种基于DWT_STA/LTA的含噪信号P波初至峰度拾取方法，包括以下几个步骤：
[0008] 步骤1 :提取微震信号X(n)
[0009] 从矿山微震信号中提取待拾取信号X(n)，n= 1，2,…，N，其中，N为所述微震信号 的采样点个数，取N= 400~600,微震信号采样频率600Hz;
[0010] 步骤2 :对微震信号采用离散小波变换DWT得到细节信号D1,i= 1，2, 3, 4 ;
[0011] 步骤3 :依据STA/LTA算法按照以下公式确定细节信号P波大致到时k1;
[0013] 其中，STA(k，D1)和LTA(k，D1)分别为短、长时窗的振幅平均值；x(n，D1)为细节信 号D1的第n个采样信号；k为采样的第k个点，k=WWA，Wm+1，…，N;WSTJPWm分别为短、 长时窗的长度，分别取10和60个米样点；
[0014] 若细节信号DJt应的短时窗和长时窗的振幅平均值的比值A(k，D1)的最大值大 于设定阀值a，则取首个触发点对应的采样点序号为该细节信号P波大致到时ki，否 贝1J，记该细节信号P波大致到时k 0。其中，i= 1，2, 3, 4 ;
[0015] 步骤4 :计算细节信号的P波大致到时Ic1= 0的个数L，i= 1，2, 3, 4 ;
[0016] 步骤5 :利用权重分析计算P波初至准确到时；
[0017] 若步骤4中L>2,则认为该信号为纯噪音信号，反之则采用以下权重公式计算P波 准确到时k' ：
[0019] 其中，1/8、3/8、3/8和1/8分别为04、03、0 2和01的权重系数洫1和1^1'分别为1?1的 最大振幅和局部峰度拾取到时，R1SD1在区间[ki-b^i+b]的信号，b为确定局部峰度拾取 区间的常量，取b= 10~15 ;
[0020]
计算X(n,Di)的滑动峰度值K(k，Di)，k为采样的第k个点， 动时窗的长度，取M= 20~30,并取K(I^D1)最大值点作为局部峰度拾取到时k/ ;
[0021] 步骤6 :依据步骤5输出拾取结果，该信号为纯噪音信号或者P波初至准确到时为 k'。
[0022]所述设定阈值a= 2~3。
[0023] 利用MalIat算法快速实现离散小波变换。
[0024] 并采用低通滤波器h和高通滤波器g获得原始信号c。（等同于本发明中的X(n)) 的小波系数{山，d2，…，七}，从而将微震信号X(n)分解得到细节信号D1 (取i= 1，2, 3, 4)。
[0027] 式中：c]+lik为第k个采样点的近似信号值c]+1，d]+1，k为第k个采样点的细节信号 值dj+1，k为采样的第k个点，1为滤波器确定指标，j为尺度参数（j= 0, 1，…，J-l)，J为 最大分解层数。为简化说明，记Cj+1= {(^+1』}和1+1= {d_j+lik}。
[0028] 设定Ici= 0时，A;= 0,k/ = 0,i= 1，2, 3, 4,简化计算过程，提高计算速度。
[0029] 有益效果
[0030] 本发明提供了一种基于DWT_STA/LTA的含噪信号P波初至峰度拾取方法，包括如 下步骤：提取矿山微震信号；离散小波变换（DWT)得到微震细节信号；计算各细节信号的 STA/LTA值；判定STA/LTA最大值是否大于阀值a，若是则取首个触发点对应的采样点序号 为该细节信号的大致到时，否则记该细节信号到时为0 ;再对到时非0的细节信号进行局部 峰态拾取；最后采用权重分析确定P波最终到时。本发明拾取纯噪音信号时，细节信号STA/ LTA触发数通常小于2,设定STA/LTA触发数大于等于2即可排除噪音信号；拾取低信噪比 信号时，STA/LTA对含噪音较大的细节信号无拾取，且细节信号无拾取的数目通常不大于 2,因此可用其余细节信号拾取初至，从而排除低信噪比的影响；通常刺突信号对STA/LTA 值影响较小，设定合适的STA/LTA触发值即可排除刺突的影响；STA/LTA通常在P波初至时 已触发，不会受尾部震荡的影响。由此，通过一种基于DWT_STA/LTA的含噪信号P波初至峰 度拾取方法有效排除了低信噪比、刺突、尾部震荡及纯噪音信号的干扰，极大地降低了错误 拾取率，增强了P波拾取的准确性。
【附图说明】
[0031] 图1是本发明所述方法流程图；
[0032] 图2是典型波形W-S/L-K法实现过程图，其中，（a)为矿山微震原始信号x(n)(n =1，2,…，400)，虚线为人工拾取P波初至点；(b)为原始信号对应的峰度值K(k)(k= 21，22,…，400)，括号中数值分别为PAI-K法拾取P波初至点及其对应的峰度值；（c)为小 波分解所得细节信号D1Q= 1，2, 3, 4) ;(d)为细节信号D1Q= 1，2, 3, 4)的STA(I^D1)/ LTA(I^D1)值A(I^D1) ;(e)为细节信号局部峰态拾取k/(i= 1，2,3,4);
[0033] 图3是PAI-K和W-S/L-K法工程拾取效果对比图，其中，（a)为PAI-K拾取误差图； (b)为W-S-L-K拾取误差图；
[0034] 图4是PAI-K和W-S/L-K法拾取典型案例图，其中，虚线、点线和实线分别对应人 工、PAI-K法和W-S/L-K法拾取到时，（a)为含刺突纯噪音信号，（b)为低信噪比信号，（c) 为起震明显信号，（d)为低信噪比且含刺突信号，（e)为极低信噪比、含刺突且尾部存在震 荡信号，（f)为起震不明显且尾部存在震荡的低信噪比信号。
【具体实施方式