CN105223614B - 一种基于dwt_sta/lta的含噪信号p波初至峰度拾取方法 - Google Patents

Abstract

本发明公开了一种基于DWT_STA/LTA的含噪信号P波初至峰度拾取方法，包括如下步骤：提取矿山微震信号；离散小波变换(DWT)得到微震细节信号；计算各细节信号的STA/LTA值；判定STA/LTA最大值是否大于阀值a，若是则取首个触发点对应的采样点序号为该细节信号的大致到时，否则记该细节信号到时为0；再对到时非0的细节信号进行局部峰态拾取；最后采用权重分析确定P波最终到时。该方法有效排除了低信噪比、刺突、尾部震荡及纯噪音信号的干扰，极大地降低了错误拾取率，增强了P波拾取的准确性。

Description

一种基于DWT_STA/LTA的含噪信号P波初至峰度拾取方法

技术领域

本发明属于信号处理技术领域，尤其是涉及一种矿山微震信号P波初至时刻联合拾取方法。

背景技术

P波初至拾取是矿山微震识别、震源定位及震源机制解释中最为基础和重要的一步。目前主要采用人工拾取P波到时，但其易受个人经验、情绪等因素影响，且拾取数据量大时非常耗时。鉴于上述原因，国内外提出了一系列P波初至自动拾取方法，但到目前为止还没有任何一种算法能够在不同的震源环境、震中距范围、噪声背景下进行微震初至的一致性检测，更不用说对后续震相初至的有效检测和拾取。

常见的P波拾取方法包括长短时窗平均值比法、分形维数法、自回归模型、人工神经网络法、高阶统计量法等。长短时窗平均值比法利用短时窗平均值与长时窗平均值之比来反映信号幅值的变化，并取比值大于某一设定阈值时为P波到时。长短时窗平均值比法具有算法简单、计算速度快、平均值特征函数多样等特点，但其拾取低信噪比信号困难，且阀值较大时，则可能拾取不到P波初至或增大拾取误差；阀值较小时，则可能过早的被噪音触发。分形维数法认为噪音与信号叠加时分形维数发生变化，并以此作为P波拾取的依据。该算法抗噪性能好，拾取精度较高，但其对插值的准确性、时窗和步长的依赖性很高，稍有不慎就会严重影响拾取结果，且其计算速度较慢。自回归模型将微震信号分为两个局部统计时段，并取自相关最小值点作为P波到时。该算法拾取精度较高，但其对低信噪比和尾部震荡信号拾取不稳定，且其采用最小值点作为P波到时，因此对纯噪音信号仍会拾取到时(Haijiang Zhang等，2003)。人工神经网络法采用峰值振幅、时窗内均方根振幅比、峰值与其前后峰值的包络斜率及噪声与信号的比值等作为神经网络输入因素，人工神经网络法是一种较为综合的P波初至拾取方法，对微震波形适应性强，但其参数计算工作量大，学习速度慢，实现复杂。高阶统计量法(PAI-S/K法)由Saragiotis等(2002,2004)提出，取峰度(或偏度)最大值点作为P波初至时刻，该算法借助于信号只包含噪音时，峰度值和偏度值趋于0，而当P波到达时峰度和偏度值开始增加，并当滑动窗口内包含部分微震信号时达到最大值。该算法拾取精度高，但未考虑尖刺和尾部震荡信号对拾取的影响，由此可能产生较大的拾取误差，且其采用最大值拾取P波到时，对纯噪音信号仍会拾取到时。

可见现有的矿山微震信号P波初至拾取方法存在很大的局限，需要研究一种拾取精度高、稳定性好的自动拾取方法。

发明内容

本发明所要解决的技术问题是提供一种改进的含噪信号P波初至峰度拾取方法，解决峰度法(PAI-K法)拾取低信噪比、刺突、尾部震荡及纯噪音信号不稳定的技术问题，该矿山微震信号P波初至拾取方法适用性强、准确性高。

针对矿山微震信号P波初至PAI-K法拾取精度较高，但对低信噪比、刺突、尾部震荡及纯噪音信号拾取不稳定的技术问题，本发明提出了W-S/L-K拾取法，该方法借助于小波分解高频细节信号能够较好的保留P波初至信息以及STA/LTA法能很好地克服刺突和尾部震荡的影响并确定P波大致到时，进而采用权重分析排除上述干扰，极大地降低了错误拾取率，增强了P波拾取的准确性。

一种基于DWT_STA/LTA的含噪信号P波初至峰度拾取方法，包括以下几个步骤：

步骤1：提取微震信号x(n)

从矿山微震信号中提取待拾取信号x(n)，n＝1,2,…,N，其中，N为所述微震信号的采样点

个数，取N＝400～600，微震信号采样频率600Hz；

步骤2：对微震信号采用离散小波变换DWT得到细节信号D_i，i＝1,2,3,4；

步骤3：依据STA/LTA算法按照以下公式确定细节信号D_i的P波大致到时k_i；

其中，STA(k,D_i)和LTA(k,D_i)分别为短、长时窗内微震信号的振幅平均值；x(n,D_i)为细节

信号D_i的第n个采样信号；k为采样的第k个点，k＝W_LTA,W_LTA+1,…,N；W_STA和W_LTA分

别为短、长时窗的长度，分别取10和60个采样点的长度；

若细节信号D_i对应的短时窗和长时窗的振幅平均值的比值λ(k,D_i)的最大值大于设定阀值a，则取首个触发点对应的采样点序号为该细节信号D_i的P波大致到时k_i，否则，记该细节信号D_i的P波大致到时k_i为0。其中，i＝1,2,3,4；

步骤4：计算细节信号的P波大致到时k_i＝0的个数L，i＝1,2,3,4；

步骤5：利用权重分析计算P波初至准确到时；

若步骤4中L>2，则认为该信号为纯噪音信号，反之则采用以下权重公式计算P波准确到时k'：

其中，1/8、3/8、3/8和1/8分别为D₄、D₃、D₂和D₁的权重系数；A_i和k_i'分别为R_i的最大振幅和局部峰度拾取到时，R_i为D_i在区间[k_i-b,k_i+b]的信号，b为确定局部峰度拾取区间的常量，取b＝10～15；

由公式计算x(n,D_i)的滑动峰度值K(k,D_i)，k为采样的第k个点，k＝k_i-b,k_i-b+1,…,k_i+b，M为滑动时窗的长度，取M＝20～30，并取K(k,D_i)最大值点作为局部峰度拾取到时k_i'；

步骤6：依据步骤5输出拾取结果，该信号为纯噪音信号或者P波初至准确到时为k'。

所述设定阈值a＝2～3。

利用Mallat算法快速实现离散小波变换。

并采用低通滤波器h和高通滤波器g获得原始信号c₀(等同于本发明中的x(n))的小波系数{d₁,d₂,…,d_J}，从而将微震信号x(n)分解得到细节信号D_i(取i＝1,2,3,4)。

式中：c_j+1,k为第k个采样点的近似信号值c_j+1，d_j+1,k为第k个采样点的细节信号值d_j+1，k为采样的第k个点，l为滤波器确定指标，j为尺度参数(j＝0,1,…,J-1)，J为最大分解层数。为简化说明，记C_j+1＝{c_j+1,k}和D_j+1＝{d_j+1,k}。

设定k_i＝0时，A_i＝0，k_i'＝0，i＝1,2,3,4，简化计算过程，提高计算速度。

有益效果

本发明提供了一种基于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波拾取的准确性。

附图说明

图1是本发明所述方法流程图；

图2是典型波形W-S/L-K法实现过程图，其中，(a)为矿山微震原始信号x(n)(n＝1,2,…,400)，虚线为人工拾取P波初至点；(b)为原始信号对应的峰度值K(k)(k＝21,22,…,400)，括号中数值分别为PAI-K法拾取P波初至点及其对应的峰度值；(c)为小波分解所得细节信号D_i(i＝1,2,3,4)；(d)为细节信号D_i(i＝1,2,3,4)的STA(k,D_i)/LTA(k,D_i)值λ(k,D_i)；(e)为细节信号局部峰态拾取k_i'(i＝1,2,3,4)；

图3是PAI-K和W-S/L-K法工程拾取效果对比图，其中，(a)为PAI-K拾取误差图；(b)为W-S-L-K拾取误差图；

图4是PAI-K和W-S/L-K法拾取典型案例图，其中，虚线、点线和实线分别对应人工、PAI-K法和W-S/L-K法拾取到时，(a)为含刺突纯噪音信号，(b)为低信噪比信号，(c)为起震明显信号，(d)为低信噪比且含刺突信号，(e)为极低信噪比、含刺突且尾部存在震荡信号，(f)为起震不明显且尾部存在震荡的低信噪比信号。

具体实施方式

下面将结合附图和实施例对本发明做进一步的说明。

如图1所示，一种基于DWT_STA/LTA的含噪信号P波初至峰度拾取方法，包括以下步骤：

步骤1：提取微震信号x(n)

从矿山微震信号中提取待拾取信号x(n)，n＝1,2,…,N，其中，N为所述微震信号的采样点个数，取N＝400～600，微震信号采样频率600Hz；

步骤2：对微震信号采用离散小波变换DWT得到细节信号D_i，i＝1,2,3,4；

利用Mallat算法快速实现离散小波变换。

并采用低通滤波器h和高通滤波器g获得原始信号c₀(等同于本发明中的x(n))的小波系数{d₁,d₂,…,d_J}，从而将微震信号x(n)分解得到细节信号D_i(取i＝1,2,3,4)。

式中：c_j+1,k为第k个采样点的近似信号值c_j+1，d_j+1,k为第k个采样点的细节信号值d_j+1，k为采样的第k个点，l为滤波器确定指标，j为尺度参数(j＝0,1,…,J-1)，J为最大分解层数。为简化说明，记C_j+1＝{c_j+1,k}和D_j+1＝{d_j+1,k}。

步骤3：依据STA/LTA算法按照以下公式确定细节信号D_i的P波大致到时k_i；

其中，STA(k,D_i)和LTA(k,D_i)分别为短、长时窗内微震信号的振幅平均值；x(n,D_i)为细节信号D_i的第n个采样信号；k为采样的第k个点，k＝W_LTA,W_LTA+1,…,N；W_STA和W_LTA分别为短、长时窗的长度，分别取10和60个采样点的长度；

若细节信号D_i对应的短时窗和长时窗的振幅平均值的比值λ(k,D_i)的最大值大于设定阀值a，则取首个触发点对应的采样点序号为该细节信号D_i的P波大致到时k_i，否则，记该细节信号D_i的P波大致到时k_i为0。其中，i＝1,2,3,4；

步骤4：计算细节信号的P波大致到时k_i＝0的个数L，i＝1,2,3,4；

步骤5：利用权重分析计算P波初至准确到时；

若步骤4中L>2，则认为该信号为纯噪音信号，反之则采用以下权重公式计算P波准确到时k'：

其中，1/8、3/8、3/8和1/8分别为D₄、D₃、D₂和D₁的权重系数；A_i和k_i'分别为R_i的最大振幅和局部峰度拾取到时，R_i为D_i在区间[k_i-b,k_i+b]的信号，b为确定局部峰度拾取区间的常量，取b＝10～15；

由公式计算x(n,D_i)的滑动峰度值K(k,D_i)，k为采样的第k个点，k＝k_i-b,k_i-b+1,…,k_i+b，M为滑动时窗的长度，取M＝20～30，并取K(k,D_i)最大值点作为局部峰度拾取到时k_i'；

为简化分析，设定k_i＝0(i＝1,2,3,4)时，A_i＝0，k_i'＝0。

步骤6：依据步骤5输出拾取结果，该信号为纯噪音信号或者P波初至准确到时为k'。

所述设定阈值a＝2～3。

实施例1：

图2为典型波形W-S/L-K法实现过程图，并将其实现过程及结果数据汇于表1。由图表知：PAI-K法拾取低信噪比信号时，P波初至附近有局部最大值，但受尾部震荡影响，全局最大值点与P波初至点相差甚大。W-S/L-K法拾取时，对信噪比低的细节信号D₄无拾取到时，记k_i'＝0、A_i＝0，进而采用式(5)计算得到P波初至。比较W-S/L-K法拾取和人工拾取，得到W-S/L-K法拾取效果较好，很好地解决了PAI-K法拾取低信噪比信号不稳定的问题。

表1 W-S/L-K法拾取过程及结果数据

注：表中k_i＝0指STA/LTA在细节信号D_i(i＝1,2,3,4)上无拾取到时，并记k_i'＝0、A_i＝0。

实施例2：

图3是PAI-K法和W-S/L-K法工程拾取效果对比图，200组微震信号从开阳磷矿用沙坝矿区IMS微震系统随机抽取得到。W-S/L-K法计算参数如下：微震信号采样频率600Hz，采样点个数N＝400，STA/LTA法阀值a＝2.5，局部峰度拾取范围常量b＝10个采样点，滑动时窗长度M＝20。

由图3知PAI-K法拾取误差为0～3、4～6、7～9及大于9个采样点的事件占总事件的比例分别为60.5％、11.5％、4％和24％；W-S/L-K法拾取误差为0～3、4～6、7～9及大于9个采样点的事件比例分别为58％、28％、11.5％和2.5％。由此可知W-S/L-K法拾取误差较大的比例很低，极大地提高了P波自动拾取的准确性。

图4是PAI-K法和W-S/L-K法拾取典型案例图，由图知PAI-K法对起震明显信号拾取效果较好(图4c)，但对纯噪音信号(图4a)、低信噪比信号(图4b、e和f)、含尖刺信号(图4d)、尾部震荡信号(图4e)拾取效果较差；而W-S/L-K法除信噪比低，且起震不明显信号(图4f)外，拾取效果均很好。由此，W-S/L-K法有效解决了PAI-K法对低信噪比、刺突、尾部震荡及纯噪音信号拾取不稳定的技术问题。

以上所述仅为本发明的实施例而已，并不用以限制本发明，凡在本发明精神和原则之内，所作任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于DWT_STA/LTA的含噪信号P波初至峰度拾取方法，其特征在于，包括以下几个步骤：

步骤1：提取微震信号x(n)

从矿山微震信号中提取待拾取信号x(n)，n＝1,2,…,N，其中，N为所述微震信号的采样点个数，取N＝400～600，微震信号采样频率600Hz；

步骤2：对微震信号采用离散小波变换DWT得到细节信号D_i，i＝1,2,3,4；

步骤3：依据STA/LTA算法按照以下公式确定细节信号D_i的P波大致到时k_i；

$\mathrm{\λ}(k,D_i)=\frac{STA(k,D_i)}{LTA(k,D_i)}=\frac{\frac{1}{W_{STA}}\underset{n=k-W_{STA}}{\overset{k}{\mathrm{\Σ}}}\left|x(n,D_i)\right|}{\frac{1}{W_{LTA}}\underset{n=k-W_{LTA}}{\overset{k}{\mathrm{\Σ}}}\left|x(n,D_i)\right|}$

其中，STA(k,D_i)和LTA(k,D_i)分别为短、长时窗内微震信号的振幅平均值；x(n,D_i)为细节信号D_i的第n个采样信号；k为采样的第k个点，k＝W_LTA,W_LTA+1,…,N；W_STA和W_LTA分别为短、长时窗的长度，分别取10和60个采样点的长度；

若细节信号D_i对应的短时窗和长时窗的振幅平均值的比值λ(k,D_i)的最大值大于设定阀值a，则取首个触发点对应的采样点序号为该细节信号D_i的P波大致到时k_i，否则，记该细节信号D_i的P波大致到时k_i为0，其中，i＝1,2,3,4；

步骤4：计算细节信号的P波大致到时k_i＝0的个数L，i＝1,2,3,4；

步骤5：利用权重分析计算P波初至准确到时；

若步骤4中L>2，则认为该信号为纯噪音信号，反之则采用以下权重公式计算P波准确到时k'：

$k^{\′}=\frac{\frac{1}{8}{A}_4\·k_4^{\′}+\frac{3}{8}{A}_3\·k_3^{\′}+\frac{3}{8}{A}_2\·k_2^{\′}+\frac{1}{8}{A}_1\·k_1^{\′}}{\frac{1}{8}{A}_4+\frac{3}{8}{A}_3+\frac{3}{8}{A}_2+\frac{1}{8}{A}_1}$

其中，1/8、3/8、3/8和1/8分别为D₄、D₃、D₂和D₁的权重系数；A_i和k_i'分别为R_i的最大振幅和局部峰度拾取到时，R_i为D_i在区间[k_i-b,k_i+b]的信号，b为确定局部峰度拾取区间的常量，取b＝10～15；

由公式计算x(n,D_i)的滑动峰度值K(k,D_i)，k为采样的第k个点，k＝k_i-b,k_i-b+1,…,k_i+b，M为滑动时窗的长度，取M＝20～30，并取K(k,D_i)最大值点作为局部峰度拾取到时k_i'；

步骤6：依据步骤5输出拾取结果，该信号为纯噪音信号或者P波初至准确到时为k'。

2.根据权利要求1所述的方法，其特征在于，所述设定阈值a＝2～3。

3.根据权利要求1或2所述的方法，其特征在于，所述步骤2中，利用Mallat算法快速实现离散小波变换。

4.根据权利要求3所述的方法，其特征在于，设定k_i＝0时，A_i＝0，k_i'＝0，i＝1,2,3,4。