0 引言

作为大震短临前兆之一的前震活动特征(陈运泰，2007)，能够有效地反映大震成核过程，是地震预测与风险评估领域长期关注的重要课题(Bouchon et al，2011；周少辉等，2016)。我国曾多次实现了短临预测，在具有明显减灾实效和社会显示度的地震短临预报中，直接前震的贡献是显见的，如1995年7月12日云南孟连7.2级和1999年11月29日辽宁岫岩5.6级地震(王林瑛等，2005)。而地震目录中所包含的前震信息受检测方法的影响，存在一定的局限性。如 图 1所示，根据中国地震正式观测报告，在2010年4月14日青海玉树7.1级地震(以下简称“玉树地震”)前约2h，发生一次4.8级前震(倪四道等，2010；陈学忠等，2012)，自这次前震发生前12h到主震发生之间，仅报告了14个前震事件(4.8级前震之前仅有1个)；2013年8月31日云南香格里拉、德钦、四川得荣交界5.9级地震(以下简称“香格里拉地震”)前约75h发生一次5.2级前震(赵小艳等，2015)，这次前震发生前12h到主震发生之间，报告前震事件为288个(5.2级前震之前为0)；2014年2月12日新疆于田7.3级地震(以下简称“于田地震”)前约31h发生一次5.4级前震(郑建常等，2015)，这次前震前12h到主震发生之间，报告前震事件为9个(5.4级前震之前为0)。而通过观测距离主震震中距最小台站的三分量波形数据，如青海玉树台(YUS，震中距30km)、云南中甸台(ZOD，震中距50km)和新疆于田台(YUT，震中距56km)，记录到的地震事件个数远不止于此。

传统人工震相检测试图尽量多地检测小震信号，所耗费的人力和时间成本较高。本文尝试人工检索玉树地震前14h玉树台的连续波形，获取78个地震事件P波初至到时的时间约为3h。随着台站密度的增大，数据量不断增加，完全采用人工震相拾取的方法已不能满足前震检测需求。近几十年来，地震科学家们凭借数字地震记录的自动处理经验，为我们提供了大量震相自动检测方法，提高了检测效率和精度，如时间域或频率域的能量瞬变方法(Withers et al，1998；Vassallo et al，2012)、自回归(autogressive，AR)方法(Leonard et al，1999)、高阶统计方法(Baillard et al，2014)、神经网络方法(Gentili et al，2006)以及小波变换方法等(Bogiatzis et al，2015)。

为探索自动拾取震相算法在前震信号探测方面的应用，补充前震活动信息，完善地震目录，从而更好地研究主震破裂前能量辐射的过程，本文采用短长时窗比、多频带滤波、深度学习和模板匹配过滤四种地震信号自动检测方法，对玉树地震、香格里拉地震和于田地震发生前数十小时的连续波形数据进行处理，并对不同方法进行对比分析，找到合适的前震自动检测方法。

1 数据和方法

按照国际时间，分别截取YUS台2010年4月13日9时40分—4月13日23时50分(约14h)、ZOD台2013年8月27日8时40分—8月31日0时10分(约87h)和YUT台2014年2月10日14时15分—2月12日9时45分(约43h)的三分量连续波形记录进行震相到时自动检测。

1.1 短长时窗比

短长时窗比方法(Short Term Average/Long Term Average，简称STA/LTA)作为以能量变化检测震相到时的经典算法之一(Withers et al，1998)，通过检测与波形振幅有关的特征值(如振幅的绝对值或平方)在短时窗内和长时窗内均值比的连续变化，当该值超过给定阈值时，认为是一次震相到达。其中，短时窗用于截取微弱的有用信号，长时窗内信号的平均值则反映背景噪声水平。本文采用Allen(1978, 1982)给出的特征方程，其对高频成分的变化更为敏锐。假设Y(i)为第i个采样点的振幅，特征函数表达式为

$ C F(i)=Y(i)^2+K(Y(i)-Y(i-1))^2 $

(1)

其中，K为权重因子。

对三次地震的最大前震事件波形进行时频分析，以玉树地震为例，图 2展示了玉树地震4.8级前震事件在YUS台垂直分量的记录，可知初至P波的主要能量集中在1~10Hz。因此，本文首先对原始数据垂直分量进行1~10Hz的带通滤波，然后分别采用STA/LTA为0.2s/5s、0.22s/20s、1s/20s和2s/80s四组参数计算短长时窗比的特征函数，K值取3，并分别取4倍均值、4倍均值、3倍均值和3倍均值作为阈值，进行了触发判定。图 3显示了YUS台BHZ分量波形片段以及相应不同STA和LTA值计算的特征函数。

1.2 多频带滤波

Lomax等(2012)提出一种基于多频带瞬时能量变化的震相自动检测方法FilterPicker，该方法对数据在多个窄带内进行滤波。用户只需给出一个最小频带值，FilterPicker即可根据数据采样率自动调整多个滤波频带，分别计算各个频带内振幅特征函数，最终合并为一个综合特征函数，触发阈值由三个参数进行计算。该方法对检测信噪比较低的地震信号有较大优势。

本文采用Chen等(2016)基于FilterPicker方法进行修改后形成的多频带滤波震相自动检测工具FBpicker。FBPicker算法的参数缺省值被应用至广泛的地震数据(Dickey et al，2019)。在对原始波形去平均和去趋势等预处理后，通过倍频程为地震记录设置多个频带，由指定的最小频带和仪器采样率共同决定滤波范围(图 4)。此外，FBpicker可由用户定义数据截断比例。

根据Lomax等(2012)的定义，FBpicker也和FilterPicker采取了相同的方法计算E_n，即

$ E_n[i]=B F_n[i]^2 $

(2)

其中，BF_n[i]为第n个频带第i个采样点的振幅值(图 5)。FBpicker通过均方根(rms)和标准差(std)的方式计算每个频带的特征函数，即

$ C F_n^{\mathrm{rms}}[i]=\frac{E_n[i]}{\operatorname{rms}\left(E_n[i-1-l: i-1]\right)} $

(3)

$ C F_n^{\mathrm{std}}[i]=\frac{E_n[i]-\operatorname{mean}\left(E_n[i-1-l: i-1]\right)}{\operatorname{std}\left(E_n[i-1-l: i-1]\right)} $

(4)

其中，l为采样点窗长。两种计算方式可由用户进行指定。FBpicker采取了和FilterPicker相同的方式计算最终的特征函数，其对每个采样点在不同频带上的最大绝对值CF_n[i]进行归一化，作为该采样点最终的CF值(图 6)。FBpicker由移动指定长度时窗内的均值倍数σ计算动态的触发阈值(图 6)。此外，FBpicker为减少误触发采取一系列措施，如设置t_up，规定在该时间段内不连续拾取。

本文采用FBpicker对三组数据的垂直分量进行处理，选择的缺省参数包括：计算CF_n的窗长t_long为5s，首个倍程频带范围的中心频率为1s，动态阈值的时窗长为10s，σ为5，t_up为2s。对CF_n的计算方式，分别采用均方根(rms)和标准差(std)两种模式进行对比。

1.3 基于深度学习的广义震相检测

随着深度学习在人工智能领域的广泛应用，越来越多的地震学家也开始尝试训练人工智能地震震相检测、到时确定、地震定位等深度学习的神经网络。如Ross等(2018)利用南加州台网人工标注的450万个地震震相波形数据(包括150万个P震相、150万个S震相和150万个噪声数据)，构建基于深度学习的广义震相检测模型(GPD)。该方法对以到时为中心的4s时间序列进行卷积、降采样和激活处理，提取数据特征，然后将该特征值作为输入层，通过4个卷积层扫描地震波形，提取局部特征，如波形的形状和时序模式。全连接层则将这些特征映射到分类结果(P波、S波或噪声)，同时使用ReLU激活函数和批量归一化来提高训练效率。

由于训练波形为震中距小于100km的3~20Hz带通滤波结果，因此本文对三组数据的三分量数据进行相同频带的滤波，滤波结果作为输入，经过处理后分别选择P波或S波概率大于95%和大于98%的震相作为检测结果(图 7)。

1.4 模板匹配过滤

利用互相关技术的模板匹配过滤方法(Matched Filter Technique，MFT)已经在余震检测(Peng et al，2009)、构造震颤与低频地震检测(Shelly et al，2007)等地震观测的多个领域得到应用。Kato等(2012, 2014)也用该方法检测了强震震前的慢滑移事件。该方法采用已知地震事件波形数据片段作为初始模板，将其与目标连续波形逐段进行互相关计算，计算公式如下

$ C C=\frac{\sum_{t 1}^{t 2}\{[X(t)-\bar{X}] \times[Y(t)-\bar{Y}]\}}{\sqrt{\sum_{t 1}^{t 2}[X(t)-\bar{X}]^2 \times \sum_{t 1}^{t 2}[Y(t)-\bar{Y}]^2}} $

(5)

其中，t1、t2为互相关时窗的起始、终止时间，X(t)为模板的时间序列，Y(t)为连续波形中的目标时间序列片段。将互相关系数达到设定阈值的片段进行叠加平均，得到信噪比更高的波形作为新的模板继续进行匹配，如此循环多次，检测到更多的事件。

本文利用中国地震正式观测报告 ^①中国地震台网中心. 中国地震正式观测报告. 中的Pg和Sg震相到时，对YUS台进行1~10Hz带通滤波后截取到时前1s至后3s、对ZOD台和YUT台1~10Hz带通滤波后截取到时前1s至后4s的波形作为初始模板，与目标连续波形相同时窗长度的片段进行模板匹配。设置最小互相关系数为0.7，即当互相关系数达到0.7时，选取时窗内互相关系数最大的点为触发时刻，最大循环次数为10次。

① 中国地震台网中心. 中国地震正式观测报告.

当初始模版具有一定相似性时，如YUS台的14个Pg震相初始模板中有10个较为相似，在匹配过滤后，其中每一个模板匹配到29~30个震相，因此看似近300个检测到的事件其实多为重复事件。为此，我们对三组数据的处理结果统一进行了去重处理。

此外，为了对比初始模板的选择对模板匹配过滤方法处理结果的影响，采用GPD概率大于95%的P波和S波震相作为初始模板，再次利用模板匹配过滤方法进行进一步检测。

2 检测结果及分析

分别按照地震信号(T)和非地震信号(F)对上述四种方法及其不同参数选择下的震相检测个数进行统计，结果如 表 1所示。

将四种方法不同参数条件下检测到的三组波形中的所有地震震相数据进行整合，共计检测到完整地震震相目录为：YUS台179个，其中包含80个P波和89个S波，平均信噪比为1.0；ZOD台701个，其中包含338个P波和363个S波，平均信噪比为1.3；YUT台317个，其中包含131个P波和186个S波，平均信噪比为1.4。为了对比不同方法所测得的震相信号特征，表 2列出了每种方法一种参数选择情况下检测到的地震震相信号的平均信噪比和平均地动速度峰值对数值。

对于每组数据的每种方法，按照以下公式分别计算其漏检率和误检率(图 8)，即

$ \text { 漏检率 }=\frac{\text { 完整震相目录个数 }- \text { 地震信号个数 }}{\text { 完整震相目录个数 }} $

(6)

$ \text { 误检率 }=\frac{\text { 其他信号个地震信号个数 }}{+ \text { 其他信号个数 }} $

(7)

此外，我们计算了3个台站检测到的P波震相4s时长的波形互相关系数矩阵，结果如 图 9所示。

较低的漏检率和误检率是考量一个震相自动检测方法是否优秀的重要参数。分析四种自动检测方法在3个地震前震信号中应用的效果，归纳总结为：

(1) STA/LTA方法的检测数量随着短窗长度由0.2s变为1s(长窗20s)而陡降。在短长窗长为0.2s/5s的情况下，3个台站的误检率均极高，即几乎将所有非地震脉冲均判为地震震相。随着短窗长增加，伴随着误检率降低，漏检率整体增加，显示出漏检率和误检率之间的平衡。对比3个台站，YUT台检测到的地震信号在信噪比和平均地动速度峰值较大的情况下，表现出较低的漏检率。

(2) FBpicker作为多频带滤波检测方法之一，其计算特征函数的两种模式中，std比rms模式漏检率有所降低，但误检率在较高水平上继续增大。总体来看，该方法对噪声抑制能力不足，误检率仍然很高，且漏检率没有明显优势。而YUT台的记录也由于其较高的信噪比和平均地动速度峰值对数值而保持了较低的漏检率。可见STA/LTA和FBpicker均以较高的误检率换取较低的漏检率，这一特点使这两种方法更适用于震相的初步检测。

(3) GPD方法在漏检和误检方面表现出较大的优势，阈值选择0.95时，漏检率较低，误检率有较好的控制；阈值选择0.98时，漏检率略升高，误检则显著下降。可见提升阈值能在可接受的漏检率增加幅度内换来较明显的误检率下降，是在漏检和误检之间比较好的折中。对比3个台站，ZOD台检测到的地震信号在较低平均地动速度峰值对数值及中等信噪比的情况下具有最低的漏检率，推测GPD方法受信号信噪比和地动速度峰值的影响较前两种方法小。

(4) 由于前震信号的相似性，模版匹配过滤方法仅利用正式观测报告中的少量震相进行匹配，也可检测出近一半或更多的震相。较高的相似性维持了较低的漏检率(如YUS台和YUT台)，而相对较低的误检率显示出该方法在前震这种特殊信号中的优越性。当使用GPD检测结果中的震相作为模板进行匹配时，将漏检率大幅度降低至几乎零漏检，同时还保持了较低的误检率。

3 结论

因前震信号可能存在一定的相似性，使得用模版匹配过滤方法对前震震相进行自动检测成为可能，但是仅使用地震目录数据又凸显了该方法的局限性。深度学习方法能够捕捉到足够多的地震信号，对短长时窗比方法和多频带滤波方法检测结果进行补充验证，将此结果作为模版进行匹配过滤，同时充分利用GPD对P波和S波震相的快速识别，可以进一步完善前震目录，为前震活动性研究提供更为充足而有效的数据。