模板匹配滤波技术在地震数据处理中的应用

李璐 王宝善 侯金欣

中国地震局地球物理研究所（地震观测与地球物理成像重点实验室），北京市民族大学南路5号 100081

0 引言

由于受背景噪声水平、记录仪器分布等多种因素的影响，许多微弱的地震信号被淹没在噪声中，无法直接观测。另一方面，大地震发生后的短时间内，由于受强震尾波的影响以及大量余震事件的发生，波形相互重叠，难以区分。如果能够自动识别出这些微弱的地震信号和遗漏的早期余震，则对于完善地震目录、全面分析地震活动性、监测地震破裂区域的震后形变等具有重要意义。

目前，在大量地震自动识别方法中，广泛应用的主要有2种：一种是工业界普遍使用的长短时窗比（Short-Term-Average／Long-Term-Average，STA／LTA）方法（Stevenson，1976；Earle et al，1994）。该方法利用包络线上短时窗内信号的平均值与长时窗内信号平均值的比值来表征信号振幅的变化，其中，短时窗用于截取微弱的有用信号，长时窗内信号的平均值则反映背景噪声水平。虽然算法简单，计算速度快，但对于信噪比过低的信号，这种方法并不适用。另一种方法是基于互相关技术的模板匹配滤波技术（Matched Filter Technique，MFT）（Gibbons et al，2006；Shelly et al，2007；Peng et al，2009）。该方法在滑动窗互相关（Sliding-Window Cross-Correlation，SCC）检测技术的基础上发展起来（Yang et al，2009），是在信噪比较低的情况下提取微弱信号的一种有效方法。其中，Gibbons等（2006）分析了1997年8月16日Kara Sea地震主震后4hr发生的1个微地震，发现传统的STA／LTA方法无法检测出此次事件，而模板匹配滤波技术在单通道记录上即可将其分辨出来。由此可见，模板匹配滤波技术在微地震的自动检测方面具有明显优势。

除了地震学，模板匹配滤波技术还在其他领域得到广泛应用。2015年曾在科学界引起轰动的重大发现——引力波的探测（Abbott et al，2016）中，科学家们利用干涉叠加提高探测器精度，于2015年9月14日09：50：45（UTC）探测到引力波信号，该信号极其微弱，应变约为1.0×10－21。为了进一步验证该微弱信号的来源，科学家们以黑洞模型模拟出理论引力波信号作为模板，在实际观测的连续波形上扫描，发现了与模板波形一致的信号，且时间上与实际观测到的引力波信号一致，从而验证了引力波信号的可靠性。Zhao等（2016）以爆炸事件波形为模板扫描连续波形，并未在2015年8月12日天津危险品仓库爆炸前发现小爆炸事件。此外，在医学、天体物理、电子通信等领域，模板匹配滤波技术还是用于造影、图像识别、语音识别等的重要手段（Hoover et al，2000；Dong et al，2008；Avadhanulu et al，2013）。

本文主要介绍模板匹配滤波技术的基本原理和数据处理流程。从寻找大地震的前震与余震、远程动态触发地震、诱发地震和非火山震颤等方面，详述该方法在地震数据处理中的应用。此外，对于模板匹配滤波技术的改进与发展也进行了简要讨论。

1 MFT基本原理

模板匹配滤波技术（MFT）利用现有地震目录中记录到的事件作为模板，扫描连续波形，识别出与模板事件波形具有一定相似度的地震事件。匹配滤波检测微地震的处理流程如图1所示。

图1 模板匹配滤波技术数据处理流程

第1步，选取合适的频段对模板事件和连续波形作同样的滤波处理。微地震检测一般为区域地震，因此，选取的频段应相对高频。

第2步，标定每个模板事件三分量波形上的P波、S波到时。标定方法可以根据已知速度模型计算理论到时，也可人为手动拾取。再对每一个模板事件的三分量分别计算P波、S波的信噪比。其中，信号项选取P波或S波到时后某一时间窗内的波形计算均方根振幅；噪声项选取在P波到时前同样长度的时间窗内的波形计算均方根振幅。分别用P波或S波的信号项除以噪声项，得到对应的P波或S波信噪比。在后面的波形扫描中，只选取信噪比高于某一阈值的模板事件波形。

第3步，对于某一模板事件，利用下式计算模板事件的波形和对应分量的连续波形的互相关系数（Cross-Correlation coefficient，CC系数））

其中，t0、t1分别为用于计算互相关系数的窗口的开始与结束时刻；X（t）、Y（t）分别为 t0到 t1时间段内的模板事件波形和连续波形。计算时每次在连续波形上移动1个采样点，最终得到该分量上的一道互相关波形。

第4步，根据模板事件上P波或S波到时与该事件发震时刻之间的到时差，将每道互相关波形平移至发震时刻，再对所有台站、所有分量的互相关波形进行叠加，取平均，得到该模板事件的平均互相关波形。

第5步，设定互相关系数的阈值，当平均互相关波形上某一时刻的互相关系数高于该阈值时，就认为这是一个新检测事件。选取阈值时，一般选择绝对中位差（Median Absolute Deviation，MAD）的 9倍（Peng et al，2009）或 12倍（Yao et al，2015）。但是，当台站分布稀疏、波形信噪比较低时，可以根据需要适当提高阈值（如15倍MAD），以减少模板事件与噪声信号间的错误检测。图2为利用MFT以2015年4月26日20：08：33发生于中国藏南地区的1次2.40级地震为模板，检测到发生于4月25日19：43：39的1次1.81级地震。如图2所示，通过扫描在连续波形上可以识别出与模板事件波形相似、没有在地震目录上出现的地震事件。基于重复地震的特性，认为新检测事件位置与对应的模板事件位置一致。根据模板事件和新检测事件在连续波形上的到时差，可以计算出新检测事件的发震时刻。根据某一时间窗内新检测事件与模板事件的振幅比，可确定新检测事件的震级。

2 M FT在地震学中的应用

在算法逐渐完善的基础上，模板匹配滤波技术已经在多种地震事件的检测中得到应用。

2.1 运用M FT寻找前震和余震

许多大地震发生前会有1个或多个前震事件发生，但前震与主震成核过程之间的联系还存在争议。建立更加完整的前震活动序列，有助于加深对地震发生物理过程的理解，为地震预测提供指导（Jones et al，1979；Ben-zion，2008）。Kato等（2012）通过波形互相关的方法，在2011年日本东北9.0级地震之前识别出2个相互分离的前震地震群，并发现它们在时间上有向主震初始破裂点缓慢移动的趋势。通过对滑动速率的分析发现，第2个前震序列产生主要的应力加载，促进主震的破裂过程。Kato等（2014a）对2014年智利8.1级地震发生前后的连续波形进行模板扫描，新检测到数量10倍于原始地震目录的地震事件，并发现前震序列中慢滑移事件同样有向着主震成核点迁移的趋势。另一方面，Wu等（2014）对1999年土耳其7.1级地震前65hr内的连续波形进行了匹配滤波分析，并未发现前震地震活动对主震破裂的促进作用，这说明前震对主震成核的促进作用以及对主震发震位置的影响并不是普遍现象。

大地震发生后的数小时内，在距主震震中几个破裂长度的范围内，地震活动性明显升高，这些被触发的地震称为余震。清楚认识早期余震的活动特性，是分析地震触发机制的重要基础（Helmstetter et al，2009）。模板匹配滤波技术有助于识别出那些淹没在主震尾波中的早期余震事件，以完善余震序列。Peng等（2009）将该方法运用到2004年帕克菲尔德6.0级地震发生后的3天内，共检测到数量11倍于地震目录的余震事件（图3）。检测结果有效补充了早期余震序列，有助于清晰地刻画出余震的迁移特性。

图2 模板匹配滤波检测地震示意图

同样运用模板匹配滤波技术，Lengliné等（2012）在2011年日本东北9.0级地震之后的12hr内检测到数量12倍于原始地震目录的事件，通过这些新检测事件在时间和空间上的分布，得到余震衰减速率（p＝1.0）和余震区沿断层走向的扩展趋势；Meng等（2013）通过匹配滤波后更加完善的地震目录发现，2003年美国圣西蒙6.5级地震后，近场的触发地震与由主震导致的静态剪切应力变化间呈正相关，说明在圣安德烈亚斯断层帕克菲尔德地区存在可能的低摩擦作用；Meng等（2014）通过对比2个距主震不同震中距的区域内的地震活动性认为，静态、动态应力的改变均对近场地震的触发有影响，其中，静态库仑应力可以作用于长期地震活动性的改变，动态应力则主要影响更大区域尺度内的瞬时地震活动性的改变；Kato等（2014b）在2007年日本能登半岛6.7级地震发生后的32天内，同样观测到新检测的余震事件沿破裂断层的走向扩展，同时发现该扩展趋势呈阶梯状，认为可能与主震发生后的抗震余滑（a seismic after slip）有关。

图3 对2004年帕克菲尔德6.0级地震后余震序列的模板匹配滤波检测结果（据Peng等（2009））

2.2 运用M FT寻找远程动态触发地震

除了在近场带来应力改变并产生触发地震之外，大地震还可以传播到几百至上千千米处，在远场触发一些地方震或慢地震（Peng et al，2010；Hill et al，2015）。这种触发多与大振幅面波引起的本地动态应力扰动有关，因此被称为动态触发。选取合适频段滤波后，人为手动拾取地方震到时的方法工作量较大，不利于系统观测。已有研究表明，火山活动区由于岩浆入侵等流体运动，容易达到应力的临界状态，是寻找动态触发现象的重点区域（Moran et al，2004；Prejean et al，2004）。运用模板匹配滤波技术，在2011年日本东北9.0级地震的面波持续期间，Yukutake等（2013）在距主震约450km的箱根火山区域检测到由主震触发的地方震。另一方面，板块交界区域由于需要密切监测地震活动，台站分布密集，有利于动态触发现象的观测，而对于内陆板块的动态触发研究则较少。Yao等（2015）通过模板匹配滤波发现，在西藏中南部由2004年苏门答腊9.1级地震和2005年尼亚斯8.6级地震动态触发的地方震，并对2次大地震动态触发的持续时间和触发区域进行了对比。

2.3 运用M FT寻找诱发地震

在远场触发地震活动中，由人为活动导致的诱发地震是一个重要组成部分。其中，非常规天然气能源——页岩气的勘探开发工作最早起源于美国，并随着水压致裂等多项技术的成熟在多个国家和地区推广使用。而在该技术迅速发展的过程中，与生产活动和提高产量有关的废水注入地下的处理方式是否直接导致了当地地震活动性的增加，一直是科学界和工业界争论的话题（Ellsworth，2013）。在实际观测中，注水带来的应力扰动十分微弱（约1MPa）（Zoback et al，1997），可触发的地震震级较小，不易观测。运用模板匹配滤波技术可以实现对当地微小地震的更加全面的探测，从而有助于判断地震活动性变化与注水活动之间的关联。van der Hilst等（2013）利用模板匹配滤波技术寻找美国中西部地区的触发地震，发现有废水注入的地区，岩石流体孔隙压增加，断层系统更接近临界状态，当有远场的大地震发生时，更容易在应力扰动下产生触发地震。Walter等（2016）同样利用模板匹配滤波技术进一步研究了美国得克萨斯州东部和路易斯安那州西北部的地震活动性，新检测到的40个地震事件在空间上可分为几个地震群，虽然没有强烈的证据表明当地地震活动增加与废水注入间有直接关系，但也并不能完全否定两者间的相关性，作者认为该地区需要更多的地震学监测。在我国，也有一些研究工作将模板匹配滤波技术应用于微小诱发地震事件的检测中。Wang等（2015）在北京房山地区于2011年日本东北9.0级地震和2012年印度洋8.6级地震发生后的短时间内，分别新检测到分布集中的大量浅部地震，并认为可能与当地采矿活动带来的地下应力变化有关。

2.4 运用M FT寻找非火山震颤

非火山震颤（Non-Volcanic Tremor，NVT）主要发生在俯冲带内，由许多相互重叠的低频事件（Low Frequency Event，LFEs）组成（Shelly et al，2007、2011），具有持续时间长、频率较低、体波震相不清晰等特点（Beroza et al，2011），其信噪比较低，不易观测。Shelly等（2009）在美国帕克菲尔德地区第1次利用波形互相关观测到LFEs，并对这一组LFEs组成的震颤进行定位，发现圣安德烈亚斯断层可能延伸至地壳底部，比震源深度最大的普通地震还要深约10km。利用模板匹配滤波在震颤中发现更多重复信号（LFEs），有助于更好地刻画俯冲带内断层的破裂方向与可能的触发机制（Aguiar et al，2009；Brown et al，2010）。

3 M FT方法的改进与发展

在模板匹配滤波技术的基本原理之上，一些研究工作对其算法进行了优化。首先，模板匹配滤波是在连续波形上逐个采样点扫描，运算量大，耗时长。因此，可以通过降采样的方式，适当缩短扫描的时间。此外，Meng等（2012）运用图形加速卡（Graphic Processing Unit，GPU）的并行计算方法，将连续波形分成几个区块，每个区块在GPU的不同线程下同时进行匹配滤波扫描。若同时启用N个线程，则计算速度可以提高到原来的N倍。其次，前人研究工作中，对模板事件波形的垂直分量和水平分量分别选取P波、S波到时附近的时间窗（Peng et al，2009），当研究区域内台站分布稀疏时，可以如图2所示，对模板事件波形的每一分量同时选取P波和S波附近的时间窗，则模板事件波形的数量翻倍，更有利于互相关波形叠加时对不相关信号的压制。最后，与重复地震的特征类似，模板匹配滤波的前提假设是，新检测到的事件与模板事件的发震位置一致。然而，在实际观测中，有些遗漏事件可能具有与模板事件高度相似的波形，但位置上略有偏差，这导致它们传至台站的用时略有不同，此时，如果将各台站的互相关波形按模板事件的到时平移至发震时刻，则会导致该遗漏事件的互相关信号没有完全对齐，此时的叠加波形上互相关系数减小，降低了该事件被检测到的可能性。Zhang等（2015a）发展了一种Match and Locate（M＆L）改进算法，在模板事件附近的1个小三维区域内进行网格搜索，将新检测事件定位在每一个网格点上，计算此时该事件和对应的模板事件之间的到时差，将各道互相关波形按对应的到时差平移后再叠加，最后比较所有网格点的叠加互相关系数和信噪比，以确定事件的位置。由于是在小范围内进行空间搜索，因此，对速度模型的依赖性较小。该方法不仅实现了对新检测事件的重定位，而且提高了被检测事件的互相关系数，从而增加了更小震级事件被探测到的可能性。Zhang等（2015c）利用此方法监测到朝鲜地区的多次核爆实验，因而将全球的探测核实验能力提高到吨的级别。运用同样的方法，Zhang等（2015b）在日本Ontake火山2007、2014年2次喷发前均探测到密集的微震活动，为火山活动的预测提供了依据。

4 结语

通过模板匹配滤波技术可以在背景噪声水平较高时利用互相关扫描的方法寻找到肉眼难以观测到的信号，该方法适用于各类微弱信号的检测。在地震学中，运用该方法监测地震活动性对于断层走向的刻画、地震破裂过程模拟、地下速度结构成像等方面的研究具有重要意义。基于模板匹配滤波技术的自动扫描、识别信号的特征，今后工作中，可将该技术应用于超密集台阵的自动化处理，以提高地震监测工作的效率。于2015年全面启动的国家地震烈度速报与预警工程，将在全国范围内建成由5000余个强震台站组成的密集观测台网用于灾情的快速分析，并指导抗震救灾工作。利用模板匹配滤波技术，将地震事件的自动识别与定位相结合将是处理这些海量地震数据的有效手段。