2017, Vol. 60
2. 国立统计数理研究所, 东京 106-8569;
3. 天津市地震局, 天津 300201
2. Institute of Statistical Mathematics, Tokyo 106-8569, Japan;
3. Tianjin Earthquake Agency, Tianjin 300201, China
作为地震预测中“可操作性”较强并可能取得明显减灾实效的余震概率预测，近年来在模型研发、应用和效能检验上取得了显著进展.包括Reasenberg-Jones(R-J)模型(Reasenberg and Jones, 1989)，“传染型余震序列”(Epidemic Type Aftershock Sequence, 简称ETAS)模型(Ogata, 1988, 1989, 2001)，“短期地震概率”(Short-Term Earthquake Probability，简称STEP)模型(Gerstenberger et al., 2005)，“地震警报短期预测”(Alarm-based Earthquake Forecast，简称EAST)模型(Shebalin et al., 2011)等短期余震概率预测模型得到显著进展.
近年来，“可操作的地震预测”(Operational Earthquake Forecasting，简称OEF)概念的提出与发展(Jordan and Jones, 2010; Jordan et al., 2011)，使得面向公众提供时间相依的地震危险性权威信息(authoritative information)、指导潜在的破坏性地震预防准备和实现有效地震减灾(Jordan et al., 2014)成为国际研究的热点.作为OEF研究的分支，“可操作的余震预测”(Operational Aftershock Forecasting，简称OAF)在美国地质调查局(USGS)和全球地震模型(GEM)项目的支持下得到快速发展.在OEF或OAF研究中，相对于长期预测模型的概率增益(probability gain，文中用Gain表示)是重要的表示方式.对STEP模型的预测效果评估结果表明，对于3~4级地震，概率增益可达Gain=10~100(Gerstenberger et al., 2007)，而对ETAS模型在美国加州和意大利地区的回溯性研究表明，其预测结果的概率增益Gain=10~100(Helmstetter et al., 2006; Console et al., 2010).结合严格的预测效果的统计检验，上述短期预测模型的发展应用以及预测策略研究(蒋长胜等, 2015; 毕金孟和蒋长胜, 2017)，将可为震后应急救援、灾民安置、次生灾害预防和震后重建提供科学决策依据.
据中国地震台网测定，北京时间2017年8月8日21时19分，在我国四川省阿坝州九寨沟县发生7.0级强震，震中位于33.20°N，103.82°E，震源深度20km.据民政部统计http://www.mca.gov.cn/article/yw/jzjz/zqkb/zqhz/201708/20170800005460.shtml，截至8月14日10时，四川九寨沟MS7.0地震已造成阿坝藏族羌族自治州九寨沟县、若尔盖县、松潘县和绵阳市平武县20.5万人受灾，25人死亡，525人因灾受伤(其中39人重伤，48人较重，438人轻伤)，8.9万人紧急转移安置；近100间房屋倒塌，5400余间严重损坏，6.7万间一般损坏.此次地震造成大量的大面积滚石、滑坡，地质灾害严重.由余震可能引起的地质次生灾害严重威胁震后应急救援和恢复重建工作.此外，有记载以来此次地震周边100km范围内曾发生20余次破坏性地震，包括1976年松潘—平武7.2级双震，对此次地震序列类型是否为双震型，以及周边断裂带的地震危险性，引起政府和社会公众的极大关注.由此，对九寨沟MS7.0地震序列的余震发生率和余震概率的“跟踪式”预测及预测策略研究，具有重要的科学价值和现实意义.本文将针对上述问题，采用国际上较为前沿的ETAS模型和得到广泛应用的R-J模型两个余震短期概率预测模型入手，通过序列参数、余震发生率和余震概率预测、余震预测的概率增益和效能检验等多角度的比较研究，试图为解决上述问题提供科学参考.
2 资料研究使用了中国地震台网中心和四川省地震局提供的全国地震编目系统《全国统一快报目录》，http://10.5.202.22/bianmu/index.jsp，查阅时间2017年8月14日.地震重新定位结果(http://www.cea-igp.ac.cn/tpxw/275881.html)显示，余震沿走向NE方向有优势分布，空间展布尺度约为30 km，垂直走向的剖面显示余震近乎垂直分布.考虑到上述余震重新定位结果，本文中根据余震的空间自然分布特征，将相对集中分布在32.8°N—33.8°N，103.6°E—104.8°E范围内(简称研究区)的事件作为余震序列进行研究.采用“震级—序号”法(见Huang, 2006; 蒋长胜和吴忠良，2011；Zhuang et al., 2012)对九寨沟MS7.0地震序列的震级完整性分析结果表明，在主震发生后的0.024天，序列的完整性震级可达Mc=ML2.0(图 1).根据上述分析结果，本文对九寨沟MS7.0地震序列的选取为，2017年8月8日21:00至8月13日14:00期间发生在研究区内的≥ML2.0地震.序列中ML2.0~2.9的地震383例、ML3.0~3.9的地震57例、ML4.0~4.9的地震13例，无5.0级以上余震.
|
图 1 研究使用的资料情况 (a)九寨沟MS7.0地震序列的空间分布，图中蓝色虚线框为地震序列选取的空间范围；(b)利用震级-序号法的九寨沟MS7.0地震序列目录完整性分析.图中横坐标为地震序列发生的先后序号，垂直虚线和水平虚线分别标出了地震序列完整的位置，对应完整性震级Mc=ML2.0，完整的时刻为主震后C0=0.024天. Fig. 1 The data used in the study (a) The spatial distribution of the Jiuzhaigou MS7.0 earthquake sequence, where the blue dashed frame is the spatial range selected for the earthquake sequence; (b) Magnitude-Rank distribution of the Jiuzhaigou MS7.0 earthquake sequence. The vertical and horizontal dashed lines indicate the completeness of ML2.0 aftershocks and the corresponding position of 0.024 day after the mainshock. |
ETAS模型是在“大森-宇津”公式(Omori, 1894; Utsu, 1961)基础上，通过定义所有的余震均可激发次级余震且震级分布独立的方式，对余震序列进行描述(Ogata, 1988, 1989).假定主震发生时刻为零时刻，对观测时间段(0, T]内的地震序列{(ti, Mi); i=1, 2, …, n}，ETAS模型的条件强度函数可表示为(Ogata, 1988)
|
(1) |
式(1) 中的t为主震发生后的时刻，一般取正值；M0为截止震级；Mi和ti分别为第i个事件的震级与发生时间；μ为震源区的背景地震发生率；pETAS表示序列衰减速率的快慢；cETAS为主震后余震频次达到峰值时的时间；常数KETAS表示余震的活跃程度；α表示序列触发次级余震的能力(Ogata, 1989)，α对于震群型序列在数值上一般小于1.0，而当被激发的次级余震较少时，α的取值一般大于1.0(Ogata, 2001).
ETAS模型参数[μ, KETAS, cETAS, α, pETAS]一般通过最大似然法(maximum likelihood estimates，简写为MLEs)进行估计.在拟合时间范围[S, T]内，ETAS模型似然函数L采用如下形式：
|
(2) |
利用上式对公式(1) 中的待估参数进行最大似然估计即可得到参数结果，相应的参数的渐进标准差，则可通过Fisher信息矩阵(即似然函数的二阶导数矩阵)的逆矩阵计算获得(Ogata, 1978).在获得Fisher信息矩阵后，5个参数的标准差即为逆矩阵各对角线元素的平方根.
利用ETAS模型进行的短期余震发生率预测，通常采用Ogata(1981)发展的修正的“瘦化算法”(thinning algorithm，Lewis and Shedler, 1979; Zhuang and Touati, 2015).“瘦化算法”基于ETAS模型的拟合参数，将地震序列转化为齐次泊松过程，并模拟未来的地震发生率.
3.2 Reasenberg-Jones(R-J)模型及其参数估计用于余震短期概率预测的Reasenberg-Jones(R-J)模型(Reasenberg and Jones, 1989)，主要基于余震衰减的“大森-宇津”公式(Omori, 1894; Utsu, 1961)和地震统计分布的G-R关系(Gutenberg and Richter, 1944)建立.其中，地震序列的频次N与时间t和震级M的关系为
|
(3) |
|
(4) |
由此，主震震级为Mm的、震级不低于M的余震，在t时刻的发生率可以表示为
|
(5) |
参数KOML、cOML、pOML可通过Ogata(1983)给出的最大似然法对大森-宇津公式拟合获得，并可求取相应的标准差.相应的对数似然函数可以表示为
|
(6) |
其中，
|
pOML≠1，对上式进行最大似然估计获得参数KOML、cOML、pOML后，还可利用最大似然法的渐进式，也即“费舍尔信息矩阵”(Fisher information matrix)估算这些参数的标准差.如用θ表示参数KOML、cOML、pOML，相应的Fisher信息矩阵J(θ; S, T)可用下式表示：
|
(7) |
对式(7) 中的J(θ; S, T)求取逆矩阵J-1，由此得到的J-1的对角线三个元素，求取平方根，即可得到参数KOML、cOML、pOML对应的标准差δKOML、δcOML和δpOML.
对参数b值可利用最大似然法(Utsu, 1965; Aki, 1965)计算获得
|
(8) |
其中的
|
(9) |
式(9) 中的n为参与计算的地震样本数量.
a值本身并不独立，需要利用KOML和b值进行计算.本文仅给出a值通过R-J模型中地震发生率λ(t, M)的两种表示形式推导，但不具体求解：
|
(10) |
其中，β=bln10.上式在M=M0的情况下可表示为
|
(11) |
对式(11) 两侧取对数并变换位置，可得
|
(12) |
除利用上式计算a值外，还可利用误差传递理论计算a值的标准差：
|
(13) |
对于采用“瘦化”算法和模型方式进行预测的ETAS模型，在进行n0次模拟时，至少发生一次预测目标地震的概率P可表示为：
|
(14) |
对概率P的误差，可利用中心极限定理(Central Limit Theorems，简称为“CLT”)来估计.CLT定理的含义是，大量相互独立的随机变量xi的均值的分布以正态分布为极限.因此在模拟预测的次数n0比较大时，预测结果的均值分布为
|
(15) |
式(15) 中的Q=1-P.对于随机变量x的方差可表示为
|
(16) |
式(16) 中的E(x)为均值.由此，在模拟时段内至少发生一次目标地震的变量xi{0, 1}的均值的分布，也即概率P的方差可由下式得到：
|
(17) |
|
(18) |
因此对概率P的标准差可由下式计算获得：
|
(19) |
而对于采用确定性的公式进行概率计算的R-J模型，相应的在时间范围[S, T]和震级范围(M1≤M≤M2)内至少发生1次地震的概率为
|
(20) |
如果确定参数KOML、cOML、pOML、b的数值，即可通过上式直接计算出给定预测时段[S, T]内的余震发生概率预测结果.余震发生率预测结果则可采用上式中右侧的对λ(t, M)的积分获得.
4 时间序列参数的稳定性分析 4.1 ETAS模型参数为考察ETAS模型对2017年九寨沟MS7.0地震序列的整体拟合情况，选定序列的起始时刻C0为主震后0.024天，对主震发生时刻直至2017年8月13日14:00(主震后4.70天)的Mc=ML2.0地震序列进行ETAS模型拟合.图 2给出了利用ETAS模型拟合的此次地震序列的条件强度曲线(单位时间内的地震发生率拟合结果)，以及获得的ETAS模型参数，其中μ=0.0000，KETAS=0.0027±0.0033，cETAS=0.0114±0.0288，α=2.0907±0.2842和pETAS=0.9912±0.0835.
|
图 2 由ETAS模型拟合给出的九寨沟MS7.0地震序列ML2.0以上地震的条件强度曲线(a)和序列M-t图(b) Fig. 2 (a) Temporal variation of the conditional intensity from fitting the ETAS model to the Jiuzhaigou MS7.0 earthquake sequence, with cutoff magnitude ML2.0; (b) M-t plot |
上述ETAS模型参数中，对描述序列衰减快慢的pETAS值，以及描述激发次级余震能力的α值对描述地震序列特征较为关键.在2010年以来中国大陆发生的与此次地震震级较为接近的一些强震相比较，例如2010年青海玉树MS7.1地震序列的α=0.9482、pETAS=1.0596(余娜等, 2016)，2013年四川芦山MS7.0地震序列的α=1.8864、pETAS=1.2205(蒋长胜等, 2013a)，2013年甘肃漳县—岷县MS6.6地震序列的α=1.7200、pETAS=1.1485(蒋长胜等, 2013b)，2014年云南鲁甸MS6.6地震序列的α=1.7177、pETAS=0.8919(蒋长胜等, 2015; Cheng et al., 2015)，2014年新疆于田MS7.3地震序列的α=1.6820、pETAS=0.6692(蒋长胜等, 2014)，2014年四川康定MS6.3地震序列的α=1.1909、pETAS=0.9781(Fang et al., 2015)，这些震例给出的平均的序列参数分别为，α=1.5242±0.3674，
在利用ETAS模型进行余震序列类型判定和短期预测时，平稳的序列参数有助于可靠地判断序列类型和进行余震发生率、余震概率预测(蒋长胜等, 2013b, 2015).为此，对九寨沟MS7.0地震序列的ETAS模型进行了连续、滑动拟合.将拟合时间窗的起始时刻固定为C0=0.024天，拟合时间窗结束时刻/序列持续时间t2自主震后0.52天开始，以0.02天步长，增加至4.70天，共进行了210个时段的滑动拟合.拟合获得的pETAS值、cETAS值、KETAS值和α值随序列持续时间的变化如图 3所示.由图可见，在序列持续时间t2=2.00天之前，九寨沟MS7.0地震序列的pETAS值、cETAS值、KETAS值和α值均变化较为剧烈，相应的标准差的数值范围也较大，显示了较不稳定的拟合状态和序列发生的复杂性.在t2=2.00天之后，上述参数均较为稳定、拟合标准差较小，其中，pETAS的均值为1.0170，平均的标准差为0.0495；cETAS的均值为0.0596，平均的标准差为0.0075；KETAS值的均值为0.0020，平均的标准差为0.0030；α值均值为2.3826，平均的标准差为0.3737，与图 2给出的最长序列持续时间t2=4.70天时的拟合值较为接近.
|
图 3 九寨沟MS7.0地震序列ETAS模型拟合参数随序列持续时间的变化 (a)参数pETAS值随序列持续时间的变化；(b)参数cETAS值随序列持续时间的变化；(c)参数KETAS值随序列持续时间的变化；(d)参数α值随序列持续时间的变化.序列截止震级Mc=ML2.0，起始时刻固定为C0=0.024天，序列持续时间t2以0.02天步长，自主震后0.52天增加至4.70天进行滑动拟合.图中灰色区域为相应的参数的标准差的范围. Fig. 3 Parameters of ETAS model against the duration time (since the mainshock occurred) in fitting the Jiuzhaigou MS7.0 earthquake sequence (a) The pETAS values; (b) The cETAS values; (c) The KETAS values; (d) The a values. The cut-off magnitude is set to Mc=ML2.0 in ETAS model fitting, the starting time fixed by C0=0.024 day and the ending time from 0.52 day to 4.70 days by sliding with 0.02 day steps. The gray area gives the standard deviation range of the parameters shown in each subgraph. |
对R-J模型参数，采用与时间序列ETAS模型相同的考察方式.首先对于九寨沟MS7.0地震序列的整体拟合中，同样选定序列的起始时刻C0为主震后0.024天，对主震发生时刻直至主震后4.70天的Mc=ML2.0地震序列，利用最大似然法进行R-J模型拟合.获得的序列持续时间t2=4.70天时的R-J模型参数分别为，pOML=0.900±0.0920，cOML=0.0810±0.0370，KOML=113.0750±9.3880和b=0.7821±0.0346.上述R-J模型参数中，同样描述序列衰减快慢的pOML值与ETAS模型拟合获得的pETAS值非常接近，按照Utsu(1965)给出的pOML的理论均值为1.0的结果，R-J模型的拟合结果也证明了此次九寨沟MS7.0地震序列的衰减较为正常.
针对九寨沟MS7.0地震序列同样考察了R-J模型参数的稳定性，进行了连续、滑动拟合.采用与ETAS模型参数稳定性研究中相同的设置方式，即对拟合时间窗的起始时刻固定为C0=0.024天，序列持续时间t2自主震后0.52天开始，以0.02天步长，增加至4.70天，共进行210个时段的滑动拟合.获得的pOML值、cOML值、KOML值和b值随序列持续时间的变化如图 4所示.由图可见，在序列持续时间t2=1.50天之前，九寨沟MS7.0地震序列的pOML值、cOML值、KOML值变化较为剧烈，相应的标准差数值范围也较大，显示了在震后早期的余震序列的复杂性.
|
图 4 九寨沟MS7.0地震序列R-J模型拟合参数随序列持续时间的变化 (a)参数pOML值随序列持续时间的变化；(b)参数cOML值随序列持续时间的变化；(c)参数KOML值随序列持续时间的变化；(d)参数b值随序列持续时间的变化.序列截止震级Mc=ML2.0，起始时刻固定为C0=0.024天，序列持续时间t2以0.02天步长，自主震后0.52天增加至4.70天进行滑动拟合.图中灰色区域为相应的参数的标准差的范围. Fig. 4 Parameters of R-J model against the duration time (since the mainshock occurred) in fitting the Jiuzhaigou MS7.0 earthquake sequence (a) The pOML values; (b) The cOML values; (c) The KOML values; (d) The b values. The cut-off magnitude is set to Mc=ML2.0 in R-J model fitting, the starting time fixed by C0=0.049 day and the ending time from 0.52 day to 4.70 days by sliding with 0.02 day steps. The gray area gives the standard deviation range of the parameters shown in each subgraph. |
在t2=1.50天之后，上述参数均较为稳定、拟合标准差较小，其中，pOML的均值为1.0123，平均的标准差为0.0359；cOML的均值为0.0873，平均的标准差为0.0081；KOML的均值为113.7176，平均的标准差为1.1288.b值展示了趋势性变化，由0.7391增加至0.7830，表明震源区应力累计水平在减小、断层愈合过程仍在持续进行中.
5 短期余震发生率和概率预测为考察ETAS模型和R-J模型对九寨沟MS7.0地震序列的余震发生率和概率预测情况，并进行比较研究，这里仅讨论对未来1天的预测情况，并采用与模型参数稳定性研究相同的设定规则.首先将拟合时间窗的起始时刻固定为C0=0.024天，序列持续时间t2自主震后0.52天开始，以0.02天步长，增加至4.70天进行连续预测.两种模型对不同的震级的“预测目标地震”(target earthquake for forecast)MT=[2.0, 3.0, …, 6.0]的预测结果如图 5所示，其中图 5a为未来1天的余震发生率预测，图 5b为相应的“至少发生一次”预测目标地震的概率.其中，在t2=4.70天时，ETAS模型对未来一天发生不低于MT=[2.0, 3.0, …, 6.0]的结果分别为22.8200次、3.3490次、0.4720次、0.0650次、0.0110次，相应的发生至少一次不低于MT=[2.0, 3.0, …, 6.0]的概率分别为1.0000、0.9590、0.3710、0.0550、0.0100；R-J模型对未来一天发生不低于MT=[2.0, 3.0, …, 6.0]的结果分别为39.7172次、6.5598次、1.0834次、0.1789次、0.0296次，相应的发生至少一次不低于MT=[2.0, 3.0, …, 6.0]的概率分别为1.000、0.9986、0.6616、0.1638、0.0291.
|
图 5 九寨沟MS7.0地震序列未来1天的余震预测结果 (a) R-J模型和ETAS模型的余震发生率预测结果；(b) R-J模型和ETAS模型的至少发生一次预测目标震级以上地震(Nfore≥1) 的概率.图中不同颜色曲线对应不同截止震级下的地震发生率预测结果，粗线和细线分别为R-J模型和ETAS模型预测结果.序列起始时刻固定为C0=0.024天，序列持续时间以0.02天步长，自震后0.52天增加至4.70天进行滑动预测. Fig. 5 Next-day′s aftershock forecasting results of the Jiuzhaigou MS7.0 earthquake sequence (a) Forecasted aftershock rate of the R-J model and ETAS mode; (b) Forecasted probability with Nfore≥1 of the R-J model and ETAS mode. The curves show the forecasting results with different cut-off magnitude, and the thick and thin lines are the ETAS model and the R-J model forecasting results respectively. The starting time fixed by C0=0.024 day and the ending time from 0.52 day to 4.70 days by sliding with 0.02 day steps. |
由图 5可见，由于ETAS模型是采用“瘦化算法”和随机模拟的方式进行余震发生率和概率预测，尽管使用了多达1000次的模拟，但余震发生率和概率随序列持续时间的预测结果仍展示了明显的涨落，尤其是在t2=2.00天之前.对于R-J模型，由于余震的发生率和概率是通过公式(10) 直接计算获得，相应的预测结果主要依赖于对序列持续时间内(C0至t2)的模型参数变化情况，因此随序列持续时间的预测结果表现出具有明显的稳定性，尤其是在t2=1.50天之后.此外，在210次连续滑动预测中，ETAS模型给出的余震发生率预测结果的数值几乎均低于R-J模型预测结果；对于余震概率预测结果，除在M≥2.0时两模型给出的预测结果均为100%外，对其他预测目标震级，ETAS模型结果几乎均低于R-J模型预测结果的数值.
6 预测结果的统计检验为进一步地检验ETAS模型和R-J模型在对九寨沟MS7.0地震序列未来1天的余震短期预测中的效能，分别采用N-test方法检验了ETAS模型和R-J模型对余震发生率预测的可靠程度，以及采用T-test方法对两个模型之间的相对预测效能的优劣进行了比较.
6.1 地震发生率预测的N-test统计检验为定量化地检验ETAS模型和R-J模型对九寨沟MS7.0地震序列未来1天余震发生率预测的可靠性，本研究使用N-test方法(Kagan and Jackson, 1995; Schorlemmer et al., 2007; Zechar, 2010)对两种方法的预测结果进行了系统的检验.N-test方法主要是用来考察预测的“目标地震”数目与实际发生地震数是否存在偏离(Kagan and Jackson, 1995; Schorlemmer et al., 2007)，多基于累积泊松分布函数生成的模拟地震目录或直接使用理论公式，并进行不同置信度的检验.当预测的地震发生数为Nfore，实际地震发生数为Nobs时，可通过两个评分量δ1和δ2分别回答发生“至少”和“不超过”Nobs个地震事件的概率(Zechar, 2010).相应的泊松累积分布函数F主要是指在给定的平均发生率λ下，发生不超过ω个事件的概率(Zechar, 2010)：
|
(21) |
与之相对应的评分变量为
|
(22) |
|
(23) |
在设定的有效显著性水平αeff=0.025的情况下，如果δ1＜αeff表明Nfore相对于Nobs过少；如果δ2＜αeff则表明Nfore相对于Nobs过多.
同样受到地震序列目录观测时间和实际发生余震数量的限制，这里仅对九寨沟MS7.0地震序列的M≥2.0级、持续时间t2=0.52~3.70天、步长0.02天的160次连续滑动预测结果进行了评价.相应的N-test检验结果如图 6a和图 6b所示.
|
图 6 对九寨沟MS7.0地震序列未来1天余震预测结果的统计检验 (a) ETAS模型的N-test统计检验结果，图中灰色垂线为δ1评分与δ2评分的连线，黑色菱形和蓝色方形分别表示δ1 < 0.025和δ2 < 0.025的结果；(b) R-J模型的N-test统计检验结果；(c) ETAS模型相对于R-J模型预测效能的T-test检验，图中分别给出了M≥2.0、M≥3.0和M≥4.0的检验结果，圆点为平均的每个地震的信息增益(information gain per earthquake)，横线为95%置信区间范围，数字标出了预测时间窗内实际发生的地震数量. Fig. 6 Statistical test of the Next-day′s aftershock forecasting of the Jiuzhaigou MS7.0 earthquake sequence (a) N-test result of the ETAS model, the vertical gray dashed lines connect the δ1 and δ2 values, the black diamonds and blue squares show the results with δ1 < 0.025 and δ2 < 0.025 respectively; (b) N-test result of the R-J model; (c) T-test applied to retrospective evaluation of the ETAS model over R-J model, with 'target magnitude' M≥2.0, M≥3.0 and M≥4.0 in aftershock forecasting. The black dots and horizontal lines represent the information gain per earthquake and 95% confidence intervals respectively of the ETAS model over R-J model, the number of observed earthquakes is shown above the T-test results. |
由图 6a可见，在ETAS模型参数较不稳定的t2=2.00天之前的震后早期阶段，出现大量的预测地震数目偏少的情况.其中，对t2=0.52~2.00天的75次预测，δ1＜αeff的高达57次，占总次数的76%；δ2＜αeff的5次，占总次数的7%；仅有12次预测的余震数目通过了检验.而在此后的t2=0.22~3.70天的85次预测中，δ1＜αeff的为0次；δ2＜αeff的为10次，占总次数的11.7%.由图 6b可见，在R-J模型参数较不稳定的t2=1.50天之前的震后早期阶段，总体预测情况较好.在t2=0.52~1.50天的50次预测中，δ1＜αeff的仅有4次，占总次数的8%；δ2＜αeff的11次，占总次数的22%；其余的35次预测的余震数目通过了检验，占总次数的70%.而在此后的t2=1.52~3.70天的110次预测中，δ1＜αeff的为0次；δ2＜αeff的高达109次，占总次数的99.1%；仅有1次预测的余震数目通过了检验.
由上述的对比可见，在ETAS模型序列参数较不稳定的震后早期阶段，高达76%的预测出现预测地震数量偏少的情况，而在序列参数稳定的时段则预测效果明显较好；对于R-J模型，在震后序列参数较不稳定的震后早期阶段，达到70%的预测可以通过N-test检验，而在此后的预测中，R-J模型预测的余震数量明显偏多.
6.2 ETAS模型相对R-J模型的预测效能评价尽管本文讨论了N-test给出的ETAS模型和R-J模型的预测效能评价，以及各自的概率增益Gain随序列持续时间的分布，但无法确定两种预测模型的相对优势.为对ETAS模型和R-J模型在九寨沟MS7.0地震序列的短期预测应用的效果进行比较，这里采用了T-test(Student′s)方法来实现.作为CSEP计划推荐应用的地震预测检验技术，T-test方法是通过计算置信区间内模型A相对于B的每个地震平均的“地震信息增益”(information gain per earthquake，IGPE)来进行比较(Imoto, 2007)：
|
(24) |
式中，N为实际发生的预测目标震级以上的地震数，lnLA和lnLB分别为模型A和模型B的似然函数，可通过下式计算获得：
|
(25) |
式(25) 中，n为离散化的空间网格个数，w为实际发生地震的序号；λ(i)为第i个网格的预测的目标震级以上地震的个数.
对于本研究的时间序列ETAS模型相对R-J模型，上式中右侧的对空间网格的计算可得到简化，直接使用相应的预测时间窗内的预测地震数量.值得注意的是，由于在实际应用中常常采用类似本文的滑动步长较小、预测时间窗相互重叠的连续预测，因此需要对相互重叠的预测时间窗进行IGPE平均值的折算.此外，T-test方法还可通过查阅Nobs与置信水平(confidence level)对照表中的系数，用于估算相应置信水平下的标准差并对预测结果进行任意置信区间的检验.按照上述的IGPE计算方法和T-test检验技术，本文对九寨沟MS7.0地震序列在持续时间t2=0.52~3.70天、步长0.02天的160次连续滑动预测，进行了ETAS模型相对R-J模型预测效能的评价，其中对M≥2.0、M≥3.0和M≥4.0目标地震在95%置信区间下的检验结果如图 6c所示.由图可见，由于获得的平均的IGPE均大于0，且95%置信区间范围也均超过0，因此在预测的有效性上ETAS模型优于R-J模型.对不同的预测目标震级，ETAS模型相对于R-J模型在M≥3.0时预测效果优势最为明显.
7 结论和讨论在对2017年四川九寨沟MS7.0地震的短期余震预测研究中，本文分别对ETAS模型和R-J模型的序列参数稳定性和余震预测效果进行了比较研究，考察了两种模型在余震发生率预测、余震概率预测，以及利用N-test和T-test检验方法的评价了预测效果.获得如下认识：
(1) 对ETAS模型和R-J模型放入序列参数稳定性研究表明，在序列持续时间t2=4.70天时的ETAS模型参数分别为μ=0.0000，KETAS=0.0027±0.0033，cETAS=0.0114±0.0288，α=2.0907±0.2842和pETAS=0.9912±0.0835；R-J模型参数分别为pOML=0.900±0.0920，cOML=0.0810±0.0370，KOML=113.0750±9.3880和b=0.7821±0.0346.其中，ETAS模型参数在t2=2.0天后较为稳定，R-J模型参数在t2=1.50天后较为稳定.两种模型给出的序列衰减速率参数均表明，此次九寨沟MS7.0地震序列的衰减较为正常.
(2) 利用ETAS模型和R-J模型分别对九寨沟MS7.0地震序列未来1天的余震发生率预测和余震概率进行了连续滑动预测，对于M≥3.0级、M≥4.0、M≥5.0级和M≥6.0级下的预测，ETAS模型给出的余震发生率和余震概率均低于R-J模型预测结果.
(3) 利用N-test方法检验了ETAS模型和R-J模型对九寨沟MS7.0地震序列未来1天余震发生率预测的效能，结果表明，对于M≥2.0级地震在序列持续时间t2=0.52~3.70天的连续预测，ETAS模型序列在参数较不稳定的震后早期阶段，预测地震数量偏少的情况高达76%，但在序列参数稳定的时段则预测效果较好.R-J模型在序列参数较不稳定的震后早期阶段预测效果较好，通过N-test检验的达到70%，但在参数稳定时段则预测的余震数明显偏多.
(4) 通过“地震信息增益”(IGPE)和计算T-test检验方法，对ETAS模型相对于R-J模型预测效能的比较研究表明，对于九寨沟MS7.0地震序列的未来1天的短期预测效果，在95%的置信区间上，ETAS模型明显优于R-J模型，尤其是对于M≥3.0余震的预测.
由上述结论可知，一方面，ETAS模型的预测效果优于R-J模型、可提供更为有效的预测信息.这与全球CSEP计划在国际上多个检验中心(Testing Center)获得的结果较为一致(Gerstenberger and Rhoades, 2010; Schorlemmer et al., 2010; Werner et al., 2010)，显示了ETAS模型在描述复杂地震序列和进行余震预测上的较大优势.另一方面，在震后早期阶段R-J模型总体上优于ETAS模型，而在序列参数稳定时段ETAS模型预测效果较好、R-J模型几乎失效.这种现象是否还适用于其他序列，以及基本原理均基于“大森-宇津”公式的两种模型，在不同的序列发展阶段预测效果差异明显的背后是否存在明确的物理机制等问题，仍需细致深入的研究.作为一种可能的解释，在震后的早期阶段，余震序列的衰减规律主要受到主震的“控制”，对于考虑激发次级余震活动的ETAS模型，如果次级余震不发育则预测效果可能有限，甚至劣于原理简单、仅考虑主震触发作用的R-J模型；而随着序列的持续发展，次级余震更为活跃、相应的次级余震的序列更为完整时，ETAS模型描述复杂余震序列的优势更为明显.此外，由于九寨沟MS7.0地震序列截至目前尚未发生5.0级及以上地震，本文中使用的ETAS模型和R-J模型在序列早期均给出了相应的偏高的概率预测结果，如何在主震后的早期阶段提高较大余震的预测准确率，可能需要进一步的研究，例如利用贝叶斯方法和区域内以往震例的平均序列参数(Reasenberg and Jones, 1989)进行余震预测等等.
如果本文的研究结果具有广泛的一般性，则可在震后早期使用R-J模型、中后期使用ETAS模型，或可取得较好的预测效果.实际上，采用不同模型优势而进行组合形成的混合模型(hybird model)，正是目前国际上地震预测研究的前沿热点问题(Rhoades and Gerstenberger, 2009; Marzocchi et al., 2014; Shebalin et al., 2014; Rhoades et al., 2016; Han et al., 2017).由于九寨沟MS7.0地震序列仍在进行中，文中使用的地震序列持续时间有限，一定程度上限制了对未来3天、7天、30天等更长时间尺度余震预测的研究.此次地震在震后早期阶段的参数不稳定性的原因，以及随着震源区应力场调整的变化规律等，都有待进一步研究.
致谢本文属中国地震局“四川九寨沟7.0级地震系统性科学考察”经费项目工作.“国际地震可预测性合作研究”(CSEP)计划中国检验中心筹备组对本研究给予了指导.研究中使用了中国地震台网中心“全国地震编目系统”提供的《统一快报目录》，四川省地震监测中心苏金蓉主任、地震预测研究中心易桂喜研究员、龙锋高级工程师、美国佐治亚理工学院彭志刚(Zhigang Peng)教授与作者进行了有益的讨论.两位评审专家提出了诸多有益建议，对稿件的提升帮助很大，在此一并表示感谢.
Aki K.
1965. Maximum likelihood estimate of b in the formula logN=a-bM and its confidence limits. Bull. Earthquake Re. Inst., Tokyo Unive, 43: 237-239.
|
|
Bi J M, Jiang C S.
2017. Evaluation on the forecasting effectiveness of short-term aftershock occurrence rate and forecasting strategies at the junction of Shanxi, Hebei and Inner Mongolia. Progress in Geophysics, 32(1): 8-17.
DOI:10.6038/pg20170102 |
|
Cheng J, Wu Z L, Liu J, et al.
2015. Preliminary report on the 3 August 2014, MW6.2/MS6.5 Ludian, Yunnan-Sichuan border, southwest China, earthquake. Seismological Research Letters, 86(3): 750-763.
DOI:10.1785/0220140208 |
|
Console R, Murru M, Falcone G.
2010. Probability gains of an epidemic-type aftershock sequence model in retrospective forecasting of M ≥ 5 earthquakes in Italy. Journal of Seismology, 14(1): 9-26.
DOI:10.1007/s10950-009-9161-3 |
|
Department of Earthquake Disaster Prevention, China Earthquake Administration. 1999.
The Catalogue of Chinese Modern Earthquakes. Beijing: China Science and Technology Press.
|
|
Department of Earthquake Disaster Prevention, State Seismological Bureau. 1995.
The Catalogue of Chinese Historical Strong Earthquakes. Beijing: Seismological Press.
|
|
Fang L H, Wu J P, Liu J, et al.
2015. Preliminary report on the 22 November 2014 MW6.1/MS6.3 Kangding earthquake, western Sichuan, China. Seismological Research Letter, 86(6): 1603-1613.
DOI:10.1785/0220150006 |
|
Gardner J K, Knopoff L.
1974. Is the sequence of earthquakes in Southern California, with aftershocks removed, Poissonian?. Bulletin of the Seismological Society of America, 64(5): 1363-1367.
|
|
Gerstenberger M C, Wiemer S, Jones L M, et al.
2005. Real-time forecasts of tomorrow's earthquakes in California. Nature, 435(7040): 328-331.
DOI:10.1038/nature03622 |
|
Gerstenberger M C, Jones L M, Wiemer S.
2007. Short-term aftershock probabilities:Case studies in California. Seismological Research Letters, 78(1): 66-77.
DOI:10.1785/gssrl.78.1.66 |
|
Gerstenberger M C, Rhoades D A.
2010. New Zealand earthquake forecast testing centre. Pure and Applied Geophysics, 167(8-9): 877-892.
DOI:10.1007/s00024-010-0082-4 |
|
Gutenberg R, Richter C F.
1944. Frequency of earthquakes in California. Bulletin of the Seismological Society of America, 34(1): 185-188.
|
|
Han P, Hattori K, Zhuang J, et al.
2017. Evaluation of ULF seismo-magnetic phenomena in Kakioka, Japan by using Molchan's error diagram. Geophysical Journal International, 208(1): 482-490.
DOI:10.1093/gji/ggw404 |
|
Helmstetter A, Kagan Y Y, Jackson D D.
2006. Comparison of short-term and time-independent earthquake forecast models for southern California. Bulletin of the Seismological Society of America, 96(1): 90-106.
DOI:10.1785/0120050067 |
|
Huang Q H.
2006. Search for reliable precursors:A case study of the seismic quiescence of the 2000 western Tottori prefecture earthquake. Journal of Geophysical Research, 111(B4): B04301.
DOI:10.1029/2005JB003982 |
|
Imoto M.
2007. Information gain of a model based on multidisciplinary observations with correlations. Journal of Geophysical Research, 112: B05306.
DOI:10.1029/2006JB004662 |
|
Jiang C S, Wu Z L.
2011. Intermediate-term medium-range Accelerating Moment Release (AMR) prior to the 2010 Yushu MS7.1 earthquake. Chinese Journal of Geophysics, 54(6): 1501-1510.
DOI:10.3969/j.issn.0001-5733.2011.06.009 |
|
Jiang C S, Zhuang J C, Long F, et al.
2013a. Statistical analysis of ETAS parameters in the early stage of the 2013 Lushan MS7.0 earthquake sequence. Acta Seismologica Sinica, 35: 661-669.
|
|
Jiang C S, Wu Z L, Han L B, et al.
2013b. Effect of cutoff magnitude Mc of earthquake catalogues on the early estimation of earthquake sequence parameters with implication for the probabilistic forecast of aftershocks:The 2013 Minxian-Zhangxian, Gansu, MS6. 6 earthquake sequence. Chinese Journal of Geophysics, 56(12): 4048-4057.
DOI:10.6038/cjg20131210 |
|
Jiang C S, Han L B, Guo L J.
2014. Parameter characteristics in the early period of three earthquake sequences in the Yutian, Xinjiang since 2008. Acta Seismologica Sinica, 36(2): 165-174.
DOI:10.3969/j.issn.0253-3782.2014.02.002 |
|
Jiang C S, Wu Z L, Yin F L, et al.
2015. Stability of early-estimation sequence parameters for continuous forecast of the aftershock rate:A case study of the 2014 Ludian, Yunnan MS6.5 earthquake. Chinese Journal of Geophysics, 58(11): 4163-4173.
DOI:10.6038/cjg20151123 |
|
Jordan T H, Jones L M.
2010. Operational earthquake forecasting:Some thoughts on why and how. Seismological Research Letters, 81(4): 571-574.
DOI:10.1785/gssrl.81.4.571 |
|
Jordan T H, Chen Y T, Gasparini P, et al.
2011. Operational earthquake forecasting:State of knowledge and guidelines for implementation. Annals of Geophysics, 54(4): 315-391.
|
|
Jordan T H, Marzocchi W, Michael A J, et al.
2014. Operational earthquake forecasting can enhance earthquake preparedness. Seismological Research Letters, 85(5): 955-959.
DOI:10.1785/0220140143 |
|
Kagan Y Y, Jackson D D.
1995. New seismic gap hypothesis:Five years after. Journal of Geophysical Research, 100(B3): 3943-3959.
DOI:10.1029/94JB03014 |
|
Lewis P A W, Shedler G S.
1979. Simulation of nonhomogeneous Poisson processes by thinning. Naval Research Logistics Quarterly, 26(3): 403-413.
DOI:10.1002/(ISSN)1931-9193 |
|
Marzocchi W, Lombardi A M, Casarotti E.
2014. The establishment of an operational earthquake forecasting system in Italy. Seismological Research Letters, 85(5): 961-969.
DOI:10.1785/0220130219 |
|
Ogata Y.
1978. The asymptotic behaviour of maximum likelihood estimators for stationary point processes. Annals of the Institute of Statistical Mathematics, 30(1): 243-261.
DOI:10.1007/BF02480216 |
|
Ogata Y.
1981. On Lewis' simulation method for point processes. IEEE Transactions on Information Theory, 27(1): 23-31.
DOI:10.1109/TIT.1981.1056305 |
|
Ogata Y.
1983. Estimation of the parameters in the modified Omori formula for aftershock frequencies by the maximum likelihood procedure. Journal of Physics of the Earth, 31(2): 115-124.
DOI:10.4294/jpe1952.31.115 |
|
Ogata Y.
1988. Statistical models for earthquake occurrences and residual analysis for point processes. Journal of the American Statistical Association, 83(401): 9-27.
DOI:10.1080/01621459.1988.10478560 |
|
Ogata Y.
1989. Statistical model for standard seismicity and detection of anomalies by residual analysis. Tectonophysics, 169(1-3): 159-174.
DOI:10.1016/0040-1951(89)90191-1 |
|
Ogata Y.
2001. increased probability of large earthquakes near aftershock regions with relative quiescence. Journal of Geophysical Research, 106(B5): 8729-8744.
DOI:10.1029/2000JB900400 |
|
Omori F.
1894. On aftershocks of earthquakes. Journal of the College of Science, Imperial University of Tokyo, 7: 11-200.
|
|
Reasenberg P A, Jones L M.
1989. Earthquake hazard after a mainshock in California. Science, 243(4895): 1173-1176.
DOI:10.1126/science.243.4895.1173 |
|
Rhoades D A, Gerstenberger M C.
2009. Mixture models for improved short-term earthquake forecasting. Bulletin of the Seismological Society of America, 99(2A): 636-646.
DOI:10.1785/0120080063 |
|
Rhoades D A, Liukis M, Christophersen A, et al.
2016. Retrospective tests of hybrid operational earthquake forecasting models for Canterbury. Geophysical Journal International, 204(1): 440-456.
DOI:10.1093/gji/ggv447 |
|
Schorlemmer D, Gerstenberger M C, Wiemer S, et al.
2007. Earthquake likelihood model testing. Seismological Research Letters, 78(1): 17-29.
DOI:10.1785/gssrl.78.1.17 |
|
Schorlemmer D, Zechar J D, Werner M J, et al.
2010. First results of the regional earthquake likelihood models experiment. Pure and Applied Geophysics, 167(8-9): 859-876.
DOI:10.1007/s00024-010-0081-5 |
|
Shebalin P, Narteau C, Holschneider M, et al.
2011. Short-term earthquake forecasting using early aftershock statistics. Bulletin of the Seismological Society of America, 101(1): 297-312.
DOI:10.1785/0120100119 |
|
Shebalin P N, Narteau C, Zechar J, et al.
2014. Combining earthquake forecasts using differential probability gains. Earth, Planets and Space, 66: 37.
DOI:10.1186/1880-5981-66-37 |
|
Shi Y L, Bolt B A.
1982. The standard error of the magnitude-frequency b-value. Bulletin of the Seismological Society of America, 72(5): 1677-1687.
|
|
Utsu T.
1961. A statistical study on the occurrence of aftershocks. Geophysical Magazine, 30: 521-605.
|
|
Utsu T.
1965. A method for determining the value of b in a formula logN=a-bM showing the magnitude-frequency relation for earthquakes. Geophysical Bulletin of the Hokkaido University, 13: 99-103.
|
|
Werner M J, Zechar J D, Marzocchi W, et al.
2010. Retrospective evaluation of the five-year and ten-year CSEP-Italy earthquake forecasts. Annals of Geophysics, 53(3): 11-30.
|
|
Yu N, Jiang C S, Ma Y H.
2016. Analysis of ETAS model parameters and associated variation characteristics for 2010 Yushu, Qinghai MS7.1 earthquake sequence. China Earthquake Engineering Journal, 38(4): 609-615, 623.
|
|
Zechar J D. 2010. Evaluating earthquake predictions and earthquake forecasts:a guide for students and new researchers.//Community Online Resource for Statistical Seismicity Analysis. doi:10.5078/corssa-77337879.Availableathttp://www.corssa.org.
|
|
Zhuang J C, Harte D, Werner M J, et al. 2012. Basic models of seismicity:temporal models.//Community Online Resource for Statistical Seismicity Analysis. doi:10.5078/corssa-79905851.Availableathttp://www.corssa.org.
|
|
Zhuang J C, Touati S. 2015. Stochastic simulation of earthquake catalogs.//Community Online Resource for Statistical Seismicity Analysis.doi:10.5078/corssa-43806322.Availableathttp://www.corssa.org.
|
|
毕金孟, 蒋长胜.
2017. 晋冀蒙交界地区余震短期发生率的预测效能评估和预测策略研究. 地球物理学进展, 32(1): 8–17.
DOI:10.6038/pg20170102 |
|
国家地震局震害防御司. 1995.
中国历史强震目录(公元前23世纪~公元1911年). 北京: 地震出版社.
|
|
蒋长胜, 吴忠良.
2011. 2010年玉树MS7. 1地震前的中长期加速矩释放(AMR)问题. 地球物理学报, 54(6): 1501–1510.
DOI:10.3969/j.issn.0001-5733.2011.06.009 |
|
蒋长胜, 庄建仓, 龙锋, 等.
2013a. 2013年芦山MS7.0地震序列参数的早期特征:传染型余震序列模型计算结果. 地震学报, 35(5): 661–669.
|
|
蒋长胜, 吴忠良, 韩立波, 等.
2013b. 地震序列早期参数估计和余震概率预测中截止震级Mc的影响:以2013年甘肃岷县-漳县6. 6级地震为例. 地球物理学报, 56(12): 4048–4057.
DOI:10.6038/cjg20131210 |
|
蒋长胜, 韩立波, 郭路杰.
2014. 新疆于田地区2008年以来3个地震序列的参数早期特征. 地震学报, 36(2): 165–174.
DOI:10.3969/j.issn.0253-3782.2014.02.002 |
|
蒋长胜, 吴忠良, 尹凤玲, 等.
2015. 余震的序列参数稳定性和余震短期发生率预测效能的连续评估——以2014年云南鲁甸6.5级地震为例. 地球物理学报, 58(11): 4163–4173.
DOI:10.6038/cjg20151123 |
|
余娜, 蒋长胜, 马玉虎.
2016. 2010年青海玉树MS7. 1地震序列ETAS模型参数及其变化特征研究. 地震工程学报, 38(4): 609–615, 623.
|
|
中国地震局震害防御司. 1999.
中国近代地震目录. 北京: 中国科学技术出版社.
|
|