引言

微震检测在地震学诸多研究中均十分重要，如大地震前震活动性（Kato等，2012；Brodsky等，2014）、地震触发（Hill等，1993；Kilb等，2000；Gomberg等，2001；Peng等，2010；Meng等，2013, 2014；Wang等，2015）、大地震成核（Kato等，2012；Brodsky等，2014；Ellsworth等，2018）、断层形态刻画（Mori等，1990；Dreger等，2000；Kuge，2003；Kao等，2007；Yang等，2009）、水力压裂诱发地震检测（House，1987；Rutledge等，2003；Eisner等，2008；Anikiev等，2014）、小当量核试验检测（Waldhauser等，2004；Wen等，2010；Zhang等，2015a）等，同时，微震检测研究具有一定挑战性。国内外学者提出了多种地震事件识别方法，如应用较多的长短时窗比检测法（Stevenson，1976；Earle等，1994；Kao等，2004, 2007；Liao等，2012；Gharti等，2010；Grigoli等，2013a, 2013b；Drew等，2013；Langet等，2014）、利用互相关技术的波形模板匹配检测方法、震源扫描叠加法和波形自相关检测技术。由于微震记录信噪比较低，长短时窗比检测法虽较简单、计算速度快，但无法有效应用于微震检测中。震源扫描叠加法检测能力较低，限制了对微震震级的检测，对震源位置进行搜索时需利用绝对走时，该过程十分依赖速度模型。因此，本文对应用较广的波形模板匹配检测方法、匹配定位技术、基于图形处理器（Graphics Processing Unit，简称GPU）加速的匹配定位技术（Liu等，2020）与波形自相关检测技术进行研究，应用基于GPU加速的匹配定位技术与双差地震定位法（Waldhauser等，2000）对北京地区19个台站记录的2015年连续地震资料进行分析，获得完备的地震目录，并进行精定位，从而为北京地区地震活动性与发震断层深部构造研究提供基础资料。

1.   波形模板匹配检测方法

存在噪声的情况下，将模板事件波形与输入的待识别连续数据进行互相关是检测已知信号的主要手段，可应用于多个台站的情况（Gibbons等，2006），以此为基础的波形模板匹配检测方法基本检测步骤为：（1）对模板事件波形与连续波形进行相同参数滤波；（2）选定滑动时间窗口与移动步长，进行模板事件波形与连续波形之间的互相关，搜索与模板事件相似的事件，并得到相应的互相关系数（Shelly等，2007；Peng等，2009；Meng等，2013）；接着对多个台站分量记录到的模板事件波形与记录的连续波形进行互相关叠加，并计算平均互相关系数，以此衡量检测到的事件与模板的相似度；（3）将检测到的高于一定阈值的事件作为结果，如果检测结果中出现重复事件，则选取窗口长度内互相关系数最高的事件。判断检测到的事件是否作为地震事件，取决于检测阈值，因此阈值的设定尤为重要。通常对每个互相关系数进行中位数绝对偏差计算（Shelly等，2007）。检测结果中的震中位置代替了相应模板检测事件的震中位置，同时，根据模板事件和检测事件之间的震级比与振幅对数比线性关系计算检测事件的震级（Peng等，2009；Wang等，2015）。在数据处理过程中，Meng等（2013）设置了数据自检环节，为使所有模板事件均实现完美自检（得到的互相关系数＞0.98），在叠加互相关系数时允许1个数据点存在偏移，即在每个时间点，其自身及2个相邻点间的最高互相关系数均用于叠加。进行偏移处理后，可叠加所有分量最佳相关窗口。

通过波形模板匹配检测方法可得到更完整的地震目录，降低地震目录的完备震级。Shelly等（2007）通过波形模板匹配检测方法，将先前记录和定位的低频地震作为模板事件，检测出了日本四国岛震颤活动，同时认为此次震颤活动为低频地震。该方法可在时间和空间上以高分辨率寻找震颤来源，表明震颤和缓慢异化是单个过程的不同表现。Peng等（2009）应用该方法，以发生在帕克菲尔德断层附近的3647个重定位地震事件作为模板，检测2004年9月28日帕克菲尔德M_w6.0地震发生后3d内漏掉的余震事件，捕捉到了比北加利福尼亚地震台网标准目录给出的余震事件多10倍的地震事件，并在此基础上研究了该次地震发生后早期余震的迁移特性。Meng等（2013）应用波形模板匹配检测方法检测了2003年圣西蒙M_w6.5地震发生前46d到发生后31d圣安德烈斯断层帕克菲尔德段的地震，得到的地震事件数量为北加利福尼亚地震台网地震目录数量的8倍，且发现沿圣安德烈斯断裂蠕动部分的地震事件数量在圣西蒙主震后呈显著下降趋势，这与主震产生的负静态应力变化（即应力阴影）密切相关。Kato等（2014）应用波形模板匹配检测方法对2014年智利伊基克M_w8.1地震震源附近记录到的连续波形数据进行了检测，得到的地震事件数量为美国地质勘探局地震台网地震目录数量的10倍，并确定了多个地震迁移序列，还识别到了重复地震，推测为前震序列中沿板块边界断层的地震滑动情况。Walter等（2015）应用波形模板匹配检测方法分析了哥斯达黎加尼科亚半岛2012年9月5日7.6级地震发生前的地震活动，大大提高了地震目录的完备性，确定了主震前约35min、距震源15km以内的小地震（震级 < 2.2），并分析了地震事件触发的可能性。Huang等（2017）通过波形模板匹配检测方法恢复2015年M_w8.0戈尔卡地震发生前80d到发生后30d的地震，检测到的余震数量是国际地震中心地震目录所列余震数量的2倍，发现2015年5月12日7.3级最大余震发生前，发生了1次沿主震破裂断层走向的余震扩展。

波形模板匹配检测方法可有效检测信噪比较低的地震信号，但进行模板事件波形与连续波形互相关时耗时较长。目前，GPU计算已从固定功能的图形设备演变为高度可编程的并行处理器。波形模板匹配检测方法在多个级别上表现出并行性：（1）存在多个模板，每个模板均需与连续波形进行匹配；（2）连续波形本身由一连串的数据点窗口组成；（3）模板和连续波形均包含多个台站记录的多个分量数据；（4）在并行结构的更底层，每个分量均包含多个数据点，这些数据点为互相关分析输入文件，可实现并行。

为降低时间成本，Meng等（2012）提出了基于GPU显著提速的检测算法。为探索这种层次结构的并行性，研究人员将计算分解为多个任务，利用分解后的任务计算模板和滑动窗口内连续波形的相关性，其中每个任务均将1个模板通过滑动窗口与相同长度的连续波形进行匹配。该检测方法的设计集中在以下算法方面：（1）将总任务分组并分配给GPU进行计算；（2）将每个任务的GPU代码优化。与串行的CPU代码相比，Meng等（2012）通过GPU处理将1张英伟达（Nvidia）显卡的速度提高了约40倍。

对于任务分配，需注意以下方面：（1）每组任务不应超过GPU的内存容量，使GPU有足够的内存容纳数据（包括输入、输出和临时处理数据）；（2）分配时应涵盖所有任务，并在各小组间保持均衡，以确保将所有任务分配给GPU，各任务GPU计算时间相对一致；（3）分配任务时应将传送和冗余的计算成本最小化。

Meng等（2012）通过使用GPU提速代码，仅用约48h完成了以下计算：使用2000个模板波形、80d的连续波形，共检测到约24000个地震事件，约为南加州地震台网地震目录数量的70倍，并根据新检测到的事件和南加州地震台网地震目录分析计算了地震活动率。Meng等（2014）利用该方法，在2010年M_w7.2El Mayor-Cucapah主震附近的索尔顿海地热田和圣哈辛托断裂带得到了更完整的地震目录，新检测到的目录揭示了2个研究区域在主震后的地震活动率均显著提高。Liu等（2019）利用波形模板匹配检测方法和双差地震定位法研究门源主震后余震的时空演变，将中国地震台网中心提供的地震目录和重定位后的地震目录中的事件作为模板，扫描了2016年连续波形数据，得到的事件较中国地震台网中心提供的地震目录和重定位后的地震目录数量分别多3倍和11倍。利用更完整的地震目录，对研究区域余震时空演变进行分析。

波形模板匹配检测方法使用的互相关技术对于微弱信号极其敏感，检测到的地震震级远小于其他方法，但仅能检测模板事件附近的地震，并假设待检测地震与模板事件发生在同一位置，不能提供检测事件的位置信息（张淼，2015）。

2.   匹配定位技术

针对波形模板匹配检测方法的局限性，Zhang等（2015a）研发了微震检测新方法，即匹配定位技术，该技术与波形模板匹配检测方法最大不同之处在于，在相关波形叠加前考虑了微震与模板之间的位置差，从而在检测的同时对检测事件进行了定位。

匹配定位技术首先根据已有台网地震目录确定扫描所需的模板事件，并确定模板事件的位置信息；然后在每个台站的所有分量上选取模板事件振幅最大的震相作为参考震相，并将其与连续波形进行滑动互相关；接着以模板事件所在位置为搜索中心，在经度、纬度和深度方向上对可能存在微震的位置进行网格化搜索，搜索过程中首先计算微震每个可能存在的格点位置与参考模板位置到同一台站的到时差，通过到时差对每个台站各分量滑动互相关值进行位置校正；最后计算各台站所有分量叠加后的滑动互相关系数平均值与信噪比，当上述2个参数大于设定的阈值时，即认为捕捉到了1个地震事件，并将该网格点的空间位置确定为相应的震源位置。捕捉到的地震事件震级依靠于所有台站分量中参考震相的振幅比中位数(张淼，2015；Peng等，2009；Meng等，2013)。

匹配定位技术将滑动互相关系数平均值与信噪比作为阈值控制输出结果，通过降低检测标准增加检测数量，降低误检概率。另外，为避免重复检测，在特定时间窗口中选取相关系数最大的地震作为检测到的地震事件。

Zhang等（2015a）应用匹配定位技术检测并定位了朝鲜于2010年5月12日进行的小当量核试验，将朝鲜2009年和2013年核试验产生的地震作为模板事件，对2010年4月1日至2010年5月31日距朝鲜核试验基地200km内的7个地震台站记录到的连续波形进行扫描，得到2010年核试验位置，进一步确认该事件的频谱比Pg/Lg，最终认定其为爆炸类型，为朝鲜2010年5月12日小当量核试验的存在提供了地震学证据，展示了利用地震数据对小当量核试验进行科学检测的能力。Zhang等（2015b）利用匹配定位技术对2007年和2014年日本Ontake火山2次活动过程中产生的连续波形进行扫描，检测到2次活动中的地震数量分别为日本气象厅提供的地震目录数量的37倍和30倍。根据新的完备地震目录，发现2014年日本Ontake火山活动中存在地震数量增多和减少交替出现的规律，该规律出现在大规模地震活动初期（张淼，2015）。

应用匹配定位技术时，首先进行相对走时校正，然后对互相关系数进行最大化叠加，随后在模板周围三维范围内搜索新识别事件的最佳位置，最后提供检测事件的位置信息。该方法虽已使用共享内存并行编程（OpenMP）进行加速计算，但由于计算机CPU内核有限，随着地震数据量的增加和搜索网格的减小，处理数据产生的时间成本仍不可忽略。计算时间主要由匹配定位技术中的波形互相关及校正时间偏移后的互相关过程决定，将任务进行合理分配是亟待解决的问题。Liu等（2020）对匹配定位技术进行改进，得到基于GPU加速的匹配定位技术，该技术为模板事件波形添加了权重因子，并进一步完善了GPU计算问题，首先将模板事件波形和连续波形加载到CPU随机存取存储器；然后计算所有模板事件位置与潜在事件位置的偏移量；接着将模板事件波形和连续波形从CPU随机存取存储器传输到GPU随机存取存储器，再传输到线程块共享内存，每个模板有2个GPU内核循环，以计算滑动相关系数，并根据偏移量将滑动相关系数叠加，将叠加后的滑动相关系数传回CPU随机存取存储器；最后合格的检测结果形成了新目录。基于GPU加速的匹配定位技术增强了检测能力，扩大了检测范围，显著提高了计算效率，且降低了可检测地震的震级，可对检测到的地震事件进行定位。Liu等（2020）利用基于GPU加速的匹配定位技术对2015年阿拉善左旗地震余震进行检测，以中国地震台网中心提供的38次余震作为模板事件波形，扫描主震后7天的连续波形，检测到的余震数量为中国地震台网地震目录数量的20倍左右，并依据余震空间分布与震源机制解，判断此次5.8级地震发生在东西向隐伏走滑断层上。

3.   波形自相关检测技术

Brown等（2008）在解决地颤动中低频地震搜索问题时，提出了波形自相关检测技术，该技术无需先验的地震信号信息，可认为是彻底的盲搜索技术。首先对连续波形进行适当的预处理，如滤波等；然后将各台站各分量连续波形数据按一定时窗长度间隔分成若干部分，由于无先验的模板信息，相对于Shelly等（2007）使用的时间窗口，波形自相关检测技术时窗长度更长；接着将所有可能成对的时窗分别进行互相关并叠加，选取超过一定阈值的事件作为候选事件；最后将所有候选事件对应在同一台站的波形进行额外互相关运算，并按互相关系数高低进行分组并叠加形成高信噪比模板，运用波形模板匹配检测方法进行补充搜索（刘翰林等，2017）。

Brown等（2008）运用波形自相关检测技术对6个台站数据进行处理，选取6s的时间窗进行计算，将大于5倍中位数绝对偏差的事件确定为候选事件，并利用这些事件进行模板匹配滤波扫描，得到287个低频地震事件，通过对低频地震事件进行分析，认为日本西南部的深部震颤是重复的低频地震群，发生在30—35km、地震发生带板块界面下倾面处。

波形自相关检测技术可在模板准备阶段提供帮助，得到更多高信噪比模板事件，但该技术依赖于大量相关计算，对于数据量庞大的研究而言，仍存在计算量大和时间成本过高等问题。

4.   基于GPU加速的匹配定位技术在2015年北京地区微震检测中的应用

4.1   台站分布与研究数据

本文收集了中国地震台网中心与北京市地震局在北京地区布设的19个固定地震台站记录到的2015年连续波形，研究区域位置和台站分布如图 1所示。在台网中心提供的地震目录中筛选245个地震事件作为模板事件，利用基于GPU加速的匹配定位技术对连续波形进行扫描检测。

图 1  研究区域台站分布图

Figure 1.  Distribution of stations in the study area

注：红色虚线表示研究区域断裂带，蓝色三角表示台站。

下载: 全尺寸图片 幻灯片

4.2   检测过程

（1）模板事件和连续波形均进行相同参数的减采样率处理，将采样率由100Hz降至20Hz，以提升计算速度。

（2）对模板事件和连续波形进行2—8Hz带通滤波，根据区域一维速度模型，使用TauP（Crotwell等，1999)计算模板事件P波和S波理论到时，并对到时进行人工核验与校正。

（3）模板扫描窗口为S波到时前1s至后5s，对模板事件每个分量波形数据质量赋予权重因子。考虑研究区域范围、模板分布与时间成本，以模板事件位置为中心，设置搜索区域为8km×8km，搜索步长为2km（图 2）。

图 2  研究区域模板事件分布与扫描范围图

Figure 2.  The distribution of templates and scanning range for each event

注：蓝色三角表示台站，黑色虚线表示断层，红点表示2015年台网地震目录中筛选出的模板事件，灰色实线框表示对应模板事件的扫描范围。

下载: 全尺寸图片 幻灯片

（4）利用双差地震定位法对地震事件进行精定位（杨智娴等，2003；朱艾斓等，2005；张广伟等，2011；胡幸平等，2013），共获得922个精确的地震位置。

4.3   结果与分析

本文利用245个模板事件，经扫描检测得到高达1229个地震事件，约为地震台网地震目录数量的5倍。模板事件为2015年1月15日2.0级地震事件，检测到的事件为2015年1月16日1.01级地震事件（图 3）。精定位后的地震分布如图 4所示，由图 4可知，地震主要集中在黄庄-高丽营断裂（图 4中白线）附近。地震数量-震级关系如图 5所示，由图 5可知，检测到的地震事件以1级以下微震为主。

图 3  利用基于GPU加速的匹配定位技术检测到的1个地震事件

Figure 3.  A seismic event detected by GPU-ML

注：图 3（a）红点表示检测到的事件，红色虚线表示检测标准为中位数绝对偏差的9倍；图 3（b）为对应图 3（a）所示窗口的互相关值分布

下载: 全尺寸图片 幻灯片

图 4  模板事件（黄色）与检测到的地震事件（红色）精定位后震中分布图

Figure 4.  The distribution of the catalog events (yellow) and detected events (red) after relocation

下载: 全尺寸图片 幻灯片

图 5  地震数量-震级关系图

Figure 5.  The magnitude–frequency relation

下载: 全尺寸图片 幻灯片

检测到地震在潭柘寺地区密集分布，震源深度主要为5—10km，分布趋势为北西向，如图 6所示，这补充了台网中心地震目录在该范围的空缺。由于该区域附近的模板事件数量较少，导致检测到的结果数量少，可见南东侧地震分布较深，北西侧地震分布较浅，提示断裂倾向为南东。在黄庄-高丽营断裂带北端，震源深度分布较广泛，以10—20km为主，震源深度较潭柘寺地区深，总体分布情况与台网地震目录相似。在垂直于断裂带走向的剖面（图 6（e）），地震集中分布情况显示出断裂倾向为南东，且断层倾角较大。检测到的地震目录完善了地震台网原始地震目录，更好地揭示了该地区断层形态。

图 6  扫描重定位后的台网地震目录事件（黄色）与检测到的地震事件（红色）震源深度沿不同剖面方向的垂直分布图

Figure 6.  The distributions of catalog events (yellow) and detected events (red) focal depth along the AB, AC, IK, and IJ sections after relocation

下载: 全尺寸图片 幻灯片

5.   总结与讨论

完备地震目录在地震学研究中十分关键且重要，近年来，微震检测与识别是地震学研究的热门课题。本文主要介绍国内外有效且应用广泛的微震识别方法，包括波形模板匹配检测方法、匹配定位技术和波形自相关检测技术，其中波形模板匹配检测方法可在信噪比较低的连续波形数据中识别出与模板相似度较高的地震事件，但不能同时进行地震定位；匹配定位技术弥补了波形模板匹配检测方法的不足，可在模板周围进行全方位检测，将识别到与模板事件相似度高的事件作为结果，且可对检测到的地震进行定位，但时间成本随着搜索网格的缩小而增加；基于GPU加速的匹配定位技术克服了匹配定位技术的缺点，通过添加权重因子，可更有效地进行地震检测和定位，更好地丰富地震目录；波形自相关检测技术虽可在模板准备阶段提供一定帮助，但由于计算量过大，时间成本过高，不便应用于数据量较大的研究。

通过对微震识别方法进行总结与对比，可知匹配定位技术适用于各类微弱信号的检测，在地震迁移过程中可提供精准的位置信息，也可较好地刻画断层走向，对于地震活动性的研究具有十分重要的意义。在今后的研究工作中，可利用该方法获取更完备的地震目录，得到更多的地下断层信息，为组织开展地震区划与各类基础工程建设提供充实的地震资料，方便进行相应的风险评估工作，提高地震防御能力。引入GPU并行计算可提升计算速度，有效提高微震识别效率与精度，不同学科的交叉与结合将成为地震学研究趋势。

利用基于GPU加速的匹配定位技术对北京地区2015年地震资料进行遗漏地震检测和识别，共识别1229个地震事件，约为地震台网地震目录数量的5倍，获得了较精确的地震位置、震级等信息。对检测到的地震事件进行精定位，可知地震在潭柘寺地区密集分布，震源深度主要为5—10km，南东侧地震分布较深，北西侧地震分布较浅，提示黄庄-高丽营断裂倾向为南东；黄庄-高丽营断裂北端震源深度以10—20km为主，揭示出断层倾向为南东，且断层倾角较大。本文研究结果为研究北京地区地震活动性与黄庄-高丽营断裂周边构造提供了重要的地震事件基础数据支撑。基于GPU加速的匹配定位技术能检测出大量台网地震目录中未记录的地震事件，遗漏事件以小于1级的微震为主，从而补充了台网地震目录。