毕金孟, 蒋长胜, 来贵娟, 宋程

1 天津市地震局， 天津 300201 2 中国地震局地球物理研究所, 北京 100081

0 引言

强震发生后随时发生的强余震严重威胁震区的经济社会和人民生命财产安全，及时发布权威可靠的强余震预测信息，对抗震救灾、应急管理决策以及灾后恢复重建等工作意义重大(Reasenberg and Jones, 1989; Marzocchi and Lombardi, 2009; Woessner et al., 2011; Nanjo et al., 2012; Ogata et al., 2013).强余震发生的时间跨度可从几秒到几年，但大多数的余震是在主震过后数天内发生的(JMA, 2009).主震过后已受损的建筑物来不及修复的情况下便再次遭受余震的继续作用，产生累计损伤效应、造成更为严重的人员伤亡和财产损失.例如2010年9月3日的新西兰7.1级地震之后的6.3级强余震造成了146人死亡、300多人失踪及建筑物的破坏(张晁军等, 2011)，1976年唐山7.8级大地震之后的滦县7.1级强余震、宁河6.9级强余震造成余震区域的房屋几乎全部倒塌及桥面铁路损毁等(吕晓健等, 2007)，2003年日本北海道8.0级地震之后的5.0级强余震造成了储油罐溢出引发火灾的次生灾害(赵金宝, 2005)，以及其他在余震作用下造成的人员财产损失和建构筑物的破坏.关于强余震造成的累计震害损失问题已在巨灾保险的灾害模型构建中引起高度关注(熊政辉, 2019).因此，作为地震预测中“可操作性”较强并可能取得明显减灾实效的早期强余震概率预测，尤其是震后3天内的强余震概率研究工作，将更为有效的提升防震减灾和公共服务水平.

对强余震的预测建模和效能评估始终是地震减灾领域的热点科学问题.例如Sornette等(1996)从物理学引入级序分析(rank-ordering analysis)方法用于预测余震序列中的最大事件，在统计地震学上基于G-R关系的b值截距法(吴开统等, 1984)和最大余震震级法(Shcherbakov and Turcotte, 2004)，以及进行概率预测的R-J模型(Reasenberg and Jones, 1989)和ETAS模型(Ogata, 1989)都在国内外得到广泛应用.Steacy等(2014)从物理模型与统计模型相结合的角度，将 Coulomb 模型和 STEP 模型结合研究了Canterbury 地震序列的余震发育情况.从聚焦到震后早期的强余震预测角度，Shcherbakov等(2018)利用地震序列中的早期事件(前震或早期余震)以及震级-频度分布和地震发生率分布，构建了极端事件震级大小的贝叶斯预测分布，可用于考虑不确定度情况下的主震或强余震的震级预测.Gentili和Di Giovambattista(2017, 2020)利用震后早期余震次数和初期几小时/几天辐射能量时空演变的统计特征发展起来的一种模式识别技术对意大利东北部和斯洛文尼亚西部的破坏性之后的强余震进行了预测，成功地检验了1976年的高破坏性震群以及2019年的三个小震群.目前国际上正在进行的全球“地震可预测性合作研究”(Collaboratory for the Study of Earthquake Predictability, CSEP)计划中，短期余震概率预测模型也得到快速发展(Schorlemmer et al., 2018).美国地质调查局(USGS)和全球地震模型(GEM)项目开展的“可操作的余震预测”(Operational Aftershock Forecasting，OAF)(Helmstetter et al., 2006; Console et al., 2010)也对规范震后早期强余震预测提供了重要的科学参考.

我国大陆地区地震强度大、频度高，当前正处于人口和财富(GPD)向大城市和城市群快速集中过程中，防范化解包括大震巨灾在内的重大风险是经济社会发展的重要保障条件，切实提升强余震的预测能力和风险处置决策的科学水平的需求强烈.但长期以来偏向于定性的强余震预测，难于满足震后地震灾区现场应急救援、疏散和灾后恢复重建，以及巨灾保险保障等客观需求.针对上述问题，本文选用当前国际上较为前沿的、可充分利用不完整地震记录技术优势的Omi-R-J模型(Omi et al., 2013, 2015, 2016)，对我国大陆地区强震开展早期余震的概率预测并量化评价预测效能，获得的余震预测模型的适用性、预测限度和制约因素，为今后开展“可操作”的强余震预测提供参考，服务于实际的地震防灾救灾工作.

1 研究所用方法和资料

Reasenberg和Jones(1989)以余震衰减关系的Omori-Utsu公式(Omori, 1894; Utsu, 1961)作为频次约束，以震级频度分布的G-R关系(Gutenberg and Richter, 1944)作为强度约束的基础上发展了余震短期概率预测的Reseanberg-Jones(R-J)模型.根据R-J模型，在地震序列中t时刻震级不低于M的余震强度函数可表示为：

(1)

式中的t为主震发生后的时刻，一般取正值.参数k控制余震活跃程度，在很大程度上取决于单个余震序列，与主震震级关系不大.参数p表示地震序列的衰减程度，p值的大小与序列的衰减快慢呈正相关关系.参数c用于调节主震发生后极短时间内余震记录的不完整性，是一个数值极小的正常数，与震源深度成负相关关系(Shebalin and Narteau, 2017).参数b表示了应力累积水平、大小上呈反比关系(Wiemer and Katsumata, 1999; Enescu et al., 2011).

Omi等(2013)在R-J模型的基础上，将地震序列早期阶段完整性震级以下的余震考虑到模型参数拟合和余震发生率预测中，提出了一种新的方法，简称为“Omi-R-J”模型.Omi等(2013)利用Ogata和Katsura(1993)给出的检测率函数q(M)的表达方式(OK1993模型)，来对地震事件记录不完整部分的检测程度进行描述.而实际记录到的地震概率密度函数可表示为：

=βe-β(M-μ)+β2σ2/2q(M|μ,σ),(2)

式中的β为bln10，μ表示检测率为50%时对应的震级、σ为相应的震级离散度，通常用μ+2σ或者μ+3σ来近似最小完整性震级Mc.余震记录的检测率会随时间动态变化，在地震序列早期阶段由于大量余震事件缺失、往往具有较高的μ值，此后逐步恢复至正常水平.

在对公式(2)的参数估计中，采用了Omi等(2013)发展的“状态-空间”模型(Omi et al., 2013)来估计随时间变化的μ(t).具体思路为，设定μ(t)为余震发震时刻序列ti≤t≤ti+1(i=1,2,…,n)对应离散的分布函数，假设一个对μ(t)平滑约束的先验分布并设定超参数V来控制其平滑程度，利用最大期望(EM)的迭代算法对参数β、σ、V进行优化和最大后验估计后，通过贝叶斯估计来获得参数μ=(μ1,μ2,…,μn)T.此后进一步结合Omori-Utsu公式和最大似然法来确定参数p，c，k及标准差.

根据参数估计结果，利用公式(1)即可计算在震级范围M>Mp和任意时间间隔[t2,te]内的余震预测数目：

(3)

为对中国大陆强震的余震序列展开研究，使用了中国地震台网中心提供的《全国统一正式编目》地震目录(1)1)全国地震编目系统.http:∥10.5.160.18/console/exit.action.[2020-10-31]..筛选出225例6.0级及以上地震事件，使用Gardner-Knopoff方法(Gardner and Knopoff, 1974)对其中的丛集地震进行删除，保留独立的主震事件134例.采用毕金孟和蒋长胜(2019)给出基于纬度-时间图、经度-时间图以及震中分布图相结合的“自然边界法”来选取计算所用的地震序列目录.为保证参数拟合的稳定性和质量，利用“震级-序号”法(蒋长胜和吴忠良, 2011)确定每个序列的最小完整性震级，保证完整性震级以上的地震数目不少于40个，由此共选取了86个地震序列开展相关的研究.相应的主震参数如表1所示，主震的震中空间分布如图1所示.

图1 中国大陆1970年以来6.0级以上地震的空间分布Fig.1 Epicentre distribution of earthquakes (M≥6.0) since 1970 in continental China

Omi-R-J模型利用对每个地震事件逐步迭代的方式，可获得研究时段内的地震检测率随时间的变化.图2给出了2020年1月19日新疆伽师MS6.4地震序列按照Ogata和Katsura(1993)模型的50%地震检测率μ(t)的分布，为考察拟合过程的稳定性，分别计算了主震发生后的0～0.05天、0～1.00天、0～2.00天、0～3.00天和0～30.00天，共计5个时段的检测率分布情况.

图2 2020年1月19日新疆伽师MS6.4地震序列的完整性分析图中各颜色的曲线分别标出了不同时段(震后0～0.05天、0～1.00天、0～2.00天、0～3.00天和0～30.00天)所对应的50%地震检测率的结果，灰色圆点代表了新疆伽师MS6.4地震事件的余震.Fig.2 Catalogue completeness analysis of the Xinjiang Jiashi MS6.4 earthquake sequence on January 19, 2020The colored curves indicate the results of the 50% detection rate calculated from data at different time periods (0～0.05 day, 0～1.00 day, 0～2.00 days, 0～3.00 days and 0～30.00 days) after the mainshock. The grey dots represent the aftershocks of Xinjiang Jiashi MS6.4 earthquake.

2 地震序列参数的拟合情况

震后快速获得较为稳定的模型序列参数，是震后短时间内提升余震预测能力关键(蒋长胜等, 2017, 2018; 毕金孟和蒋长胜, 2017)，而震后早期小震波形的“淹没”所导致的余震序列缺失，可能是影响模型参数稳定性的主要因素之一.震后早期阶段余震大量缺失、传统的严重依赖于完备性震级的ETAS、R-J模型难以给出准确的余震概率预测结果，而Omi-R-J模型充分利用了震后全部记录的余震事件，并在震后短期内(0.05天)快速开展参数拟合，以2020年1月19日新疆伽师MS6.4地震序列为例，选取了0～0.05天、0～1.00天、0～3.00天、0～30.00天四个早期时间段进行了说明.震后早期序列持续时间t2=0.05天时的模型参数为β=1.733±0.292，k=0.011±0.075，p=1.020±0.140，c=0.016±0.065；序列持续时间t2=1.00天时的模型参数为β=1.729±0.106，k=0.023±0.011，p=0.899±0.119，c=0.094±0.049；序列持续时间t2=3.00天时的模型参数为β=1.859±0.083，k=0.012±0.006，p=0.853±0.0.072，c=0.071±0.037，序列持续时间t2=30.00天时的模型参数为β=1.853±0.054，k=0.015±0.005，p=1.024±0.0.033，c=0.192±0.047.

表1 中国大陆部分强震序列的Omi-R-J模型预测的N-test检验结果Table 1 A list of N-test parameters for some strong earthquake sequences in continental China based on the forecast results of Omi-R-J model

续表1

续表1

表2 中国大陆地区地震序列Omi-R-J模型参数拟合的整体统计特征Table 2 The overall statistical characteristics of the fitting parameters for the earthquake sequences in continental China based on the Omi-R-J model

图3 基于震后1天和30天数据的86个地震序列的Omi-R-J模型参数对比图Fig.3 Comparison of aftershock sequence parameters estimated from the first 1-day span data against those from the first 30-days span data after the 86 mainshocks based on the Omi-R-J model

3 强余震概率预测效能评估

为更好的验证Omi-R-J模型预测的地震数目与实际发生的地震数目的偏离程度，选用了CSEP计划中针对地震数检验的N-test方法(Kagan and Jackson, 1995; Schorlemmer et al., 2007; Zechar, 2010)，并选用有效显著性水平αeff=0.025来对N-test检验结果进行单边检验，当δ1<αeff时，表示存在预测余震数目“过少”的情况；当δ2<αeff时则表明存在预测余震数目“过多”的情况.作为实例，图4和图5给出了Omi-R-J模型基于不同时间段的数据、不同时间尺度、不同震级档分别在震级-频次、时间-频次上的对比结果，可以看出，分级分段分档的实际发生地震数均在预测数的95%置信区间内，显示了较好的预测效能.

图4 利用Omi-R-J模型对2020年1月19日新疆伽师MS6.4地震序列在震级-频次上的强余震概率预测结果示例图中分别是基于震后[0.05，1.00，3.00]天的数据预测未来1天(a、d和g)、3天(b、e和h)和7天(c、f和i)的强余震发生次数，其中黑色圆点为实际观测结果，红色线段为预测结果，粉色区域为预测结果的95%置信区间.Fig.4 Strong aftershock probability forecasting of magnitude-cumulative number of aftershocks result for the Xinjiang Jiashi MS6.4 earthquake sequence using the Omi-R-J model on January 19, 2020Based on the data of [0.05, 1.00, 3.00] days after the mainshock, the number of strong aftershocks in the next 1 day (a, d and g), 3 days (b, e and h) and 7 days (c, f and i), respectively. The black dots represent the actual observed aftershocks, the red lines represent the forecasting result, and the pink areas represent the 95% confidence interval of the forecasting result.

图5 利用Omi-R-J模型对2020年1月19日新疆伽师MS6.4地震序列在时间-频次上的强余震概率预测示例图中分别是基于震后1天的数据预测未来1、3、7天在震级档M-3.0(a、d和g)、M-2.5(b、e和h)和M-2.0(c、f和i)的强余震发生次数，其中黑色曲线为实际观测结果，红色曲线为预测结果，红色虚线标出了预测结果的95%置信区间，黑色垂直虚线标出了预测起始时刻.Fig.5 Strong aftershock probability forecasting of time-cumulative number of aftershocks for the Xinjiang Jiashi MS6.4 earthquake sequence using the Omi-R-J model on January 19, 2020Based on the 1-day after earthquake, the number of strong aftershock forecasting with different magnitude M-3.0 (a, d and g), M-2.5 (b, e and h) and M-2.0 (c, f and i) in the next 1 day, 3 days, 7 days, respectively. The black curves represent the actual observation results, the red curves represent the forecasting results, the red dotted lines indicate the 95% confidence interval of the forecasting results,and the black vertical dotted lines indicate the forecasting start time.

为了对1970年以来中国大陆地震序列的早期余震(震后3.00天内)的预测效能进行系统性评估，因每个地震序列受到区域地震监测能力等因素的制约，Omi-R-J模型需要一定的地震数目才可以进行强余震预测，在0.05～3.00天的时间内尽可能的去“逼近”发震时刻，来进行分级分段分档的概率预测，并对预测结果进行了“量化”的效能评估.因实际未发生地震不能进行N-test检验，为了更加客观的进行评价，因此本文将预测有一定概率的地震发生，而实际没发生地震的情况，并未考虑在统计范围内.

面对复杂的地震序列，即使利用充分考虑小震事件的、公认的预测效果较好的Omi-R-J模型仍出现一些预测“失效”的情况.在所有参与计算的时间段内，其预测效能评估结果如图6所示，从统计结果来看，预测过多占优的序列(39/86)大于预测过少占优的序列(19/86)；预测过少的比例大于预测过多的比例.从所有序列预测时段整体上来看，对于Omi-R-J模型，δ1<0.025的391次，即存在预测余震数目“过少”比例占12.14%；δ2<0.025的216次，即存在预测余震数目“过多”的比例占6.70%；总的预测失效的次数为607，预测失效的比例为18.84%，所有序列预测评估结果如表1所示.平均预测有效率指的是满足不同条件下序列有效率的平均值，来探究不同情况下序列的整体预测效能.对所选单个序列总的平均预测有效率为83.82%；对所选序列1天、3天和7天的预测平均有效率分别为88.90%、82.76%、80.82%；对不同震级档[M-3.0,M-2.5,M-2.0]预测的平均有效率为78.89%、86.14%、92.35%；基于震后早期不同时段([0.05, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0])预测的平均有效率为84.31%、92.06%、91.12%、90.74%、92.49%、91.87%、92.90%，具体整体的分级分段分档结果如表3所示.由于预测次数的差异，整体的预测失效情况和单个序列预测平均失效结果存在稍小的差异.

4 讨论

快速获得震后较为稳定的序列参数，是提高强余震概率预测有效性的重要因素，Omi-R-J模型震后初期获得的b、c、k值与稳定时段的序列参数具有较小的差异，这与利用匹配滤波技术拾取被常规地震台网检测遗漏的余震事件(Peng et al., 2006)、Zhuang等(2017)利用余震序列重新构建技术给出的地震序列参数稳定性所获得的效果一致.为进一步讨论震后早期预测失效的制约因素，统计了强震震级、序列参数和震后早期强余震数目(M-3)之间的关系，如图7所示.由图可知，震级的大小与震后早期的序列参数k、强余震数目(M-3)之间没有明显的相关性，而序列参数k与余震数目成一定的正相关关系，其控制着余震的数目.即使是在相同的构造条件和相同的震级条件下，其参数也存在一定的差异，这在很大程度上取决于单个余震序列.包括日本对2004年Chuetsu 地震和2007 年的Chuetsu-oki 地震研究也获得了类似的结论(JMA, 2009)，地震序列参数k是影响强余震数目的关键因素.

表3 中国大陆地区地震序列Omi-R-J模型效能评估统计结果Table 3 Statistical results of effectiveness evaluation of the Omi-R-J model in continental China

图6 Omi-R-J模型强余震发生率预测的N-test统计检验结果(a)、(c)、(e)分别为Omi-R-J模型在不同震级档下对未来1天、3天、7天的δ1评分结果；(b)、(d)、(f)分别为Omi-R-J模型在不同震级档下对未来1天、3天、7天的δ2评分结果.图中蓝色色调代表在95%置信水平预测有效的时段，红色色调代表δ1<0.025或δ2<0.025的结果，即预测失效的时段.空白区域为预测时段没有发生一次预测震级以上的地震事件，斜线区域代表因数据不足原因无法参与计算的时段.Fig.6 N-test results for the strong aftershock forecasting based on the Omi-R-J model(a), (c) and (e) show δ1-score for future 1 day, 3 days and 7 days aftershock forecasting in different magnitudes, respectively. (b), (d) and (f) show δ2-score for future 1 day, 3 days and 7 days aftershock forecasting in different magnitudes, respectively. The blue block represents the period when the prediction is valid at 95% confidence level. The red block represents the result of δ1<0.025 or δ2<0.025, i.e. the forecasting period of failure. The blank area indicates that no earthquake events of or above forecasting magnitude have occurred in the forecasting period. And the slash area represents the time period that cannot participate in the calculation because of insufficient data.

图7 地震序列参数与震后早期强余震数目的关系颜色的深浅代表了主震震级的大小.Fig.7 The relationship between the sequence parameters and the number of strong aftershocks in the early stage after the mainshockThe color represents the magnitude of the mainshock.

在基于不同时间尺度、不同震级档的条件下，随着预测时长和震级档的下调，预测有效率反而降低，这和之前的一些研究(蒋长胜等，2015; 毕金孟和蒋长胜，2017)相反，主要可能是因为预测时段次数的增加导致了失效率的增加.除基于0.05天的预测有效率低于80%外，其他预测时段的有效率均大于80%，从0.05天到0.5天预测失效率由30.96%下降到16.23%，随着更多余震事件参与计算，地震的预测有效率得到显著提升，早期有效率的偏低可能与参与计算的余震数目过少，即序列的发育程度有关.具体的基于不同时段、不同震级档以及不同时段的预测失效次数和失效比例如表3所示.监测能力的显著提升是提高余震概率预测有效性的重要因素之一(Bi and Jiang, 2020)，且对Omi-R-J模型比对包括ETAS、R-J模型依赖于完备性震级的影响稍小一些.蒋长胜等(2018)利用OK1993模型和事件“删除”概率函数构建的随机地震序列目录测试结果表明，Omi-R-J模型对比于其他模型在余震监测能力有限的震后早期阶段逐步恢复到震前水平过程中，具有更好的适用性.

5 结论

造成震害叠加效应和进一步加重次生灾害的强余震，往往发生在震后的早期.为系统评估中国大陆地区强震的早期余震概率预测效能及制约因素，本文利用可充分考虑震后早期小震信息的Omi-R-J模型，对中国大陆的86个地震序列进行了系统性的预测研究，并利用N-test方法对预测结果进行了统计检验.获得如下结论：

(1)利用Omi-R-J模型对中国大陆的强震序列进行分震级、分时段的预测结果显示，Omi-R-J模型对中国大陆早期的强余震具有较好的预测能力，总的预测有效率为81.16%，平均预测有效率为83.82%，预测过少的比例大于预测过多的比例；对所选序列1天、3天和7天的预测平均有效率分别为88.90%、82.76%、80.82%；对不同震级档[M-3.0,M-2.5,M-2.0]预测的平均有效率为78.89%、86.14%、92.35%；基于震后早期不同时段([0.05, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0])预测的平均有效率为84.31%、92.06%、91.12%、90.74%、92.49%、91.87%、92.90%.

(2)通过分析震后早期1天和30天数据拟合参数的相关性发现，b、k、c在早期与此后更长的序列持续时间具有较强的相关性，差异较小，表明Omi-R-J模型可在震后早期记录相对不完整的阶段更早的获得较为稳定的序列参数. 更重要的是，这也意味着在震后1天内即可对震源区应力状态、余震触发能力获得可靠认识，这对强余震危险性研判和震后应急处置有重要意义；而p值的相关性弱，可能与p值表示的是余震活动的长期衰减特性、需要较长时间的数据来精确估计有关.由此表明在提升早期余震概率预测的准确性上，需要更为关注主震后p值的动态变化.

(3)地震序列参数k是影响地震数目的关键因素，而震后余震序列的发育程度可能是制约震后0.05天地震预测有效率重要因素之一.由此也表明除了利用类似于带有克服早期余震序列缺失问题的技术优势的Omi-R-J模型外，继续从利用模板匹配滤波、深度学习等技术提高早期余震的检测率也是提高余震预测效能的重要途径.

致谢研究中使用了中国地震台网中心“全国地震编目系统”提供的“统一正式目录”，日本东京大学生产技术研究所Takahiro Omi博士为本研究提供了程序和技术支持，两位评审专家提出了诸多建设性修改建议、对稿件质量提升帮助很大，在此一并表示感谢.