地球物理学报  2019, Vol. 62 Issue (4): 1405-1412   PDF    
引用本文
付继华, 王旭, 李智涛, 等. 2019. 强噪声环境下基于信噪比的地震P波到时自动提取方法. 地球物理学报, 62(4): 1405-1412, doi: 10.6038/cjg2019L0745.  
Fu J H, Wang X, Li Z T, et al. 2019. Automatic picking up earthquake's P waves using signal-to-noise ratio under a strong noise environment. Chinese J. Geophys. (in Chinese), 62(4): 1405-1412, doi: 10.6038/cjg2019L0745.  
强噪声环境下基于信噪比的地震P波到时自动提取方法
付继华, 王旭, 李智涛, 谭巧, 王建军     
中国地震局地壳应力研究所, 北京 100085
收稿日期 2017-12-06, 2019-01-22 收修定稿
基金项目: 国家自然科学基金（41631073，41874019）和中央级公益性科研院所基本科研业务专项资助（JDZ2017-19）共同资助
第一作者简介: 付继华, 男, 博士, 副研究员, 主要从事精密测量与智慧系统研究.E-mail:fujh@email.eq-icd.cn
摘要：大数据量、强噪声环境给地震P波到时的自动提取带来很大挑战.针对此问题，本文通过构建特殊的特征函数，建立SNR与STA/LTA的内在联系，提出两种基于SNR的地震P波到时自动提取方法，即基于SNR的STA/LTA方法与基于SNR的综合方法.这两种方法分别是运用SNR概念对传统STA/LTA方法和STA/LTA与AIC综合方法的改进.仿真分析结果表明：对于弱噪声环境（10 dB）和一般噪声环境（6 dB），本文方法较传统STA/LTA方法对地震P波到时提取的准确度更高；而对于强噪声环境（3 dB），本文方法仍能准确提取地震P波到时，而传统STA/LTA方法则出现了较大的误判率（10%）与漏判率（65%）.本文方法为STA/LTA赋予了明确的物理意义，使其阈值的选取建立在严密的数学推导之上.另外，本文方法在进行地震P波到时自动提取的同时，兼具数据预处理功能，无需额外的基线校正或高通滤波，因而具有较好的实时性.
关键词: P波到时      地震预警      信噪比      STA/LTA      AIC     
Automatic picking up earthquake's P waves using signal-to-noise ratio under a strong noise environment
FU JiHua, WANG Xu, LI ZhiTao, TAN Qiao, WANG JianJun     
Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China
Abstract: Big data and strong noise environments bring great challenges to pick up P waves automatically. In order to solve this problem, a special characteristic function is constructed with the conception of Signal-to-Noise Ratios (SNR). By using this function, an internal relationship between the SNR and the Short-Term Average and Long-Term Average ratio (STA/LTA) is built. And two novel SNR-based P waves' picking methods are proposed, namely the SNR-based STA/LTA method and the SNR-based comprehensive method which are respectively the improvements of the traditional STA/LTA method and the comprehensive method of STA/LTA and Akaike Information Criteria (AIC) by using the SNR conception. The simulation analysis shows that under a weak noise circumstance (10 dB) and normal noise circumstance (6 dB) the two proposed methods have higher accuracy than the traditional STA/LTA method. And under a strong noise circumstance (3 dB), both the methods can accurately pick up the seismic P waves without any regulations, whereas the traditional STA/LTA method has a big error ratio (10%) and a large missing ratio (65%). The proposed methods give STA/LTA a clear physical meaning, and their thresholds can be obtained based on rigorous mathematical derivation. In addition, the proposed methods have favorable real-time performances because they can process data without requirement of additional baseline correction or high-pass filtering.
Keywords: P wave's arrival    Earthquake early warning    Signal-to-noise ratios    Short-Term Average and Long-Term Average ratio (STA/LTA)    Akaike Information Criteria (AIC)    
0 引言

地震预警利用地震实时监测台网，在地震发生后，破坏性的S波和面波到达前，提前几秒至几十秒发出警报，以便采取必要的紧急措施，以达到减轻地震灾害，防止次生灾害的目的(Kanamori et al., 1997；Allen and Kanamori, 2004；张红才等, 2013).日本、墨西哥、美国先后建立了各自的地震预警系统，并取得了一定的减灾效益；国内地震预警领域的相关研究与应用也在逐步推进(陈会忠等, 2011; 殷海涛等, 2012).地震预警系统性能很大程度上受限于地震监测网络的密度.为了降低地震预警系统的建设成本，增加地震监测网络的密度，一些相对廉价的地震动传感器也被应用到地震动监测中，如微机电系统(Micro-Electromechanical System, MEMS)加速度传感器，即美国国家地震监测台网系统(Advanced National Seismic System, ANSS)定义的C类传感器(Evans et al., 2014).C类传感器在0~50 Hz的频带范围内其动态范围内不低于60 dB.由于C类传感器的引入，地震预警系统中地震动信号的噪声水平被显著地增强了.然而强噪声环境给地震P波到时的自动提取，乃至地震预警三要素(时间、地点与震级)的确定带来了很大的困难.因此在大数据量和强噪声环境下，如何快速、准确地实现地震P波的自动提取是地震预警系统亟待解决的问题.

Allen提出的通过构建特征函数来计算长短时比值(Short-Term Average and Long-Term Average ratio, STA/LTA)的方法是目前地震动P波自动提取应用最广泛的一种算法(Allen, 1978; Munro, 2004).其优点是计算简便，实时性强，能够适合大数据量的应用需求.但是在强噪声环境下，如信噪比(Signal-to-Noise Ratios, SNR)小于3.5 dB，STA/LTA的误判率和漏判率则显著增加(Han et al., 2009).另外，STA/LTA方法的阈值、长/短时窗等关键参数的选择没有理论公式可以直接遵循，需要反复通过识别的误报率与漏报率进行多目标优化.Akaike提出的一种统计学的方法(Akaike Information Criteria, AIC)被应用到地震动P波中(Akaike, 1974).AIC方法假设背景噪声与地震动信号分别是两个不同的平稳随机过程，通过平稳随机过程的状态变化判定地震事件，自动提取.其优点是地震动P波的提取精度较高，能够抑制强背景噪声的干扰.但是AIC是一种基于极值计算的地震动P波方法，不适合连续、大数据量的应用需求.马强等提出了一套综合运用STA/LTA与AIC的多步骤地震动P波自动提取技术(马强等, 2013).该方法通过STA/LTA初步获得地震动P波的位置，然后再使用AIC方法提高提取的精度.其初衷是充分利用了STA/LTA与AIC各自的优点，但是AIC方法只能提高提取的准确度，而不能降低STA/LTA方法的漏报率.在大数据量、强噪声干扰条件下，其他的方法或技术手段，如分形法、偏振分析法等，也存在一定的局限性(Bataille and Chiu, 1991; Tosi et al., 1999；刘建华等, 2006; Wang et al., 2018).

本文从SNR的概念入手，通过构建特殊的特征函数，从而建立了SNR与STA/LTA的内在联系，提出了两种基于SNR的地震P波到时自动提取方法.这两种方法，为STA/LTA赋予了明确的物理意义，使其阈值的选取建立在严密的数学推导之上.另外，在进行地震P波到时自动提取的同时，兼具数据预处理功能，无需额外的基线校正或高通滤波，因而具有较好的实时性.

1 基于SNR的地震P波自动提取原理 1.1 SNR

在噪声环境下，SNR经常被用来描述有效信号与背景噪声的相对强度.SNR通常可以按照公式(1)或(2)获得，即：

(1)

(2)

式中，SNR为信噪比，Pe为有效信号的功率，Pn为背景噪声的功率，而dSNR为以线性表示的信噪比.

Pe与Pn通常采用相对功率加以描述，相对功率定义为

(3)

式中，P为信号的相对功率，xi(i=1, 2, …, N)是信号值.

1.2 STA/LTA方法

图 1所示为STA/LTA方法的示意图.左侧矩形为计算LTA的长窗，而右侧矩形为计算STA的短窗，两个窗口的右端重合，并且长窗包含短窗.

图 1 STA/LTA方法示意图 Figure 1 Schematic diagram of the STA/LTA method

STA/LTA方法分别在长、短窗内计算地震动信号特征函数的均值，并求得二者之比，作为地震P波是否到达的判据.长窗均值LTA表征了背景噪声的能量，而短窗均值STA代表了地震动信号P波的能量.当STA与LTA之比超过一定阈值时，说明P波的能量显著超出了背景噪声的能量，则可以认为地震P波到达了.

STA/LTA的计算公式为

(4)

式中，STAi和LTAi分别是短窗与长窗的特征函数均值，CF为特征函数，而S和L分别是短窗与长窗的宽度.

1.3 基于SNR的STA/LTA方法

通过比较公式(2)、公式(3)和公式(4)，可以发现SNR与STA/LTA类似，也是一种均值比，即功率比.而SNR与STA/LTA的不同之处在于，SNR是窗长度相等的功率比，而STA/LTA是窗长度不同的特征函数均值比.因此可以考虑选择适当的特征函数CF建立SNR与STA/LTA的内在联系.文中构建特征函数为

(5)

将公式(5)代入公式(4)可得新的STA/LTA公式为

(6)

根据公式(3)，公式(6)可以改为

(7)

式中，PS和PL分别是短窗与长窗内地震动信号的相对功率.

PS和PL跟Pe与Pn的关系如图 2所示.并且可由公式(8)与公式(9)给出，即：

图 2 信噪比与长短窗相对功率比的关系框图 Figure 2 Relation block diagram of the SNR and the relative power ratio of long and short sliding windows

(8)

(9)

将公式(8)、(9)代入公式(7)，并经过推导可得STA/LTA与SNR的关系式为

(10)

因此，当给定一个地震P波可能达到的SNR值时，则可以通过公式(10)求得地震P波到时的判定阈值.至此基于SNR的STA/LTA阈值、长/短时窗等关键参数便建立在严密的数学推导之上了，不再需要反复通过识别的误报率与漏报率进行多目标优化.

另外当某一地震P波到达时，设其STA/LTA值为R，那么此地震P波的SNR值的计算公式为

(11)

通过式(11)便可以研究地震P波到时的瞬时SNR特征.

1.4 基于SNR的综合方法

为了进一步增强上述STA/LTA方法的提取精度，可以将STA/LTA方法与AIC方法加以综合，多步骤地实现地震P波的自动提取.

首先对地震动信号实施基于SNR的STA/LTA方法.当该STA/LTA方法提取到地震P波时，算法自提取结果(C点)向前回溯一个长窗宽度，并向后延时一个短窗宽度，从而构建如图 3所示的AIC计算窗口.根据公式(12)计算该计算窗口内的AIC值，计算窗口内AIC的极小值代表着地震P波到达的精确时刻，公式(12)为

图 3 AIC计算窗口示意图 Figure 3 AIC′s calculation window

(12)

式中，xk(k=1, 2, …, N)是计算窗口中的地震动信号，N是AIC窗口的长度，var(x[1, k])和var(x[k+1, N])分别是数据集合x[1, k]与x[k, N]的方差.

2 算例分析 2.1 噪声分析

为了研究各种算法在强噪声环境下的实际性能，设计了一个如公式(13)所示的分段函数，来代表一个理想的地震动信号，即：

(13)

式中，当t < T0时，地震动信号f(t)为理想的零值；当t≥T0时，地震动信号f(t)为特定频率fc的正弦信号；T0则代表地震P波到时.

然后根据公式(14)，将理想的地震动信号f(t)与白噪声信号nw(t)叠加，构成叠加信号s(t)以模拟含有噪声的实际的地震动信号.白噪声信号nw(t)的强度可由其均方根值ARSM描述.这样则可以通过改变白噪声信号nw(t)的强度，来研究不同背景噪声下各种地震P波自动提取方法的实际性能.公式(15)所示为叠加信号的SNR.公式(14)和公式(15)为

(14)

(15)

之后假设地震P波到时T0为10 s，正弦信号幅值A为5 mg，而正弦信号特征频率fc为5 Hz.并且分别选择10 dB、6 dB和3 dB三种不同的SNR来代表弱噪声环境、一般噪声环境和强噪声环境.在这三种不同的背景噪声下，分别采用传统的STA/LTA方法、基于SNR的STA/LTA方法，以及基于SNR的综合方法对叠加信号s(t)进行地震P波自动提取，并根据提取结果分析各种方法在不同背景噪声下的实际性能.

传统STA/LTA方法的特征函数为

(16)

另外传统STA/LTA的短窗宽度S设为0.2 s，长窗宽度L设为3 s，P波到时自动提取的阈值设为2.4.为了更好的与传统STA/LTA方法比较，基于SNR的STA/LTA方法与综合方法也选择同样的长、短窗宽度.基于SNR的STA/LTA方法中给定的SNR设为3 dB，则dSNR等于1.9953.由公式(10)可求得地震P波提取的阈值为2.644.

在上述给定的条件下，每一个不同的SNR(10 dB、6 dB和3 dB)，分别随机产生20组不同的白噪声信号nw(t)，而地震P波到时的自动提取结果如表 1所示.表中误判率是指地震P波自动提取误差大于1 s的事件所占的比例；而漏判率是指未提取到地震P波的事件所占的比例.

表 1 不同背景噪声下地震P波自动提取统计表 Table 1 Automatic picking results of seismic P waves under different background noise

如表 1所示，在10 dB的弱噪声环境下，传统STA/LTA方法、基于SNR的STA/LTA方法与基于SNR的综合方法均能够无误判、无漏判的自动提取地震P波.基于SNR的综合方法其提取结果平均偏差和标准偏差分别为0.02 s和0.007 s，优于传统STA/LTA方法与基于SNR的STA/LTA方法.在6 dB的一般噪声环境下，传统的STA/LTA方法出现了5%的误判率、15%的漏判率，正确率为80%.而基于SNR的STA/LTA方法与基于SNR的综合方法均可以无误判、无漏判地提取地震P波到时.并且基于SNR的综合方法(平均偏差0.012 s，标准偏差0.020 s)明显地改善了基于SNR的STA/LTA方法的提取结果(平均偏差0.138 s，标准偏差0.036 s).在3 dB的强噪声环境下，传统STA/LTA方法误判率和漏判率分别上升至10%和65%，其正确率仅为25%.相比之下，基于SNR的STA/LTA方法与基于SNR的综合方法仍可以无误判、无漏判地提取地震P波到时.基于SNR的STA/LTA方法的提取结果比在6 dB一般噪声环境下的提取结果其准确度有所下降.而基于SNR的综合方法仍然保持着较高的提取准确度，其平均偏差0.036 s，标准偏差0.061 s.

图 4为在3 dB强噪声环境下一组地震动模拟信号及其上述三种不同方法的提取结果.其中图 4a为理想的不包含噪声的地震动模拟信号.图 4b为叠加了强噪声的地震动模拟信号.此地震动模拟信号中包含了较强的噪声信号，使得地震P波的到时变得模糊.如图 4c所示，传统的STA/LTA方法，求得的STA/LTA值在地震P波到时处并无显著优势.因此该方法很难通过某一阈值自动地识别地震P波到时，从而导致该方法的误判率和漏判率很高.而基于SNR的STA/LTA方法，如图 4d所示，其求得的STA/LTA值在地震P波到时处优势显著，通过设定的阈值其自动提取结果为10.15 s.如图 4e所示，在基于SNR的STA/LTA方法提取结果的基础上，向前回溯3 s、向后延时0.2 s，计算AIC值.搜索AIC在计算窗口内的极小值，此极小值代表着地震P波到达的精确时刻(10.03 s).

图 4 一组强噪声环境下地震P波自动提取的结果 (a)理想地震动信号；(b)叠加白噪声的地震动信号；(c)传统STA/LTA方法的提取结果；(d)基于SNR的STA/LTA方法的提取结果；(e)基于SNR的综合方法的提取结果. Figure 4 Automatic picking of seismic P waves under strong background noise (a) Ideal ground motion signal; (b) Ground motion signal with white noise superimposed; (c) Picking result of traditional STA/LTA; (d) Picking result of SNR based STA/LTA; (e) Picking result of SNR based comprehensive method.
2.2 汶川余震的P波自动提取

为了进一步验证本文方法对于实际地震P波提取的性能，将基于SNR的STA/LTA方法与基于SNR的综合方法应用于MS8.0汶川地震的30条余震记录，并将两种算法的提取结果与人工的提取结果进行了比较.余震P波到时的自动提取，选取的是余震记录的竖直分量.另外，短窗宽度S设为0.2 s，长窗宽度L设为2 s，地震P波到时的SNR阈值设为5 dB，由公式(10)可求得地震P波到时提取的阈值为4.两种方法的提取结果如表 2所示.

表 2 汶川余震P波自动提取结果 Table 2 Automatic P waves′ picking results for the Wenchuan aftershocks

由表 2可知，当地震P波到时的SNR阈值、长/短时窗等关键参数给定后，不需要反复权衡误报率与漏报率便可以获得较为理想的地震P波到时自动提取结果.在给定参数下，基于SNR的STA/LTA方法地震P波到时自动提取的平均偏差为0.060 s，而标准偏差为0.078 s.而基于SNR的综合方法，其提取的平均偏差更缩小为-0.017 s，而标准偏差也优于基于SNR的STA/LTA方法，具体为0.075 s.

传统的STA/LTA方法在进行地震P波到时自动提取前，通常需要对地震记录进行基线校正或高通滤波等预处理，以消除基线漂移与低频信号干扰.

如图 5a所示，上述序号为13的余震记录就包含了一定的低频干扰.由图 5b可知，采用传统的STA/LTA算法求得余震记录13的STA/LTA序列.该STA/LTA序列在真实的地震P波到时前，出现了一个较大的尖峰.由于此尖峰与真实地震P到时的尖峰大小相当，很难从阈值的角度加以剔除.因此若不进行地震信号的预处理，将导致地震P波到时的自动提取发生错误.而由图 5c可知，对于基于SNR的STA/LTA方法，不仅求得的STA/LTA序列未在真实的地震P波到时前出现较大的尖峰干扰，而且地震P波自动提取的结果5.07 s与人工提取结果5.003 s也非常接近.这说明基于SNR的STA/LTA方法，不会受到基线漂移或低频信号的干扰，具有较强的鲁棒性.

图 5 余震记录13的原始记录与两种STA/LTA值 (a)序号为13的余震记录；(b)传统的STA/LTA值；(c)基于SNR的STA/LTA值. Figure 5 Aftershock No.13 and its STA/LTAs (a) Original seismic data of aftershock No.13; (b) STA/LTAs obtained by traditional STA/LTA method; (c) STA/LTAs obtained by SNR based method.
2.3 地震P波的SNR估计

仍然以上述30组汶川余震记录为例，运用基于SNR的STA/LTA方法，求得其STA/LTA序列在地震P波到时附近的最大值.并且根据公式(11)估算此地震P波的SNR最大值.

由表 3可知，这30组汶川余震记录中P波的SNR大致分布在9.9 dB到43.5 dB范围内.基于SNR的STA/LTA方法在提取地震P波时，选择的SNR阈值5 dB小于30组汶川余震记录P波SNR估计值的最小值9.9 dB.因此5 dB的SNR阈值能够准确地提取到这30组汶川余震的P波到时.

表 3 汶川余震P波的SNR估计 Table 3 Estimated SNRs for Wenchuan aftershocks′ P waves
3 结论和讨论

在大数据量和强噪声环境下，针对地震P波到时自动、快速、准确提取的问题，本文对比了SNR与STA/LTA，通过构建特殊的特征函数，建立了SNR与STA/LTA的内在联系，并且提出了基于SNR的STA/LTA方法，以及改进的STA/LTA与AIC综合的地震P波自动提取方法.

通过理论分析，以及仿真与算例分析，基于SNR的地震P波到时自动提取方法继承了传统STA/LTA方法的优点，同时它们还具有以下新特征：

(1) 基于SNR的STA/LTA方法，建立了SNR与STA/LTA的内在联系，给STA/LTA赋予了明确的物理意义.进而其STA/LTA阈值的选取也建立在严密的数学推导之上，而无需像传统的STA/LTA方法那样通过算例反复权衡误报率与漏报率来获得.

(2) 基于SNR的STA/LTA方法，具有较强的鲁棒性.它能够适应弱噪声、一般噪声，以及强噪声环境，并且误判率和漏判率均较低.而传统的STA/LTA方法在弱噪声环境下能够正确的提取地震P波到时，而对于一般噪声和强噪声环境，其提取结果的误判率和漏判率显著增加.另外，基于SNR的STA/LTA方法，在进行地震P波到时自动提取时，可以不进行基线校正或高通滤波，以提高算法的实时性.而传统的STA/LTA方法必须进行基线校正或高通滤波.

(3) 基于SNR的STA/LTA方法，其求得的STA/LTA可以用于估算地震P波的SNR值.这对研究地震动信号特征以及确定SNR阈值有很大帮助.

(4) 基于SNR的综合方法，提高了基于SNR的STA/LTA方法的提取准确度.它充分利用了基于SNR的STA/LTA方法的鲁棒性，以及AIC方法的准确性.

References
Akaike H. 1974. Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes. Annals of the Institute of Statistical Mathematics, 26(1): 363-387. DOI:10.1007/BF02479833
Allen R M, Kanamori H. 2004. The potential for earthquake early warning in southern California. Science, 300(5620): 786-789.
Allen R V. 1978. Automatic earthquake recognition and timing from single trace. Bulletin of the Seismological Society of America, 68(5): 1521-1532.
Bataille K, Chiu J M. 1991. Polarization analysis of high-frequency, three-component seismic data. Bulletin of the Seismological Society of America, 81(2): 622-642.
Chen H Z, Hou Y Y, He J Y, et al. 2011. The well development of earthquake warning system of Japan. Recent Developments in World Seismology (in Chinese), (4): 10-15.
Evans J R, Allen R M, Chung A I, et al. 2014. Performance of several low-Cost accelerometers. Seismological Research Letters, 85(1): 147-158. DOI:10.1785/0220130091
Han L J, Wong J, Bancroft J C. 2009. Time picking on noisy microseismograms. CREWES Research Report, 21: 30.1-30.13.
Kanamori H, Hauksson E, Heaton T. 1997. Real-time seismology and earthquake hazard mitigation. Nature, 390(6659): 461-464. DOI:10.1038/37280
Liu J H, Liu F T, Xu Y. 2006. Polarization Analysis of Three-component Seismic Data. Process in Geophysics (in Chinese), 21(1): 6-10. DOI:10.3969/j.issn.1004-2903.2006.01.002
Ma Q, Jin X, Li Y S, et al. 2013. Automatic P-arrival detection for earthquake early warning. Chinese Journal of Geophysics (in Chinese), 56(7): 2313-2321. DOI:10.6038/cjg20130718
Munro K A. 2004. Automatic event detection and picking of P-wave arrivals. CREWES Research Report, 18: 12.1-12.10.
Tosi P, Barba S, De Rubeis V, et al. 1999. Seismic signal detection by fractal dimension analysis. Bulletin of the Seismological Society of America, 89(4): 970-977.
Wang X, Fu J H, Tang C P, et al. 2018. A P waves' automatic picking by detecting the changes of seismic signals' stationary random process through similarity analysis. Soil Dynamics and Earthquake Engineering, 115: 225-231. DOI:10.1016/j.soildyn.2018.08.006
Yin H T, Liu X Q, Li J, et al. 2012. Current Earthquake Early Warning technology and it's development in China. Earthquake Research in China (in Chinese), 28(1): 1-9.
Zhang H C, Jin X, Li J, et al. 2013. Progress of research and application of earthquake early warning system (EEWs). Progress in Geophysics (in Chinese), 28(2): 706-719. DOI:10.6038/pg20130219
陈会忠, 侯燕燕, 何加勇, 等. 2011. 日本地震预警系统日趋完善. 国际地震动态, 4: 10-15.
刘建华, 刘福田, 胥頤. 2006. 三分量地震资料的偏振分析. 地球物理学进展, 21(1): 6-10. DOI:10.3969/j.issn.1004-2903.2006.01.002
马强, 金星, 李山有, 等. 2013. 用于地震预警的P波震相到时自动拾取. 地球物理学报, 56(7): 2313-2321. DOI:10.6038/cjg20130718
殷海涛, 刘希强, 李杰, 等. 2012. 现今地震预警技术及其在国内发展状况的探讨. 中国地震, 28(1): 1-9.
张红才, 金星, 李军, 等. 2013. 地震预警系统研究及应用进展. 地球物理学进展, 28(2): 706-719. DOI:10.6038/pg20130219