0 引言

一般情况下, 浅部大地震之后会产生相当多的余震, 随着时间延续, 余震的发生率呈现幂次性衰减(Omori, 1894；Gutenberg et al, 1941；Aki, 1965), 在空间上, 余震由主震破裂区域向外迁移, 而余震的时空演化及触发机制对研究地震断裂活动性和评估地震灾害至关重要。地震的触发机制可能有多种, 如静态应力的变化(Stein, 1999；Meng et al, 2013)、慢滑移(Hearn et al, 2002；Lohman et al, 2007；Peng et al, 2009；Borghi et al, 2016)、动态应力触发机制(Brodsky et al, 2000)及壳内流体(Shapiro et al, 2006)等。大地震之后, 会产生数量巨大的余震事件, 但大地震的尾波信号往往掩盖了震级较小的地震信号, 而且很多小地震的震相重合, 这使得常规方法得到的地震目录中缺失很多余震事件(Peng et al, 2009)。因此, 提取出更多的余震事件, 获得较完整的地震目录, 对分析大地震后早期余震分布和迁移规律至关重要。

检测地震事件的常用方法有STA/LTA(长短时间窗平均值比)、MFT方法。STA/LTA(Stevenson, 1976)方法主要基于长时间窗信号平均值与短时间窗信号平均值的比值来判断地震事件。但这种方法只能检测出高信噪比事件, 仍会遗漏很多早期余震。MFT(the matched filter technique, 模板匹配滤波)(Gibbons et al, 2006)方法通过模板地震事件在连续波形上滑动求取互相关值来检测地震事件。此方法对信噪比要求较低, 更依赖于事件波形之间的相关性, 相比STA/LTA方法, 能够行之有效地检测出微小的地震事件, 模板匹配滤波技术在微地震检测方面具有明显优势(李璐等, 2017)。基于该性质, MFT方法不仅可用来检测天然地震的前震和余震事件(Peng et al, 2009；Kato et al, 2012；Sugan et al, 2014；谭毅培等, 2014；Wu et al, 2017；Huang et al, 2017；侯金欣等, 2017), 还可被用于非火山颤动和低频地震(Shelly et al, 2007)、油气勘探中水压致裂诱发地震(Bao et al, 2016)、爆炸监测(Zhao et al, 2016)等研究领域。

2016年8月24日1时36分(国际标准时), 意大利中部亚平宁山脉附近发生M_W6.0地震, 震源深度约8km, 地震对当地居民和建筑造成了较大的伤害。2个月之后, 2016年10月26日在M_W6.0地震震中的北部, 2h内相继发生了M_W5.5、M_W5.9地震。紧接着4天以后, 在M_W6.0、M_W5.9地震震中的中间区域又发生了M_W6.5地震, 形成了序列地震。在过去20年里, 该区域曾发生2次M_W≥6.0地震(即1997年Colfiorito M_W6.0、2009年L'Aquila M_W6.3地震)(图 1(a))。

本文利用意大利中部密集台网的地震记录, 采用模板匹配滤波方法检测2016年8月24日M_W6.0地震发生后80天内的余震事件, 研究此地震序列的余震时空分布规律, 探索余震迁移和地震触发机制。

1 地质背景

亚平宁山脉中部地区位于意大利中部, 处于阿德里亚微板块和欧洲板块的挤压和俯冲作用带上(图 1(b)), 是意大利境内活跃的地震活动带之一。前人对该研究区的地质研究显示, 亚平宁中部地区现今主要的地震活动是由一系列NW-SE走向的正断层构成的断裂系引起的, 该断裂系形成的地质时期为第四纪(Boncio et al, 2004；Galli et al, 2008；Pizzi et al, 2009；Lavecchia et al, 2016), 大多数都是倾向SW, 以雁列式分布在山脉附近。主要的断裂有Colforito断裂(CF)、Bove Faul断裂(BF)、Vettore断裂(VF)、Norcia断裂(NF)和Gorzano断裂(GF)等(图 1(a))。该断裂系又被2组先存构造带分割(图 1(a)中的OAST和MCT), 这2组逆冲断裂带形成于第四纪之前的挤压构造活动(Pizzi et al, 2009)。

亚平宁中部地区以中强度的正断层地震活动为主(Chiarabba et al, 2005), 图 1(a)中几次大地震的震源机制解均显示了拉张为主的活动性质(Chiaraluce et al, 2005；http://cnt.rm.ingv.it), 主张应力轴呈NE-SW向(图 1(a))。研究区的大地测量结果显示为向山脉两侧扩张, 形变速率约3mm/a, 方向主要为NE-SW向(Serpelloni et al, 2005；Petricca et al, 2015), 与震源机制显示的拉张主应力轴方向一致。

此外, 研究区域内的地震活动常常以序列地震的形式发生。1997年、2009年序列地震都包含了数个5级以上的地震(Amato et al, 1998；Chiarabba et al, 2009；Chiaraluce et al, 2011), 1997年主要的地震分布在2016年地震的NNW方向, 2009年L'Aquila地震分布在SSE方向, 且震源机制解均显示以拉张型的断裂活动为主(Ekström et al, 1998; Chiaraluce et al, 2011)。2016年8月底至10月底的地震活动主要发生在1997年、2009年地震空区处(Liu et al, 2017；Chiaraluce et al, 2017), 余震主要分布在BF、VF与NF两组平行的断裂带附近。

2 数据与方法 2.1 地震数据来源

本文使用的地震数据是由欧洲台网中心(http://www.seismicportal.eu)获取的三分量宽频带连续波形数据, 分别选自意大利CAMP和AQU等13个台站(图 1(c)), 连续波形采样率为100Hz, 记录时间从2016年8月24日0时~11月12日0时(UTC时间), 总共80天。

地震目录来源于国际台网中心ISC(http://www.isc.ac.uk), 我们共获取2016年8月24日~11月1日发生的1426个地震事件作为后续匹配滤波扫描的模板。

2.2 地震模板截取

在模板地震事件截取之前, 首先对地震数据进行了预处理, 为了削弱远震事件的影响, 采取滤波范围1~8Hz, 连续波形数据重采样至25Hz, 模板数据从预处理后的连续波形数据中截取。

对于微小地震来说, 由于S波能量强于P波, 更容易从噪声中检测出来, 因此模板选自S波段。根据地震目录, 选取S波理论到时的前2s、后6s为时间窗截取模板信号, 选用P波到时前12~4s作为背景噪声信号计算模板的信噪比, 取信噪比大于5作为初选模板。为保证最终模板的可信度, 需足够多的模板数据用于叠加, 因此要求选取的每个模板数据有4个以上台站的记录(即至少12个分量的记录数据)。其中, P、S波理论到时由TAUP软件(Crotwell et al, 1999)计算得到, 速度模型采用当地CIA模型(Herrmann et al, 2011)。最终, 选出了1381个模板事件(图 1(c))。

2.3 互相关计算

第1步, 从每天记录的零时开始, 计算某一分量模板与相同分量、时长连续波形数据的互相关系数(cross-correlation coefficient, CC), 然后将模板在连续波形上后移1个采样点, 每移动1次计算1个互相关系数, 对1天的记录得到该分量的互相关系数序列。第2步, 将不同台站同一模板事件的S波按照到时对齐, 对不同台站所有分量计算平均互相关值(mean CC)。第3步, 假定每一天的信号序列是符合高斯分布的随机过程, 信号序列的标准差为1.4826倍绝对中位差(median absolute deviation, MAD), 相关值大于12倍MAD(8.1倍标准偏差)的概率为3.3×10^-16。由于1天记录的数据点为2 160 000个, 选取1天的平均互相关值组成序列的12倍MAD作为阈值, 每天的错误检测率小于1个。当序列上某个时刻的平均互相关系数高于阈值, 则认定是检测出来的事件(图 2)。

2.4 检测事件的确认

对于同一时间窗内出现多个模板匹配的情况, 保留平均CC值最高的模板。对于新检测出的事件, 由于事件波形和模板事件相似度最高, 其震源位置选用模板事件的位置；新事件的发震时刻由CC值最高的时刻相对模板S波到时的差值, 再加上模板的发震时刻计算得到；检测事件的震级, 通过与模板在所有分量上最大绝对振幅的比值的中位值计算得到, 依据振幅每增加10倍, 震级增加1级(Peng et al, 2009)。

3 结果与分析 3.1 地震目录完备性分析

利用选取的1381个模板事件, 从2016年8月24日~11月12日总共80天的连续波形中新检测出71525个余震(12倍MAD), 数量是选用模板的数十倍。新检测出的余震震级范围为-0.8~5.3级, 与模板事件的震级范围(0.8~6.5级)相比, 检测的余震事件震级更小(图 3)。对新检测的余震与模板事件作地震频度与震级大小间关系的分析, 采取Wiemer等(2000)提出的GOF(goodness of fit)方法, 当曲线拟合度达到90%时, 得到完备震级M_c, 计算结果如图 4所示。由检测结果和模板数据计算得到的M_c约为1.0, 相比模板目录的M_c(2.5), 有了大幅度的降低。检测的余震事件中震级大于1.0的地震事件共有29035个, 加上1381个模板事件, 组成了新的事件目录, 新目录共包含30416个地震事件, 由此提高了地震目录的完备性。我们利用新的地震目录, 分析意大利中部地区2016年8月24日~11月12日余震的时空分布规律。

3.2 余震分布的总体特征

对于本文研究中2016年意大利中部发生的序列地震, 依据其发震时刻和所在位置, 将其分为Event 1(2016年8月24日M_W6.0地震)、Event 2(2016年10月26日M_W5.5、M_W5.9地震)以及Event 3(2016年10月30日M_W6.5地震), 以进一步分析序列地震中不同事件余震的时空演化过程与差异。

图 5为2016年8月24日~11月12日余震分布, 纵轴为沿测线A-A′(图 1(c))的距离, 图 5一定程度上反映了该序列地震的余震分布规律。

为了更清晰地分析余震的分布特征, 按照4个主事件位置沿测线分成4段进行余震的发震频度统计(即每小时发生地震的个数)(图 5中彩色曲线)：橙色线表示Event 1(测线0位置)的东南段、投影距离为-20~0km(测线A′至Event1之间)的余震频度曲线；黄色线表示沿测线0~17km(Event1至Event3之间)的频度曲线；红色线为沿测线17~25km(Event3至Event2-2之间)的频度曲线；绿色线显示25km西北侧的频度曲线。由图 5可见, Event 1的余震事件主要发生在测线-15~20km的位置, 主震东南侧较快地趋于平静(图 5橙色频度曲线), 而西北侧的余震迁移至Event 2震中附近, 并且持续到Event 2之前仍然较活跃(图 5蓝色点分布)；Event 2(包括Event 2-1和Event 2-2)的余震事件主要分布在测线15~45km之间, 余震活动以Event 2-1西北段的为主(图 5红色、绿色频度曲线)；Event 3的余震事件分布范围较广, 大致为测线-10~45km的范围内, 其东南侧的活动频率明显高于西北侧(图 5黄色、红色频度曲线)。

3.3 早期余震的空间分布

为了进一步分析序列地震中几次主要地震的早期余震空间分布特征, 分别从新余震事件目录中提取Event 1、Event 2、Event 3发生后3天内的余震, 绘制平面分布图(图 6)。由图 6可见, Event 1发生后10⁴ s(~3h)内余震(图 6(a)蓝色点)主要集中在NF断裂东侧, 并沿断裂延伸方向朝主震的两侧扩散, 呈现沿断层走向向南北两侧对称发生迁移的特征；在10⁴s之后的时间内(图 6(a)红色点), 余震向南扩散的范围未见明显变化, 而向北的范围明显扩大；余震活动主要集中分布在NF、GF、VF断裂之间(图 6(a))。

Event 2的早期余震分布特征(图 6(b))表现为：Event 2-1发生后约10³s(~0.5h)内余震主要集中在震中附近, 随后余震沿断裂走向逐渐向北迁移至Event 2-2附近, 2h后发生M_W5.9地震；在10⁴ s(~3h)左右逐渐向东北部转移到BF断裂附近(图 6(b)蓝色点), 此时向南迁移的余震较少；相反, 10⁴ s之后(图 6(b)黄、红色点)余震活动呈现向南迁移的特征, NF、VF断裂之间逐渐被余震覆盖, 但基本上没有超过OAST断裂, 表明受到先存断裂OAST的限制(Pizzi et al, 2009)。Event 2的余震活动向南迁移的过程明显较向北的迁移更晚, 呈现“南北不对称”的余震迁移规律。

Event 3的余震分布范围更广。在10⁴s(~ 3h)以内, 余震主要沿着断裂走向朝南迁移, 10⁴s之后, 余震在断裂北端活动加强(图 6(c)红色点), 余震迁移表现出“南北不对称”的规律。Event 3的大部分余震在平面上分布于NF、BF、VF断裂附近, 仅少部分分布在GF断裂附近, 总体上以分布在OAST断裂北部为主, 西北部基本上没有超过MCT断裂, 表明先存的MCT逆冲断裂对VF、BF断裂的活动起到一定限制作用(Pizzi et al, 2009)。

3.4 余震活动性与迁移速率

我们通过沿A-A′方向余震随对数时间的分布来估算余震向周围迁移的速率(图 7), 以进一步分析余震的迁移活动规律。鉴于M_W6.0地震的余震在2个月后已经相当微弱(图 5), 我们把Event 2-1之后时段分离出来分析Event 2、Event 3余震活动特点与差异(图 7(b))。

在8月24日M_W6.0地震之后1天内, 余震迁移前端沿断裂带向南、北以2~3km/decade速率扩散到15km附近, 之后向南部的扩展逐渐减弱, 相比之下, 北端余震仍不断扩展, 约1周后出现迁移速率明显加快的现象(图 7(a))。以沿投影距离每2km余震数目达30个以上来计算余震活动曲线(Kato et al, 2014)(图 7红色虚线)可更细致反映余震在迁移过程中的表现, 并显示不同时段、不同位置余震的活跃程度。由图 7可见, 在Event 2、Event 3发震前出现了几处迁移速率明显加快的现象, 余震活动曲线呈阶梯状。Event 1发生后的1周内, 余震活动较平稳地向北、向南扩散；1周后, 北部活动曲线明显变化, 迁移速率加大, 显示向北的活动出现了加快的现象, 扩散到Event 2-1震中附近, 而向南的余震活动在测线-15 km附近不再向南扩展(图 7(a))。在Event 2-1发生后1~2 h之间, 余震活动曲线显示向北的迁移速率有小幅度加快, 向Even 2-2震中迁移, 2h后发生了Even 2-2。Even 2-2发生后, 余震活动曲线上向南迁移的余震出现迁移略有加快的现象, 达到Event 3震中附近。这可能反映了M_W6.0地震的余震迁移在Event 2、Event 3震中附近造成了应力加载, 促使了其发生。

3.5 对余震成因机制的讨论

余震的触发、迁移机制可能受多种因素共同影响, 如静态应力、动态应力触发、慢滑动和流体迁移等。Liu等(2017)、Convertito等(2017)通过计算静态库仑应力变化认为, 2016年意大利中部M_W6.0地震的余震活动主要由同震应力加载所致；Vuan等(2017)拟合余震迁移前端, 计算不同深度地壳流体扩散因子D(Shapiro et al, 2006)时发现, M_W6.0地震的余震活动反映的D值(浅部D值为25~50m²/s, 深部D值为160m²/s)过高, 排除了粘弹性松弛原因, 并认为可能是慢滑动所致。我们通过拟合余震迁移前端发现, M_W6.0地震早期余震的两侧迁移较符合lgt的特征, 表明可能与慢滑动的迁移有关。除此以外, M_W5.5地震(2016年10月26日17时)之后朝南侧的迁移也符合lgt特征, 说明早期余震迁移活动与慢滑动的迁移有关。我们推测慢滑动迁移可能引起了相邻区域的应力加载, 对后续M_W5.9、M_W6.5地震的发生起到一定程度的触发作用。但要严格论证慢滑动造成的应力加载, 还需进一步的研究工作。

4 结论

我们利用模板匹配滤波方法(MFT)有效地检测出意大利中部2016年8月24日序列地震的余震, 获得了数十倍于原模板数量的地震检测结果, 得到检测余震的完备震级M_c由2.5降至1.0, 提高了地震目录的完备性。

对余震时空分布特征的分析显示, 序列地震中不同事件余震的扩散方式和迁移速率存在一定差异。M_W6.0地震的早期余震分布呈现向南北两侧对称扩散的特征, 迁移速率2~3km/decade。M_W5.5~5.9、M_W6.5地震的早期余震分布表现出“南北不对称”的余震迁移现象：M_W5.5~5.9地震的早期余震(3h内)主要向北迁移, 然后再向南部扩散；而M_W6.5地震早期余震表现出先向南、后向北迁移的现象, 且整体的活动范围更广。

M_W6.0地震向北迁移的余震, 在主震发生1周后迁移加快, 可能反映了M_W6.0地震的余震活动造成后续大地震震中附近的应力加载。M_W6.0地震向两侧迁移的早期余震前端和10月26日M_W5.5地震向南迁移的早期余震前端符合lgt的特征, 说明早期余震迁移活动与慢滑动的迁移有关。