CN105487114B

CN105487114B - 一种微震信号p波初至点综合拾取方法

Info

Publication number: CN105487114B
Application number: CN201510899985.0A
Authority: CN
Inventors: 尚雪义; 李夕兵; 董陇军; 王泽伟; 张楚旋; 刘栋; 周勇勇
Original assignee: Central South University
Current assignee: Central South University
Priority date: 2015-12-08
Filing date: 2015-12-08
Publication date: 2018-02-02
Anticipated expiration: 2035-12-08
Also published as: CN105487114A

Abstract

本发明公开了一种微震信号P波初至点综合拾取方法，包括如下步骤：导入微震时间序列数据；拾取STA/LTA法首个阀值触发点k₁；拾取PAI‑K法峰度最大值点k₂；拾取AIC法自相关最小值点k₃；综合法确定P波初至点k：(1)k₁无拾取时，取k＝k₃；(2)k₁有拾取时：①|k₁‑k₃|≤a，取k＝k₃；②|k₁‑k₂|≤a，且|k₁‑k₃|>a，取k＝k₂；③|k₁‑k₂|>a，且|k₁‑k₃|>a，取k＝k₁。通过综合拾取法有效解决了STA/LTA法无拾取事件较多、拾取精度较低，PAI‑K法和AIC法拾取稳定性较差的技术问题。此方法具有处理简便、适用性强、准确性高等特点。

Description

一种微震信号P波初至点综合拾取方法

技术领域

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

背景技术

微震监测技术作为一种先进且行之有效的区域灾害监测和预测手段，在国内外矿山工程、油气开采、边坡稳定和隧道工程等领域得到了广泛应用。微震监测技术主要包括：监测台网布置、数据处理和微震事件定位(Ge，2005)，而微震信号P波初至点拾取是数据处理的核心工作之一。目前，主要采用人工拾取P波初至点，但人工拾取效率低，且易受个人因素影响。为此，研究一种精度高、稳定性好的P波初至点拾取方法非常重要。

Stevenson(1976)提出了长短时窗均值比(STA/LTA)法，Allen(1978，1982)、Baer&Kradolfer(1987)、Earle&Shearer(1994)和刘晗&张建中(2014)等提出了不同的STA/LTA特征函数，增强了STA/LTA值对微震信号幅值和频率变化的响应；Saragiotis等(2002,2004)提出了基于高阶统计量-偏度和峰度的PAI-S/K法，Galiana-Merino等(2008)、等(2010)和Baillard等(2014)结合小波分解或改进拾取的功能函数，实现PAI-S/K法的提升拾取；Maeda(1985)提出了无需使用AR系数的VAR-AIC法，Sleeman&van Eck(1999)和Leonard(1999,2000)提出了自回归池赤准则(AR-AIC)法。STA/LTA法利用了P波初至点微震信号幅值和频率的变化，该算法简单易懂、计算速度快，但其采用阀值触发点作为P波初至点：阀值较小时，微震信号容易被噪音提前触发；阀值较大时，微震信号可能不被触发或触发较晚，致使误差较大，由此很难找到一个适合不同波形的阀值。PAI-K法和AIC法充分利用了P波初至前后微震信号的奇异性，拾取精度较高(图3和表2)，但易受低信噪比、尖刺和尾部振荡信号影响(图2)，拾取稳定性较差。

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

发明内容

本发明所要解决的技术问题是提供一种微震信号P波初至点综合拾取方法，该微震信号P波初至点综合拾取方法计算简便、适用性强、准确性高。

发明的技术解决方案如下：

一种微震信号P波初至点综合拾取方法，包括以下步骤：

步骤1：导入微震时间序列x(n)

导入微震时间序列x(n)，n＝1,2,…,N，其中N为微震信号的总采样点数，取N＝4000～7000，微震信号采样频率f＝4000～7000Hz；

步骤2：拾取STA/LTA法首个阀值触发点k₁

由公式计算x(n)的长短时窗均值比值STA(k')/LTA(k')，k'为第k'个采样点，k'＝W_LTA,W_LTA+1,…,N，取STA/LTA首个阀值触发点作为第一P波初至点k₁；

其中，STA/LTA阀值λ取2～3；

其中，W_STA、W_LTA分别为短、长时窗采样点个数，W_STA＝50～70，W_LTA＝500～700；

步骤3：拾取PAI-K法峰度最大值点k₂

由公式计算x(n)的峰度值K(k″)，k″＝M,M+1,…,N，取峰度最大值点作为第二P波初至点k₂。

其中，M为滑动时窗的长度，M＝200～300；

步骤4：拾取AIC法自相关最小值点k₃

由公式AIC(k″′)＝k″′·log{var(x[1,k″′])}+(N-k″′-1)·log{var(x[k″′+1,N])}计算x(n)的自相关值AIC(k″′)，k″′＝1,2,…,N，取自相关最小值点作为第三P波初至点k₃。

其中，var(x[1,k″′])为第1点至第k″′点的微震时间序列方差，var(x[k″′+1,N])为第k″′+1点至第N点的微震时间序列方差；

步骤5：综合法确定P波初至点k

综合法确定P波初至点k：(1)k₁无拾取时，取k＝k₃；(2)k₁有拾取时：①|k₁-k₃|≤a，取k＝k₃；②|k₁-k₂|≤a，且|k₁-k₃|>a，取k＝k₂；③|k₁-k₂|>a，且|k₁-k₃|>a，取k＝k₁。

其中，a＝180～300。

误差范围阈值a由STA/LTA法与PAI-K法、AIC法拾取结果合理差值确定：误差范围阈值过小，则会减少PAI-K法和AIC法拾取结果的比例，最终拾取结果误差较大；误差范围阈值过大，则可能取PAI-K法和AIC法得到的不稳定结果。经过试验得到a＝180～300为合理值。

有益效果：

本发明提供的一种微震信号P波初至点综合拾取方法，主要解决STA/LTA法无拾取事件较多、拾取精度较低，PAI-K法和AIC法拾取稳定性较差的问题。本方法包括如下步骤：导入微震时间序列数据；拾取STA/LTA法首个阀值触发点k₁；拾取PAI-K法峰度最大值点k₂；拾取AIC法自相关最小值点k₃；综合法确定P波初至点k：(1)k₁无拾取时，k＝k₃；(2)k₁有拾取时：①|k₁-k₃|≤a，取k＝k₃；②|k₁-k₂|≤a，且|k₁-k₃|>a，取k＝k₂；③|k₁-k₂|>a，且|k₁-k₃|>a，取k＝k₁。本发明借助于STA/LTA法有拾取时稳定性较好，PAI-K法和AIC法拾取精度较高，且AIC法的稳定性和精度均优于PAI-K法的特性(图3和表2)：STA/LTA法的k₁无拾取时，取拾取结果较优的k₃作为P波初至点；STA/LTA法有拾取时，k₁即为P波大致初至点，再设定k₁与k₂、k₃的允许误差范围a，排除PAI-K法和AIC法稳定性较差的问题，从而提高P波初至点拾取的稳定性和精度。此方法具有处理简便、适用性强、准确性高等特点。

附图说明

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

图2是STA/LTA、PAI-K和AIC法拾取典型案例图。其中，(a)、(b)为起振明显信号，(c)为信噪比低，起振不明显信号，(d)为尾部振荡严重信号，(e)为含尖刺信号。

图3是STA/LTA、PAI-K和AIC法工程拾取误差图。

图4是综合法工程拾取误差图。

具体实施方式

下面将结合附图1～4，对本发明提出的一种综合拾取方法作进一步说明。

本发明算法思想的描述如下：针对微震信号P波初至点STA/LTA法无拾取事件较多、拾取精度较低，PAI-K法和AIC法拾取稳定性较差的问题，提出一种微震信号P波初至点综合拾取法，该方法借助于STA/LTA法有拾取时稳定性较好，PAI-K法和AIC法拾取精度较高，且AIC法的稳定性和精度均优于PAI-K法的特性(图3)：STA/LTA法无拾取时，取拾取结果较优的k₃作为P波初至点；STA/LTA法有拾取时，k₁即为P波大致初至，再设定k₁与k₂、k₃的允许误差范围a，排除PAI-K法和AIC法稳定性较差的问题，从而提高P波初至点拾取的稳定性和精度。

如图1所示，一种微震信号P波初至点综合拾取法，包括以下步骤：

步骤1：导入微震时间序列x(n)

导入微震时间序列x(n)(n＝1,2,…,N)，其中N为微震信号的总采样点数，取N＝4000～7000，微震信号采样频率f＝4000～7000Hz；

步骤2：拾取STA/LTA法首个阀值触发点k₁

由公式计算x(n)的长短时窗均值比值STA(k')/LTA(k')，k'为第k'个采样点，k'＝W_LTA,W_LTA+1,…,N，W_STA、W_LTA分别为短、长时窗采样点个数，W_STA＝50～70，W_LTA＝500～700，取STA/LTA首个阀值触发点作为第一P波初至点k₁；其中，STA/LTA阀值λ取2.5；

步骤3：拾取PAI-K法峰度最大值点k₂

由公式计算x(n)的峰度值K(k″)，k＝M,M+1,…,N，M为滑动时窗的长度，M＝200～300，取峰度最大值点作为第二P波初至点k₂；

步骤4：拾取AIC法自相关最小值点k₃

由公式AIC(k″′)＝k″′·log{var(x[1,k″′])}+(N-k″′-1)·log{var(x[k″′+1,N])}计算x(n)的自相关值AIC(k″′)，k″′＝1,2,…,N，var(x[1,k″′])为第1点至第k″′点的微震时间序列方差，var(x[k″′+1,N])为第k″′+1点至第N点的微震时间序列方差，取自相关最小值点作为第三P波初至点k₃；

步骤5：综合法确定P波初至点k

综合法确定P波初至点k：(1)k₁无拾取时，取k＝k₃；(2)k₁有拾取时：①|k₁-k₃|≤a，取k＝k₃；②|k₁-k₂|≤a，且|k₁-k₃|>a，取k＝k₂；③|k₁-k₂|>a，且|k₁-k₃|>a，取k＝k₁。其中，a＝180～300。

实施例1：

如图2(a)～(e)所示，分别为起振明显信号(图2a、b)、低信噪比，起振不明显信号(图2c)、尾部振荡严重信号(图2d)和含尖刺信号(图2e)的STA/LTA、PAI-K和AIC法拾取图，现将STA/LTA、PAI-K、AIC和综合法拾取结果汇于表1。其中，参数设定如下：传感器采样频率6000Hz，微震时间序列长度为4000个采样点，STA/LTA法的阀值为2.5，长、短时窗长度分别为600和50个采样点，PAI-K法的滑动时窗长度为200个采样点，综合法允许误差范围阈值a为240个采样点。由表知：综合法拾取的稳定性和精度较高，有效解决STA/LTA法拾取精度较低，PAI-K法和AIC法拾取稳定性较差的问题。

表1 STA/LTA、PAI-K、AIC和综合法典型波形拾取结果

实施例2：

图3是STA/LTA、PAI-K和AIC法工程拾取误差图，并将STA/LTA、PAI-K和AIC法拾取误差汇于表2，图中580组微震信号为开阳磷矿用沙坝矿区IMS微震系统随机抽取得到。参数设定如下：传感器采样频率6000Hz，微震时间序列长度为4000个采样点，STA/LTA法的阀值为2.5，长、短时窗长度分别为600和50个采样点，PAI-K法的滑动时窗长度为200个采样点。由图表知：STA/LTA法有拾取时稳定性较好，但其无拾取事件较多，且拾取精度较低；PAI-K法和AIC法拾取精度较高，但稳定性较STA/LTA法有拾取时差，且AIC法的稳定性和精度均优于PAI-K法。

图4是综合法工程拾取误差图，并将综合法拾取误差汇于表2，其中综合法允许误差范围阈值a为240个采样点。由图表知：综合法的拾取精度和稳定性均很好，误差在10ms内事件比例为92.76％，误差在30ms内事件达到了97.41％。综合法有效解决了STA/LTA法无拾取事件较多、拾取精度较低，PAI-K法和AIC法拾取稳定性较差的技术问题。

表2 STA/LTA、PAI-K、AIC和综合法工程拾取误差统计

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

Claims

1.一种微震信号P波初至点综合拾取方法，其特征在于，包括以下步骤：

步骤1：导入微震时间序列x(n)

步骤2：拾取STA/LTA法首个阀值触发点k₁

由公式计算x(n)的长短时窗均值比值STA(k')/LTA(k')，k'为第k'个采样点，k'＝W_LTA,W_LTA+1,…,N，W_STA、W_LTA分别为短、长时窗采样点个数，W_STA＝50～70，W_LTA＝500～700，取STA/LTA首个大于阀值的触发点作为第一P波初至点k₁；

其中，STA/LTA阀值为λ取2～3；

步骤3：拾取PAI-K法峰度最大值点k₂

由公式计算x(n)的峰度值K(k”)，M为滑动时窗的长度，M＝200～300，k”＝M,M+1,…,N，取峰度最大值点作为第二P波初至点k₂；

步骤4：拾取AIC法自相关最小值点k₃

由公式AIC(k″′)＝k″′·log{var(x[1,k″′])}+(N-k″′-1)·log{var(x[k″′+1,N])}计算x(n)的自相关值AIC(k”')，k”'＝1,2,…,N，var(x[1,k”'])指第1点至第k”'点的微震时间序列方差，var(x[k”'+1,N])指第k”'+1点至第N点的微震时间序列方差，取自相关最小值点作为第三P波初至点k₃；

步骤5：综合法确定P波初至点k

按照以下准则确定P波初至点k：

(1)k₁无拾取时，取k＝k₃；

(2)k₁有拾取时：①|k₁-k₃|≤a，取k＝k₃；②|k₁-k₂|≤a，且|k₁-k₃|>a，取k＝k₂；③|k₁-k₂|>a，且|k₁-k₃|>a，取k＝k₁，a表示误差范围阈值。

2.根据权利要求1所述的微震信号P波初至点综合拾取方法，其特征在于，取N＝4000～7000，f＝4000～7000Hz，W_STA＝50～70，W_LTA＝500～700，λ＝2～3，M＝200～300，a＝180～300。