余震自动识别技术研究进展
李力 张建国 李盛乐
摘要：基于余震识别自动化对地震监测、灾后应急、确定发震构造、分析预报和科学研究的重要意义，总结和分析了余震自动识别技术的现状和未来发展。首先，总结单台震相到时自动拾取方法及其特征函数，包括单特征方法、互相关方法和机器学习方法。其次，分别总结基于震相走时的常规识别方法和基于特征函数的偏移叠加方法。最后，对余震识别技术的未来发展进行分析和展望，认为偏移叠加方法的性能改进将成为未来余震识别的重点发展方向，深度学习方法作为一种新技术，将在包括震相拾取和余震识别在内的地震数据处理领域发挥重要作用。
关键词：余震；自动识别；震相自动拾取；走时定位；偏移叠加
中图分类号：P315.61 文献标识码：A 文章编号：1000-0666（2018）01-0001-13
0 引言
破坏性地震发生之后，随之而来的是大量的余震，从连续波形中识别余震是开展余震震级测定、震源机制解反演等后续监测工作的基础。余震信息是评估灾情和指导救援力量部署（张苏平等，2013；郑韵，2015）、确定发震构造和揭示隐伏断层（黄媛等，2008；徐锡伟等，2008，2013a，b）、判定震后趋势（庄建仓，马丽，2000；宋金，蒋海昆，2009；余娜等，2016）的重要依据，也是震后地下速度结构变化、余震迁移等科研工作的基础数据。人工识别余震是一项十分耗时且枯燥的工作，结果中不可避免地会引人主观因素，随着观测台网的建设，数据的产出量越来越大，完全由人工处理变得越来越不现实。大震之后，识别高频度的余震需要持续投入大量的时间和人力，这让有限的监测资源更加捉襟见肘，数据产出时效也受到限制。
余震识别的自动化有重要的应用和科研价值，但是现有技术仍然有很多问题亟待解决和完善，地震监测领域也没有专门的余震自动识别系统。为了提高对余震自动识别技术的关注和投入，2017年5月中国地震局地球物理研究所和阿里云联合举办余震捕捉AI大赛，吸引来自国内外的近千支队伍报名参赛，一时间余震自动识别成为地震领域的热门话题和技术热点。
余震发生后，其能量以地震波的形式向外传播，当震相到达各个观测台站时，会引起波形记录的变化。自动识别余震时，首先需要从连续观测波形中提取特征，并自动拾取震相到时，然后联合多个台站的数据确定余震信息。对于震相清晰的独立地震，可以使用单台波形进行识别，并初步确定震源信息。但是余震在时间上高度密集，并且震级普遍较小，记录到的信号受主震尾波和其他余震的强烈干扰，多台联合是提高余震识别可靠性的必要技术。同时，余震在构造上与主震有密切联系，在空间上沿发震构造密集展布，发震机理有较高的相似性，这对于震相识别和震相联合是可以利用的有利条件。传统震相拾取方法所使用的特征函数在偏移叠加和人工智能等新技术中也有非常广泛的应用，因此本文从单特征方法出发，总结几种重要的震相自动拾取方法及其特征函数。然后介绍基于震相走时定位的常规方法的发展，并分析其在余震自动识别中的适用性，最后对基于特征函数的偏移叠加方法进行综述，并在此基础上展望余震自动识别技术的未来发展方向，以期对开展相关研究工作有所帮助。震相特征的提取是走时定位和偏移叠加方法的基础，为便于理解，对STA/LTA、VAR-AIC、峰度、极性分析等单台特征提取算法进行了简单实现，相关示例代码发布于github版本库https：//github.com/dualli/trigger。
1 震相自动拾取方法及其特征函数
当地震震相到达观测台站时，会改变记录信号的能量、频谱、极性和整体形态等特征，这些变化可以通过不同的特征函数提取，当特征函数达到设定阈值时拾取相应的震相到时。多个台站的特征函数和震相信息成为后续地震识别和定位的基本数据。目前，震相自动拾取方法种类繁多，前人对这些方法按照时间域、频率域、时频域、综合方法等进行过分类综述（Withers el al，1998；周彦文，刘希强，2007； Nippress el al，2010；王彩霞等，2013：Ross，Ben-Zion，2014；刘翰林，吴庆举，2017）。本文选取常用的长短窗、赤池信息准则、峰态和偏态、偏振分析等单特征方法和基于整体形态的互相关方法、基于多特征的机器学习方法，分别阐述基本原理和特征函数演化，并比较分析其特点。
1.1 长短窗均值比方法
长短窗均值比方法（Short-Term-Average/Long-Term-Average，简称STA/LTA）基于震相到达时特征值的平均值在短窗口比长窗口变化快的特点拾取震相。STA/LTA方法的发展历史久远，应用广泛，形成了反映波形不同变化特点的众多特征函数。Stevenson（1976）首次提出长短窗方法时，使用振幅的绝对值|Y（i）|作为特征函数，其中Y（i）为第i个采样点的振幅值。但是使用|Y（i）|或者Y（i）2作为特征函数时，长短窗比值只对振幅的变化敏感，因此Allen（1978，1982）引入了特征函数CF=Y（i）2+K（Y（i）-Y（i-1））2，其中K为根据采样率和台站噪声特征设置的权重常数，该特征函数的引入使得长短窗均值比同时随振幅和频率变化而变化（图1）。Bear和Kradolfer（1987）使用基于包络函数平方的特征函数，同时引入动态阈值，既避免了瞬变噪声引起的误触发，也适用于拾取地方震到远震的不同震相，同时对低信噪比微震信号也变得更敏感。为了更好地拾取各种微弱的次生震相，Earle和Shearer（1994）定义了基于Hilbert变换的包络函数作为特征函数。高淑芳等（2008）也提出了反映振幅和频率变化的新特征函数。余建华等（2011）引入加权系数法减少震相漏捡率。
1.2 赤池信息准则方法
赤池信息准则方法（Akaike Information Criteri-on，简称AIC）沿每個采样点将信号分为2段，当分割方法与2个独立过程对应时，对应的AIC值最小，即全局最小值对应震相到时。AIC方法也有非常多的变种，以下介绍几种较为常见的AIC函数。AR-AIC（Sleeman，Van Eek，1999；Leonard，Kennett，1999）使用自回归方法计算AIC值，即分割点对应2个自回归过程，自回归阶数需要反复试验得到。Maeda（1985）提出的基于二阶统计量的VAR-AIC方法，可以直接从地震信号计算AIC值（图2）.0刘希强等（2009）提出TOC-AIC方法，使用三阶累积量代替方差来压制高斯有色噪声，在低信噪比情况下有更高的精度。赵大鹏等（2012）提出峰度一赤池信息准则方法（Kurtosis-AIC），使用峰度作为特征函数替代VAR一AIC中的方差，有效减少了误报率和漏报率。
1.3 偏态/峰态方法
偏态/峰态方法（P Arrival Identification-Skewness/Kurtosis，简称PAI-S/K）由Saragiotis等（2002）提出，当震相到达观测台站时，记录数据的分布由具有高斯性的随机噪声向非高斯性转变。峰态（Kurtosis）和偏态（Skewness）是反映数据分布形态的高级统计量（Higher-Order Statis-tics，HOS），数据服从高斯分布时峰态等于3，分布变宽时小于3，变窄时大于3，同理偏态在数据对称分布时等于0，左偏时小于0，右偏时大于0（图3）。刘劲松等（2013）引入了新的到时判定算法，改进了PAI-SIK方法的拾取精度。Baillard等（2014）使用滑动窗峰态函数和极化分析识别P、S波。峰态和偏态值方法有较好的低信噪比拾取能力和精度，属于常用的震相拾取方法之一（Panagiotakis et al，2008；Kiiperkoch et al，2010；Nippress et al，2010；Hibert et al，2014；Langet etal，2014；Ross，2016）。
1.4 偏振分析方法
常用的2种偏振分析方法是协方差矩阵方法和奇异值分解方法，在多数情况下二者效果接近。偏振分析的思想最早由Flinn（1965）提出，他使用协方差矩阵的本征值分解来识别S波偏振，进而约束震源机制解。Montalbetti和Kanasewieh（1970）首次详细阐述如何用对称本征问题来进行偏振分析。分解有限长度三分量数据的协方差矩阵，得到本征值和本征向量矩阵，可以用来计算质点运动的多种偏振特征，沿着时间轴滑动窗口，重复计算就可以得到偏振特性随时间的变化（图4），进而完成P或S波能量压制、波形分解等操作，是提升S震相拾取率的常用方法（Jurkevies，1988；Robertset al，1989；Bai，Kennett，2000：Diallo et al，2006：Amoroso et al，2012；Ross et al，2014，2016）。Jackson等（1991）首次引入奇异值分解方法（Singular Value Decomposition，简称SVD）进行偏振分析，一步即可完成数据的分解操作。在沿时间轴滑动时，传统的偏振分析方法要使用窗口内数据重复求解协方差矩阵，Rosenberger（2010）引入迭代奇异值分解法，通过追加数据矢量来更新解，提升了奇异值分解方法的计算效率，使其适用于实时数据处理。Kurzon等（2014）基于迭代奇异值分解方法将地震波形实时分解成P分量和S分量，然后使用视出射角、奇异值、信噪比拾取器拾取P、S震相，并且尝试多种指示变量来增强P、S波的拾取。
1.5 互相关方法
互相关系数是描述2段信号相似程度的量，地震波形的相似性主要由震源机制和经过的路径决定，因此互相关方法可以将事件波形作为模板，从连续波形中寻找相似地震（Poupinet et al，1984；Rubin，2002；Schaff，Richards，2004，2011），该方法对于拾取低信噪比震相和震颤等特殊信号非常有用（Shelly et al，2007；Peng，Gomberg，2010），能检测到比正式目录震级低1级的地震事件（Schaff，2010；Schaff，Waldhauser，2010；Gib-blon，Ringdal，2006）。Yang等（2009）提出聯合三分量的滑动窗互相关方法（Sliding-WindowCross-correlation，简称SCC），比单独的三分量互相关方法拾取到更多低信噪比震相。
互相关方法必须先制作好模板，并且灵敏度高度依赖模板和待检测信号的相似性，Brown等（2008）提出自相关方法来，能够从连续波形中提取模板，但是计算十分耗时。Yoon等（2015）基于数据挖掘算法提出（Fingerprint And SimilarityThresholding，简称FAST），用关键特征代替波形，速度是自相关方法的140倍，能够在2h内处理单台一周的连续波形数据。FAST方法有效地提升了模板的制作速度，但是运算量依然会随数据长度的增加而快速增长，处理长达数月甚至数年的余震活动数据时，运算效率依然是需要解决的问题。Shelly等（2007）为了拾取更多的低频地震，尝试降低时间精度来延拓互相关系数的时间范围，提高对空间位置变化的适应性。
1.6 机器学习方法
机器学习（Machine Learning，简称ML）是一类从数据中寻找潜在模式，进而提升性能的庞大算法集合的总称，大致分为经典机器学习算法和人工神经网络算法（Artificial Neural Networks，简称ANN），这2种算法都包括样本数据训练模型和模型预测新数据的2个阶段。ANN是ML的子集，用简单数学函数构建神经元网络，模拟高度复杂的非线性关系，孔令文（2014）按照网络结构进行了分类。ML在地震数据处理的很多领域均有应用，但是主要集中在处理分类问题。Rawles和Thurber（2015）在微震定位中使用邻域法测定潜在震相的类别和到时。Kortstrom等（2016）使用支持向量机方法，基于多频带不同震相段的特征函数值，从监测系统自动识别的地震中分离出爆破和误触发事件。Maggi等（2017）使用随机森林方法实现1套自动地震事件分类器，使用多个台站的波形和频谱特征组合，对留尼旺火山岛的地震事件进行分类。曲均浩等（2012）使用反向传播神经网络分类近震和远震。Esposito等（2013）使用多层感知机分类意大利南部斯特龙博利火山的滑坡和其他事件。Horstmann等（2013）使用自组织映射神经网络检测和分类地震。Giudicepietro等（2017）使用多层感知机分类地方震、区域震和远震。
震相自动拾取主要采用反向传播神经网络（Back Propagation Neural Network，简称BPNN），Wang和Teng（1995，1997）设计使用长短窗序列和滑动窗频率谱作为输入的神经网络算法，与STA/LTA相比，在拾取精度和低信噪比适应性上取得一定优势。Dai和Mocbeth（1995）对ANN在震相到时拾取方面的性能进行了测试。张范民等（1998）使用三分量地震数据的矢量模作为ANN输入，确定震相类型和拾取到时。Zhao和Takano（1999）使用BPNN从IRIS数据中拾取初至，用于地幔结构层析成像研究。Gentili和Michelini（2006）用具有高泛化能力的树形神经网络IU-ANT2拾取P和S震相。
1.7 总结与分析
STA/LTA、AIC、PAI-S/K、偏振分析屬于单特征方法，所使用的特征函数反映波形记录的某一项或几项特征，有相对明确的物理含义，都发展出了大量的衍生变种和组合算法。其中STA/LTA方法的突出优点是计算简单快速，不需要先验信息，只使用震相监测点之后少量的采样点（图1a），因此适用于对实时性和处理效率要求较高的地震监测、地震预警和大数据量处理。STA/LTA方法的缺点之一是低信噪比时，到时拾取精度会受到影响，因此通常配合其他方法进行精确拾取。而Lomax等（2012）基于多频带滤波和Bear和Kradolfer（1987）的方法提出FilterPicker方法，较好地解决了余震震相受其他地震尾波干扰的问题。STA/LTA方法的另一个缺点是需要根据不同应用调整参数和选择特征函数，Vassallo等（2012）对自动拾取方法的优化策略进行了研究，刘晗和张建中（2014）对各种特征函数和改进算法进行了分析。
AIC方法的突出优点是到时拾取精度高。其计算效率相对较低，并且假设背景噪声和震相信号是2个独立的过程，例如，当数据窗口内的震相数量大于一时便不符合这一假设，拾取结果就会受到影响。因此，AIC方祛的窗口选择十分重要，通常使用STA/LTA提供的先验信息进行精确拾取（马强等，2013；Bastin et al，2013；牟磊育，吴建平，2015）
PAI-S/K方法与STA/LTA一样，是常用的初至拾取方法，使用统计信息代替信噪比，在背景噪声接近高斯分布时有很好的拾取精度，计算简单高效，低信噪比环境下表现更稳定。缺点是在背景噪声不符合高斯分布时，拾取精度会受到影响。
偏振分析方法充分利用三分量信号的运动特征，可以得到多种偏振特征。震相的类型是前述方法无法或者很难直接获取的。由于与P波尾波重叠，出现转换震相，或者可能出现剪切波分裂，拾取S震相变得很困难。通过偏振分析方法能够完成震相类型识别、P波能量压制，进而提高S震相拾取率，这对于识别余震很有帮助。
余震与主震在构造上有密切联系，发震机理相似，空间上高度密集，大多数情况下波形相似度非常高，前述单特征方法都没有利用这些特点。由于噪声分布有随机性，互相关操作时会被有效的压制。互相关方法比对事件波形，是目前低信噪比拾取能力最强，拾取精度最高的方法。当然，事件波形形态受震源机制和传播路径等与具体事件相关的因素影响，这也导致波形互相关方法缺乏通用性。
机器学习算法基于多特征值和概率预测，最终的效果主要由具体算法和参数、特征选择和训练数据决定，震相拾取最常用的是经典机器学习算法和BPNN。经典机器学习算法需要人工提取和选择特征值。BPNN实际上比SVM、随机森林等经典机器学习方法更早出现，理论上通过训练可以自行提取特征值，但是层数增加时需要训练的权重系数很多，容易产生过拟合，受限于当时的硬件条件通常只采用三层网络，依然依赖人为选择特征，并且反复调参试验，实际使用需要一定的理论基础和经验。经典机器学习和全连接神经网络取得了优于传统算法的低信噪比适应性和拾取精度，通用性也要优于互相关方法，但是在实际监测工作中并没有得到广泛应用。
余震震相识别最大的难度来源j二震相信号受到主震和其他余震尾波的强烈干扰，从单台拾取震相到时的角度，FilterPicker和互相关方法的适用性和稳定性最好。但是随着偏移叠加技术的兴起，传统单特征方法所发展出的众多特征函数也可以用于低信噪比环境下的余震识别，本文将在第三章中进行介绍。
2 走时定位法研究现状
基于走时的常规方法被地震监测和预警所普遍采用，可以分为2部分：从连续观测波形中拾取震相；使用震相信息确定震源位置和发震时刻（图5）。一次余震有唯一的初始破裂位置和发震时刻，如果已知台站的空间位置，使用合适的理论模型可以预估震相到时。因此，利用多个台站拾取到的震相，通过走时定位方法，能够寻找到合适的震源位置和发震时刻，使得走时残差小于设定阈值，即可识别出余震。
走时定位方法通过最小化观测与理论走时的残差确定震源位置和发震时刻，其发展已经较为成熟，目前的难点在于深度定位（郑勇，谢祖军，2017），
常用的走时定位方法依据求解方式可分为传统几何定位方法、线性定位法和非线性方法，依据残差计算方式可分为绝对定位方法和相对定位方法。
几何定位法使用均匀介质走时方程和几何知识，在平面上作图确定震源点，是早期地震定位中的主要方法，常见的有方位角法、石川法、交切法。虽然精度有限，几何方法在余震识别中仍然是有效的约束手段。
线性定位方法最早由Geiger（1912）提出，其基本原理是将关于震源时空参数的非线性方程组作泰勒展开，简化为线性方程组，用最小二乘法求解，通过取初值迭代得到震源参数。由于理论速度模型与真实速度模型存在偏差，在Geiger方法基础上衍生出众多改进型方法。Douglas（1967）提出的JED方法和Dewey（1972）提出的JHD方法，认为残差由台站下方速度结构的差异引起，引入台站校正系数校正；Crosson（1976）提出的SSH方法，认为走时残差由速度结构偏差引起，同时进行地震定位和地下速度结构反演。
当2次地震发生在同一个震源区，并且震源区的范围远小于到某个台站的距离，相同震相到达这个台站所经过的路径绝大部分是相同的，也就是说速度模型引起的走时残差近似相等，因此可以假设2个震相间的走时差是由震源间的距离引起的，利用该思想的定位方法称之为相对定位方法。主事件方法以一个定位比较准确的地震作为参考事件，对周围的其他地震进行相对定位（Evernden，1969a，b； Spence，1980），由于被定位事件必须距离参考事件足够近，主事件方法的定位范围受到一定限制。Waldhauser（2000）提出双差定位方法（DD），将震源区距离较近的地震构成地震对，利用事件对的走时差约束相对位置，由于不需要主事件，适用范围比主事件更广。
地震定位是一个非线性反演问题，可以使用非线性最优化方法求解，国内也有很多研究，常见的有Powell方法（唐国兴，1979）、遗传算法（万永革，李鸿吉，1995；周民都等，1999；王文萍，王庆良，1999）、模拟退火法（高星等，2002；孙丽，2016）、网格搜索法（赵仲和，牟磊育，2005；宋维琪，杨晓东，2011）、单纯形法（赵珠等，1994）。非线性方法的优势是可以较好地避免结果陷人局部极小值，但是运算量通常要高于线性方法。
研究发震构造的空间展布是余震自动识别的重要应用，对余震的定位精度提出较高的要求，由于绝对定位会引入理论模型和真实地下结构在整个射线路径上的偏差，因此无法满足精度需求。余震序列在空间上高度密集，这一特征是采用主事件和DD等相对定位方法的有利条件。DD方法目前是震后重定位余震，确定发震构造最常用的方法（黄媛等，2008；房立华等，2011；苏金蓉等，2013：张广伟等，2014）。
3 偏移叠加法研究现状
基于走时定位识别余震只使用震相到时和类型，识别能力受制于单台震相拾取算法，如果震相引起的特征变化无法达到触发阈值，显然无法通过联合走时来定位和识别余震。一次余震会引发多个台站的波形记录变化，在不同台站记录中独立进行震相的识别，没有利用这种多台站震相信息间的一致性。偏移叠加方法不需要拾取到时或者识别震相类型，直接对齐各个台站的波形或者特征值，加权平均，能有效压制噪声，提高地震识别准确率和定位精度。由于石油、页岩气、水库监测等领域对诱发微震监测的投入，偏移叠加方法发展十分迅速，可以将其分为2类：一类是基于波动方程的波场反向传播聚焦方法，这类方法完全通过数值计算，因而计算量很大，也容易受强噪声和地下结构不均匀的影响（Mcmechan et al，1982，1985；Gajewski，Tessmer，2005）；另一类是直接叠加特征函数的震源扫描方法（Source Scanning Algorithm，简称SSA）及其衍生方法。
SSA由Kao和Shan（2004，2007）提出，在大致确定震源参数的情况下，搜索潜在震源位置和发震时间，叠加各个台站归一化的记录振幅绝对值，得到的亮度函数最大值点即认为是最佳解。Baker等（2005）用P波包络进行叠加，并引入基尔霍夫重建法（Kirchhoff reconstruction method）。Gharti等（2010）根据理论预测到时，截取数据并旋转至射线路径，最后同时使用P波和S波包络进行叠加。Grigoli等（2013，2014，2016）在SSA方法基础的上使用STA/LTA作为特征函数进行叠加，并引入主事件方法和震源台站修正项，进一步提高定位精度。Langet等（2014）则提出基于连续峰度函数的偏移叠加识别微震。
模板匹配技术（Matched Filter Technique，简称MFT）结合波形互相关和台阵叠加技术，将多通道的互相关波形偏移到模板事件的发震时刻，叠加取平均值得到平均互相关波形，最大值超过设定阈值就认为检测到一次相似地震（Gibbons，Ringdal，2006；Shelly et al，2007；Peng，Zhao，2009），模板的位置视为被匹配地震的位置，模板发震时刻偏移后确定发震时刻。MFT被广泛运用于完备目录中遗漏的早期余震（Peng，Zhao，2009；Lenglinee etal，2012：Meng et al，2013：Rastin，2013：Kato，2012，2014a，2014b； Wu et al，2014）和研究远场动态触发地震（Yukutake，2013； Meng，Peng，2014）0在MFT和SSA基础上，Zhang和Wen（2015a）提出匹配定位方法（Match And Locate，简称M&L;），通过搜索模板事件附近的三维空间网格，将各个台站的互相关波形对齐叠加，互相关系数最高的点即为被检测事件的空间位置。M&L;方法不仅进一步提高对低信噪比事件的检测能力，同时也给出精确的空间位置，被用于精确定位朝鲜核试验（Zhang，Wen，2015b）和监测火山喷发前的微震活动（Zhang，Wen，2015c）。
偏移叠加方法直接叠加特征函数，比传统走时定位方法有更高的地震识别能力，其函数既可以选择传统震相拾取方法所用特征函数，也可以是互相关波形，不同偏移叠加方法的识别能力和特点与所选叠加函数有关。余震震源机制和震中位置重复性好，这对基于互相关方法的MFT和M&L;十分有利，MFT是目前检测能力最强的余震自动识别方法之一，可以识别出10倍于正式目录的地震，常被应用于完备余震目录。M&L;采用相对定位方式，非常适合余震的高精度定位，是已知事件大致位置和时间的情况下，目前识别能力和定位精度最高的方法。但是二者对模板有非常高的依赖性，需要大量的先验信息，通用性较差，在识别余震时需要准备大量的模板。因此，选用其他单特征震相拾取方法所用的特征函数，更利于提升性能和通用性，但是检测能力会有所降低。余震在空间上密集分布，能將搜索区域限定在相对小的范围内，这对于计算量庞大的偏移叠加方法也是十分有利的条件。
4 讨论与展望
如前所述，余震自动识别仍然有许多技术难点亟待解决，震相自动识别技术是余震自动识别技术的基础，也是最大的难点所在，多台震相联合方法是增强余震识别能力的有效手段。本文以传统震相拾取方法和互相关方法为出发点，首先总结现有震相识别技术及其特征函数，然后对基于走时定位的识别方法进行综述，同时总结了近些年兴起的偏移叠加方法的发展，并分别分析了各项技术在余震自动识别中的优缺点和适用性。下面将从偏移叠加和震相识别方法2方面，分析存在的问题和原因，并简述对余震自动识别技术未来工作的设想。
4.1 优化偏移叠加方法
偏移叠加方法通过叠加特征函数，不需要拾取震相，能识别单台震相信号无法达到触发阈值的余震，识别能力明显强于基于走时定位的常规方法，各种衍生方法都可以用于余震的识别。常规走时定位方法的优点是计算简单快速，容易并行执行，不需要任何先验信息，适用于日常监测和地震预警等对时效要求比较高的应用。偏移叠加方法有非常高的识别能力，对潜在震源时间和空间进行搜索需要耗费大量的计算资源，虽然余震和人工诱发地震类似，震源的密集空间分布缩小了搜索范围，但是时效上仍然需要提升。因此，从实际应用的角度出发，偏移叠加方法的性能优化将成为未来余震识别的重点发展方向。
4.2 探索深度学习在余震识别领域的应用
波形特征函数是震相到时拾取和偏移叠加方法的基础，因此单台震相拾取方法仍然是提升余震自动识别能力的重点。特征值是对观测数据的抽象，传统方法选取人为设计的数种特征值，通常不需要先验信息，具有计算简单快速的特点，但是单一的特征值无法描述震相信号的所有变化，参数设置也很难适应复杂多变的余震波形。互相关方法利用波形的整体形态特征，压制噪声的能力非常强，震相拾取能力远超传统方法。波形形态中包含与具体事件相关的信息，如震源机制和空间位置，余震的震源机制和空间位置虽然多数非常接近，但也存在偏差较大的情况。因此，互相关方法通常需要准备大量模板，如Peng和Zhao（2009）在識别余震时使用3647个地震事件的波形作为模板，这导致其通用性和时效性较差，自相关方法和FAST方法的出现也没有从根本上解决这一问题。
笔者认为介于传统方法和互相关方法之间的综合特征方法是未来发展的方向，即从波形数据中尽可能多地抽象出与具体余震事件无关的特征，综合这些特征进行震相拾取，应当注意的是，这里的综合特征并不等同于多特征。多特征方法中最典型的属经典机器学习方法，这些方法依赖于人工从数据中提取出的数种特征值。BPNN等全连接神经网络方法理论上可以自动提取特征值，但是层数增加时需要训练的权重过多，实际应用中以3层网络为主，一般通过增加神经元数量来提升性能，与经典机器学习方法一样，比较依赖人为选择。
深度学习方法是近年兴起的ANN子集，优点是拥有更多隐含层，通过多层次的抽象比简单的增加神经元数量具有更强的描述能力，可以从原始数据中抽象出所需特征，这些特征也许我们从未注意到过，人工智能现今的优异成绩几乎都是由深度学习方法所取得。Perol等（2017）使用CNN方法识别和分类美国中部的诱发地震，在取得比匹配滤波方法更高识别能力的同时，计算耗时和存储使用量比自相关方法和FAST方法明显降低。卷积神经元网络（Convolutional Neural Network，简称CNN）是深度神经网络的一种，使用卷积核作为中介，识别数据中存在的固有局部结构，降低权重参数数量，可以有效避免过拟合，十分适合处理图像和声学数据。地震连续波形数据与声学数据有很多共同点，都是记录介质的震动信息，并且存在固有的局部结构，卷积神经网络方法识别震相可能会取得与处理声学数据相近的效果。
深度学习方法的运算量集中在模型训练阶段，使用模型进行预测的速度则非常快，模型的泛化可以提升适应性，因此具备同时兼顾运行效率和通用性的能力，在余震识别领域尚处于起步阶段。深度学习方法，属于黑盒算法，即无法从其训练模型中分析震相特征的抽象过程，以及如何组合特征识别震相，因此无法提升对地震信号的认识，但是优异的性能对于解决震相识别这类实际问题很有意义，未来将成为震相识别，乃至地震监测技术的重要发展方向。
如果深度学习能够准确拾取震相到时，那么在识别余震时，常规走时定位法所面临的低信噪比震相无法拾取和精确性不足问题就迎刃而解，而偏移叠加方法也能够更好地约束搜索范围，提高定位精度。当然，立足更长远的未来，深度学习技术也有可能替代偏移叠加，以更复杂和更抽象的方式实现多台震相联合，甚至进一步由余震识别这个小的领域向地震监测和地震数据处理发展。
本人关于地震定位、震相拾取的理论学习和算法实现都是在中国地震台网中心马延路博士的指导下完成，审稿专家对本文也提出了宝贵意见，在此一并感谢。
参考文献：
房立华，吴建平，张天中，等.2011.2011年云南盈江M55.8地震及其余震序列重定位[J].地震学报，33（2）：262-267.
高淑芳，李山有，武东坡，等.2008.一种改进的STA/LTA震相自动识别方法[J].世界地震工程，24（2）：37-41.
高星，王卫民，姚振兴.2002.用于地震定位的SAMS方法[J].地球物理学报，45（4）：525-532.
黄媛，吴建平，张天中，等，2008.汶川8.0级大地震及其余震序列重定位研究[J].中国科学：地球科学，38（10）：1242-1249.
孔令文.2014.神经网络的分类及其应用[J].电子技术与软件工程，（17）：25.
刘晗，张建中.2014.微震信号自动检测的STA/LTA算法及其改进分析[J].地球物理学进展，29（4）：1708-1714.
刘翰林，吴庆举.2017.地震自动识别及震相自动拾取方法研究进展[J].地球物理学进展，32（3）：1000-1007.
刘劲松，王赟，姚振兴.2013.微地震信号到时自动拾取方法[J].地球物理学报，56（5）：1660-1666.
刘希强，周彦文，曲均浩，等.2009.应用单台垂向记录进行区域地震事件实时检测和直达P波初动自动识别[J].地震学报，31（3）：260-271.
马强，金星，李山有，等.2013.用于地震预警的P波震相到时自动拾取[J].地球物理学报，56（7）；2313-2321.
牟磊育，吴建平.2015.利用单台触发检测台阵数据库事件[J].地球物理学进展，30（6）；2511-2516.
曲均浩，刘希强，吴丹彤，等.2012.基于神经网络的近震与远震识别[J].地震研究，30（3）：360-366.
宋金，蒋海昆.2009.序列衰减与余震激发研究进展及应用成果[J].地震地质，31（3）：559-571.
宋维琪，杨晓东.2011.解域约束下的微地震事件网格搜索法、遗传算法聯合反演[J].石油地球物理勘探，46（2）：259-266，160.
苏金蓉，郑钰，杨建思，等.2013.2013年4月20日四川芦山7.0级地震与余震精确定位及发震构造初探[J].地球物理学报，56（8）：2636-2644.
孙丽.2016.国家地震台网地震定位方法的改进[J].国际地震动态，（11）：12-16.
唐国兴.1979.用计算机确定地震参数的一个通用方法[dl.地震学报，1（2）：186-196.
万永革，李鸿吉.1995.遗传算法在确定震源位置中的应用[J].地震地磁观测与研究，16（6）：1-7.
王彩霞，白超英，王馨.2013.地震震相初至自动检测技术综述[J].地球物理学进展，28（5）：2363-2375.
王文萍，王庆良.1999.利用遗传算法和最小二乘联合反演共和地震位错参数[J].地震学报，21（3）；62-67.
徐锡伟，陈桂华，于贵华，等.2013a.芦山地震发震构造及其与汶川地震关系词论[J].地学前缘，20（3）：11-20.
徐锡伟，闻学泽，韩竹军，等.2013b.四川芦山7.0级强震：一次典型的盲逆断层型地震[J].科学通报，58（20）：1887-1893.
徐锡伟，闻学泽，叶建青，等.2008.汶川MS8.0地震地表破裂带及其发震构造[J].地震地质，30（3）；597-629.
余建华，李丹丹，韩国栋.2011.特征函数响应特性分析及STA/LTA方法的改进[J].国外电子测量技术，30（7）：17-19，23，27.
余娜，张晓清，袁伏全.2016.2016年青海门源MS6.4地震余震序列参数稳定性分析[J].地震研究，39（S1）；62-68，133.
张范民，李清河，张元生，等.1998.利用人工神经网络理论对地震信号及地震震相进行识别[J].西北地震学报，20（4）；43-49.
张广伟，雷建设，梁姗姗，等.2014.2014年8月3日云南鲁甸MS6.5级地震序列重定位与震源机制研究[J].地球物理学报，57（9）：3018-3027.
张苏平，孙艳萍，陈文凯.2013.余震信息在岷县漳县6.6级地震震后应急救援重灾区快速判定时作用探讨[J].地震工程学报，35（3）：465-470.
赵大鹏，刘希强，李红，等.2012.峰度和AIC方法在区域地震事件和直达P波初动自动识别方面的应用[J].地震研究，35（2）；220-225，295.
赵仲和，牟磊育.2005.单台地震自动定位网格搜索法及其MATLAB试验[J]地震地磁观测与研究，26（4）：1-12.
赵珠，丁志峰，易桂喜，等.1994.西藏地震定位——一种使用单纯形优化的非线性方法[J].地震学报，16（2）；212-219.
郑勇，谢祖军.2017.地震震源深度定位研究的现状与展望[J].地震研究，40（2）：167-175，333.
郑韵，姜立新，杨天青，等.2015：利用余震频度分布进行宏观震中快速判定[J].地震，35（2）：121-132.
周民都，张元生，张树勋.1999.遗传算法在地震定位中的应用[J].西北地震学报，21（2）：167-171.
周彦文，刘希强.2007.初至震相自动识别方法研究与发展趋势[J].华北地震科学，25（4）：18-22.
庄建仓，马丽.2000.主震和余震——从大森公式到ETAS模型[J].国际地震动态，（5）：12-18.
Allen R V.1978.Automatic earthquake recognition and timing from singletraces[J].Bulletin of the Seismological Society of America，68（5）：1521-1532.
Allen R.1982.Automatic phase pickers；Their present use and futureprospects[J].Bulletin of the Seismological Society of America，72（6）：225-242.
Amoroso O，Maercklin N，Zollo A.2012.S-wave identification by polari-zation filtering and waveform coherence analyses[J].Bulletin of theSeismological Society of America，102（2）：854-861.
Baer M，Kradolfer U.1987.An automatic phase picker for local andteleseismic events[J].Bulletin of the Seismological Society of Amer-ica，77（4）：1437-1445.
Bai C Y，Kennett B L N.2000.Automatic phase-detection and identifi-cation by full use of a single three-component broadband seismo-gram[J].Bulletin of the Seismological Society of America，90（1）：187-198.
Baillard C，Crawford W C，Ballu V，et al.2014.An automatic kurtosis-based P-and S-phase picker designed for local seismic networks[J].Bulletin of the Seismological Society of America，104（1）：394-409.
Baker T，Granat R，Clayton R W.2005.Real-time earthquake locationusing Kirchhoff reconstruction[J].Bulletin of the Seismological So-ciety of America，95（2）：699-707.
Brown J R，Beroza G C，Shelly D R.2008.An autocorrelation method todetect low frequency earthquakes within tremor[J].Geophysical Re-search Letters，35（16）：1-5.
Crosson R S.1976.Crustal structure modeling of earthquake data，1，si-multameous least squares estimation of hypocenter and velocity pa-rameters[J].Geophys Res，81（17）：3036-3046.
Dai H，Macbeth C.1995.Automatic picking of seismic arrivals in localearthquake data using an artificial neural network[J].GeophysicalJournal International，120（3）：758-774.
Dewey James W.1972.Seismicity and tectonics of western venezuela[J].Bulletin of the Seismological Society of America，62（6）：1711-1751.
Diallo M S，Kulesh M，Holschneider M，et al.2006.Instantaneous polari-zation attributes based on an adaptive approximate covariance method[J].Geophysics，71（5）：V99.
Douglas A.1967.Joint epicentre determination[J].Nature，215（5096）：47-48.
Earle P S，Shearer P M.1994.Characterization of global seismograms u-sing an automatic-picking algorithm[J].Bulletin of the Seismologi-cal Society of America，84（2）：366-376.
Esposito A M，D Auria L，Giudicepietro F，et al.2013.Automatic recogni-tion of landslides based on neural network analysis of seismic sig-nals：an application to the monitoring of stromboli volcano（southernItaly）[J].Pure and Applied Geophysics，170（11）：1821-1832.
Evernden J F.1969a.Identification of earthquakes and explosions by useof teleseismic data[J].Journal of Geophysical Research，74（15）：3828-3856.
Evernden J F.1969b.Precision of epicenters obtained by small numbersof world-wide stations[J].Bulletin of the Seismological Society ofAmerica，59（3）：1365-1398.
Flinn E A.1965.Signal analysis using rectilinearity and direction of parti-cle motion[J].Proceedings of the IEEE，53（12）：1874-1876.
Gajewski D，Tessmer E.2005.Reverse modelling for seismic event charac-terization[J].Geophysical Journal International，163（1）：276-284.
Geiger L.1912.Probability method for the determination of earthquakeepicenters from arrival time only[J].Bull St Louis Univ，8；60-71.
Gentili S，Michelini A.2006.Automatic picking of P and S phases using aneural tree[J].Journal of Seismology，10（1）：39-63.
Gharti H N，Oye V，roth M，Kiihn D.2010.Automated microearthquakelocation using envelope stacking and robust global optimization[J].Geophysics，75（4）：27-46.
Gibbons S J，Ringdal F.2006.The detection of low magnitude seismic e-vents using array-based waveform correlation[J].GeophysicalJournal International，165（1）：149-166.
Giudicepietro F，Esposito A M，Ricciolino P.2017.Fast discrimination oflocal earthquakes using a neural approach[J].Seismological Re-search Letters，88（4）：1089-1096.
Grigoli F，Cesca S，Amoroso O，et al.2014.Automated seismic event loca-tion by waveform coherence analysis[J].Geophysical Journal Inter-national，196（3）：1742-1753.
Grigoli f，Cesca S，Krieger L，et al.2016.Automated microseismic eventlocation using Master-Event Waveform Stacking[J].Scientific Re-ports，6：1-13.
Grigoli F，Cesca S，Vassallo M，et al.2013.Automated Seismic Event Lo-cation by Travel-Time Stacking；An Application to Mining InducedSeismicity[J].Seismological Research Letters，84（4）：666-677.
Hibert C，Mangeney A，Grandjean G，et al.2014.Automated identifica-tion，location，and volume estimation of rockfalls at Piton de la Four-naise volcano[J].Journal of Geophysical Research：Earth Surface，119（5）：1082-1105.
Horstmann T，Harrington R M，Cochran E S.2013.Semiautomated tremordetection using a combined cross-correlation and neural networkapproach[J].Journal of Geophysical Research；Solid Earth，118（9）：4827-4846.
Jackson G M，Mason I M，Greenhalgh S A.1991.Principal componenttransforms of triaxial recordings by singular value decomposition[J].GEOPHYSICS，56（4）：528-533.
Jurkevics A.1988.Polarization analysis of three-component array data[J].Bulletin of the Seismological Society of America，78（5）：1725-1743.
Kao H，Shan S J.2004.The Source-Scanning Algorithm；Mapping thedistribution of seismic sources in time and space[J].GeophysicalJournal International，157（2）：589-594.
Kao H，Shan S J.2007.Rapid identification of earthquake rupture planeusing Source-Scanning Algorithm[J].Geophysical Journal Interna-tional，168（3）：1011-1020.
Kato A，Nakagawa S.2014a Multiple slow-slip events during a fore-shock sequence of the 2014 Iquique，Chile Mw 8.1 earthquake[J].Geophysical Research Letters，41（15）：5420-5427.
Kato A，Obara K，Igarashi T，et al.2012.Propagation of slow slip leadingup to the 2011 MW9.0 Tohoku-Oki earthquake[J].Science，335（6069）：705-708.
Kato A，Obara K.2014b.Step-like migration of early aftershocks follow-ing the 2007 Mw 6.7 Noto-Hanto earthquake，Japan[J].Geophysi-cal Research Letters，41（11）：3864-3869.
Kortstrom J，Uski M，Tiira T.2016.Automatic classification of seismic e-vents within a regional seismograph network[J].Computers and Ge-osciences，87：22-30.
Kiiperkoch L，Meier T，Lee]，et al.2010.Automated determination of P-phase arrival times at regional and local distances using higher orderstatistics[J].Geophysical Journal International，181：1159-1170.
Kurzon I，Vernon F L，Rosenberger A，et al.2014.Real-time automaticdetectors of P and S waves using singular value decomposition[J].Bulletin of the Seismological Society of America，104（4）：1696-1708.
Langet N，Maggi A，Michelini A，et al.2014.Continuous kurtosis-basedmigration for seismic event detection and location，with application toPiton de la Fournaise volcano，La Reunion[J].Bulletin of the Seis-mological Society of America，104（1）：229-246.
Lenglinee O，Enescu B，Peng Z，et al.2012.Decay and expansion of theearly aftershock activity following the 2011，MW9.0 Tohoku earth-quake[J].Geophysical Research Letters，39（18）：L18309.
Leonard M，Kennett B L N.1999.Multi-component autoregressive tech-niques for the analysis of seismograms[J].Physics of the Earth andPlanetary Interiors，113（1-4）：247-263.
Lomax A，Satriano C，Vassallo M.2012.Automatic picker developmentsand optimization：FilterPicker-a robust，broadband picker for real-time seismic monitoring and earthquake early warning[J].Seismo-logical Research Letters，83（3）：531-540.
Maeda N.1985.A method for reading and checking phase time in auto-processing system of seismic wave data[J].Zisin，38（3）：365-379.
Maggi A，Ferrazzini V，Hibert C，et al.2017.Implementation of a multista-tion approach for automated event classification at piton de la Four-naise volcano[J].Seismological Research Letters，88（3）：878-891.
Mcmechan G A，Luetgert J H，Mooney W D.1985.Imaging of earthquakesources in Long Valley Caldera，alifornia，1983 [J].Bulletin of theSeismological Society of America，75（4）：1005-1020.
Mcmechan G A.1982.Determination of source parameters by wavefieldextrapolation[J].Geophysical Journal of the Royal Astronomical So-ciety，71（3）：613-628.
摘要：基于余震识别自动化对地震监测、灾后应急、确定发震构造、分析预报和科学研究的重要意义，总结和分析了余震自动识别技术的现状和未来发展。首先，总结单台震相到时自动拾取方法及其特征函数，包括单特征方法、互相关方法和机器学习方法。其次，分别总结基于震相走时的常规识别方法和基于特征函数的偏移叠加方法。最后，对余震识别技术的未来发展进行分析和展望，认为偏移叠加方法的性能改进将成为未来余震识别的重点发展方向，深度学习方法作为一种新技术，将在包括震相拾取和余震识别在内的地震数据处理领域发挥重要作用。
关键词：余震；自动识别；震相自动拾取；走时定位；偏移叠加
中图分类号：P315.61 文献标识码：A 文章编号：1000-0666（2018）01-0001-13
0 引言
破坏性地震发生之后，随之而来的是大量的余震，从连续波形中识别余震是开展余震震级测定、震源机制解反演等后续监测工作的基础。余震信息是评估灾情和指导救援力量部署（张苏平等，2013；郑韵，2015）、确定发震构造和揭示隐伏断层（黄媛等，2008；徐锡伟等，2008，2013a，b）、判定震后趋势（庄建仓，马丽，2000；宋金，蒋海昆，2009；余娜等，2016）的重要依据，也是震后地下速度结构变化、余震迁移等科研工作的基础数据。人工识别余震是一项十分耗时且枯燥的工作，结果中不可避免地会引人主观因素，随着观测台网的建设，数据的产出量越来越大，完全由人工处理变得越来越不现实。大震之后，识别高频度的余震需要持续投入大量的时间和人力，这让有限的监测资源更加捉襟见肘，数据产出时效也受到限制。
余震识别的自动化有重要的应用和科研价值，但是现有技术仍然有很多问题亟待解决和完善，地震监测领域也没有专门的余震自动识别系统。为了提高对余震自动识别技术的关注和投入，2017年5月中国地震局地球物理研究所和阿里云联合举办余震捕捉AI大赛，吸引来自国内外的近千支队伍报名参赛，一时间余震自动识别成为地震领域的热门话题和技术热点。
余震发生后，其能量以地震波的形式向外传播，当震相到达各个观测台站时，会引起波形记录的变化。自动识别余震时，首先需要从连续观测波形中提取特征，并自动拾取震相到时，然后联合多个台站的数据确定余震信息。对于震相清晰的独立地震，可以使用单台波形进行识别，并初步确定震源信息。但是余震在时间上高度密集，并且震级普遍较小，记录到的信号受主震尾波和其他余震的强烈干扰，多台联合是提高余震识别可靠性的必要技术。同时，余震在构造上与主震有密切联系，在空间上沿发震构造密集展布，发震机理有较高的相似性，这对于震相识别和震相联合是可以利用的有利条件。传统震相拾取方法所使用的特征函数在偏移叠加和人工智能等新技术中也有非常广泛的应用，因此本文从单特征方法出发，总结几种重要的震相自动拾取方法及其特征函数。然后介绍基于震相走时定位的常规方法的发展，并分析其在余震自动识别中的适用性，最后对基于特征函数的偏移叠加方法进行综述，并在此基础上展望余震自动识别技术的未来发展方向，以期对开展相关研究工作有所帮助。震相特征的提取是走时定位和偏移叠加方法的基础，为便于理解，对STA/LTA、VAR-AIC、峰度、极性分析等单台特征提取算法进行了简单实现，相关示例代码发布于github版本库https：//github.com/dualli/trigger。
1 震相自动拾取方法及其特征函数
当地震震相到达观测台站时，会改变记录信号的能量、频谱、极性和整体形态等特征，这些变化可以通过不同的特征函数提取，当特征函数达到设定阈值时拾取相应的震相到时。多个台站的特征函数和震相信息成为后续地震识别和定位的基本数据。目前，震相自动拾取方法种类繁多，前人对这些方法按照时间域、频率域、时频域、综合方法等进行过分类综述（Withers el al，1998；周彦文，刘希强，2007； Nippress el al，2010；王彩霞等，2013：Ross，Ben-Zion，2014；刘翰林，吴庆举，2017）。本文选取常用的长短窗、赤池信息准则、峰态和偏态、偏振分析等单特征方法和基于整体形态的互相关方法、基于多特征的机器学习方法，分别阐述基本原理和特征函数演化，并比较分析其特点。
1.1 长短窗均值比方法
长短窗均值比方法（Short-Term-Average/Long-Term-Average，简称STA/LTA）基于震相到达时特征值的平均值在短窗口比长窗口变化快的特点拾取震相。STA/LTA方法的发展历史久远，应用广泛，形成了反映波形不同变化特点的众多特征函数。Stevenson（1976）首次提出长短窗方法时，使用振幅的绝对值|Y（i）|作为特征函数，其中Y（i）为第i个采样点的振幅值。但是使用|Y（i）|或者Y（i）2作为特征函数时，长短窗比值只对振幅的变化敏感，因此Allen（1978，1982）引入了特征函数CF=Y（i）2+K（Y（i）-Y（i-1））2，其中K为根据采样率和台站噪声特征设置的权重常数，该特征函数的引入使得长短窗均值比同时随振幅和频率变化而变化（图1）。Bear和Kradolfer（1987）使用基于包络函数平方的特征函数，同时引入动态阈值，既避免了瞬变噪声引起的误触发，也适用于拾取地方震到远震的不同震相，同时对低信噪比微震信号也变得更敏感。为了更好地拾取各种微弱的次生震相，Earle和Shearer（1994）定义了基于Hilbert变换的包络函数作为特征函数。高淑芳等（2008）也提出了反映振幅和频率变化的新特征函数。余建华等（2011）引入加权系数法减少震相漏捡率。
1.2 赤池信息准则方法
赤池信息准则方法（Akaike Information Criteri-on，简称AIC）沿每個采样点将信号分为2段，当分割方法与2个独立过程对应时，对应的AIC值最小，即全局最小值对应震相到时。AIC方法也有非常多的变种，以下介绍几种较为常见的AIC函数。AR-AIC（Sleeman，Van Eek，1999；Leonard，Kennett，1999）使用自回归方法计算AIC值，即分割点对应2个自回归过程，自回归阶数需要反复试验得到。Maeda（1985）提出的基于二阶统计量的VAR-AIC方法，可以直接从地震信号计算AIC值（图2）.0刘希强等（2009）提出TOC-AIC方法，使用三阶累积量代替方差来压制高斯有色噪声，在低信噪比情况下有更高的精度。赵大鹏等（2012）提出峰度一赤池信息准则方法（Kurtosis-AIC），使用峰度作为特征函数替代VAR一AIC中的方差，有效减少了误报率和漏报率。
1.3 偏态/峰态方法
偏态/峰态方法（P Arrival Identification-Skewness/Kurtosis，简称PAI-S/K）由Saragiotis等（2002）提出，当震相到达观测台站时，记录数据的分布由具有高斯性的随机噪声向非高斯性转变。峰态（Kurtosis）和偏态（Skewness）是反映数据分布形态的高级统计量（Higher-Order Statis-tics，HOS），数据服从高斯分布时峰态等于3，分布变宽时小于3，变窄时大于3，同理偏态在数据对称分布时等于0，左偏时小于0，右偏时大于0（图3）。刘劲松等（2013）引入了新的到时判定算法，改进了PAI-SIK方法的拾取精度。Baillard等（2014）使用滑动窗峰态函数和极化分析识别P、S波。峰态和偏态值方法有较好的低信噪比拾取能力和精度，属于常用的震相拾取方法之一（Panagiotakis et al，2008；Kiiperkoch et al，2010；Nippress et al，2010；Hibert et al，2014；Langet etal，2014；Ross，2016）。
1.4 偏振分析方法
常用的2种偏振分析方法是协方差矩阵方法和奇异值分解方法，在多数情况下二者效果接近。偏振分析的思想最早由Flinn（1965）提出，他使用协方差矩阵的本征值分解来识别S波偏振，进而约束震源机制解。Montalbetti和Kanasewieh（1970）首次详细阐述如何用对称本征问题来进行偏振分析。分解有限长度三分量数据的协方差矩阵，得到本征值和本征向量矩阵，可以用来计算质点运动的多种偏振特征，沿着时间轴滑动窗口，重复计算就可以得到偏振特性随时间的变化（图4），进而完成P或S波能量压制、波形分解等操作，是提升S震相拾取率的常用方法（Jurkevies，1988；Robertset al，1989；Bai，Kennett，2000：Diallo et al，2006：Amoroso et al，2012；Ross et al，2014，2016）。Jackson等（1991）首次引入奇异值分解方法（Singular Value Decomposition，简称SVD）进行偏振分析，一步即可完成数据的分解操作。在沿时间轴滑动时，传统的偏振分析方法要使用窗口内数据重复求解协方差矩阵，Rosenberger（2010）引入迭代奇异值分解法，通过追加数据矢量来更新解，提升了奇异值分解方法的计算效率，使其适用于实时数据处理。Kurzon等（2014）基于迭代奇异值分解方法将地震波形实时分解成P分量和S分量，然后使用视出射角、奇异值、信噪比拾取器拾取P、S震相，并且尝试多种指示变量来增强P、S波的拾取。
1.5 互相关方法
互相关系数是描述2段信号相似程度的量，地震波形的相似性主要由震源机制和经过的路径决定，因此互相关方法可以将事件波形作为模板，从连续波形中寻找相似地震（Poupinet et al，1984；Rubin，2002；Schaff，Richards，2004，2011），该方法对于拾取低信噪比震相和震颤等特殊信号非常有用（Shelly et al，2007；Peng，Gomberg，2010），能检测到比正式目录震级低1级的地震事件（Schaff，2010；Schaff，Waldhauser，2010；Gib-blon，Ringdal，2006）。Yang等（2009）提出聯合三分量的滑动窗互相关方法（Sliding-WindowCross-correlation，简称SCC），比单独的三分量互相关方法拾取到更多低信噪比震相。
互相关方法必须先制作好模板，并且灵敏度高度依赖模板和待检测信号的相似性，Brown等（2008）提出自相关方法来，能够从连续波形中提取模板，但是计算十分耗时。Yoon等（2015）基于数据挖掘算法提出（Fingerprint And SimilarityThresholding，简称FAST），用关键特征代替波形，速度是自相关方法的140倍，能够在2h内处理单台一周的连续波形数据。FAST方法有效地提升了模板的制作速度，但是运算量依然会随数据长度的增加而快速增长，处理长达数月甚至数年的余震活动数据时，运算效率依然是需要解决的问题。Shelly等（2007）为了拾取更多的低频地震，尝试降低时间精度来延拓互相关系数的时间范围，提高对空间位置变化的适应性。
1.6 机器学习方法
机器学习（Machine Learning，简称ML）是一类从数据中寻找潜在模式，进而提升性能的庞大算法集合的总称，大致分为经典机器学习算法和人工神经网络算法（Artificial Neural Networks，简称ANN），这2种算法都包括样本数据训练模型和模型预测新数据的2个阶段。ANN是ML的子集，用简单数学函数构建神经元网络，模拟高度复杂的非线性关系，孔令文（2014）按照网络结构进行了分类。ML在地震数据处理的很多领域均有应用，但是主要集中在处理分类问题。Rawles和Thurber（2015）在微震定位中使用邻域法测定潜在震相的类别和到时。Kortstrom等（2016）使用支持向量机方法，基于多频带不同震相段的特征函数值，从监测系统自动识别的地震中分离出爆破和误触发事件。Maggi等（2017）使用随机森林方法实现1套自动地震事件分类器，使用多个台站的波形和频谱特征组合，对留尼旺火山岛的地震事件进行分类。曲均浩等（2012）使用反向传播神经网络分类近震和远震。Esposito等（2013）使用多层感知机分类意大利南部斯特龙博利火山的滑坡和其他事件。Horstmann等（2013）使用自组织映射神经网络检测和分类地震。Giudicepietro等（2017）使用多层感知机分类地方震、区域震和远震。
震相自动拾取主要采用反向传播神经网络（Back Propagation Neural Network，简称BPNN），Wang和Teng（1995，1997）设计使用长短窗序列和滑动窗频率谱作为输入的神经网络算法，与STA/LTA相比，在拾取精度和低信噪比适应性上取得一定优势。Dai和Mocbeth（1995）对ANN在震相到时拾取方面的性能进行了测试。张范民等（1998）使用三分量地震数据的矢量模作为ANN输入，确定震相类型和拾取到时。Zhao和Takano（1999）使用BPNN从IRIS数据中拾取初至，用于地幔结构层析成像研究。Gentili和Michelini（2006）用具有高泛化能力的树形神经网络IU-ANT2拾取P和S震相。
1.7 总结与分析
STA/LTA、AIC、PAI-S/K、偏振分析屬于单特征方法，所使用的特征函数反映波形记录的某一项或几项特征，有相对明确的物理含义，都发展出了大量的衍生变种和组合算法。其中STA/LTA方法的突出优点是计算简单快速，不需要先验信息，只使用震相监测点之后少量的采样点（图1a），因此适用于对实时性和处理效率要求较高的地震监测、地震预警和大数据量处理。STA/LTA方法的缺点之一是低信噪比时，到时拾取精度会受到影响，因此通常配合其他方法进行精确拾取。而Lomax等（2012）基于多频带滤波和Bear和Kradolfer（1987）的方法提出FilterPicker方法，较好地解决了余震震相受其他地震尾波干扰的问题。STA/LTA方法的另一个缺点是需要根据不同应用调整参数和选择特征函数，Vassallo等（2012）对自动拾取方法的优化策略进行了研究，刘晗和张建中（2014）对各种特征函数和改进算法进行了分析。
AIC方法的突出优点是到时拾取精度高。其计算效率相对较低，并且假设背景噪声和震相信号是2个独立的过程，例如，当数据窗口内的震相数量大于一时便不符合这一假设，拾取结果就会受到影响。因此，AIC方祛的窗口选择十分重要，通常使用STA/LTA提供的先验信息进行精确拾取（马强等，2013；Bastin et al，2013；牟磊育，吴建平，2015）
PAI-S/K方法与STA/LTA一样，是常用的初至拾取方法，使用统计信息代替信噪比，在背景噪声接近高斯分布时有很好的拾取精度，计算简单高效，低信噪比环境下表现更稳定。缺点是在背景噪声不符合高斯分布时，拾取精度会受到影响。
偏振分析方法充分利用三分量信号的运动特征，可以得到多种偏振特征。震相的类型是前述方法无法或者很难直接获取的。由于与P波尾波重叠，出现转换震相，或者可能出现剪切波分裂，拾取S震相变得很困难。通过偏振分析方法能够完成震相类型识别、P波能量压制，进而提高S震相拾取率，这对于识别余震很有帮助。
余震与主震在构造上有密切联系，发震机理相似，空间上高度密集，大多数情况下波形相似度非常高，前述单特征方法都没有利用这些特点。由于噪声分布有随机性，互相关操作时会被有效的压制。互相关方法比对事件波形，是目前低信噪比拾取能力最强，拾取精度最高的方法。当然，事件波形形态受震源机制和传播路径等与具体事件相关的因素影响，这也导致波形互相关方法缺乏通用性。
机器学习算法基于多特征值和概率预测，最终的效果主要由具体算法和参数、特征选择和训练数据决定，震相拾取最常用的是经典机器学习算法和BPNN。经典机器学习算法需要人工提取和选择特征值。BPNN实际上比SVM、随机森林等经典机器学习方法更早出现，理论上通过训练可以自行提取特征值，但是层数增加时需要训练的权重系数很多，容易产生过拟合，受限于当时的硬件条件通常只采用三层网络，依然依赖人为选择特征，并且反复调参试验，实际使用需要一定的理论基础和经验。经典机器学习和全连接神经网络取得了优于传统算法的低信噪比适应性和拾取精度，通用性也要优于互相关方法，但是在实际监测工作中并没有得到广泛应用。
余震震相识别最大的难度来源j二震相信号受到主震和其他余震尾波的强烈干扰，从单台拾取震相到时的角度，FilterPicker和互相关方法的适用性和稳定性最好。但是随着偏移叠加技术的兴起，传统单特征方法所发展出的众多特征函数也可以用于低信噪比环境下的余震识别，本文将在第三章中进行介绍。
2 走时定位法研究现状
基于走时的常规方法被地震监测和预警所普遍采用，可以分为2部分：从连续观测波形中拾取震相；使用震相信息确定震源位置和发震时刻（图5）。一次余震有唯一的初始破裂位置和发震时刻，如果已知台站的空间位置，使用合适的理论模型可以预估震相到时。因此，利用多个台站拾取到的震相，通过走时定位方法，能够寻找到合适的震源位置和发震时刻，使得走时残差小于设定阈值，即可识别出余震。
走时定位方法通过最小化观测与理论走时的残差确定震源位置和发震时刻，其发展已经较为成熟，目前的难点在于深度定位（郑勇，谢祖军，2017），
常用的走时定位方法依据求解方式可分为传统几何定位方法、线性定位法和非线性方法，依据残差计算方式可分为绝对定位方法和相对定位方法。
几何定位法使用均匀介质走时方程和几何知识，在平面上作图确定震源点，是早期地震定位中的主要方法，常见的有方位角法、石川法、交切法。虽然精度有限，几何方法在余震识别中仍然是有效的约束手段。
线性定位方法最早由Geiger（1912）提出，其基本原理是将关于震源时空参数的非线性方程组作泰勒展开，简化为线性方程组，用最小二乘法求解，通过取初值迭代得到震源参数。由于理论速度模型与真实速度模型存在偏差，在Geiger方法基础上衍生出众多改进型方法。Douglas（1967）提出的JED方法和Dewey（1972）提出的JHD方法，认为残差由台站下方速度结构的差异引起，引入台站校正系数校正；Crosson（1976）提出的SSH方法，认为走时残差由速度结构偏差引起，同时进行地震定位和地下速度结构反演。
当2次地震发生在同一个震源区，并且震源区的范围远小于到某个台站的距离，相同震相到达这个台站所经过的路径绝大部分是相同的，也就是说速度模型引起的走时残差近似相等，因此可以假设2个震相间的走时差是由震源间的距离引起的，利用该思想的定位方法称之为相对定位方法。主事件方法以一个定位比较准确的地震作为参考事件，对周围的其他地震进行相对定位（Evernden，1969a，b； Spence，1980），由于被定位事件必须距离参考事件足够近，主事件方法的定位范围受到一定限制。Waldhauser（2000）提出双差定位方法（DD），将震源区距离较近的地震构成地震对，利用事件对的走时差约束相对位置，由于不需要主事件，适用范围比主事件更广。
地震定位是一个非线性反演问题，可以使用非线性最优化方法求解，国内也有很多研究，常见的有Powell方法（唐国兴，1979）、遗传算法（万永革，李鸿吉，1995；周民都等，1999；王文萍，王庆良，1999）、模拟退火法（高星等，2002；孙丽，2016）、网格搜索法（赵仲和，牟磊育，2005；宋维琪，杨晓东，2011）、单纯形法（赵珠等，1994）。非线性方法的优势是可以较好地避免结果陷人局部极小值，但是运算量通常要高于线性方法。
研究发震构造的空间展布是余震自动识别的重要应用，对余震的定位精度提出较高的要求，由于绝对定位会引入理论模型和真实地下结构在整个射线路径上的偏差，因此无法满足精度需求。余震序列在空间上高度密集，这一特征是采用主事件和DD等相对定位方法的有利条件。DD方法目前是震后重定位余震，确定发震构造最常用的方法（黄媛等，2008；房立华等，2011；苏金蓉等，2013：张广伟等，2014）。
3 偏移叠加法研究现状
基于走时定位识别余震只使用震相到时和类型，识别能力受制于单台震相拾取算法，如果震相引起的特征变化无法达到触发阈值，显然无法通过联合走时来定位和识别余震。一次余震会引发多个台站的波形记录变化，在不同台站记录中独立进行震相的识别，没有利用这种多台站震相信息间的一致性。偏移叠加方法不需要拾取到时或者识别震相类型，直接对齐各个台站的波形或者特征值，加权平均，能有效压制噪声，提高地震识别准确率和定位精度。由于石油、页岩气、水库监测等领域对诱发微震监测的投入，偏移叠加方法发展十分迅速，可以将其分为2类：一类是基于波动方程的波场反向传播聚焦方法，这类方法完全通过数值计算，因而计算量很大，也容易受强噪声和地下结构不均匀的影响（Mcmechan et al，1982，1985；Gajewski，Tessmer，2005）；另一类是直接叠加特征函数的震源扫描方法（Source Scanning Algorithm，简称SSA）及其衍生方法。
SSA由Kao和Shan（2004，2007）提出，在大致确定震源参数的情况下，搜索潜在震源位置和发震时间，叠加各个台站归一化的记录振幅绝对值，得到的亮度函数最大值点即认为是最佳解。Baker等（2005）用P波包络进行叠加，并引入基尔霍夫重建法（Kirchhoff reconstruction method）。Gharti等（2010）根据理论预测到时，截取数据并旋转至射线路径，最后同时使用P波和S波包络进行叠加。Grigoli等（2013，2014，2016）在SSA方法基础的上使用STA/LTA作为特征函数进行叠加，并引入主事件方法和震源台站修正项，进一步提高定位精度。Langet等（2014）则提出基于连续峰度函数的偏移叠加识别微震。
模板匹配技术（Matched Filter Technique，简称MFT）结合波形互相关和台阵叠加技术，将多通道的互相关波形偏移到模板事件的发震时刻，叠加取平均值得到平均互相关波形，最大值超过设定阈值就认为检测到一次相似地震（Gibbons，Ringdal，2006；Shelly et al，2007；Peng，Zhao，2009），模板的位置视为被匹配地震的位置，模板发震时刻偏移后确定发震时刻。MFT被广泛运用于完备目录中遗漏的早期余震（Peng，Zhao，2009；Lenglinee etal，2012：Meng et al，2013：Rastin，2013：Kato，2012，2014a，2014b； Wu et al，2014）和研究远场动态触发地震（Yukutake，2013； Meng，Peng，2014）0在MFT和SSA基础上，Zhang和Wen（2015a）提出匹配定位方法（Match And Locate，简称M&L;），通过搜索模板事件附近的三维空间网格，将各个台站的互相关波形对齐叠加，互相关系数最高的点即为被检测事件的空间位置。M&L;方法不仅进一步提高对低信噪比事件的检测能力，同时也给出精确的空间位置，被用于精确定位朝鲜核试验（Zhang，Wen，2015b）和监测火山喷发前的微震活动（Zhang，Wen，2015c）。
偏移叠加方法直接叠加特征函数，比传统走时定位方法有更高的地震识别能力，其函数既可以选择传统震相拾取方法所用特征函数，也可以是互相关波形，不同偏移叠加方法的识别能力和特点与所选叠加函数有关。余震震源机制和震中位置重复性好，这对基于互相关方法的MFT和M&L;十分有利，MFT是目前检测能力最强的余震自动识别方法之一，可以识别出10倍于正式目录的地震，常被应用于完备余震目录。M&L;采用相对定位方式，非常适合余震的高精度定位，是已知事件大致位置和时间的情况下，目前识别能力和定位精度最高的方法。但是二者对模板有非常高的依赖性，需要大量的先验信息，通用性较差，在识别余震时需要准备大量的模板。因此，选用其他单特征震相拾取方法所用的特征函数，更利于提升性能和通用性，但是检测能力会有所降低。余震在空间上密集分布，能將搜索区域限定在相对小的范围内，这对于计算量庞大的偏移叠加方法也是十分有利的条件。
4 讨论与展望
如前所述，余震自动识别仍然有许多技术难点亟待解决，震相自动识别技术是余震自动识别技术的基础，也是最大的难点所在，多台震相联合方法是增强余震识别能力的有效手段。本文以传统震相拾取方法和互相关方法为出发点，首先总结现有震相识别技术及其特征函数，然后对基于走时定位的识别方法进行综述，同时总结了近些年兴起的偏移叠加方法的发展，并分别分析了各项技术在余震自动识别中的优缺点和适用性。下面将从偏移叠加和震相识别方法2方面，分析存在的问题和原因，并简述对余震自动识别技术未来工作的设想。
4.1 优化偏移叠加方法
偏移叠加方法通过叠加特征函数，不需要拾取震相，能识别单台震相信号无法达到触发阈值的余震，识别能力明显强于基于走时定位的常规方法，各种衍生方法都可以用于余震的识别。常规走时定位方法的优点是计算简单快速，容易并行执行，不需要任何先验信息，适用于日常监测和地震预警等对时效要求比较高的应用。偏移叠加方法有非常高的识别能力，对潜在震源时间和空间进行搜索需要耗费大量的计算资源，虽然余震和人工诱发地震类似，震源的密集空间分布缩小了搜索范围，但是时效上仍然需要提升。因此，从实际应用的角度出发，偏移叠加方法的性能优化将成为未来余震识别的重点发展方向。
4.2 探索深度学习在余震识别领域的应用
波形特征函数是震相到时拾取和偏移叠加方法的基础，因此单台震相拾取方法仍然是提升余震自动识别能力的重点。特征值是对观测数据的抽象，传统方法选取人为设计的数种特征值，通常不需要先验信息，具有计算简单快速的特点，但是单一的特征值无法描述震相信号的所有变化，参数设置也很难适应复杂多变的余震波形。互相关方法利用波形的整体形态特征，压制噪声的能力非常强，震相拾取能力远超传统方法。波形形态中包含与具体事件相关的信息，如震源机制和空间位置，余震的震源机制和空间位置虽然多数非常接近，但也存在偏差较大的情况。因此，互相关方法通常需要准备大量模板，如Peng和Zhao（2009）在識别余震时使用3647个地震事件的波形作为模板，这导致其通用性和时效性较差，自相关方法和FAST方法的出现也没有从根本上解决这一问题。
笔者认为介于传统方法和互相关方法之间的综合特征方法是未来发展的方向，即从波形数据中尽可能多地抽象出与具体余震事件无关的特征，综合这些特征进行震相拾取，应当注意的是，这里的综合特征并不等同于多特征。多特征方法中最典型的属经典机器学习方法，这些方法依赖于人工从数据中提取出的数种特征值。BPNN等全连接神经网络方法理论上可以自动提取特征值，但是层数增加时需要训练的权重过多，实际应用中以3层网络为主，一般通过增加神经元数量来提升性能，与经典机器学习方法一样，比较依赖人为选择。
深度学习方法是近年兴起的ANN子集，优点是拥有更多隐含层，通过多层次的抽象比简单的增加神经元数量具有更强的描述能力，可以从原始数据中抽象出所需特征，这些特征也许我们从未注意到过，人工智能现今的优异成绩几乎都是由深度学习方法所取得。Perol等（2017）使用CNN方法识别和分类美国中部的诱发地震，在取得比匹配滤波方法更高识别能力的同时，计算耗时和存储使用量比自相关方法和FAST方法明显降低。卷积神经元网络（Convolutional Neural Network，简称CNN）是深度神经网络的一种，使用卷积核作为中介，识别数据中存在的固有局部结构，降低权重参数数量，可以有效避免过拟合，十分适合处理图像和声学数据。地震连续波形数据与声学数据有很多共同点，都是记录介质的震动信息，并且存在固有的局部结构，卷积神经网络方法识别震相可能会取得与处理声学数据相近的效果。
深度学习方法的运算量集中在模型训练阶段，使用模型进行预测的速度则非常快，模型的泛化可以提升适应性，因此具备同时兼顾运行效率和通用性的能力，在余震识别领域尚处于起步阶段。深度学习方法，属于黑盒算法，即无法从其训练模型中分析震相特征的抽象过程，以及如何组合特征识别震相，因此无法提升对地震信号的认识，但是优异的性能对于解决震相识别这类实际问题很有意义，未来将成为震相识别，乃至地震监测技术的重要发展方向。
如果深度学习能够准确拾取震相到时，那么在识别余震时，常规走时定位法所面临的低信噪比震相无法拾取和精确性不足问题就迎刃而解，而偏移叠加方法也能够更好地约束搜索范围，提高定位精度。当然，立足更长远的未来，深度学习技术也有可能替代偏移叠加，以更复杂和更抽象的方式实现多台震相联合，甚至进一步由余震识别这个小的领域向地震监测和地震数据处理发展。
本人关于地震定位、震相拾取的理论学习和算法实现都是在中国地震台网中心马延路博士的指导下完成，审稿专家对本文也提出了宝贵意见，在此一并感谢。
参考文献：
房立华，吴建平，张天中，等.2011.2011年云南盈江M55.8地震及其余震序列重定位[J].地震学报，33（2）：262-267.
高淑芳，李山有，武东坡，等.2008.一种改进的STA/LTA震相自动识别方法[J].世界地震工程，24（2）：37-41.
高星，王卫民，姚振兴.2002.用于地震定位的SAMS方法[J].地球物理学报，45（4）：525-532.
黄媛，吴建平，张天中，等，2008.汶川8.0级大地震及其余震序列重定位研究[J].中国科学：地球科学，38（10）：1242-1249.
孔令文.2014.神经网络的分类及其应用[J].电子技术与软件工程，（17）：25.
刘晗，张建中.2014.微震信号自动检测的STA/LTA算法及其改进分析[J].地球物理学进展，29（4）：1708-1714.
刘翰林，吴庆举.2017.地震自动识别及震相自动拾取方法研究进展[J].地球物理学进展，32（3）：1000-1007.
刘劲松，王赟，姚振兴.2013.微地震信号到时自动拾取方法[J].地球物理学报，56（5）：1660-1666.
刘希强，周彦文，曲均浩，等.2009.应用单台垂向记录进行区域地震事件实时检测和直达P波初动自动识别[J].地震学报，31（3）：260-271.
马强，金星，李山有，等.2013.用于地震预警的P波震相到时自动拾取[J].地球物理学报，56（7）；2313-2321.
牟磊育，吴建平.2015.利用单台触发检测台阵数据库事件[J].地球物理学进展，30（6）；2511-2516.
曲均浩，刘希强，吴丹彤，等.2012.基于神经网络的近震与远震识别[J].地震研究，30（3）：360-366.
宋金，蒋海昆.2009.序列衰减与余震激发研究进展及应用成果[J].地震地质，31（3）：559-571.
宋维琪，杨晓东.2011.解域约束下的微地震事件网格搜索法、遗传算法聯合反演[J].石油地球物理勘探，46（2）：259-266，160.
苏金蓉，郑钰，杨建思，等.2013.2013年4月20日四川芦山7.0级地震与余震精确定位及发震构造初探[J].地球物理学报，56（8）：2636-2644.
孙丽.2016.国家地震台网地震定位方法的改进[J].国际地震动态，（11）：12-16.
唐国兴.1979.用计算机确定地震参数的一个通用方法[dl.地震学报，1（2）：186-196.
万永革，李鸿吉.1995.遗传算法在确定震源位置中的应用[J].地震地磁观测与研究，16（6）：1-7.
王彩霞，白超英，王馨.2013.地震震相初至自动检测技术综述[J].地球物理学进展，28（5）：2363-2375.
王文萍，王庆良.1999.利用遗传算法和最小二乘联合反演共和地震位错参数[J].地震学报，21（3）；62-67.
徐锡伟，陈桂华，于贵华，等.2013a.芦山地震发震构造及其与汶川地震关系词论[J].地学前缘，20（3）：11-20.
徐锡伟，闻学泽，韩竹军，等.2013b.四川芦山7.0级强震：一次典型的盲逆断层型地震[J].科学通报，58（20）：1887-1893.
徐锡伟，闻学泽，叶建青，等.2008.汶川MS8.0地震地表破裂带及其发震构造[J].地震地质，30（3）；597-629.
余建华，李丹丹，韩国栋.2011.特征函数响应特性分析及STA/LTA方法的改进[J].国外电子测量技术，30（7）：17-19，23，27.
余娜，张晓清，袁伏全.2016.2016年青海门源MS6.4地震余震序列参数稳定性分析[J].地震研究，39（S1）；62-68，133.
张范民，李清河，张元生，等.1998.利用人工神经网络理论对地震信号及地震震相进行识别[J].西北地震学报，20（4）；43-49.
张广伟，雷建设，梁姗姗，等.2014.2014年8月3日云南鲁甸MS6.5级地震序列重定位与震源机制研究[J].地球物理学报，57（9）：3018-3027.
张苏平，孙艳萍，陈文凯.2013.余震信息在岷县漳县6.6级地震震后应急救援重灾区快速判定时作用探讨[J].地震工程学报，35（3）：465-470.
赵大鹏，刘希强，李红，等.2012.峰度和AIC方法在区域地震事件和直达P波初动自动识别方面的应用[J].地震研究，35（2）；220-225，295.
赵仲和，牟磊育.2005.单台地震自动定位网格搜索法及其MATLAB试验[J]地震地磁观测与研究，26（4）：1-12.
赵珠，丁志峰，易桂喜，等.1994.西藏地震定位——一种使用单纯形优化的非线性方法[J].地震学报，16（2）；212-219.
郑勇，谢祖军.2017.地震震源深度定位研究的现状与展望[J].地震研究，40（2）：167-175，333.
郑韵，姜立新，杨天青，等.2015：利用余震频度分布进行宏观震中快速判定[J].地震，35（2）：121-132.
周民都，张元生，张树勋.1999.遗传算法在地震定位中的应用[J].西北地震学报，21（2）：167-171.
周彦文，刘希强.2007.初至震相自动识别方法研究与发展趋势[J].华北地震科学，25（4）：18-22.
庄建仓，马丽.2000.主震和余震——从大森公式到ETAS模型[J].国际地震动态，（5）：12-18.
Allen R V.1978.Automatic earthquake recognition and timing from singletraces[J].Bulletin of the Seismological Society of America，68（5）：1521-1532.
Allen R.1982.Automatic phase pickers；Their present use and futureprospects[J].Bulletin of the Seismological Society of America，72（6）：225-242.
Amoroso O，Maercklin N，Zollo A.2012.S-wave identification by polari-zation filtering and waveform coherence analyses[J].Bulletin of theSeismological Society of America，102（2）：854-861.
Baer M，Kradolfer U.1987.An automatic phase picker for local andteleseismic events[J].Bulletin of the Seismological Society of Amer-ica，77（4）：1437-1445.
Bai C Y，Kennett B L N.2000.Automatic phase-detection and identifi-cation by full use of a single three-component broadband seismo-gram[J].Bulletin of the Seismological Society of America，90（1）：187-198.
Baillard C，Crawford W C，Ballu V，et al.2014.An automatic kurtosis-based P-and S-phase picker designed for local seismic networks[J].Bulletin of the Seismological Society of America，104（1）：394-409.
Baker T，Granat R，Clayton R W.2005.Real-time earthquake locationusing Kirchhoff reconstruction[J].Bulletin of the Seismological So-ciety of America，95（2）：699-707.
Brown J R，Beroza G C，Shelly D R.2008.An autocorrelation method todetect low frequency earthquakes within tremor[J].Geophysical Re-search Letters，35（16）：1-5.
Crosson R S.1976.Crustal structure modeling of earthquake data，1，si-multameous least squares estimation of hypocenter and velocity pa-rameters[J].Geophys Res，81（17）：3036-3046.
Dai H，Macbeth C.1995.Automatic picking of seismic arrivals in localearthquake data using an artificial neural network[J].GeophysicalJournal International，120（3）：758-774.
Dewey James W.1972.Seismicity and tectonics of western venezuela[J].Bulletin of the Seismological Society of America，62（6）：1711-1751.
Diallo M S，Kulesh M，Holschneider M，et al.2006.Instantaneous polari-zation attributes based on an adaptive approximate covariance method[J].Geophysics，71（5）：V99.
Douglas A.1967.Joint epicentre determination[J].Nature，215（5096）：47-48.
Earle P S，Shearer P M.1994.Characterization of global seismograms u-sing an automatic-picking algorithm[J].Bulletin of the Seismologi-cal Society of America，84（2）：366-376.
Esposito A M，D Auria L，Giudicepietro F，et al.2013.Automatic recogni-tion of landslides based on neural network analysis of seismic sig-nals：an application to the monitoring of stromboli volcano（southernItaly）[J].Pure and Applied Geophysics，170（11）：1821-1832.
Evernden J F.1969a.Identification of earthquakes and explosions by useof teleseismic data[J].Journal of Geophysical Research，74（15）：3828-3856.
Evernden J F.1969b.Precision of epicenters obtained by small numbersof world-wide stations[J].Bulletin of the Seismological Society ofAmerica，59（3）：1365-1398.
Flinn E A.1965.Signal analysis using rectilinearity and direction of parti-cle motion[J].Proceedings of the IEEE，53（12）：1874-1876.
Gajewski D，Tessmer E.2005.Reverse modelling for seismic event charac-terization[J].Geophysical Journal International，163（1）：276-284.
Geiger L.1912.Probability method for the determination of earthquakeepicenters from arrival time only[J].Bull St Louis Univ，8；60-71.
Gentili S，Michelini A.2006.Automatic picking of P and S phases using aneural tree[J].Journal of Seismology，10（1）：39-63.
Gharti H N，Oye V，roth M，Kiihn D.2010.Automated microearthquakelocation using envelope stacking and robust global optimization[J].Geophysics，75（4）：27-46.
Gibbons S J，Ringdal F.2006.The detection of low magnitude seismic e-vents using array-based waveform correlation[J].GeophysicalJournal International，165（1）：149-166.
Giudicepietro F，Esposito A M，Ricciolino P.2017.Fast discrimination oflocal earthquakes using a neural approach[J].Seismological Re-search Letters，88（4）：1089-1096.
Grigoli F，Cesca S，Amoroso O，et al.2014.Automated seismic event loca-tion by waveform coherence analysis[J].Geophysical Journal Inter-national，196（3）：1742-1753.
Grigoli f，Cesca S，Krieger L，et al.2016.Automated microseismic eventlocation using Master-Event Waveform Stacking[J].Scientific Re-ports，6：1-13.
Grigoli F，Cesca S，Vassallo M，et al.2013.Automated Seismic Event Lo-cation by Travel-Time Stacking；An Application to Mining InducedSeismicity[J].Seismological Research Letters，84（4）：666-677.
Hibert C，Mangeney A，Grandjean G，et al.2014.Automated identifica-tion，location，and volume estimation of rockfalls at Piton de la Four-naise volcano[J].Journal of Geophysical Research：Earth Surface，119（5）：1082-1105.
Horstmann T，Harrington R M，Cochran E S.2013.Semiautomated tremordetection using a combined cross-correlation and neural networkapproach[J].Journal of Geophysical Research；Solid Earth，118（9）：4827-4846.
Jackson G M，Mason I M，Greenhalgh S A.1991.Principal componenttransforms of triaxial recordings by singular value decomposition[J].GEOPHYSICS，56（4）：528-533.
Jurkevics A.1988.Polarization analysis of three-component array data[J].Bulletin of the Seismological Society of America，78（5）：1725-1743.
Kao H，Shan S J.2004.The Source-Scanning Algorithm；Mapping thedistribution of seismic sources in time and space[J].GeophysicalJournal International，157（2）：589-594.
Kao H，Shan S J.2007.Rapid identification of earthquake rupture planeusing Source-Scanning Algorithm[J].Geophysical Journal Interna-tional，168（3）：1011-1020.
Kato A，Nakagawa S.2014a Multiple slow-slip events during a fore-shock sequence of the 2014 Iquique，Chile Mw 8.1 earthquake[J].Geophysical Research Letters，41（15）：5420-5427.
Kato A，Obara K，Igarashi T，et al.2012.Propagation of slow slip leadingup to the 2011 MW9.0 Tohoku-Oki earthquake[J].Science，335（6069）：705-708.
Kato A，Obara K.2014b.Step-like migration of early aftershocks follow-ing the 2007 Mw 6.7 Noto-Hanto earthquake，Japan[J].Geophysi-cal Research Letters，41（11）：3864-3869.
Kortstrom J，Uski M，Tiira T.2016.Automatic classification of seismic e-vents within a regional seismograph network[J].Computers and Ge-osciences，87：22-30.
Kiiperkoch L，Meier T，Lee]，et al.2010.Automated determination of P-phase arrival times at regional and local distances using higher orderstatistics[J].Geophysical Journal International，181：1159-1170.
Kurzon I，Vernon F L，Rosenberger A，et al.2014.Real-time automaticdetectors of P and S waves using singular value decomposition[J].Bulletin of the Seismological Society of America，104（4）：1696-1708.
Langet N，Maggi A，Michelini A，et al.2014.Continuous kurtosis-basedmigration for seismic event detection and location，with application toPiton de la Fournaise volcano，La Reunion[J].Bulletin of the Seis-mological Society of America，104（1）：229-246.
Lenglinee O，Enescu B，Peng Z，et al.2012.Decay and expansion of theearly aftershock activity following the 2011，MW9.0 Tohoku earth-quake[J].Geophysical Research Letters，39（18）：L18309.
Leonard M，Kennett B L N.1999.Multi-component autoregressive tech-niques for the analysis of seismograms[J].Physics of the Earth andPlanetary Interiors，113（1-4）：247-263.
Lomax A，Satriano C，Vassallo M.2012.Automatic picker developmentsand optimization：FilterPicker-a robust，broadband picker for real-time seismic monitoring and earthquake early warning[J].Seismo-logical Research Letters，83（3）：531-540.
Maeda N.1985.A method for reading and checking phase time in auto-processing system of seismic wave data[J].Zisin，38（3）：365-379.
Maggi A，Ferrazzini V，Hibert C，et al.2017.Implementation of a multista-tion approach for automated event classification at piton de la Four-naise volcano[J].Seismological Research Letters，88（3）：878-891.
Mcmechan G A，Luetgert J H，Mooney W D.1985.Imaging of earthquakesources in Long Valley Caldera，alifornia，1983 [J].Bulletin of theSeismological Society of America，75（4）：1005-1020.
Mcmechan G A.1982.Determination of source parameters by wavefieldextrapolation[J].Geophysical Journal of the Royal Astronomical So-ciety，71（3）：613-628.
