地震目录是区域地震活动性研究、地震危险性分析和地震地质研究的重要基础资料。地震目录是否完整、震源参数是否准确，直接影响后续研究的可靠性与科学性。强震发生后或震群活动时，短时间内会发生大量地震，不同的地震事件波形相互重叠，仅基于观测现象识别震级较小的地震难度较大，此时信号自动识别方法里的模板匹配方法可以很好地解决该问题。模板匹配方法是以互相关技术为基础，它的前提假设是：如果被同一台站记录到的两个波形呈现出很高的相似性，我们可以认为这两个震源相距很近，且震源机制相同。国内外利用模板匹配方法研究遗漏地震以及地震活动性等均获得了比较好的效果，如Peng等^[1]利用模板匹配方法，对2004年Parkfield地震的余震目录进行完备，并根据完备的地震目录研究余震的迁移; Meng等^[2]利用模板匹配方法研究完备的地震目录发现，2003年San Simeon M_W6.5地震之后，Parkfield地区的地震活动性与静剪切应力间存在比较好的对应关系; Zhang等^[3]利用匹配滤波的方法实现了对朝鲜地区多次核爆的检测，且给出了几次核爆的相对位置。

2016-01-21 01：13青海门源发生M_S6.4地震，这是甘青交界地区自2013-07-22甘肃门源地震后发生的最大地震^[4]。主震后3 h内国家测震台网总共统计余震序列332次，其中单台共57次，单台事件地震观测报告仅给出其震级，而未给出震中位置。本文利用模板匹配方法检测遗漏地震震相到时，如果能在4个以上台站检测到震相，则利用测震台网常用的HYPOSAT方法估计其震中位置^[5]。速度模型选择甘肃-青海区域甘青一维平均速度模型^[6]，并利用多个台站记录波形与模板地震的振幅比估计其震级。之后计算主震发生后不同时间的最小完备性震级，并通过线性拟合得到最小完备性震级随时间变化的表达式，以分析此地震余震序列的目录完备性。

1 遗漏地震检测与单台震相补充 1.1 模板匹配方法简介

大震发生后地震目录往往会遗漏较多的余震事件^[1]，造成余震遗漏的主要原因在于主震及较大余震的面波等后续震相及尾波的干扰，其中所含低频成分较多^[7]。为了抑制低频干扰的影响，本研究使用4阶零相移Butterworth滤波器对连续波形进行截止频率为15 Hz的高通滤波，选取M_L2.0以上地震事件作为模板地震，然后进行2~8 Hz的滤波。模板波形每个分量的信噪比(SNR)都要大于3，噪声能量水平由P波到时前6~2 s计算得到。根据观测报告以直达S波到时为中心，截取其前2 s~后2 s的波形作为模板波形，总共挑选出满足SNR条件且M_L>2.0模板地震事件112个，扫描时长为主震后3 h。计算模板与连续波形之间互相关系数(correlation coefficient，CC)计算公式为：

$ {\rm{CC}} = \frac{{\sum\limits_{{t_0}}^{{t_1}} {\left( {X\left( t \right) - \bar X{\rm{ }}} \right)\left( {Y\left( t \right) - \bar Y{\rm{ }}} \right)){\rm{ }}} }}{{\sqrt {\sum\limits_{{t_0}}^{{t_1}} {{{\left( {X\left( t \right) - \bar X{\rm{ }}} \right)}^2}\cdot\sum\limits_{{t_0}}^{{t_1}} {{{\left( {Y\left( t \right) - \bar Y{\rm{ }}} \right)}^2}} } } }} $

(1)

式中，t₀、t₁为模板窗的起始时间，X(t)、Y(t)分别为模板波形与连续波形的时间序列，${\bar X} $、${\bar Y} $分别是模板波形与连续波形的时间序列的平均值。

选取两个水平向S波到时前1 s~后3 s、垂直向P波到时前1 s~后3 s，将主震发震时刻作为0点，用4 s窗与连续波形进行滑动扫描，扫描间隔一个采样点，取三分向互相关系数平均值。通过计算序列的绝对离差中位数(median absolute deviation，MAD)检测微小地震。绝对离差中位数表达式为：

$ {\rm{MAD}} = {\rm{median}}(|{X_i}| - \bar X{\rm{ }}){\rm{ }} $

(2)

式中，X_i为第i个互相关系数序列，$ {\bar X}$为其平均值。

在利用模板匹配方法进行遗漏事件检测过程中，阈值的选取需要权衡，阈值太大可能会造成事件遗漏，阈值太小可能会包含太多错误事件，阈值的选取必须通过多次的实验来检验。Peng等^[1]、谭毅培等^[7]都将9倍的MAD值作为阈值，本文参考前人研究结果，取9倍MAD值作为阈值，相关系数大于阈值的事件则被认为是一个地震事件，当相邻两个大于阈值的点过于接近(小于3 s)，判断为同一个事件。图 1(波形模板为Eq01210318门源台(MEY)记录，时间段为2016-01-21 03:00~04:00)展示了模板地震Eq01210318在01-21 03:00~04:00连续波形上进行互相关扫描的结果。计算MAD得到阈值为0.283 1，大于阈值的点共82个，其中余震目录上对应的有53个，剩下的29个点为疑似地震事件。图 2为模板匹配实际检测地震示意图，模板匹配检测地震判定依据为至少需要3个台站平均CC参数都大于0.3，图中3个CC值为3个台站各分向CC值的平均值，左边黑色波形部分为目录事件，红色虚线波形为模板波形，模板与蓝色连续波形重合之处即为检测遗漏事件。

在得到疑似地震事件之后，使用原始波形(每s 100个采样点)和模板波形进行互相关系数扫描，确认疑似地震，搜索遗漏事件P波和S波的到时。图 3(蓝色实线为连续波形，红色实线分别为P、S波震相模板)给出了离主震最近台门源台(MEY)波形互相关提取震相到时示意图，两分向均有明显的CC峰值，其中垂直向CC_max为0.72，水平分向的东南向CC_max为0.91，最大CC的点即为P波、S波的到时，本文遗漏地震震相到时与单台震相到时补全均使用该方法。垂直和两个水平向波形分别对应检测P波和S波。鉴于区域测震台网手动拾取震相到时存在一定误差，为截取比较完整的P波、S波波形，在连续波形垂直向截取P波到时前4 s~到时后4 s，模板垂直向截取P波到时前2 s~到时后2 s波形互相关扫描; 在连续波形两个水平方向分别截取S波到时前4 s~到时后4 s，模板水平向截取S波到时前2 s至到时后2 s波形进行互相关扫描。扫描得到的互相关序列中互相关系数最大值即为遗漏事件P波、S波到时。

1.2 震级计算方法

Peng等^[8]通过研究日本多个地震认为，以1.5级为界包络振幅对应的震级与台网测量震级间存在两组线性关系。谭毅培等^[9]研究甘肃门源M_S6.6地震的时候给出了包络差峰值振幅与震级的线性拟合关系式。本文针对门源地震余震震级特性，余震震级偏小，M_L4级以上余震个数为1，因此本文设定包络差峰值振幅与震级的线性关系式为：

$ M = a{\rm{lg}}({\rm{Env}}) + b $

(3)

式中，M为震级，Env为高频包络差峰值振幅。本文一共选取了主震后01-23全天130个余震事件，震级下限取M_L0.1对包络差峰值振幅与目录震级进行了拟合计算。计算结果为:a=0.89, b=1.54，代入式(3)得到:M=0.89lg(Env)+1.54。

1.3 检测结果与分析

本文选取主震后M_L2.0以上余震112个作为模板事件，通过匹配滤波的方式扫描出遗漏地震123个，并利用波形互相关提取震相到时的方法对单台事件补充其他台站震相到时，对单台事件重新定位，最后得出41个单台事件的震中位置，重新定位后的余震空间分布见图 4(a)。震级计算依据震级与包络差峰值振幅的线性关系，计算结果显示，遗漏地震震级偏低，震级范围为M_L-0.4~1.6，且大部分集中在M_L0.0~1.0区域，如图 4(b)所示。根据模板匹配结果与震级估算，从时间域来分析，遗漏地震频次随时间推移有逐渐衰减的现象，主震后短时间内遗漏地震最多，原因之一是通常主震应力释放能量最大引起断层移动导致震源附近地下较脆结构发生级联破裂，短时间内余震频次最大; 其次是受主震尾波干扰、波形叠加影响，使人工分析地震难度加大导致。

为确保检测质量，通过波形对比及经验分析对123个检测地震事件再研究，结果发现，遗漏地震波形直观辨认上有一定程度的包络特征，但震相到时不清晰，不符合台网中心地震编目震相标注要求——高于背景信号3倍以上。同时发现，遗漏地震事件可以大体分为两类：一类属于无波形叠加情况下震级太小、能量太弱、波形特征不明显或是噪声干扰较大导致，如图 5(a)为检测事件Eq01210135，此次事件之后20 s左右为目录事件; 另一类是波形叠加干扰情况，湮没在震级较大地震事件的尾波中或多个地震叠加一起难以分辨，如图 5(b)为检测事件Eq01210150，此次事件与目录事件尾波叠加，因此无法分辨准确到时。

2 地震目录完备性与b值分析

完备震级(magnitude of completeness, M_C)是国际上普遍采用的评估地震台网检测能力的一个定量标准，同时也是进行地震活动性研究的重要参数。地震目录的完备震级主要取决于台站分布，同时还受到处理方法的影响。Schaff等^[10]研究表明，利用模板匹配方法能有效降低目录的完备震级。为进一步研究地震活动性特征，本文首先在检测事件的基础上分析获得检测目录的完备震级，采取的方法是最大曲率法^[11]，该方法认为震级-频度曲线一阶导数最大值对应的震级为完备震级。利用最大曲率法分析了主震后3 h内青海台网目录和检测目录，得到的完备震级由台网目录的M_L1.1降低到检测目录的M_L0.6(图 6(a))，这说明模板匹配方法能有效减轻由于强震尾波干扰而造成的目录缺失问题，提高地震检测能力。

地震数量与震级通常遵循G-R关系LogN=a－bM(N表示大于震级M的地震累积频度，a、b为拟合常数)^[12]。研究表明，b值与材料的不均匀性等因素有关，求解b值的方法包括最小二乘法、最大似然法等。模板匹配检测遗漏事件震级一般是小震级事件，因此对所有地震赋予同样权重的最大似然法更适合本文研究。利用最大似然法计算震后3 h内检测前后b值的相对大小，起始震级取最大拐角震级。从G-R关系图分布上可以得到，检测前后的最大拐角震级分别与检测前后的完备震级M_L0.6和M_L1.1一致(图 6(b)虚线处)，因此起始震级分别取检测前后的完备震级，b值最终结果为检测前b=0.74±0.03，检测后b=0.67±0.02(图 6(b))。检测后地震活动性b值略微减小，结果更佐证了强震后短时间内b值明显降低的普遍现象^[5]。

3 结语

本文使用模板匹配方法对2016-01-21青海门源M_S6.4主震震后3 h的连续波形进行遗漏地震检测，总共检测出目录以外的123个地震事件，约为青海测震台网目录给出余震数量332的40%。针对目录中57次单台地震事件没有定位结果的缺陷，本文利用波形互相关方法对57次单台事件补充了其他台站到时参数，最终对41个单台事件震相到时补全了至少4个台站的数据。根据地震波振幅与震级线性关系对震级进行估算，最后将检测前后余震目录最小完整性震级进行对比，使最小完整性震级M_C由M_L1.1减小到M_L0.6，地震活动性b值由检测前的0.74±0.03减小到0.67±0.02。遗漏地震通常是较小地震，加入检测地震后，余震目录得到进一步完善。

本文工作获得的地震目录与震相到时为门源地震余震序列增加了更多的样本，为其他研究如门源地区地下三维速度结构、余震衰减规律、未来地震活动性等提供更多的资料和参考依据。完备震级的变化说明，区域台网常规产出的地震目录反映的监测能力较实际结果偏低，而使用模板匹配方法能有效减轻由于强震尾波干扰而造成的目录缺失问题，提高地震检测能力，而b值的进一步精确计算对我们判断强震余震活动性以及震后跟踪具有重要意义。