地震P波和S波自动识别方法的研究与应用

2019-01-09冯红武颜文华张炜超王建昌

防灾科技学院学报 2018年4期

冯红武，颜文华，张炜超，王建昌

(陕西省地震局，陕西西安 710068)

0 引言

震相识别是现代地震学研究中的重要课题之一，也是地震定位、地球内部结构、震源机制等一系列研究的基础。随着数字化地震台网和地震科学台阵大规模建立和发展，仅靠人工捡拾震相已满足不了实际需要。震相的自动捡拾不但可以提高地震观测数据的处理效率，而且提升了地震台网的快速响应能力，也是地震预警技术的一个关键环节。因此，研究更加有效和准确的震相自动识别方法是我们亟需解决的一个科学问题。

一般来说，地震事件的自动检测是震相自动识别的基础。目前，长短时平均比方法(STA/LTA方法)是一种既可判断有无地震事件，又能自动捡拾初至震相到时的方法。这种方法因算法简单、捡拾效率高、适应性强等优点，已被广泛地应用于地震波初至到时的自动捡拾[1-7]。但该方法缺点是当地震波信噪比较低或初动不清晰时，自动捡拾震相的效果会较差。除此之外，主要的震相自动识别方法还有分形分维法[8-9]、基于小波变化的主成分分析法[10-11]、人工神经网络法[12-14]、基于幅值和频率的P波识别方法[15-16]等，都取得了一定的效果。

STA/LTA方法和AIC准则相结合的多步骤捡拾方法是目前P波震相自动捡拾的主要方法，该方法实质是综合了几种单一算法的优点，避免单一识别方法受使用范围和条件限制的缺点，提高了自动识别的精度，是一种综合分析法[17]。本文在STA/LTA方法和AIC准则相结合的多步骤捡拾方法的基础上，针对地震波的变化特征，将地震信号的瞬时振幅和瞬时频率作为基本参数，重新构建识别P波和S波的STA/LTA方法的新特征函数，提出了基于新特征函数的STA/LTA和VAR-AIC多步骤捡拾P波和S波的方法，以期提高震相捡拾的准确性和可靠性，最后对流动地震科学台阵记录的地震震相进行离线拾取，来检验该方法的识别效果。

1 STA/LTA的特征函数构建

地震波到达后，台站记录到的地震波初动振幅和频率明显区别于台站记录的背景噪声，而且后至震相的振幅或频率较初至波也有明显的变化。为了准确捡拾P波和S波，STA/LTA特征函数应能充分体现地震波到达后的振幅增大和频率改变(变大或变小)的变化特征，同时考虑到P波和S波的变化特征一般分别在垂直向和水平向突出。为此，本文将瞬时振幅和瞬时频率作为基本参量，构建了拾取P波和S波的特征函数。

由于地震波是具有时变特性的非线性、非平稳信号，我们利用美籍华裔科学家Huang等人提出的希尔伯特-黄变换方法(HHT)[18]计算地震波的瞬时频率。HHT方法可将复杂信号分解为有限个固有模态分量，从而赋予了瞬时频率合理定义和物理意义。

P波特征函数：

(1)

S波特征函数：

(2)

2 震相自动捡拾方法

2.1 STA/LTA法

STA/LTA方法是一种能量方法，其原理是用信号的短时平均值STA和长时平均值LTA的比值来反映信号水平或能量的变化。当地震波到来时，STA要比LTA变化的快，STA/LTA值相应会有一个显著的增加，若该值大于设定的触发阈值，即可判定地震信号的到来。本文采用递归STA/LTA算法，计算公式如式(3)、式(4)：

(3)

(4)

式中STAi和LTAi分别为地震信号在时刻点i的短时平均值和长时平均值，CF(i)为地震信号在i时刻点的特征函数值，Nsta和Nlta分别为短时平均值和长时平均值的时间窗内所包含的记录点数。本文选取新构建的P波和S波的特征函数。

2.2 基于AIC准则的震相识别方法

由于STA/LTA方法仅能粗略捡拾震相的到时，而且识别出的震相到时往往滞后于实际到时，即该方法得到的结果只是震相到达的大致位置。

1973年日本学者Akaike提出一个基本信息量的定阶准则，即AIC准则[19]。震相到时检测的自回归(AR)技术是假设震相到时前后的地震记录是两个不同的稳态过程，基于这个假设，Sleeman和orild于1999年提出了AR-AIC方法[20]用于对地震震相进行识别。不同于AR-AIC方法，Maeda[21]建议直接根据地震波形数据而不需要计算AR系数来得到AIC函数，该方法被称为VAR-AIC方法。对于地震记录x(i)(i=1,2,…,L)，将AIC检测器定义为：

AIC(k)=klog{var(x[1,k])}+

(N-k-1)log{var(x[k+1,N])}

(5)

式中，k的范围是地震图某窗长内所有的采样点，AIC函数局部最小值对应的时刻即为震相到时。

不同于在整个地震波信号上直接应用AIC准则，本文是利用STA/LTA方法捡拾得到的粗捡拾点，在固定窗内对该点进行AIC精确捡拾。

2.3 基于新特征函数的STA/LTA和VAR-AIC多步骤捡拾方法

本文对P波捡拾主要分为两步：(1)粗略捡拾：利用STA/LTA方法对P波粗捡拾，特征函数选取本文新构建的P波特征函数，为了避免不同阈值造成的误触发或漏触发现象，直接取STA/LTA的最大值，最大值点出现时刻即为P波的粗捡拾点(即P波到达的大致位置)，记为tP；(2)精确捡拾：以P波粗捡拾点tP为中心，对其前推和后推一定时间窗，前后时间窗各取Δt时间长度，对时间窗长[tP-Δt,tP+Δt]的数据应用VAR-AIC方法对P波到时进行精确捡拾。

由于S波容易受到P波尾波以及各种反射波的干扰，与P波捡拾相比，S波捡拾过程通常较复杂而且可靠性低。为了提高对S波捡拾的可靠性，对S波的捡拾主要分为两步：(1)粗略捡拾：首先应用STA/LTA方法对S波进行粗略捡拾，特征函数选取本文新构建的S波特征函数。长、短时间窗从捡拾到的P波到时点开始，同步向前逐点滑动对S波特征函数曲线进行处理，将每个点的计算结果作为该点的STA/LTA值。同时为了解决阈值选择的问题，直接取STA/LTA值的最大值，最大值点出现的时刻即为S波的粗捡拾到时(即S波到达的大致位置)；(2)精确捡拾：我们将P波段和S波段认为是两个不同的稳态过程，应用AIC准则对S波到时进行精确捡拾。同样以S波粗捡拾的到时tS为中心，对其前推和后推一定时间窗，前后时间窗各取Δt时间长度，对时间窗长[tS-Δt,tS+Δt]的数据应用VAR-AIC方法对S波到时进行精确捡拾。

3 实际应用

为验证本文方法的有效性和准确性，我们将以上识别方法应用于甘东南流动地震科学台阵的地震记录。由于文中的震相识别方法旨在用计算机自动识别来代替人工处理地震科学台阵大量地震数据，我们感兴趣的是对近震Pg和Sg震相的自动捡拾，按照甘青区域地震波走时表[22]，选择震中距在360km以内台站记录的波形数据。以台阵记录的2011年11月2日甘肃岷县ML3.6地震为例，对其进行震相识别，记录地震较清晰台站数为29个，P波和S波有效数分别为29个和28个。

3.1 P波自动捡拾结果及分析

在应用STA/LTA方法对P波粗略捡拾时，选取时间窗长，主要基于以下考虑：长时平均窗LTA描述的是相对于待捡拾信号的背景噪声水平的平均大小，其取值应能反应背景噪声水平；短时平均窗STA描述的主要是信号幅值的瞬时变化，一般取长于待捡拾信号的几个周期左右，时间窗太短则对短周期的干扰更敏感，容易产生误触发，如果时间窗太长则显示不出待捡拾信号的瞬时特征，容易产生漏触发。因此，根据已有研究结果[7,16,23]，并且多次反复试验，取长时间窗30s，短时间窗0.5s。对于触发阈值的选取，若阈值过大，就有可能会出现漏触发，若阈值太小，就对很多干扰会误触发，所以在选取阈值时，要考虑台站的背景噪声、地震震级、地震的远近以及仪器频带范围等，经过试验，设定的触发阈值为5，且基本满足要求。

精确捡拾时，考虑到地震数据原始采样率为40Hz，以及多次反复试验，以P波粗捡拾点tP为中心，对其前后各取2s，对时间窗[tP-2s,tP+2s]的数据应用VAR-AIC方法对P波到时进行精确识别。对该地震事件P波进行识别，P波自动捡拾结果与人工识别结果见表1和图1，P波震相检测曲线见图2。可以看出，该方法对P波的识别效果较好，人工识别和自动识别误差的平均值为0.034s，误差很小。

3.2 S波自动捡拾结果及分析

在应用STA/LTA方法对S波粗略捡拾时，同样基于以上考虑以及相关研究成果[23]，我们取长时间窗4s，短时间窗取0.5s，特征函数的触发阈值设定为5。精确捡拾时，考虑到地震数据原始采样率为40Hz，以及多次反复试验，以S波粗捡拾点tS为中心，对其前后各取3s，对时间窗[tS-3s,tS+3s]的数据应用VAR-AIC方法对S波到时进行精确识别。对该地震事件的S波进行自动识别，S波自动捡拾结果与人工识别结果见表2和图3，S波震相检测曲线见图4。可以看出，除308台站识别误差较大，该方法对S波的识别效果较好，人工识别和自动识别误差的平均值为0.105s，误差很小。按照S波识别误差小于0.2s的精度要求，识别有效率高达96.43%，有效识别平均误差0.091s。

表1 P波到时识别结果比较

图1 P波自动捡拾与人工捡拾的结果比较误差图Fig.1 The error graph of comparison of results between manual and automatic P-wave pickups

图2 112台站的P波捡拾曲线Fig.2 P-wave picking curve of Seismic station 112

台站代码人工识别/s自动识别/s误差台站代码人工识别/s自动识别/s误差101 33.150 33.1500.000210 16.825 17.000-0.175102 30.850 30.8250.025211 14.200 14.325-0.125 106 30.075 30.0750.000212 14.275 14.175 -0.100107 27.200 27.325-0.125213 15.325 15.500-0.175109 22.550 22.625-0.075214 19.025 19.125 -0.100 111 22.150 22.250-0.100215 22.050 22.150 -0.100 112 21.600 21.675-0.075217 29.700 29.850 -0.150113 22.675 22.700-0.025308 38.325 37.825 0.500 114 30.450 30.400 0.050312 27.700 27.825-0.125115 32.325 32.425 -0.100313 28.675 28.825 -0.150201 47.275 47.225 0.050314 29.875 30.075 -0.200 205 30.875 30.950-0.075701 14.375 14.475 -0.100 208 21.075 21.150-0.075702 14.675 14.800-0.125 209 16.750 16.7250.025703 18.425 18.450-0.025