基于改进的STA/LTA方法的微地震P波自动拾取技术

引用本文 [复制中英文]

刘晓明, 赵君杰, 王运敏, 彭平安. 基于改进的STA/LTA方法的微地震P波自动拾取技术[J]. 东北大学学报:自然科学版, 2017, 38(5): 740-745.

LIU Xiao-ming, ZHAO Jun-jie, WANG Yun-min, PENG Ping-an. Automatic Picking of Microseismic Events P-wave Arrivals Based on Improved Method of STA/LTA[J]. Journal of Northeastern University Nature Science, 2017, 38(5): 740-745. DOI: 10.3969/j.issn.1005-3026.2017.05.027.

[复制英文]

基金项目

国家自然科学基金资助项目 (51604301, 51374242)；金属矿山安全与健康国家重点实验室开发基金资助项目 (2015-JSKSSYS-01)；中央高校基本科研业务费专项资金资助项目 (2015zzts073)

作者简介

刘晓明 (1982-)，男，江西宜春人，中南大学讲师，博士;
王运敏 (1956-)，男，安徽砀山人，中钢集团马鞍山矿山研究院教授级高级工程师。

文章历史

收稿日期：2015-12-23

Contents Abstract Full text Figures/Tables PDF

基于改进的STA/LTA方法的微地震P波自动拾取技术

刘晓明^1,2, 赵君杰^1,2, 王运敏³, 彭平安¹

1. 中南大学资源与安全工程学院, 湖南长沙 410083;
2. 金属矿山安全与健康国家重点实验室, 安徽马鞍山 243004;
3. 中钢集团马鞍山矿山研究院有限公司, 安徽马鞍山 243004

收稿日期：2015-12-23

基金项目：国家自然科学基金资助项目 (51604301, 51374242)；金属矿山安全与健康国家重点实验室开发基金资助项目 (2015-JSKSSYS-01)；中央高校基本科研业务费专项资金资助项目 (2015zzts073)。

作者简介：刘晓明 (1982-)，男，江西宜春人，中南大学讲师，博士;
王运敏 (1956-)，男，安徽砀山人，中钢集团马鞍山矿山研究院教授级高级工程师。

摘要：基于STA/LTA方法, 提出了一种改进的自动拾取P波到时算法.引入权重因子K, 构建了一种新的特征函数, 提出寻找STA/LTA曲线上的最大值对应时刻为P波到时的方法, 并运用Matlab对冬瓜山铜矿微震事件进行初至拾取分析.结果表明：该方法能够准确地拾取微地震P波初至, 减少了传统的STA/LTA算法需要调整触发阈值 (Thr) 的过程, 克服了固定阈值不能适应复杂多变信号到时拾取的缺点, 具有抗噪性强、稳定性强、拾取精确度高等优势.

关键词：微震监测 P波初至自动拾取 STA/LTA方法 Matlab

Automatic Picking of Microseismic Events P-wave Arrivals Based on Improved Method of STA/LTA

LIU Xiao-ming^1,2, ZHAO Jun-jie^1,2, WANG Yun-min³, PENG Ping-an¹

1. School of Resources and Safety Engineering, Central South University, Changsha 410083, China;
2. State Key Laboratory of Safety and Health for Metal Mines, Maanshan 243004, China;
3. Maanshan Institute of Mining Research, SinoSteel Group, Maanshan 243004, China

Corresponding author: ZHAO Jun-jie, E-mail:812108841@qq.com

Abstract: An improved method of automatic picking P-wave first arrival time based on the method of STA/LTA was put forward. A new kind of characteristic function was established by introducing weighting factor K, the global maximum on the curve of STA/LTA was found, then the corresponding moment is identified as the arrival of P-wave, and Matlab was used to analyze seismic events. Results show that the method can accurately pick up P-wave arrivals of seismic event, and reduce the process of adjusting the trigger threshold (Thr) to pick up the P-wave of the traditional algorithm (STA/LTA), and overcome the shortcoming that a fixed threshold can not pick the arrival time of complex signals. And the method has advantages of strong noise resistance, good stability, etc.

Key Words: microseismic monitoring automatic picking of first arrival time of P-wave method of STA/LTA Matlab

微震监测技术作为一种地压稳定性监测和管理的重要手段, 目前已在石油页岩气非常规开采、二氧化碳和放射性核废料的封存、大坝和边坡及硐室的稳定性监测、以及矿山地压安全稳定性监测等多个工程领域取得了广泛的应用和发展.随着我国浅地表矿产资源逐渐减少, 深部开采已是一种常态, 而矿床深部开采高应力诱发的岩爆灾害严重制约着矿山的安全生产, 给社会以及人们的生命财产造成了巨大的损失, 微震监测技术必将在我国矿山安全生产中发挥重要作用.

微地震数据处理中, 准确拾取P波的初至对微震事件的识别、震源定位以及震源机制分析具有非常重要的意义, 通过人工观察震动图在一定程度上能够较为准确地拾取P波初至, 但是, 对于地震台网记录的庞大微地震事件, 若采用人工手动拾取P波到时不仅耗时耗力, 效率低下而且经常由于主观因素的介入, 不同专业人员拾取到的P波初至不尽相同, 对后续的分析也将带来不同程度的影响.因此, 人们都在不断地寻求能够自动、快速、精确、稳定拾取P波初至的方法.

国内外学者对微震P波到时自动拾取做了详细深入的研究.Stevenson^[1]研究了地震自动分析理论并提出了长短时平均能量比法 (STA/LTA)；Sabbione等对单通道地震记录的地震自动监测和初至拾取进行研究和分析, 并总结了震相自动识别的现状, 提出了以特征函数作为输入对象的STA/LTA算法, 该思想广为流传, 经过不断地发展和探索, 相继提出了许多改进特征函数的STA/LTA算法^[2-4]；张唤兰等提出了AIC信息准则法, 以及随后基于自回归模型的AR-AIC算法^[5]；分形分维法和人工神经网络法^[6-7], 一定程度上实现困难且效率低；Saragiotis提出了基于高阶统计量偏斜度和峰度的PAI-S/K法^[8], 该方法是结合数学统计发展起来的应用于信号分析与处理的技术, 具有抗干扰性强的特点, 得到了广泛的应用^[9-10]；以及小波变换在信号处理分析方面发展日益成熟^[11-12].

STA/LTA方法是目前应用最广泛的自动拾取方法, 但其拾取结果的准确性和稳定性严重依赖特征函数、移动时窗长度和触发阈值的选取, 因此一定程度上无法实现对P波到时的自动拾取.针对以上问题, 本文提出一种稳定性强、运行速度快、计算准确的P波自动拾取算法, 该方法能够适应于各种复杂多变的微震信号, 具有广泛的工程领域应用价值.

1 STA/LTA拾取P波初至的算法分析

长短时平均值比法 (STA/LTA), 最早由Stevenson^[1]提出, 用来准确拾取微震事件的到达时刻, 其中STA (短时窗平均值) 表示微震事件信号振幅或能量的变化趋势, LTA (长时窗平均值) 表示背景噪声信号振幅或能量的变化趋势.

采用该种方法首先需要定义特征函数, 然后以特征函数为基础在滑动时窗内分别计算STA与LTA的值.鉴于信号本身特征, 当信号到达时, STA的值比LTA的值变化快, LTA相对比较稳定, 相应STA/LTA值在信号到达时会出现突然增大.在人为设定好的阈值 (Thr) 条件下, 当STA/LTA大于该阈值 (Thr) 时, 认为此点就是P波的初至时刻, 从而便可以实现微震事件的自动检测和到时拾取.STA/LTA算法的基本公式为

(1)

(2)

(3)

式中：i表示采样时刻；long表示长时窗长度；short表示短时窗长度；λ表示设定的触发阀值；CF (i) 表示微震信号在i时刻对应的特征函数值, 表征微震数据的振幅、能量的变化情况.

STA/LTA算法主要影响因素体现在以下3个方面.

1) 时窗长度.特征函数确定情况下, 不同的岩石破坏类型产生的微震事件, 其到时拾取结果受时窗长度影响程度不同, 如硬岩剪切破坏类型的事件受时窗长度的影响较小；而硬岩拉伸破坏的事件, 时窗长度对拾取结果的影响较为显著^[12].

2) 特征函数.特征函数反映信号振幅与能量变化特征, 特征函数的选取对准确拾取微震到时极为重要, 不同的特征函数对信号变化的灵敏度不同, 常用的特征函数有

(4)

(5)

(6)

(7)

以某原始信号1为例, 分别将式 (4)~式 (7)4种不同的特征函数代入到式 (1) 和式 (2) 中计算相应的STA, LTA及两者的比值STA/LTA, 得到STA/LTA变化曲线, 分析其对信号的响应特性如图 1所示.从图中可以看出, 图 1b~图 1e对应的4种不同特征函数对信号的响应表现各有差异, 图 1b, 图 1c, 图 1e所对应的特征函数F₁, F₂, F₄能够反映出P波到时的特征 (即STA/LTA比值突然增大), 但是同样极易受到噪声的影响, 图 1d所对应的特征函数F₃效果很不理想.

图 1 STA/LTA变化曲线 Fig.1 Corresponding curves of STA/LTA (a)—原始信号1；(b)—F₁; (c)—F₂; (d)—F₃；(e)—F₄.

3) 触发阈值.触发阈值的选取在拾取过程中也至关重要, 不同信噪比的信号所选阈值也不同.当信噪比较高时, 为了确保拾取的准确性, 阈值应设置高些, 反之如果阈值设置较低, 见图 1c, 当阈值设为2时, 任何异常的噪声扰动都可能被误认为是微震事件的初至, 就无法准确完成到时拾取；因此针对不同信号如何选取合理的触发阈值对准确拾取微震事件的初至至关重要.

2 改进的STA/LTA算法

STA/LTA到时拾取算法, 需要根据实际情况选取特征函数、设置长短时窗长度和触发阈值.然而, 触发阈值的选取具有很大波动性, 不同信号阈值的选取不尽相同, 对于高信噪比和低信噪比的信号其触发阈值选取均在一定的范围内.阈值的不确定性因素对准确自动拾取初至带来诸多误差.因此, 在STA/LTA算法基础上, 引入权重因子K, 提出一种新的特征函数F, 再通过移动时窗求出长短时窗平均值之比, 提出在STA/LTA曲线上寻找全局最大值点的方法来确定P波的初至.

2.1 特征函数的改进

特征函数应尽可能简洁可行, 微震事件到时拾取结果很大程度上取决于特征函数的选取.当P波到达时, 在整个微震信号序列上, 此刻信号主要体现为频率或者振幅, 或者两者都出现较大的变化, 因此, 所选特征函数在P波到达时必须能够尽可能理想地反映出这些变化.式 (4)~式 (7) 的4种特征函数在一定程度上能够实现到时的自动拾取, 但是对于不同岩石破坏类型、不同信噪比和能量大小的信号其稳定性较差, 造成拾取上的较大误差.特征函数F₁=|X(i)|, F₂=X(i)²一定程度上可以反映出信号在P波到达时振幅上的变化, 但是却不能体现频率上的变化, 因此提出一种新的特征函数：

1) 引进一个权重因子K.对于实数的微震信号序列X(i), 求取信号的绝对值之和, 再对信号求一阶导数, 取一阶导数的绝对值之和, 两者之比为权重因子K.

(8)

2) 构建特征函数.信号中包含的有效信息主要有能量、频率变化、振幅等, 根据信号的该特性, 考虑权重因子获得构造特征函数F：

(9)

其中:权重因子K是根据信号采样频率和台站固有的噪声性质对信号振幅和频率的权重进行分配, X(i+1)²表示信号的振幅特征; K(X(i+1)-X(i))²表示信号的频率特征, 因此该特征函数包含了详细的信号振幅和频率参数特性, 当微震事件发生时, 信号的振幅和频率两者之中的任意一个或者两者均会出现突变, F恰能反应该变化特性, 迅速做出反应.

分别采用特征函数式 (9) 与特征函数式 (4)~式 (7) 计算原始信号2的STA与LTA, 以及两者之比, 得到STA/LTA变化曲线, 原始信号1采用特征函数 (4)~(7) 所对应的STA/LTA变化曲线如图 1所示, 图 2为采用特征函数 (9) 对应的STA/LTA变化曲线.

图 2 STA/LTA变化曲线 Fig.2 Corresponding curves of STA/LTA (a)—原始信号1；(b)—式 (9).

由图 1和图 2可以看出, 根据图 1b, 图 1c, 图 1e, 所对应的特征函数 (即上文所述的特征函数 (4)~(5), (7)) 得到的STA/LTA曲线对噪声与实际微震信号P波到达时的敏感程度相差不大, 容易受到噪声的干扰, 不能清晰地显示P波到时；而根据图 2f所对应的特征函数 (即本文提出的特征函数 (9)) 得到的STA/LTA曲线能明显反应出微震信号P波到时, 即当P波到达时刻对应的STA/LTA的值相对于P波到达前所有的微震信号序列对应的STA/LTA的值要大得多, 完全不受噪声信号的影响.

2.2 P波到时拾取方法

由于井下信号 (凿岩、爆破、出矿、卡车运行等信号) 和微震信号并行存在, 设置固定触发阈值显然不能满足实际工作需要.在改进的特征函数基础上, 提出采用全局最大值拾取P波到时的方法.基本流程如图 3所示, P波到时拾取结果如图 4所示.

图 3 流程图 Fig.3 Flow diagram

图 4 P波到时拾取结果图 Fig.4 Results of P-wave arrival time

从图 4中可以看出, 改进后的算法能够快速准确地拾取到P波的初至, 省去了传统方法需要多次调整阈值的繁琐过程.

3 案例分析

冬瓜山铜矿是我国首批迈入千米深井大规模开采的有岩爆倾向的硬岩金属矿山, 也是我国实施微震监测较早的矿山之一.运用Matlab对冬瓜山铜矿记录的1 536个微震事件进行初至拾取研究分析：

1) 采用传统的STA/LTA算法, 选取F₂=X(i)²作为特征函数, 触发阈值取Thr=2, 3, 4, 4.1, 4.2, 4.3, 4.4, 分别对22个原始微震信号进行P波到时拾取, 结果见表 1.

表 1 P波到时拾取结果 Table 1 Results of P-wave arrival time

s
信号名称	触发阈值Thr							人工判别拾取到时	改进算法拾取到时	误差/%
信号名称	2	3	4	4.1	4.2	4.3	4.4	人工判别拾取到时	改进算法拾取到时	误差/%
原始信号1	0.006 5	0.016 6	0.023 6	0.023 6	0.023 6	0.023 6	阈值太高	0.023 6	0.023 6	0
原始信号2	0.005 4	0.014 5	0.014 5	0.021 4	0.021 4	0.02 1 4	阈值太高	0.021 4	0.021 4	0
原始信号3	0.002 4	0.016 5	0.016 5	0.023 5	0.023 5	阈值太高	—	0.023 5	0.023 6	0.4
原始信号4	0.004 9	0.010 7	0.021 5	0.021 5	0.021 5	0.021 5	0.021 5	0.021 5	0.021 5	0
原始信号5	0.012 2	0.016 6	0.016 6	0.016 6	0.016 6	0.021 5	阈值太高	0.021 5	0.021 6	0.4
原始信号6	0.002 5	0.016 9	0.023 8	0.023 8	0.023 8	0.023 8	0.023 8	0.023 8	0.023 9	0.4
原始信号7	0.000 4	0.015 1	0.017 2	0.017 2	0.017 2	0.017 2	0.021 8	0.021 8	0.021 9	0.4
原始信号8	0.016 9	0.016 9	0.016 9	0.016 9	0.016 9	0.023 9	0.023 9	0.023 9	0.023 8	0.4
原始信号9	0.003 5	0.003 4	0.015 9	0.015 9	0.018 9	0.022 9	0.022 9	0.022 9	0.022 9	0
原始信号10	0.005 5	0.013 0	0.016 0	0.016 0	0.016 0	0.022 9	0.022 9	0.022 9	0.023 0	0.4
原始信号11	0.009 4	0.016 5	0.020 0	0.023 4	0.023 4	0.023 4	0.023 4	0.023 4	0.023 4	0
原始信号12	0.005 3	0.005 3	0.023 3	0.023 3	0.023 3	0.023 3	阈值太高	0.023 3	0.023 5	0.8
原始信号13	0.003 5	0.015 2	0.023 1	0.023 1	0.023 1	0.023 1	0.024 4	0.024 4	0.024 3	0.1
原始信号14	0.002 3	0.021 3	0.023 0	0.023 0	0.023 0	0.024 3	0.024 3	0.024 3	0.024 3	0
原始信号15	0.017 7	0.017 6	0.024 3	0.024 3	0.024 3	0.024 3	0.024 3	0.024 3	0.024 3	0
原始信号16	0.005 5	0.006 7	0.022 7	0.022 7	0.022 7	0.022 7	0.023 9	0.023 9	0.023 9	0
原始信号17	0.017 8	0.023 5	0.023 5	0.023 5	0.024 8	0.024 8	0.024 8	0.024 8	0.024 8	0
原始信号18	0.000 7	0.022 5	0.024 4	0.024 4	0.024 4	0.024 4	0.024 4	0.024 4	0.024 4	0
原始信号19	0.005 0	0.022 6	0.022 6	0.022 6	0.022 6	0.022 6	0.022 6	0.022 6	0.022 6	0
原始信号20	0.008 7	0.023 2	0.024 4	0.024 4	0.024 4	0.024 4	0.024 4	0.024 4	0.024 4	0
原始信号21	0.006 1	0.023 1	0.023 1	0.023 1	阈值太高	—	—	0.023 1	0.023 1	0
原始信号22	0.005 1	0.023 0	0.024 2	0.024 2	0.024 2	0.024 2	0.024 2	0.024 2	0.024 2	0
注：表 1中标记的“阈值太高”表示STA/LTA曲线上没有高于该触发阈值 (Thr) 的点.

表 1 P波到时拾取结果 Table 1 Results of P-wave arrival time

2) 采用改进的STA/LTA算法, 选取特征函数式 (9) 对22个原始微震信号进行P波到时拾取, 并与传统的STA/LTA算法计算的结果进行对比分析, 结果见表 1.

图 5为部分信号传统的STA/LTA算法与文中改进的STA/LTA算法拾取过程结果图.

图 5 不同阈值条件下及本文算法P波到时拾取结果 Fig.5 Results of P-wave arrival time under different thresholds and proposed method (a)—原始信号1；(b)—Thr=3; (c)—Thr=4; (d)—Thr=4; (e)—本文算法.

从表 1和图 5中可以看出：

1) 传统的STA/LTA算法对于某一固定的触发阈值Thr不能适应所有信号, 需要多次调整触发阈值才能够准确拾取P波的初至, 文中提出的拾取STA/LTA曲线上全局最大值方法能快速拾取到P波的到时；

2) 文中所采用的改进算法其拾取到时结果与人工识别到时结果两者的误差：

|改进拾取到时-传统拾取到时|/传统拾取到时×100% 低于1%, 很大程度上保证了拾取结果的准确性.

4 结论

1) 传统的STA/LTA算法拾取结果准确性和稳定性严重依赖特征函数、移动时窗长度和触发阈值选取, 通过引入权重因子K, 构建了一种新的特征函数, 并提出了在STA/LTA曲线上寻找全局最大值点的方法来确定P波的初至.

2) 以冬瓜山铜矿实际采集的微震信号为研究对象, 利用Matlab软件对其进行到时拾取分析, 结果表明, 该方法能快速、准确、稳定地拾取微震P波到时, 拾取误差控制在1%以内, 避免了传统算法多次调节阈值的冗余步骤, 缩短了微震信号处理的时间, 极大地提高了工作效率.

3) 采用本文方法对于低信噪比或是高信噪比的信号均能取得很好的拾取结果, 因此可以应用于矿山工程领域微震监测.

参考文献

[1]	Stevenson R. Microearthquakes at Flathead Lake, Montana:a study using automatic earthquake processing[J]. Bulletin of the Seismological Society of America, 1976, 66: 61–79.
[2]	Sabbione J I, Velis D R. A robust method for microseismic event detection based on automatic phase pickers[J]. Journal of Applied Geophysics, 2013, 99: 42–50. DOI:10.1016/j.jappgeo.2013.07.011
[3]	刘晗, 张建中. 微震信号自动检测的STA/LTA算法及其改进分析[J]. 地球物理学进展, 2014, 29(4): 1708–1714. ( Liu Han, Zhang Jian-zhong. STA/LTA algorithm analysis and improvement of microseismic signal automatic detection[J]. Progress in Geophysics, 2014, 29(4): 1708–1714. DOI:10.6038/pg20140429 )
[4]	马强, 金星, 李山有, 等. 用于地震预警的P波震相到时自动拾取[J]. 地球物理学报, 2013, 56(7): 2313–2321. ( Ma Qiang, Jin Xing, Li Shan-you, et al. Automatic P-arrival detection for earthquake early warning[J]. Chinese Journal of Geophys, 2013, 56(7): 2313–2321. DOI:10.6038/cjg20130718 )
[5]	张唤兰, 朱光明, 王云宏. 基于时窗能量比和AIC的两步法微震初至自动拾取[J]. 物探与化探, 2013, 37(2): 269–273. ( Zhang Huan-lan, Zhu Guang-ming, Wang Yun-hong. Automatic microseismic event detection and picking method[J]. Geophysical & Geochemical Exploration, 2013, 37(2): 269–273. )
[6]	王馨, 王彩霞, 白超英. 地震震相初至自动检测技术综述[J]. 地球物理学进展, 2013, 28(5): 2363–2375. ( Wang Xin, Wang Cai-xia, Bai Chao-ying. Review of automatic onset time picking for seismic arrivals[J]. Progress in Geophysics, 2013, 28(5): 2363–2375. DOI:10.6038/pg20130518 )
[7]	武晔, 李信富, 李小凡. 分形理论在地震学中的应用研究[J]. 地球物理学进展, 2007, 22(2): 411–417. ( Wu Ye, Li Xin-fu, Li Xiao-fan. Application of fractal theory in seismology[J]. Progress in Geophysics, 2007, 22(2): 411–417. )
[8]	Saragiotis C D. PAI-S/K:a robust automatic seismic P phase arrival identification scheme[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(6): 1395–1404. DOI:10.1109/TGRS.2002.800438
[9]	Küperkoch L, Meier T, Lee J, et al. Automated determination of P-phase arrival times at regional and local distances using higher order statistics[J]. Geophysical Journal International, 2010, 181(2): 1159–1170.
[10]	刘劲松, 王赟, 姚振兴. 微地震信号到时自动拾取方法[J]. 地球物理学报, 2013, 56(5): 1660–1666. ( Liu Jin-song, Wang Yun, Yao Zhen-xing. On microseismic first arrival identification:a case study[J]. Chinese Journal of Geophys, 2013, 56(5): 1660–1666. DOI:10.6038/cjg20130523 )
[11]	李仕雄, 吴治涛, 骆循. 联合小波变换与偏振分析自动拾取微地震P波到时[J]. 地球物理学进展, 2012, 27(1): 131–136. ( Li Shi-xiong, Wu Zhi-tao, Luo Xun. United wavelet transform and polarization analysis automatically identify microseismic P-arrival[J]. Progress in Geophysics, 2012, 27(1): 131–136. DOI:10.6038/j.issn.1004-2903.2012.01.015 )
[12]	宋维琪, 吕世超. 基于小波分解与Akaike信息准则的微地震初至拾取方法[J]. 石油物探, 2011, 50(1): 14–21. ( Song Wei-qi, Lyu Shi-chao. Microseismic event picking-up method based on wavelet transform and Akaike criterion[J]. Geophysical Prospecting for Petroleum, 2011, 50(1): 14–21. )