据中国地震台网测定,2014年8月3日16时30分在云南省昭通市鲁甸县(27.1°N,103.3°E)发生M_S6.5地震,并造成人员重大伤亡。地震发生后,关于此次地震序列类型、序列参数和是否发生强余震等问题引起广泛关注。

由于地震序列参数尤其是参数的早期特征,对序列类型快速判定(Guo,Ogata,1997; 蒋海昆等,2007)、研究震源区构造应力水平和区域大地热流值(Kagan et al.,2010)等有重要参考意义。因此,在强震后获得科学、准确的地震序列参数具有重要的现实意义。近年来,基于点分支过程的“传染型余震序列”(ETAS)模型得到持续发展,由于可接近真实地描述余震序列活动,ETAS模型在地震序列参数早期特征研究中发挥了重要作用。

为研究云南鲁甸M_S6.5地震序列参数早期特征,本文考虑了不同的截止震级、序列持续时间的影响,获得了云南鲁甸M_S6.5地震序列的ETAS模型参数,以期对此次地震的序列类型判定等提供参考。

作为修正的“大森—宇津”公式(Omori,1894; Utsu,1961)的推广,ETAS模型假设所有的余震均可激发自己的余震,且相互独立,被用于对接近真实地震序列的模拟(Ogata,1988):

λ(t)=μ+K∑_ti<t(e^α(M_i-M₀))/((t-t_i+c)^p), M_i>M₀.(1)

式中,M_i和t_i分别为事件i的震级和发生时刻,M₀为参考震级,一般取为截止震级。p值表示序列衰减的快慢,p值越大衰减越快,反之越慢。α值用于描述触发次级余震的能力(Ogata,1989,1992)。

ETAS模型参数可用最大似然法进行估计,利用在时间段[t_c,T]内的似然函数:

lgL=∑_i:S≤ti<Tlgλ(t_i)-∫^T∫_Sλ(t)dt.(2)

即可对ETAS模型参数[μ,K,c,α,p]进行估计。为对ETAS模型的拟合效果进行评估,这里使用了“转换时间”域的“残差分析”方法(Ogata,1988; Daley,Vere-Jones,2003)。“残差分析”方法的基本原理是:条件强度曲线λ(t)逐点积分后即为单调递增函数:

t|→τ=Λ(t)=∫^t∫₀λ(u)du.(3)

根据随机点过程理论,这一积分结果可成为稳态泊松过程。由此,可通过积分转换,将时间序列{t_i}逐点变换为“转换时间”{τ_i},并通过判断“转换时间”域{τ_i}下的实际地震序列累积地震数是否接近直线判断模型的拟合效果(Zhuang et al.,2012)。

本研究使用的鲁甸M_S6.5地震序列目录来自中国地震台网中心“全国地震编目系统”提供的《云南地震台网速报目录》 全国统一地震编目(日报)查询系统,http://www.csndmc.ac.cn/newweb/catalog_direct_link.htm#,查阅时间2014年8月16日.,选定2014年8月3日16时30分11秒至8月12日16时17分4.8秒,空间范围为(26.8°～27.4°N,103°～103.4°E)的地震目录作为序列目录。在云南鲁甸M_S6.5地震震后9天内,目录中共包括M_L≥4.0地震5次,M_L3.0～3.9地震41次,M_L2.0～2.9地震263次,M_L1.5～1.9地震340次。

由于大地震之后,受主震面波振幅较强、余震大量集中发生等因素影响,余震区的地震监测能力会出现短暂的显著下降(Iwata,2008)。为确保地震序列参数计算的可靠性,这里使用了震级—序号方法(Huang,2006; 蒋长胜,吴忠良,2011; Zhuang et al.,2012)对鲁甸M_S6.5地震序列目录进行完整性分析,结果如图1所示,图中红色五角星代表主震的位置,两条垂直虚线分别标出了截止震级为M_L2.0和M_L1.5时地震目录完整的位置,即分别为C₀=0.031 7 d和C₀=0.251 9 d。在利用震级—序号图考察完备震级时,一方面需考察图中地震事件分布的下包络线,即记录的最小震级; 另一方面,需确保所选定的“完整震级”起始位置之后的地震,在震级—序号图中连续、无缺失地分布,由此方可获得地震序列目录完整的起始位置C₀。由图1可见,完整性震级由主震发生后的M_L2.5逐渐降低至M_L0.8左右,反映了余震序列监测能力逐渐提升的过程。当截止震级M_C=M_L2.0和M_C=M_L1.5时,C₀可分别选定为主震后0.031 7 d和0.251 9 d。

笔者为从不同角度考察鲁甸M_S6.5地震序列参数的早期特征,分别考察ETAS模型的参数。

图1 云南鲁甸MS6.5地震序列目录的震级—序号图
Fig.1 Magnitude-order diagram of the Ludian MS6.5 earthquake sequence

在利用ETAS模型对地震序列参数进行分析时,需在“转换时间”域进行“残差分析”。如在残差分析中,转换时间域的累积频次近似为直线,并与理论曲线较为接近,则拟合效果较好。为测试鲁甸M_S6.5地震序列目录对ETAS模型的拟合情况,图2和图3分别给出了主震之后0～9 d、C₀=0.031 7 d、截止震级M_C=M_L2.0时的残差分析和条件强度曲线情况。

图2a给出了鲁甸M_S6.5地震序列累积地震数与ETAS拟合曲线在“转换时间”域{τ}的比较。由图可见,在“转换时间”域{τ}下,实际地震数目的累积曲线接近一条直线,与理论的稳态泊松过程(图中黑色虚线)较为接近。由此表明,在相应的计算参数设置下,鲁甸M_S6.5地震序列可被ETAS模型较好地描述。在图2b给出的M-τ图可见,在经过转换后,地震序列的时间{t}已被逐个变换为“转换时间”域的时间序列{τ},从{τ}的时间分布特征上,原有的地震序列{t}从震后密集发生的丛集特征已转换为间隔随机分布的稳态泊松过程。

图3a给了ETAS模型的条件强度曲线,即单位时间内地震发生率随时间的变化情况。由于ETAS模型考虑了余震激发次级余震的情况,图中显示了在接近震后5 d和7 d时,次级余震活动所引起的地震发生率的显著变化。图3b显示了序列原始的M-t图。上述图2a和图3a中还标出了在

图2 利用ETAS模型对云南鲁甸MS6.5地震序列的“残差分析”(a)在“转换时间”(τ)域的ETAS拟合曲线与累积地震数的比较;(b)M-τ图
Fig.2 Residual analysis of the Yunnan Ludian MS6.5 earthquake sequence by ETAS model (a)cumulative numbers of earthquakes(the deep blue curve) against ETAS fitting curve in transformed time τn(the thick dashed line);(b)M-τ plot

图3 利用ETAS模型拟合给出的云南鲁甸MS6.5地震序列条件强度曲线(a)地震序列条件强度曲线;(b)M-t图
Fig.3 Temporal evolution of the conditional intensity curve fitted by ETAS model to the Ludian MS6.5 earthquake sequence(a)conditional intensity curve of earthquake sequence;(b)M-t plot

地震序列参数与序列持续时间(蒋长胜等,2013a,b,2014)和截止震级M_C的选取有关(Sornette,Werner,2005; Schoenberg et al.,2010; Wang et al.,2010),为研究云南鲁甸M_S6.5地震序列的早期参数特征,本文分别选用M_C=M_L2.0与M_C=M_L1.5以及相应的C₀=0.031 7 d与C₀=0.251 9 d,并对序列持续时间t={1.5,2.0,2.5,…9.0}进行ETAS模型拟合,计算获得的p值和α值分别如图4和图5所示,图中蓝色和绿色分别给出了M_C=M_L1.5和M_C=M_L2.0的结果,图中误差棒同时给出了p值的误差。

由图4可见,当M_C=M_L2.0时,p值随序列持续时间的变化较为稳定、误差较小,分布在0.89～1.10之间,总体较低,显示了序列较慢的衰减特征。当M_C=M_L1.5时,p值变化范围较大,在0.78～1.43之间。尽管截止震级较低、序列中的地震事件较多,但误差相比M_C=M_L2.0时较大,表明了ETAS模型拟合的相对不稳定性。

由图5给出的α值随序列持续时间的变化可见,在M_C=M_L1.5和M_L2.0时,α值均显示了随序列持续时间的增加而减小的变化趋势,尤其是在震后5 d后出现大幅下降。对于M_C=M_L2.0,α值由3.76减小至1.67,在M_C=M_L1.5时,α值由3.10减小至0.46。α值在震后5 d的剧烈变化,显示了云南鲁甸M_S6.5地震序列激发次级余震能力的显著变化,即由较弱变为相对较强。

图4 ETAS模型p值随着不同截止震级和序列持续时间的变化
Fig.4 The p values of ETAS model varied with the different cutoff magnitude and duration time of sequence

图5 ETAS模型α值随着不同截止震级和序列持续时间的变化
Fig.5 The α values of ETAS model varied with the different cutoff magnitude and duration time of sequence

本文利用ETAS模型,对云南鲁甸M_S6.5地震震后0～9 d的序列参数早期特征进行了考察。分别考虑截止震级M_C=M_L1.5和M_C=M_L2.0以及相应的序列起始时间C₀,对震后序列持续时间t={1.5,2.0,…9.0}天的p值和α值的演化进行了研究,获得了如下方面的认识:

(1)当选取截止震级M_C=M_L2.0时,云南鲁甸M_S6.5地震序列可获得相对稳定、可靠的p值和α值等序列参数。

(2)当M_C=M_L2.0时,p值在震后1.5～9.0天内随序列持续时间的变化较为稳定,且数值较低,p值分布在0.89～1.10之间,显示出较慢的衰减特征。

(3)当M_C=M_L2.0时,α值在震后5.0 d后出现大幅下降,由主震后较高的3.76减小至1.67,表明激发次级余震的能力出现显著变化。

为比较云南鲁甸M_S6.5地震序列参数早期特征,本文对同样采用ETAS模型获得的2013年四川芦山M_S7.0地震(蒋长胜等,2013b)、甘肃岷县—漳县M_S6.6地震(蒋长胜等,2013a)、2014年新疆于田地区3次地震序列(蒋长胜等,2014)参数进行了对比。其中,2013年四川芦山7.0级地震序列在主震后10 d内,p=0.93～1.24,在震后24 d稳定在1.22左右; 2013年甘肃岷县—漳县6.6级地震在震后1.678 d,p=1.15; 2008年新疆于田M_S7.3地震序列在主震后20天内p=1.15～1.30; 2012年新疆于田M_S6.3地震序列在主震后20 d内,p=1.16～1.19; 2014年新疆于田M_S7.3地震序列在主震后3.5 d内p=0.67～0.80。对于α值,2013年芦山7.0级地震序列在主震后10 d内,α=2.47～2.10; 2013年甘肃岷县漳县6.6级地震序列在主震后1.678 d,α=1.70; 2008年新疆于田M_S7.3地震序列在主震后20 d内,α=1.01～1.30; 2012年新疆于田M_S6.3地震序列在主震后20 d内,α=1.30～1.64; 2014年新疆于田M_S7.3地震序列在主震后3.5 d内,α=1.46～1.73。

由上述比较可见,云南鲁甸M_S6.5地震序列的p值仅高于2014年新疆于田M_S7.3地震序列,显示了序列较慢的衰减变化,而这种低p值、慢衰减的特征尽管不能排除后续发生强余震的可能,但受p值可能的区域特征所限,仍需对此次地震周边地区地震序列的比较来判断这种慢衰减是否属于余震活动的异常现象。α值未呈现随序列持续时间平稳变化的趋势,在震后5 d显著减小,表明激发次级余震的能力显著增强,这一点也可从震后5 d发生的4.0级余震相应的次级余震较为丰富上看出。α值的突然变化,可能表明在震后早期,余震区以主震的剩余能量释放为主、序列相对简单,而5天后余震区复杂的应力调整过程激发稍强余震及次级余震。目前,云南鲁甸M_S6.5地震序列仍在持续发展中,目前呈现的p值较低、α值有变化的早期特征尚需进一步跟踪考察。

研究中使用的地震目录由中国地震台网中心“全国地震编目系统”提供,日本统计数理研究所庄建仓教授提供了在访问地球物理研究所期间编制的ETAS模型程序,并予以了指导。本项工作属中国地震局“云南鲁甸6.5级地震专题研究”项目工作,两位审者提出的建设性意见对本文帮助很大。在此一并表示衷心感谢。