0 引言

2017年8月8日21时19分，四川省阿坝州九寨沟县发生M_S7.0地震，震源深度20km。此次地震是继2008年汶川M_S8.0地震和2013年芦山M_S7.0地震之后在青藏高原东部边缘发生的又一次破坏性地震，地震波及到四川、甘肃、青海、宁夏、陕西等省份，共造成24人死亡，近500人受伤(郭安宁等，2017)。地震发生在东昆仑断裂带东端、构造转换地区、呈帚状分叉的次级走滑树正断裂带上，向南东为雄黄山隆起区，其边界断裂为岷江断裂、塔藏断裂、雪山断裂和虎牙断裂等(袁道阳等，2016)。此次地震为走滑型地震(USGS)，发震断层—树正断裂位于岷江断裂带与塔藏断裂中间，地震时断层向两端破裂扩展，北段向塔藏断裂方向扩展，南段与虎牙断裂贯通。该区历史上地震活动较为活跃，震中附近200km范围内共发生26次5级以上地震(高曙德等，2017)。

当一次强震发生后，人们最关心的问题莫过于震后是否还会有强余震或者更大的破坏性地震发生。对余震序列衰减特征以及早期序列参数特征的研究一直是地震学和地球物理学研究的热点问题(蒋海昆等，2007)。余震序列衰减特征研究是震后趋势判定、强余震预测和发震构造分析等研究的重要组成部分，对于理解主震发震过程、分析地震危险性和抗震救灾等工作具有重要作用(宋金等，2009)。余震序列的衰减特征和激发余震的能力可通过序列的统计参数来表征，由于震源区的孕震条件和物理状态的不同，不同的构造区域、序列类型、主震的震源机制类型、震源区局部构造应力水平及大地热流等地球物理特征都会有不同的序列参数特征(Kagan et al，2010)。因此，深入分析强震的余震序列活动特征具有重要的意义。

余震序列活动特征主要包括序列时间演化、空间展布及强度分布特征等3个方面。其中，序列时间演化，一方面可基于b值、h值、k值等序列活动参数判定早期发展趋势(刘正荣等，1986；刘正荣，1984；朱传镇等，1989；吴开统等，1996；马玉虎等，2012；宋春燕等，2014)；另一方面可采用具有一定物理背景的统计模型，如修正的大森公式(Utzu，1961、1965；谭毅培等，2015a、2015b)、ETAS模型(Ogata，1988、1989、2001；蒋海昆等，2007；蒋长胜等，2013)以及BASS模型(Sheherbakov et al，2004；Turcotte et al，2007)等来判定。

地震目录是研究余震序列衰减的基础资料，目录的完整性直接关系着研究结果的科学性和可信性。强震发生后短时间内往往伴随着大量的余震事件，其地震波形受主震和强余震的尾波干扰，波形相互叠加，一些地震事件的震相难以清晰识别，因此，测震台网的地震目录往往会遗漏大量的地震事件。研究显示，主震后短时间内遗漏的余震可达到地震目录给出余震数量的数倍甚至十几倍(Peng et al，2007、2009；Lengliné et al，2012；谭毅培等，2015c)，如此大量的遗漏地震势必会影响余震序列地震目录的完整性，进而影响以地震目录为基础的地震活动性、发震构造和发震机理等研究的可靠性。

基于波形互相关的模板匹配技术近期发展迅速且应用广泛。Gibbons等(2006)证实了模板匹配方法在微震检测方面更有优势；Shelly等(2007)研究了非火山地脉动与低频地震间的关系；Schaff(2008)及Schaff等(2010)进一步证实了在地震监测中模板匹配方法的识别能力；Peng等(2009)完备了2004年Parkfield地震的余震目录并分析了余震的迁移特征；Yang等(2009)确认了2008年Illinois Monut Carmel M_W5.8地震的断层面；Meng等(2013)发现了2003年San Simeon M6.5地震之后Parkfield地区的地震活动性与静剪切应力间的对应关系；谭毅培等(2016)研究了滦县震群的发震构造；马腾飞(2015)研究了2008年汶川地震序列中的“重复地震”；Zhang等(2015)实现了对朝鲜地区多次核爆位置的精确检测。虽然模板匹配技术已成功应用在地震学研究的多个领域，但由于其计算量过于庞大，传统的CPU计算会花费大量的时间，Meng等(2012)利用GPU并行计算实现了模板匹配技术算法加速，并完备了2010年Mayor-Cucapah M_W7.2地震之后Salton Sea地区的地震目录。

本文利用基于GPU加速的模板匹配技术来检测九寨沟M_S7.0地震前后的遗漏地震，再基于补充遗漏地震后的完备地震目录，利用b值和M_c值对地震序列的时间演化特征进行分析。

1 遗漏地震检测

强震序列震后短时间内波形相互交叠，不易标识，容易造成大量余震的缺失。研究表明，西南地区的强震序列在震后5天以内衰减的能量最为显著，此后，激发次级余震的能力明显减弱(郭路杰等，2014)。2017年8月8日九寨沟M_S7.0地震发生之前，仅在8月7日发生1次1.5级地震，前震序列偏少。因此，为了考察完整的地震序列的活动特征，本研究选取了主震前7天和主震后5天的连续波形进行分析计算，所用数据为地震序列周边的九寨沟台(JZG)、若尔盖台(REG)、松潘台(SPA)、平武台(PWU)、舟曲台(ZHQ)和文县台(WXT)等6个固定台站记录的三分量连续波形(图 1)，数据采样率为100Hz。

模板匹配技术是在滑动窗互相关(sliding-window cross-correlation，SCC)检测技术的基础上发展起来的(Yang et al，2009)，是在低信噪比和稀疏台网情况下提取信号的一种有效方法。模板匹配技术是以已知的震相清晰的地震事件作为模板，利用波形相关性，扫描连续波形数据，寻找与模板事件波形相似度高的一段连续波形作为新检测事件，以实现对遗漏地震的自动识别和检测。

利用研究时段内台网记录到的1336个M_L≥0地震事件作为候选模板，其中，最大震级为M_S7.0。对数据进行去均值、去线性趋势等预处理后，采用2~12Hz的Butterworth带通滤波器滤波；利用震相报告中的P、S波到时计算模板事件的信噪比，舍弃信噪比小于5的模板事件，最终选取了满足条件的1033个地震事件作为模板。利用P波、S波到时的前1~5s作为扫描窗长，分别与被检测连续波形的垂直分量和水平分量作互相关计算

$ {\rm{CC = }}\frac{{\sum\nolimits_{{t_0}}^{{t_1}} {\left\{ {\left[ {x\left(t \right) - \bar x} \right] \times \left[ {y\left(t \right) - \bar y} \right]} \right\}} }}{{\sqrt {\sum\nolimits_{{t_0}}^{{t_1}} {\left\{ {{{\left[ {x\left(t \right) - \bar x} \right]}^2} \times {{\left[ {y\left(t \right) - \bar y} \right]}^2}} \right\}} } }} $

(1)

其中，x、y分别为x(t)、y(t)的均值；t₀、t₁为扫描时窗的开始、结束时间；CC为互相关系数(cross-correlation coefficient，CC)。

将各台站垂直、水平分量的CC值所对应的时间分别减去模板对应的P、S波到时与连续波形起始时间之差，将各台站、各分量的互相关系数序列进行偏移后的叠加处理，最后，除以所用分量的有效个数来计算平均CC值，得到平均互相关系数波列。同时，计算CC值波列的绝对离差中位数(median absolute deviation，MAD)，若平均CC值大于其阈值，则判定为1次地震事件。本文采用12倍MAD作为阈值(Peng et al，2009；谭毅培等，2014；侯金欣等，2017)，考虑到采样率为100Hz时，0.01s的移动扫描间距会造成大量的误检测事件，因此，取2s内相关系数最大的地震为检测事件。

基于模板匹配方法的假设，检测事件与模板事件具有相同的走时，由此得到检测事件的时间。检测事件的震级可以通过检测事件与其对应的模板事件在各台站的水平分量上S波到时前2s至后2s内的最大振幅比的平均值求得(图 2)。

通过计算，共检测到九寨沟地震序列2017年8月1~12日间的5887个地震事件，其中，包括1033个模板事件的自检测，因此，共4854个检测事件。中国地震台网中心的地震目录共记录到地震事件3057个，其中，可定位地震事件1460个，单台事件1597个。由此可见，本文的检测事件个数为台网可定位地震个数的3倍多，共检测到1797个遗漏地震事件，该数量远大于可定位地震和单台事件的数量。单台事件为震后仅在震中距最近的1个台站能检测到P、S波震相的地震事件，多出现在主震发生后较短的时间内，是受强震后续震相干扰、余震波形相互重叠等因素的影响，因记录的余震波形信噪比较低而造成的。此次地震序列的单台事件数略大于可定位事件数，约占台网目录的52.2%，与谭毅培对岷县漳县M_S6.6地震序列主震后1h内单台事件所占的比例一致(谭毅培等，2015b)。单台事件震中缺失且震级可靠性较差，而检测事件是基于多个台站检测到的地震事件，对震级的估计更为准确，因此，基于补充了遗漏地震的地震目录对地震序列的活动特征进行分析具有更强的科学性和可靠性。

将检测到的遗漏地震补充到台网目录中(台网目录截至2017年9月30日)，得到更加完整的地震序列目录，这使得8月8日主震之前的前震序列更加丰富完整，而不仅是台网目录中记录的1次地震(图 3)。从前震序列的演化可见，8月1日以来存在着地震活动的平静—增强过程，特别是主震前地震活动短时间内的增强，可能是应力的积累导致震源区发生微小破裂所致，这一阶段可能是深部更大的断层发生破裂之前加速蠕动的临震阶段(牛志仁，1978)。同时，主震后5天内的余震序列目录也更加完整，主要体现在0.0~1.0级震级档的补充，遗漏地震事件大多集中在这一震级区间内(图 3)，而从震后5天以后的台网目录记录地震的分布可见，0.0~1.0级震级区间的地震较震后5天内的数量明显增多，台网目录的完整性较好，因此，仅对震后5天作了遗漏地震检测。从补充后的整个地震序列来看，震级主要集中在1.0~2.0级区间内，地震频次的统计也呈现较好的正态分布。

2 地震目录完备性分析

地震目录是进行地震活动性分析、活动构造研究、地震预测与地震危险性评估的重要基础资料(徐伟进等，2014)，而地震目录的质量则依赖于研究区域地震台网的检测能力。完备震级M_c是国际上普遍采用的评估地震台网检测能力的定量标准，通常被定义为在一个时空范围内，地震能被台网100%监测到的最小震级(Rydelek et al，1989)。假设地震的发生是自相似过程，则观测到的地震累积发生频率应满足G-R定律。在震级较大处，地震发生频率偏离线性G-R定律，其原因可能是因大地震数目太少而引起的随机波动或特征地震现象(黄亦磊等，2016)。而小震级处的偏离，则被解释为台网监测能力的不足造成震级小于完备震级的地震事件可能未被台网监测而在地震台网目录中被遗漏，导致地震目录的不完备，从而影响到基于地震目录的相关研究结果的可靠性(李智超等，2014)。而通过检测遗漏地震来补充地震台网目录，则会使得地震目录更加完备(Schaff，2008)，更有利于提高分析结果的科学性和可靠性。

本文利用最大曲率法MAXC(The Maximum Curvature technique)(Wiemer et al，2000)和拟合优度测试法GFT(The Goodness of Fit Test)(Wiemer et al，2000)分析地震序列的完备震级，其中，最大曲率法选取震级-频率曲线中斜率最大值所对应的震级作为M_c。在实际中，这个震级往往对应非累积震级-频率分布中拥有最多地震数目的震级。而拟合优度法则是通过计算不同假定起始震级与G-R关系间的拟合情况来确定起始震级的方法。

通过最大曲率法对8月1~12日的检测目录和台网目录进行分析(图 4)。由图 4可见，检测目录的频次明显大于台网目录，但非累计频次与震级的关系均显示2种目录在1.2级处记录的地震频次最高，因此，完备震级都是1.2级。利用最大似然法对累计频度与震级的拟合发现，完备震级也相差不大，检测目录为1.2级，台网目录为1.3级。检测目录拟合得到的b值为0 8558，略高于台网目录的0 7682(图 4)。由于拟合优度法对于每一个震级下限都采用最似然估计计算，因此，在处理台网探测地震能力随震级变化较快和有不均匀性的地震目录时更加有优势。在拟合优度测试法的结果中，检测目录的完备震级为1.4级，小于台网目录的1.6级(图 5)。研究表明，主震后短时间内记录的单台事件给出的震级明显偏低(谭毅培等，2005b)。2种目录的完备震级相差不大，可能是因为大量单台事件的加入使得台网目录在小震级区间内的地震数增多，使得完备震级被低估。

3 序列活动特征分析

基于补充了遗漏地震的完整地震目录(截至2017年9月30日)，以天为单位，逐渐增大统计时窗，分别利用拟合优度法和最大似然法计算统计时段的完备震级和b值，分析余震序列的完备震级和b值随时间的演化特征。主震发生后短时间内的完备震级较高，出现3 5级的突跳，这可能与主震尾波的严重干扰有关，随着时间的推移，完备震级缓慢降低并趋于稳定，在1.2级左右波动，且呈现准周期的波动(图 6)，这可能与晚间噪声水平较低有关(侯金欣等，2017)。

研究表明，b值与研究区域的应力状态(Chan et al，2012；易桂喜等，2011)等因素有关，序列的低b值反映了震源区的高应力状态。主震前的低b值反映了震源区应力积累到了较高的水平，伴随着主震前地震的密集增强后而发震，此时可能正处于主震前应力加速释放的临震阶段。主震发生后，随着震源区应力的释放，b值快速升高，8月15日以后逐渐平缓，在0 5左右波动，并仍有缓慢升高的趋势(图 6)，表明余震序列的衰减过程减缓，而该区的背景b值为1.0左右(易桂喜等，2011)，说明余震序列仍会在未来较长一段时间内处于衰减的状态。

4 讨论与结论

本文利用模板匹配技术，结合GPU并行运算加速，识别了2017年8月8日九寨沟M_S7.0地震序列在主震前7天和震后5天的遗漏地震，检测地震事件数量是台网可定位地震事件的3倍多，将余震的完备震级从1.6下降到1.4级，并基于补充了遗漏地震的完整地震目录，分析了余震序列的活动特征。

以九寨沟M_S7.0地震为中心，对半径为100km的范围进行统计，2008年汶川M_S8.0地震后地震活动开始明显增强，尤其是2013年芦山M_S7.0地震以来地震活动进一步增强，一直处于较活跃状态。但2017年以来，该区域地震活动明显减弱，特别是3月之后，地震活动处于低活动背景，这与传统意义上强震前地震活动由增强到平静的过程相符。

从完备目录的震级-时间分布来看，此次地震序列为前震-主震-余震型。前震序列存在平静—增强的过程，特别是主震前的地震活动短时间内增强，明显高于该区域2017年以来每天1次的地震日次均值，这种高于背景地震活动水平的增强可能是应力的积累导致震源区发生微小破裂所致，这一阶段可能是深部更大的断层发生破裂之前加速蠕动的临震阶段。与之对应，主震前序列的b值处于低值水平，表明了震源区处于较高的应力状态，而主震发生后b值随时间的推移快速回升，维持在0 5左右，仍存在缓慢升高的趋势。与之相反，序列的完备震级在主震发生后短时间内达到最高值，这可能与主震尾波的严重干扰有关，随后，完备震级缓慢降低并逐渐稳定在1.2级左右。与该区域1.0的背景b值相比，目前序列的b值仍处于低值水平，从b值缓慢升高的趋势来看，余震序列仍会在未来较长一段时间处于衰减的状态。

本文仅对九寨沟M_S7.0地震序列震前、震后几天的连续波形进行扫描，并基于发震前后几天的完备目录结合b值进行初步分析。对于序列衰减系数和其他地震参数的分析，以及对于前震序列的相关推断的论证，将在后续工作中进一步研究。