CN112526602B - 一种基于长短时窗和ar模型方差激增效应的p波到时拾取方法 - Google Patents

Abstract

发明提供一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法。该方法包括提取微震信号波形数据、初步判断P波到时所处区间和拾取P波到时等步骤。本方法借助STA/LTA法初步判断P波到时所处区间以及信号与噪声水平差异导致的自回归模型方差激增效应，极大地降低了拾取错误率，提高了拾取效率，同时也增强了P波拾取的稳定性。此方法具有拾取精度高、速度快、稳定性好、适用性强等特点。

Description

一种基于长短时窗和AR模型方差激增效应的P波到时拾取 方法

技术领域

本发明涉及信号处理技术领域，特别涉及一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法。

背景技术

随着我国经济建设的快速发展，能源开采及存储工程、隧道工程、核废料深埋工程不断兴建，矿震、岩爆等地压灾害剧增，严重威胁了国民生命财产安全，因此有效预测岩爆等灾害已成为刻不容缓的研究课题。微震监测作为一种有效监测技术，已广泛应用于矿震、岩爆等灾害的预警中。微震监测主要包括波形识别、震相拾取、震源定位、震源机制分析和微震活动预测等，其中震相拾取是震源定位以及震源机制分析中最为关键的一步。因S波常叠加于P波尾部，难以识别，故震相拾取一般选用P波到时拾取，其准确性直接影响到震源定位和震源机制分析的准确性，其处理速度直接影响微震监测的效率。

P波到时拾取最初采用人工拾取，虽保证了拾取的精度，但效率太低。因此，国内外学者相继提出了一系列P波到时自动拾取方法，提高了拾取的效率，但它们都存在各自的局限：STA/LTA算法中长短时间窗口的长度、阈值的选择都直接影响了识别效果，且低信噪比或初动不明显时，识别效果较差；小波变换易陷入Fourier分析谱泄露的局限；高阶统计法对波形清晰程度要求较高；PAI-S/K方法引入微震波形峰度和偏斜度函数，计算简单、速度快，但抗噪不佳；分形维数方法和神经网络法计算耗时长。基于自回归模型(AR模型)的AR-AIC算法方法拾取精度较高，但采用双序列AR模型导致拾取效率不高，同时对于低信噪比、刺突信号拾取效果较不稳定。

因此，急需研究出一种拾取精度高、速度快、稳定性好的自动拾取方法。

发明内容

本发明的目的是提供一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法，以解决现有技术中存在的问题。

为实现本发明目的而采用的技术方案是这样的，一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法，包括以下步骤：

1)从微震信号中截取待拾取信号波形数据x(n)。其中，n＝1，2，…，N。N为微震信号的采样点个数。N＝1000～2000。

2)采用STA/LTA法定位包含P波到时的区间X′。

3)拾取X′的VIC自相关最小值点作为P波到时。

进一步，步骤2)具体包括以下步骤：

2.1)给定一个滑动的长时间窗，在长时间窗内取一个短时间窗。所述长时间窗和短时间窗的终点重合。短时窗的长度为ST。长时窗的长度为LT。用短时窗信号平均值STA和长时窗信号平均值LTA之比反映信号振幅的变化。STA/LTA算法公式为：

S(k)＝STA(k)/LTA(k) (1)

式中，i为采样时刻，X(i)为在i时刻的关于微震信号的特征函数值。k为采样的第k个点。X(i)＝|xi|。

2.2)拾取STA/LTA法首个阈值触发点k₁。STA/LTA阀值为3。

2.3)取K＝k₁+LT-1作为初步P波作为初步P波初至点。

2.4)利用初步P波初至点左右各推移250采样点，得到P波到时所处区间X′＝[X(K-250)，X(K+250)]。

进一步，步骤3)中由公式(2)计算x(n)的自相关值VIC(k)。选取使VIC最小的k作为P波到时。

式中，k＝1,2,…,n。为信号[X′(1)，X′(k)]的自回归模型利用最小二乘法估计得到的标准差估计值。a为方差激增系数。

进一步，步骤3)之后，还具有对拾取的P波到时进行评价和后续分析处理的相关步骤。

本发明还提供一种P波到时拾取设备，包括：

消息交换模块，用于接收并截取微震信号。

存储器，用于存储计算机程序。

处理器，用于执行所述存储器中存储的计算机程序，实现如权利要求1～4任一所述的P波到时拾取方法。

数据共享模块，用于对P波到时进行报送。

本发明的技术效果是毋庸置疑的：

A.解决了AR-AIC法拾取精度较高，但计算效率低下且对低信噪比、刺突信号拾取不稳定的技术问题；

B.借助STA/LTA法初步判断P波到时所处区间以及信号与噪声水平差异导致的自回归模型方差激增效应，极大地降低了拾取错误率，提高了拾取效率，同时也增强了P波拾取的稳定性；

C.拾取精度高、速度快、稳定性好且适用性强。

附图说明

图1为方法流程图；

图2为和nlnN/N随N变化趋势图；

图3为AR-VIC法P波到时拾取示意图；

图4为常规信号AR-AIC拾取图；

图5为刺突信号AR-AIC拾取图；

图6为低信噪比信号AR-AIC拾取图。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

实施例1：

参见图1，本实施例公开一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法，包括以下步骤：

1)从微震信号中截取待拾取信号波形数据x(n)。其中，n＝1，2，…，N。N为微震信号的采样点个数。N＝1000～2000。

2)采用STA/LTA(short term averaging/long term averaging，长短时窗能量比)法定位包含P波(primary wave)到时的区间X′。

2.1)给定一个滑动的长时间窗，在长时间窗内取一个短时间窗。所述长时间窗和短时间窗的终点重合。短时窗的长度为ST,通常取8-15，长时窗的长度为LT，通常取(5-10)ST。用短时窗信号平均值STA和长时窗信号平均值LTA之比反映信号振幅的变化。STA/LTA算法公式为：

S(k)＝STA(k)/LTA(k) (1)

式中，i为采样时刻，X(i)为在i时刻的关于微震信号的特征函数值，k为采样的第k个点，其中X(i)＝|xi|。

2.2)拾取STA/LTA法首个阈值触发点k₁。STA/LTA阀值为3。

2.3)取K＝k₁+LT-1作为初步P波初至点。本实施例中ST取8，LT取50。

2.4)利用初步P波初至点左右各推移250采样点，得到P波到时所处区间X′＝[X(K-250)，X(K+250)]。

3)拾取X′的VIC自相关最小值点作为P波到时。由公式(2)计算x(n)的自相关值VIC(k)。选取使VIC最小的k作为P波到时。

式中，k＝1,2,…,n。为信号[X′(1)，X′(k)]的自回归模型利用最小二乘法估计得到的标准差估计值。a为方差激增系数。

参见图2，波形数据其AR模型方差的对数lnσ随着N的增加呈现出稳定-激增-稳定的变化趋势，那么考虑方差激增效应，并引入一个减小因子nlnN/N，从而构造VIC值，使其第一阶段的稳定态变为下降状态，前期下降快而后期下降趋缓，同时保留第三阶段的稳定态，则第一阶段的VIC的极小值点有理由认为是P波到时点。为定阶AR模型利用最小二乘法估计得到的标准差估计值，a为方差激增系数，目的在于放大方差激增幅度，保证第三阶段VIC最终值始终小于第一阶段的VIC极小值。噪声与信号幅值数量级相差越小，其值越大，但不能过大，否则/>波动将变得很大。a可以考虑取8/log(b),b为信号与噪声最大幅值之比，此时/>大于20，只要噪声长度x₁不过小(x₁<100)，信号长度x₂不过长(x₂>2000)，nlnN/N肯定小于20。n为减小因子的放大系数，n取太小会造成后期减小因子减幅过小，使得整个VIC值受方差控制，方差波动较大时，一旦方差的最小值出现在VIC拐点前，此点VIC值可能最小，从而造成识别错误。因此噪声波动较大，n应取大值。N为AR模型的容量，经大量试算，n取整个序列容量时AR-VIC法识别效果最佳。

值得说明的是，AR-AIC法采用[X′(1)，X′(k)]和[X′(k+1)，X′(n)]两个序列的AR模型的AIC值之和的最小值所对应的k作为P波到时，而本发明的STA/LTA_VIC法仅需[X′(1)，X′(k)]序列便可完成P波到时拾取，故本发明方法效率更佳。

图4～6中，(a)微震信号原始波形和STA/LTA_VIC法拾取示意图，虚线为人工和STA/LTA_VIC法拾取P波到时(两者相同)；(b)整个微震信号对应的AIC值，括号中数字分别对应可能拾取点和该点AIC值；(c)整个微震信号对应的STA/LTA值，括号中数字分别对应可能拾取点和该点STA/LTA值；(d)P波到时所在区间的信号的AIC值，括号中数字分别对应可能拾取点和该点AIC值。

图4至图6表明：相比于AR-AIC法、STA/LTA法和STA/LTA-AIC，STA/LTA-VIC法拾取常规信号、刺突信号和低信噪比信号的P波到时精度更高。将其结果汇于下表1。实际微震信号拾取效果(表1)表明STA/LTA-VIC拾取效果在这几种方法中最优。表1为三种方法拾取结果。从表1可看出，STA/LTA_VIC对常规信号拾取精度高，且有效解决了低信噪比和刺突对AIC拾取结果稳定的影响。

表1

表1中，微震序号1～3分别对应常规信号、刺突信号和低信噪比信号，Manual表示手动拾取，SL表示STA/LTA法，SL-A表示STA/LTA_AIC法，SL-V表示STA/LTA_VIC法。

表2是STA/LTA法、STA/LTA_AIC法和STA/LTA_VIC法拾取误差表，图中50组微震信号为加入一定信噪比的高斯噪声的声发射信号。由表2知STA/LTA_VIC法拾取精度最高，其次是STA/LTA_AIC法，STA/LTA法最差，且在低信噪比下，STA/LTA_VIC法仍有极佳的拾取精度，此外随着信噪比的降低，三种方法拾取结果的差异逐渐增大。

表2

4)对拾取的P波到时进行评价和后续分析处理。微震时差定位法利用拾取的P波到时进行微震定位。通过识别P波到时从而提取出P波初动振幅，继而进行微震类型的判别。

实施例2：

本实施例公开一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法，包括以下步骤：

1)从微震信号中截取待拾取信号波形数据x(n)。其中，n＝1，2，…，N。N为微震信号的采样点个数。N＝1000～2000。

2)采用STA/LTA(short term averaging/long term averaging，长短时窗能量比)法定位包含P波(primary wave)到时的区间X′。

3)拾取X′的VIC自相关最小值点作为P波到时。

实施例3：

本实施例公开一种P波到时拾取设备包括：

消息交换模块，用于接收并截取微震信号。

存储器，用于存储计算机程序。

处理器，用于执行所述存储器中存储的计算机程序，实现如实施例1或2任一所述的P波到时拾取方法。

数据共享模块，用于对P波到时进行报送。

Claims (3)

1.一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法，其特征在于，包括以下步骤：

1)从微震信号中截取待拾取信号波形数据x(n)；其中，n＝1，2，…，N；N为微震信号的采样点个数；N＝1000～2000；

2)采用STA/LTA法定位包含P波到时的区间X′；步骤2)具体包括以下子步骤：

2.1)给定一个滑动的长时间窗，在长时间窗内取一个短时间窗；所述长时间窗和短时间窗的终点重合；短时窗的长度为ST；长时窗的长度为LT；用短时窗信号平均值STA和长时窗信号平均值LTA之比反映信号振幅的变化；STA/LTA算法公式为：

S(k)＝STA(k)/LTA(k) (1)

式中，i为采样时刻，X(i)为在i时刻的关于微震信号的特征函数值；k为采样的第k个点；X(i)＝|xi|；

2.2)拾取STA/LTA法首个阈值触发点k₁；STA/LTA阀值为3；

2.3)取K＝k₁+LT-1作为初步P波作为初步P波初至点；

2.4)利用初步P波初至点左右各推移250采样点，得到P波到时所处区间X'＝[X(K-250),X(K+250)]；

3)拾取X′的VIC自相关最小值点作为P波到时；其中，由公式(2)计算x(n)的自相关值VIC(k)；选取使VIC最小的k作为P波到时；

式中，k＝1,2,…,n；为信号[X′(1)，X′(k)]的自回归模型利用最小二乘法估计得到的标准差估计值；a为方差激增系数。

2.根据权利要求1所述的一种基于长短时窗和AR模型方差激增效应的P波到时拾取方法，其特征在于：步骤3)之后，还具有对拾取的P波到时进行评价和后续分析处理的相关步骤。

3.一种P波到时拾取设备，其特征在于，包括：

消息交换模块，用于接收并截取微震信号；

存储器，用于存储计算机程序；

处理器，用于执行所述存储器中存储的计算机程序，实现如权利要求1或2所述的P波到时拾取方法；

数据共享模块，用于对P波到时进行报送。