引言

紫坪铺水库位于四川省都江堰市岷江上游的麻溪乡，于2005年9月开始蓄水，其正常蓄水位为877 m，当水库水位达到850 m的汛限水位时，水库水深最大约为100 m。大地构造位置分布显示，紫坪铺水库库区及邻近区域（103.3°E—103.7°E，30.85°N—31.15°N）位于青藏高原东缘的松潘—甘孜地块与扬子陆块过渡带上的龙门山断裂带内中南段，处于中央北川—映秀断裂与前山安县—灌县断裂之间（周宏伟等，2007）。

紫坪铺水库距离2008年汶川M_S8.0地震震中仅十几千米，地震发生后，国内外许多学者对紫坪铺水库与汶川地震之间可能的关系展开了研究（Kerr，Stone，2009，2010；Deng et al，2010；Ruan et al，2017）。先前的研究大多使用四川区域地震台网（程万正等，2010）或当地紫坪铺水库地震台网（雷兴林等，2008；马文涛等，2011）数据资料，但是一些震级较小的地震由于信噪比较低往往会遗失在原始地震目录中，使地震目录的完整性降低。地震目录是区域地震活动性研究重要的基础资料，编目的完整性直接影响后续研究结果的可信性（Woessner，Wiemer，2005）。

地震的识别一般利用两类方法（Grigoli et al，2018；Zhang et al，2019）：基于震相拾取（Allen，1982；Withers et al，1998）方法和基于波形（Kao，Shan，2007；Shelly et al，2007；Zhang，Wen，2015）识别方法。在检测小震事件时，基于波形的方法通常比基于震相拾取的方法表现更好（Pesicek et al，2014；Grigoli et al，2018），其中基于波形互相关的匹配滤波技术（Match Filter）通过模板波形与连续波形中潜在小震信号之间的叠加互相关来检测小震事件，已得到了广泛的应用（Shelly et al，2007；Peng，Zhao，2009）。然而，该方法只能检测那些紧邻参考模板的地震，不能对检测到的地震事件进行定位。Zhang和Wen （2015）提出了匹配定位方法（Match and Locate），该方法采用了与匹配滤波技术中波形互相关检测类似的概念，但考虑了小震事件与模板可能存在的位置差异，同时提供了高精度的事件位置信息。但由于该方法需要对模板周围的三维空间进行搜索，因此运行起来非常耗时。为了提高匹配定位的计算效率和探测能力，Liu等（2020）发展了基于图像处理器加速的匹配定位（graphics processing unit-based match and locate，缩写为GPU-M&L）方法，综合考虑了几何扩散和地震波形的信噪比对模板事件各个分量赋予不同的权重系数，在通过图像处理器加速的同时也提高了地震事件的探测能力。

本文基于紫坪铺水库地震台网所提供的2005—2008年连续波形数据，首先以四川省地震局提供的地震目录重定位后的地震事件作为模板，使用GPU-M&L （Liu et al，2020）对紫坪铺水库及邻区进行遗漏地震事件的检测。然后，利用基于深度学习算法的去噪技术Deep Denoiser （Zhu et al，2019）对检测出的事件进行去噪处理，并设计了实现噪声和事件分类的卷积神经网络模型对检测出来的事件进行筛选。最后，再次使用相对定位方法hypoDD （Waldhauser，Ellsworth，2000）对筛选后的检测事件精定位，依据得到的地震目录和重定位结果进行地震时空分布特征、地震活动性与深部地质构造等分析，以期为后续的研究工作提供一定的数据支撑。

1.   数据和方法

本文使用四川省地震局提供的地震目录和紫坪铺水库地震台网2005—2008年的连续波形资料，研究区范围为（30.8°N—31.2°N，103.3°E—103.7°E），涉及7个地震台站，台站分布如图1所示。地震目录中记录研究区域内的地震事件1 241个，连续波形数据记录到了1 017次地震事件，其中最大震级为M_L5.2，M＞4地震一共有17次。

图  1  研究区断裂、台站及模板事件重定位前后的震中分布图

F₁：汶川—茂汶断裂；F₂北川—映秀断裂；F₃安县—灌县断裂

Figure  1.  Distribution of the faults，stations，and epicenters of template events before and after relocation in the study area

F₁：Wenchuan-Maowen fault；F₂：Beichuan-Yingxiu fault；F₃：Anxian-Guanxian fault

下载: 全尺寸图片 幻灯片

数据处理主要分为以下几个步骤：

1） 利用深度学习算法PhaseNet对地震目录中事件进行P波、S波震相拾取，依据拾取结果挑选模板事件，并利用相对定位法hypoDD对其进行精定位；

2） 利用GPU-M&L方法进行微震检测；

3） 利用基于深度学习的去噪技术对检测到的事件进行去噪处理，并依据去噪后波形各分量的时频图，训练卷积神经网络模型对检测事件进行筛选；

4） 再次利用PhaseNet和hypoDD对筛选后的检测事件进行震相拾取和重定位。

1.1   模板事件选取与精定位

Zhu和Beroza （2019）发展了基于深度学习的震相拾取算法PhaseNet，该算法使用三分量地震波作为输入，输出P波、S波和噪声的概率分布。与现有方法相比，该算法在震相拾取的准确度方面有显著的改进，特别是对于S波。本文利用该算法对四川省地震局所提供地震目录中的事件进行P波、S波初动拾取，根据经验设置P波和S波的拾取可信度阈值为0.5，并舍去拾取台站数小于3的PhaseNet拾取结果。当在同一事件中拾取到多个相同震相时，保留可信度较高的震相到时，并利用TauP计算理论到时标记未被PhaseNet拾取到的震相，最终通过人工读取波形进行校正。基于上述步骤拾取的P波、S波初动，本文以S波到时前1 s为起始的6 s时间窗作为信号窗口，P波到时前6 s为起始的6 s时间窗作为噪声窗口，计算了地震事件每个分量的信噪比，并依据信噪比大于5的分量至少为9个的原则，共挑选了718个高质量的地震事件。

为提高模板事件震源位置的分辨率，本文基于PhaseNet得到的P波、S波到时记录，利用相对定位方法hypoDD对挑选的718个高质量事件进行了精定位。在该过程中，利用波形互相关技术对所有的地震对走时差进行校正，用于校正的波形窗口长度为P波到时前0.2 s至后1 s，滑动窗口为0.6 s；S波到时前0.5 s至后1.5 s，滑动窗口为1 s。地震对之间的最大震源位置差为20 km，最少震相连接个数为4，并选用前人层析成像结果作为初始速度模型。最终得到了578个模板事件的高分辨率震源位置（图1）。

1.2   基于图像处理器加速的匹配定位技术

基于上述重定位后的模板事件，本文利用基于图像处理器加速的匹配定位法（GPU-M&L）进行微震识别，以提高匹配定位的计算效率和探测能力，其程序与步骤主要遵循先前的研究（Liu et al，2020），仅在此简要描述。首先，对连续波形和模板事件都进行2—8 Hz带通滤波，来放大有用的地震信号，减轻背景噪声和大型事件尾波的影响，同时为减小运行内存提升运行速度，将数据从100 Hz减采样到20 Hz。然后，扫描时间窗口为模板事件S波到时前1 s至后5 s，并设置当在6 s的时间窗口内检测到多个事件时，只保留互相关系数（CC）值最高的事件。在扫描过程中，匹配定位法会将每一个网格点内各个台站分量的滑动互相关系数进行叠加，当平均滑动互相关系数大于所设阈值时，这一时刻被设置为一个可能事件的发震时刻，这一空间位置被设置为一个可能事件的震源位置，并根据所有分量中检测事件与模板事件之间峰值幅度比的中值计算其震级大小。

在以往的相关研究中，常常使用9—15倍的中值绝对偏差（MAD）作为检测阈值（Meng et al，2013；Walter et al，2015；Meng，Peng，2016；Ruan et al，2017）。为确定适用于本研究的阈值条件，表1统计了分别以9，12和15倍的MAD作为阈值条件下的检测事件数量和低CC （＜0.2）检测事件所占的比例。综合考虑检测数量与错误检测率，本文选取了12倍MAD值作为检测阈值。同时，由于低CC （＜0.2）值表明其可能是一个错误的检测（Meng et al，2013；Huang et al，2017），本文以CC＝0.2作为另一个条件对结果进行约束，最终得到2万零746个检测事件。图2显示了当地时间2008年6月6日02:21:33检测到的事件，这一事件未在当地地震目录中列出。该检测模板为2008年6月6日M_L1.4地震事件，检测的平均CC值为0.792 1，推断震级为M_L2.10。

表  1  不同阈值下的检测事件

Table  1.  The detected events under different thresholds

阈值	检测事件个数	互相关系数＜0.2的事件占比
9MAD	109 822	21.15%
12MAD	21 213	2.20%
15MAD	8 643	2.13%
注：MAD为中值绝对偏差。

下载: 导出CSV 

| 显示表格

图  2  基于图像处理器加速的匹配定位技术检测结果示例 （起始时间为2008−06−06 00:00:00）

（a） 平均互相关系数CC值，红色虚线表示检测阈值为0.3，红色圆点表示检测到的地震事件；（b） 对应图（a）中所示窗口互相关系数的分布统计图；（c） 连续波形（灰色）与模板事件（红色）在时间窗口中的比较，台网台站和各分量CC值分别标记于波形左右两侧

Figure  2.  An example of detection based on the GPU-M&L technique

（a） Mean CC value trace，where the red dashed line marks the detection threshold of 0.3，and the positive detection is shown by the red dot；（b） Histogram of the CC value distribution corresponding to the window shown in Fig. （a）；（c） Comparisons of continuous waveforms （gray） with templates （red） in a zoom-in time window. The stations of seismic networks and CC values are given to the left and right，respectively

下载: 全尺寸图片 幻灯片

1.3   检测事件筛选

为进一步提高检测事件的真实性，本文对上述获得的检测事件进行了如下的筛选。首先，通过利用基于深度学习的去噪方法DeepDenoiser （Zhu et al，2019）对检测事件进行了去噪处理，并依据短时傅里叶变换对去噪后各分量时间域波形绘制了发震时刻后30 s的时频谱图，共得到36万7079幅时频谱图。然后，结合去噪后的时间域波形图，将时频谱图人工分类标记为“事件”和“噪声”两类标签，得到了分类后的事件和噪声的时频图各5 000幅，图3为去噪后事件和去噪后噪声的时频谱图示例。

图  3  去噪后地震事件（a，b）和去噪后噪声（c，d）的时频谱图

Figure  3.  Time-frequency spectra of events （a，b） and noises （c，d） after denoising

下载: 全尺寸图片 幻灯片

将上述1万幅标记好的时频谱图以3 ∶ 1 ∶ 1的比例分为训练集、验证集和测试集，训练卷积神经网络模型，实现噪声和事件分类。卷积神经网络结构如图4所示，具体包括3个卷积层，卷积核个数分别为32，64，128，卷积核个数的增加保证了对图像特征的充分提取；卷积核大小均为3×3，移动步长均为1；每一个卷积层后增加Relu层进行非线性激活，并在卷积层与非线性部分之间使用批量归一化层来加速网络训练并降低对网络初始化的敏感度（Ioffe，Szegedy，2015）。3层池化层均采用最大池化的方法，池化窗口大小均为2×2，滑动窗口步长均为2；两个全连接层的隐藏单元数分别为512和2，分别使用Relu和Softmax激活函数进行激活，并在第一个全连接层前加入Dropout层防止过拟合，提升模型的泛化能力（Srivastava et al，2014）。最终得到的模型在测试集上的准确率为96%。将该模型应用于上述得到的检测事件的各分量时频谱图，并将被分类为“事件”的时频谱图总数大于9的检测事件筛选为一个真实的地震事件，最终得到了1万6 836个筛选后的地震事件。

图  4  卷积神经网络结构

Figure  4.  Convolutional neural network structure

下载: 全尺寸图片 幻灯片

1.4   检测事件精定位

小震活动包含了丰富的构造信息，其位置的精确程度直接关系到对研究区构造活动的认识。对于水库地震，震源深度也是判断水库触发地震的重要基础（姜金钟等，2016）。因此，本文再次按照与1.1节中相似的步骤和参数，利用深度学习算法PhaseNet对筛选后的1万6836个检测事件进行震相拾取，并利用波形互相关技术得到了互相关系数大于0.5的96万1899条P波和33万1331条S波互相关走时差数据。双差定位法可以同时使用绝对走时和互相关走时差数据进行定位，但当一个地震对在同一台站的绝对走时差和互相关走时差相差较大时，容易造成定位结果不稳定。因此，为了避免这种情况，本文以精度更高的互相关走时差数据为主进行双差定位，最终得到了1万5745个事件的高分辨率震源位置。

2.   结果分析

2.1   地震事件检测结果

本文以台网目录中经过重定位的578个事件作为模板事件，利用基于图像处理器加速的匹配定位技术进行事件检测，筛选后得到的事件高达1万6836个，约为台网地震目录的地震记录数量的13倍。检测目录的M-t图 （图5）显示检测出来的地震事件在2005年至2008年4月期间分布较为零散，集中分布在2008年5月之后的时段内。图6给出了台网地震目录与检测目录的震级-频率关系，可以看出，GPU-M&L识别出了大量台网地震目录中遗漏的地震事件，且为震级集中于M_L2.0以下的微震事件。同时，本文使用地震活动性分析软件ZMAP （Wiemer，2001），基于最大曲率法来估算地震目录的完备震级M_c，相较于台网震级目录，基于GPU-M&L获得的检测目录使得完备震级从1.4降低到−0.1。

图  5  检测目录的M-t图

Figure  5.  The M-t plot of the detected catalog

下载: 全尺寸图片 幻灯片

图  6  台网地震目录和检测目录的震级-频率关系图

Figure  6.  The magnitude-frequency relationship for the local catalog and detected catalog

下载: 全尺寸图片 幻灯片

根据古登堡-里希特（Gutenberg-Richter，缩写为G-R）的震级-频次关系式（Gutenberg，Richter，1944），b值直接反映了大小地震之间的比例关系，也是用于分析水库触发地震的重要参数。实验室及大量的实际观测（Scholz，1968；Wiemer，Benoit，1996；Murru et al，2007）都充分表明了b值与应力呈负相关关系，b值越低，表明应力越高。与同震级天然构造地震相比，由于水库蓄水造成地下介质孔隙压力增大或者水的润滑作用使断层摩擦力减小等影响，水库触发地震的应力降值明显偏低，即b值较高（do Nascimento et al，2004；Mekkawi et al，2004）。小震级事件的遗漏可能会导致计算的b值偏低，大地震后的低b值也可能指示了地震台网无法检测到许多小震级的地震事件（Chang et al，2007）。因此，本文利用MPU-M&L补充了更多小震级事件的检测目录，统计自2005年9月蓄水以来研究区大于完备震级（M_c−0.1）的震级-频度关系（图7），并通过最大似然法（Aki，1965）计算得到b值及拟合误差为0.632±0.006，这一结果与Zhang和Zhou （2016）利用1985—2012年中国地震台网中心地震目录及刘雁冰和裴顺平（2017）利用1992—2014年龙门山地区全部M＞2.0地震目录计算所得的研究区所属区域b值相比较小，指示较高的区域应力，即研究区蓄水后的地震活动可能主要是构造应力积累导致的天然地震活动。

图  7  水库蓄水后的检测目录震级-频率关系图

Figure  7.  The magnitude-frequency relationship for the detected catalog after reservoir impoundment

下载: 全尺寸图片 幻灯片

2.2   检测事件精定位结果

本文再次利用PhaseNet对筛选后的检测目录中地震事件进行P波、S波初动拾取，联合波形互相关技术使用相对定位法hypoDD对其进行精定位，定位结果水平和垂直方向的平均误差分别为42 m，51 m和38 m，发震时刻误差为7 ms，最终得到了1万5745个事件的高分辨率震源位置（图8），其水平方向误差分布如图9所示。重定位后的地震事件集中分布在北川—映秀断裂和汶川—茂汶断裂之间，延展分布趋势与区域内三条大断裂的北东走向基本一致。将重定位后的事件分别沿平行断裂走向方向（剖面AB）和垂直断裂走向方向（剖面CD）进行投影。剖面AB （图8a）显示，重定位后的事件震源集中分布在0—20 km深度范围内，且沿剖面方向距剖面A端约20—30 km位置处存在一个地震空区，与三维空间分布中（图10）绿色圆圈所指示的位置一致，对应于水库库体覆盖区以下的区域。剖面CD （图8b）显示，重定位后震源深度在0—5 km和10—15 km存在优势分布，两层大致平行，呈水平状。在深度7 km左右存在一地震空区，与龙门山地区三叠系雷口坡组和嘉陵江组主滑脱面位置一致（周斌等，2010）。重定位后检测目录的P波和S波走时分布如图11所示，其数据离散程度和整体趋势较为合理，但考虑到研究区范围过小且深度剖面上存在的空区较小，在走时分布图中未能呈现明显的与震源深度分布相对应的“两支现象”。

图  8  检测事件重定位后的震中分布（a）及沿剖面AB和CD方向的事件投影分布（b，c）

F₁：汶川—茂汶断裂； F₂：北川—映秀断裂；F₃：安县—灌县断裂

Figure  8.  The distribution of relocated earthquakes （a） and projection of the detected events along AB and CD （b，c）

F₁：Wenchuan-Maowen fault；F₂：Beichuan-Yingxiu fault；F₃：Anxian-Guanxian fault

下载: 全尺寸图片 幻灯片

图  9  检测事件的水平方向重定位误差分布

Figure  9.  The distribution of horizontal relocation errors of the detected events

下载: 全尺寸图片 幻灯片

图  10  检测事件重定位后的三维分布图

Figure  10.  3-D distribution of the detected events after relocation

下载: 全尺寸图片 幻灯片

图  11  检测目录重定位后的P波、S波走时分布图

Figure  11.  Travel time distribution of P wave and S wave of the detected events after relocation

下载: 全尺寸图片 幻灯片

震源深度也是判断水库触发地震的重要参数，水库触发地震的震源深度一般不超10 km （马文涛等，2013）。杨智娴等（2004）对中国中西部地区（21°N—36°N，98°E—112°E）1992—1999年间的6 459次地震进行了重定位，精定位后的震源深度在0—10 km，0—15 km和0—20 km范围内的地震事件比例分别为52%，71%和82%。而Yang等（2005）利用双差定位法对中国中西部地区（21°N—36°N，98°E—112°E）1992—1999年间的1万零57次地震进行重定位的结果显示，平均震源深度为11 km，震源深度在0—15 km和0—20 km范围内的地震事件占地震总数的比例分别为77%和91%。统计本文双差定位结果显示，水库蓄水后库区平均震源深度为5.14 km，震源深度小于10 km的地震数占地震总数的80.9%。将其与天然构造地震深度进行对比（杨智娴等，2004；Yang et al，2005），指示库区内存在着水库触发地震活动。

3.   讨论与结论

汶川地震之后，为探究紫坪铺水库与汶川地震之间可能的关系，国内外许多研究者开展了有关紫坪铺水库蓄水后地震活动性的研究。Deng等（2010）认为根据库仑应力改变，紫坪铺水库的蓄水仅能造成水库附近地区浅层地震（＜5 km）的频率增加。Ruan等（2017）通过对汶川地震前5个月微震活动的研究，未发现当地地震活动性与紫坪铺水库水位之间可能的相关性。张致伟等（2009）和程万正等（2010）则认为，紫坪铺水库蓄水以来增加的地震活动从表面上看与水库的加卸载过程有关，但从长期序列来看，其频次和强度仍在区域地震活动的正常起伏范围之内。

本文基于紫坪铺水库地震台网提供的2005—2008年连续波形数据，利用GPU-M&L技术进行了地震检测。相较于前人的研究（Ruan et al，2017；Huang et al，2020；Liu et al，2020）， 本文为了确保检测事件的真实性，增加了对检测事件进行筛选的步骤。通过使用深度学习算法DeepDenoiser对检测事件进行去噪处理，并基于去噪后波形绘制的时频谱图，设计了一个卷积神经网络模型对其进行筛选，最终得到的地震目录事件增加到1万6836个，极大地丰富了原始地震目录，提高了地震活动性分析的可靠性。

将最终得到的检测目录与台网地震目录给出的重定位前后的震源分布进行对比，可以发现，台网地震目录重定位前（图12中灰色圆点）在平面和深度剖面上的分布较为分散，无明显的优势震源分布；重定位后的台网地震目录（图12中红色圆点）在平面分布上的特征与检测目录的结论类似，在沿断层走向和垂直断层走向两个方向的深度剖面中存在一个深度处于5—10 km范围的地震空区。检测目录中大量地震事件重定位后的结果将该地震空区的范围缩小，在沿断层走向深度剖面上勾勒出了更明显的水库库体覆盖区下方的地震空区，在垂直断层走向的深度剖面上将上述地震空区的范围缩小到7 km深度左右。

图  12  台网目录重定位前后震源分布（a）及沿剖面AB和CD方向的事件投影分布（b，c）

F₁：汶川—茂汶断裂； F₂：北川—映秀断裂；F₃：安县—灌县断裂

Figure  12.  The distribution of earthquakes before and after relocation （a） and projections of the events along AB and CD profiles，respectively （b，c）

F₁：Wenchuan-Maowen fault；F₂：Beichuan-Yingxiu fault；F₃：Anxian-Guanxian fault

下载: 全尺寸图片 幻灯片

龙门山断裂带发育多层次的滑脱层，历史地震的优势震源深度可能反映龙门山断裂带的推覆逆冲面的铲式变化及其下伏的滑脱层位置（李勇等，2013）。滑脱层作为一个塑性软弱区域，易发生形变，且较难积累应变能，应力容易在其上部的脆性地层中积累（李大虎等，2021）。检测目录重定位后的结果显示，地震活动集中分布在0—20 km深度范围内脆性程度高、渗透稳定性低的岩层中，而7 km左右深度处为岩性较软弱、渗透稳定性高的三叠系雷口坡组和嘉陵江组白云岩和膏盐，少有地震发生，表现为一个明显的地震空区。依据本文的震源分布的优势深度，推测在研究区壳内20 km左右的深度存在一滑脱层顶面，与朱介寿（2008）对跨龙门山爆破地震剖面研究得出的松潘−甘孜地块及龙门山推覆体的中地壳内（20—30 km深）普遍存在一个厚约10 km的滑脱层的结果相一致。这一结果也证实了本文定位结果的可靠性。

前人研究总结出水库触发地震的震源深度一般不超过10 km，震中分布一般距库区不超过10 km，b值相较于区域构造地震的b值较高（胡毓良，陈献程，1979；丁原章等，1983；华卫等，2012；马文涛等，2013）。本文的重定位结果显示紫坪铺水库库体覆盖区下方的地震活动非常少，并在水库覆盖区下方呈现出明显的地震空区，与周斌等（2010）的研究结果一致。利用地震检测目录统计自2005年9月蓄水后研究区的震级−频度关系计算得到的b值结果显示，与前人计算的研究区所属区域b值相比本研究结果较小，指示了在研究时段内区域整体的平均应力水平较高。然而重定位结果显示，80%左右的地震震源深度不超过10 km，结合水库蓄水后研究区地震活动与水库加卸载过程有一定程度吻合的结果（程万正等，2010；张致伟等，2009；周斌等，2010），推测研究区蓄水后的地震主要是天然构造地震活动，同时混杂发生水库触发地震。

总之，本文利用紫坪铺水库地震台网7个台站2005—2008年的波形记录，采用基于图像处理器加速的匹配定位技术（GPU-M&L）对其遗漏的地震事件进行检测和筛选，最终得到的检测目录事件高达1万6 836个，约为台网目录记录的地震数量的13倍，并使完备震级从台网目录的1.4降低到−0.1。然后，通过再次使用结合波形互相关的双差定位法进行重定位，结果显示地震事件沿北东方向展布于北川—映秀断裂和汶川—茂汶断裂之间，水库覆盖区下方存在一个明显的地震空区，且地震活动集中分布在0—20 km深度范围内脆性程度高、渗透稳定性低的岩层中。结合研究区域地质构造背景，优势震源深度较好地指示了研究区地壳内滑脱层位置。最后，根据检测目录计算的b值和震源深度分析结果推测研究区蓄水后的地震主要是构造应力积累导致的天然地震活动，并伴随水库触发地震混杂发生。

四川省地震局为本文提供了数据支持，香港中文大学刘敏博士为本文提供了GPU-M&L代码，加州理工学院朱尉强博士为本文提供了DeepDenoiser以及PhaseNet代码，审稿专家为本文提供了宝贵的修改意见及建议，本文大部分图件采用 GMT 绘制，作者在此一并表示感谢。