, 王运敏3, 彭平安11. 中南大学 资源与安全工程学院, 湖南 长沙 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-ming1,2, ZHAO Jun-jie1,2
, WANG Yun-min3, PENG Ping-an11. 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 monitoringautomatic picking of first arrival time of P-wavemethod of STA/LTAMatlab
微震监测技术作为一种地压稳定性监测和管理的重要手段, 目前已在石油页岩气非常规开采、二氧化碳和放射性核废料的封存、大坝和边坡及硐室的稳定性监测、以及矿山地压安全稳定性监测等多个工程领域取得了广泛的应用和发展.随着我国浅地表矿产资源逐渐减少, 深部开采已是一种常态, 而矿床深部开采高应力诱发的岩爆灾害严重制约着矿山的安全生产, 给社会以及人们的生命财产造成了巨大的损失, 微震监测技术必将在我国矿山安全生产中发挥重要作用.
微地震数据处理中, 准确拾取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) |
STA/LTA算法主要影响因素体现在以下3个方面.
1) 时窗长度.特征函数确定情况下, 不同的岩石破坏类型产生的微震事件, 其到时拾取结果受时窗长度影响程度不同, 如硬岩剪切破坏类型的事件受时窗长度的影响较小；而硬岩拉伸破坏的事件, 时窗长度对拾取结果的影响较为显著[12].
2) 特征函数.特征函数反映信号振幅与能量变化特征, 特征函数的选取对准确拾取微震到时极为重要, 不同的特征函数对信号变化的灵敏度不同, 常用的特征函数有
 | (4) |
 | (5) |
 | (6) |
 | (7) |
图 1(Fig. 1)
 | 图 1 STA/LTA变化曲线Fig.1 Corresponding curves of STA/LTA (a)—原始信号1；(b)—F1; (c)—F2; (d)—F3；(e)—F4. |
3) 触发阈值.触发阈值的选取在拾取过程中也至关重要, 不同信噪比的信号所选阈值也不同.当信噪比较高时, 为了确保拾取的准确性, 阈值应设置高些, 反之如果阈值设置较低, 见图 1c, 当阈值设为2时, 任何异常的噪声扰动都可能被误认为是微震事件的初至, 就无法准确完成到时拾取；因此针对不同信号如何选取合理的触发阈值对准确拾取微震事件的初至至关重要.
2 改进的STA/LTA算法STA/LTA到时拾取算法, 需要根据实际情况选取特征函数、设置长短时窗长度和触发阈值.然而, 触发阈值的选取具有很大波动性, 不同信号阈值的选取不尽相同, 对于高信噪比和低信噪比的信号其触发阈值选取均在一定的范围内.阈值的不确定性因素对准确自动拾取初至带来诸多误差.因此, 在STA/LTA算法基础上, 引入权重因子K, 提出一种新的特征函数F, 再通过移动时窗求出长短时窗平均值之比, 提出在STA/LTA曲线上寻找全局最大值点的方法来确定P波的初至.
2.1 特征函数的改进特征函数应尽可能简洁可行, 微震事件到时拾取结果很大程度上取决于特征函数的选取.当P波到达时, 在整个微震信号序列上, 此刻信号主要体现为频率或者振幅, 或者两者都出现较大的变化, 因此, 所选特征函数在P波到达时必须能够尽可能理想地反映出这些变化.式 (4)~式 (7) 的4种特征函数在一定程度上能够实现到时的自动拾取, 但是对于不同岩石破坏类型、不同信噪比和能量大小的信号其稳定性较差, 造成拾取上的较大误差.特征函数F1=|X(i)|, F2=X(i)2一定程度上可以反映出信号在P波到达时振幅上的变化, 但是却不能体现频率上的变化, 因此提出一种新的特征函数：
1) 引进一个权重因子K.对于实数的微震信号序列X(i), 求取信号的绝对值之和, 再对信号求一阶导数, 取一阶导数的绝对值之和, 两者之比为权重因子K.
 | (8) |
 | (9) |
分别采用特征函数式 (9) 与特征函数式 (4)~式 (7) 计算原始信号2的STA与LTA, 以及两者之比, 得到STA/LTA变化曲线, 原始信号1采用特征函数 (4)~(7) 所对应的STA/LTA变化曲线如图 1所示, 图 2为采用特征函数 (9) 对应的STA/LTA变化曲线.
图 2(Fig. 2)
 | 图 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)
 | 图 3 流程图Fig.3 Flow diagram |
图 4(Fig. 4)
 | 图 4 P波到时拾取结果图Fig.4 Results of P-wave arrival time |
从图 4中可以看出, 改进后的算法能够快速准确地拾取到P波的初至, 省去了传统方法需要多次调整阈值的繁琐过程.
3 案例分析冬瓜山铜矿是我国首批迈入千米深井大规模开采的有岩爆倾向的硬岩金属矿山, 也是我国实施微震监测较早的矿山之一.运用Matlab对冬瓜山铜矿记录的1 536个微震事件进行初至拾取研究分析：
1) 采用传统的STA/LTA算法, 选取F2=X(i)2作为特征函数, 触发阈值取Thr=2, 3, 4, 4.1, 4.2, 4.3, 4.4, 分别对22个原始微震信号进行P波到时拾取, 结果见表 1.
表 1(Table 1)
 表 1 P波到时拾取结果Table 1 Results of P-wave arrival time
| 表 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(Fig. 5)
 | 图 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.) |
