2010年青海玉树MS7.1地震序列ETAS模型参数及其变化特征研究①

余　娜,蒋长胜,马玉虎

(1.青海省地震局,青海 西宁 810000; 2.中国地震局地球物理研究所,北京 100081)



2010年青海玉树MS7.1地震序列ETAS模型参数及其变化特征研究①

余娜1,蒋长胜2,马玉虎1

(1.青海省地震局,青海 西宁 810000; 2.中国地震局地球物理研究所,北京 100081)

采用时间序列的“传染型余震序列”(ETAS)模型和最大似然法对2010年4月14日青海玉树MS7.1地震序列的参数进行了估计。为考察其结果的稳定性,设定不同的截止震级和不同的拟合截止时间分别进行序列参数估算。选用截止震级MC=ML1.5对玉树MS7.1地震序列整体的参数进行拟合,计算获得α=0.948 2,p=1.059 6和b=0.817 3。与中国大陆M>7.0地震序列的平均序列参数相比较表明,玉树MS7.1地震序列表现为触发次级余震的能力较弱,序列衰减速率较快。研究表明,截止震级MC对α、k和p值有一定的影响。随截止震级逐渐增大,参数k值逐渐减小,α值总体上有增加的趋势,而对p值的影响较小。在地震序列的早期阶段,ETAS 模型参数中的p值和α值在震后14天内随时间的变化幅度较为明显,其后各参数变化相对平稳。

玉树MS7.1地震; 余震; 地震序列; ETAS模型

0　引言

据中国地震台网中心测定,2010年4月14日7时49分(北京时间)在青海省玉树藏族自治州玉树县发生MS7.1地震。此次地震造成重大人员伤亡,是继2008年汶川8.0级地震后我国大陆发生的又一次灾难性地震。由于玉树MS7.1地震发生在玉树活动断裂带上,也是1997年以来中国大陆7级以上地震集中发生的巴颜喀拉地块的南边界[1],此次地震的发生引起了社会各界的高度关注。

不同的构造区域、地震序列类型、震源区构造应力水平和区域大地热流值等均可能表现为地震序列参数的差异[2],获得准确可靠的地震序列参数对震后序列类型快速判定、强余震预测、断层愈合过程等研究具有重要的参考价值和现实意义[3]。目前国际上对地震序列参数计算主要采用时间序列的“传染型余震序列”(epidemic type aftershock sequence,简称ETAS)模型[4-6]。蒋海昆等[7]系统研究了中国大陆5级以上地震序列的统计特征,并从物理上讨论了不同条件下的序列衰减和余震激发问题。蒋长胜等[8-10]利用ETAS模型对1976年唐山地震序列、2013年甘肃岷县漳县MS6.6地震序列和四川芦山MS7.0地震序列进行了研究。

由于地震序列的复杂性,截止震级的选取、研究区构造的一致性、震后早期阶段震源区应力场的调整等均可能影响到地震序列参数并引起其变化[8-11],这给震后快速、准确地获取地震序列参数带来了较大挑战。从现实意义上积累区域内更多的、更为可靠的地震序列参数,并充分研究其受各种影响因素的变化特征是解决上述问题的重要基础。我国地震观测系统在2009年进入“十五”之后得到飞跃式发展,地震观测能力显著提升,也为获得科学可靠的地震序列参数提供了数据基础。针对上述问题和新情况,本研究以2010年青海玉树MS7.1地震为例,利用ETAS模型考察其序列参数,并探讨其变化特征和影响因素。

1　时间序列ETAS模型

用于拟合地震序列的ETAS模型[4-5,12]假设所有的余震均可按照大森-宇津公式[13-14]激发自己的余震:

f(t)=μ+k/(t+c)p

(1)

式中:f(t)为t时刻对应的地震序列发生率;t为主震发生后的离逝时间;μ为背景地震发生率;p表示余震序列衰减的快慢;c为主震后余震频次达到峰值时所对应的时间;k为常数,用于描述余震的活跃程度。ETAS模型假定主震的发生为初始零时刻,在其后的一个观测时间段[0,T]内的地震序列{(ti,Mi);i=1,2,…,N}的强度函数可表示为[4]:

(2)

其中:Mi和ti分别表示第i个事件的震级和发生时间;M0为参考震级,一般可取截止震级。p与式(1)中p的物理含义相同,p越大序列衰减越快;α表示触发次级余震的能力,对于震群型序列一般情况α<1,而当地震序列中无明显的被激发的次级余震时α一般大于1[12]。

(3)

为考察模型对数据的适用性情况,常使用Akaike信息准侧(AIC)[15]进行判断:

AIC=-2lgL+2k

(4)

AIC值较小表示适用效果较好。

2　研究区域和序列目录的完整性分析

对于青海玉树MS7.1地震序列的选取,在空间范围上设定31.80°～34.45°N,95.10°～98.15°E内的地震事件作为此次地震序列进行研究(图1)。序列的地震事件基本参数来自中国地震台网中心提供的《全国统一地震编目目录》,选用2014年3月10日—5月27日期间1.0级以上的地震用于研究。

图1　玉树地区活动断裂与地震分布Fig.1　Distribution of active faults and earthquakes in the Yushu region

在实际的地震序列选取中,为降低不同算法的人为主观性,参照蒋长胜等[10]的“自然边界法”,采用纬度-时间图、经度-时间图及震中分布图相结合的方式,根据地震时空的自然边界选取地震序列,选取过程如图2所示。此外,由于在一些强震发生后的短期内,主震的波形振幅较大,面波等波列持续时间较长,可出现将随后大量发生的震级较小的余震“淹没”、余震区甚至更大范围内的地震监测能力显著降低[15]的现象。为确保地震序列的完整性,研究中采用“震级-序号”法[17-19]确定玉树地震序列目录的完整性震级MC。震级-序号法按地震发生时间的先后顺序排序,地震密度较大的区域连线大致为MC的时序变化。玉树地震发生后的短时间内,MC为ML2.9左右,其后余震的监测能力逐渐提升,MC逐渐减小。通过设定ETAS模型拟合的起始时间C0=0.026 7 d,避开序列的不完整时段,综合考虑,选取玉树地震序列截止震级MC=ML1.5。

3　地震序列参数估算结果

利用最大似然法对玉树MS7.1地震序列进行ETAS模型参数估计,得到参数μ=0,k=0.040 7,c=0.005 6,α=0.948 2和p=1.059 6。图3(a)为ETAS模型拟合得到的地震序列条件强度曲线,即单位时间内的地震发生率。由图3可知,在主震后4天内序列衰减较为平稳,其后由较大余震引起的条件强度曲线的变化较为明显。

为考察ETAS模型的拟合效果,采用通常使用的“残差分析”[4]方法,将地震序列转换为在“转换时间(τ)”域的分布,并考察实际地震序列和理论值的拟合情况。对条件强度函数λ(t)可采用下式对时间序列{ti}进行转换:

t

(5)

采用式(5)即可将{ti}转化为服从单位速率的稳态泊松分布的{τi}[19]。若ETAS模型对数据拟合较好,在转换时间域的地震累积数目则表现为线性,接近标准稳态泊松过程的理论直线。图4(a)给出了玉树MS7.1地震序列累积地震数与ETAS拟合曲线在“转换时间”域的比较情况,从图中可以看出模型与数据拟合得较好,表现为线性。

图2　玉树MS7.1地震序列目录的完整性分析(图中红色五角星代表MS7.1地震, 蓝色圆点代表用自然边界法获得的地震序列,紫色圆点为背景地震)Fig.2　Catalogue completeness of the Yushu MS7.1 earthquake sequence (The red pentagram indicates the MS7.1 earthquake,blue dots indicates the earthquake sequence obtained by natural boundary method,and purple dots indicates the background earthquake)

图3　ETAS模型给出的玉树MS7.1地震序列ML1.5以上地震的条件强度曲线和M-t图Fig.3　Temporal variation of the conditional intensity from fitting the ETAS model to the Yushu MS7.1 earthquake sequence with cutoff magnitude ML1.5 and M-t plot

由图3和图4可知,对于玉树MS7.1地震序列,在设定的MC=ML1.5和C0=0.026 7 d条件下获得的ETAS模型参数表明,玉树MS7.1地震序列p=1.059 6表现为相对正常的衰减过程。对于α=0.948 2的结果,由于激发高阶余震能力较强的震群型序列一般α<1,而被激发次级余震不明显时α一般大于1[12],此外,相较于2013年芦山7.0级地震序列在主震后10天内α=2.47～2.10[9]、2013年甘肃岷县漳县6.6级地震序列在主震后1.678天α=1.70[8]以及新疆于田地区3次6级以上地震α值也显著大于1.0的情况[10],表明此次玉树MS7.1地震序列展示了较强的激发次级余震的能力。

图4　利用ETAS模型对玉树MS7.1地震序列的拟合Fig.4　Residual analysis results from fitting the ETAS model to the Yushu MS7.1 earthquake sequence

4　序列参数稳定性分析

截止震级的设置对ETAS模型参数估计有一定的影响[8-9,20],而且所用地震序列的长度或拟合截止时间不同获得的ETAS模型参数也可能不同[9]。为考察地震序列参数的稳定性及受影响程度,本研究设定不同的截止震级,即MC=ML1.5,1.6,…,2.1,分别进行最大似然估算。获得的ETAS模型参数如表1所列。从表1可知,当截止震级MC从ML1.5到2.1逐渐增大时,参数k值逐渐减小;α值总体上有增加的趋势,表明随着截止震级的增加触发次级余震的能力逐渐变弱;p值随截止震级的增大而减小,但是变化幅度不大,约为0.11。

表1　不同截止震级下ETAS模型拟合参数

研究中还考察了序列拟合的截止时间tend对序列参数的影响。设定tend=[1.0,2.0,3.0,…,42.0],利用最大似然法分别进行参数估计。同时考察了b值的变化。当截止震级为ML1.5和ML2.0时,估算获得的α值、p值和b值分别如图5所示。由图可见,在震后14天内各参数均有明显的变化,标准差幅值较大,其中α值变化较为剧烈,第14天有一个突变。当截止震级MC=1.5时[图5(a)],震后14天内α值由4.37减小至0.88,其后逐渐稳定在0.89左右;p值也同样出现突跳变化,且在震后18天内变化幅度较大,约为0.85～1.14,其后稳定在1.1;b值无明显突跳,逐渐由0.65增加至0.77,此后稳定在0.79左右。当截止震级MC=2.0时[图5(b)],震后14天内α值由4.07减小至0.89,p值变化范围为0.80～1.14,b值逐渐由0.60增加至0.74;其后α值逐渐减小至0.95,p值稳定在0.98左右,b值相对稳定,保持在0.78左右。

对上述两个不同截止震级下ETAS模型参数随序列拟合截止时间的变化来看,对于玉树MS7.1地震序列,ETAS模型参数α值和p值均在震后14天出现突跳变化,尤其是激发次级余震能力的α值变化最为剧烈,p值在震后18天内有一定幅度波动;b值变化相对稳定,但在震后14天内有逐渐增加的变化。对于序列参数“突然”变化的原因,可从图3(b)的序列活动进行简单分析,发现在主震后的第14天内发生了一个较大的4级以上余震,而此后的序列在余震发生率、次级余震的激发程度上相比震后早期均有较大变化。然而这种序列活动仅仅是一种表象分析,对于物理的震源区应力场调整或者余震破裂机制的调整变化还需要深入研究。

图5　ETAS模型参数α、p和b随拟合截止时间的变化Fig.5　Parameters α,p and b against the end time tend in ETAS fitting

5　讨论与结论

为考察2010年青海玉树MS7.1地震序列参数并研究影响其变化的可能因素,采用ETAS模型进行研究;采用“自然边界法”选取地震序列,并利用最大似然法对不同的地震序列截止震级、序列拟合截止时间下的地震序列参数进行拟合研究,获得了如下认识:

(1)选用截止震级MC=ML1.5、拟合起始时间C0=0.026 7 d的拟合结果表明,玉树MS7.1地震序列ETAS模型参数μ=0.000,k=0.040 7,b=0.817 3,α=0.948 2和p=1.059 6。比较蒋海昆等[3,7]给出的中国大陆M>7.0地震序列的平均序列参数α=0.509,b=0.956和p=0.684,玉树MS7.1地震序列的α值明显高于上述结果,表明序列触发次级余震的能力较弱;b值与上述结果比较接近;p值高于上述结果,表明序列衰减速率较快。

(2)为检测地震序列参数特征的稳定性,设定截止震级MC=ML1.5,1.6,1.7,…,2.1,分别进行参数拟合和标准差估计。研究表明,截止震级MC对α、k和p值有一定的影响。随截止震级逐渐升高,参数k值逐渐减小,α值总体上有增加的趋势,p值减小,但变化幅度不大。

(3)为考察玉树MS7.1地震序列ETAS模型参数的拟合结果随时间变化情况,分别估算了MC=1.5和2.0时的参数变化情况。结果表明,选用MC=ML1.5时α值在震后早期14天由较高的4.37下降至0.88,其后稳定在0.89左右;p值由0.91逐渐增加至1.09,b值由0.65增加至0.77,其后稳定在0.79左右。当截止震级MC=ML2.0时,在震后14天内α值由4.07下降至0.89,p值由0.80增加至1.14,b值由0.60逐渐增加至0.74,其后保持在0.78左右。上述分析表明,玉树MS7.1地震序列在震后早期阶段触发次级余震的能力较弱,在震后14天内各参数均有明显的变化,其后参数变化较平稳。

本研究在考察2010年玉树MS7.1地震序列ETAS模型参数过程中发现,一方面玉树地震序列的参数受到截止震级的影响,但变化幅度不大;另一方面在震后约14天内的早期阶段,参数变化较为剧烈,尤其是描述激发次级余震能力的α值。这些分析结果与已有研究获得的2013年芦山7.0级地震[8]、甘肃岷县漳县6.6级地震[9]和2014年新疆于田7.3级地震[11]的结果较为接近,但具体参数的变化方式不同。上述结果,尤其是ETAS模型参数在震后早期阶段的可能变化,对于序列参数快速获取并用于序列类型和震后趋势判定等工作都带来重要警示,即需要谨慎使用序列参数并动态跟踪其变化。

致谢:本文是在中国地震局地球物理研究所访问学习期间完成,“国际地震可预测性合作研究”(CSEP)计划中国检验中心筹备组对本研究给予了指导。研究中使用了中国地震台网中心“全国地震编目系统”提供的“统一正式目录”,日本统计数理研究所庄建仓教授在访问地球物理研究所期间编制、提供了ETAS模型和余震短期概率预测程序,并予以了指导,在此表示感谢。

References)

[1]温燕林,宋治平,赵文舟,等.东昆仑断裂带玛曲—玛沁段大震危险性分析和探讨[J].地震工程学报,2015,37(1):175-180.WEN Yan-lin,SONG Zhi-ping,ZHAO Wen-zhou,et al.Analysis and Discussion of Large Earthquake Risk along the Maqu—Maqin Segment of the East Kunlun Fault Zone[J].China Earthquake Engineering Journal,2015,37(1):175-180.(in Chinese)[2]Kagan Y Y,Bird P,Jackson D D.Earthquake Patterns in Diverse Tectonic Zones of the Globe[J].Pure Appl Geophys,2010,167(6/7):721-741.

[3]蒋海昆,郑建常,吴琼,等.传染型余震序列模型震后早期参数特征及其地震学意义[J].地球物理学报,2007,50(6):1778-1786.JIANG Hai-kun,ZHENG Jian-chang,WU Qiong,et al.Earlier Statistical Features of ETAS Model Parameters and Seismological Meanings[J].Chinese J Geophys,2007,50(6):1778-1786.(in Chinese)

[4]Ogata Y.Statistical Models for Earthquake Occurrences and Residual Analysis for Pointprocesses[J].J Amer Statist Assoc,1988,83(401):9-27.

[5]Ogata Y.Statistical Model for Standard Seismicity and Detection of Anomalies by Residualanalysis[J].Tectonophysics,1989,169(1/2/3):159-174.

[6]Zhuang J C.Next-day Earthquake Forecasts for the Japan Region Generated by the ETAS Model[J].Earth Planets Space,2011,63(3):207-216.

[7]蒋海昆,曲延军,李永莉,等.中国大陆中强地震余震序列的部分统计特征[J].地球物理学报,2006,49(4):1110-1117.

JIANG Hai-kun,QU Yan-jun,LI Yong-li,et al.Some Statistic Features of Aftershock Sequences in Chinese Mainland[J].Chinese J Geophys,2006,49(4):1110-1117.(in Chinese)

[8]蒋长胜,庄建仓,龙锋,等.2013年芦山MS7.0地震序列参数的早期特征:传染型余震序列模型计算结果[J].地震学报,2013,35(5):661-669.

JIANG Chang-sheng,ZHUANG Jian-cang,LONG Feng,et al.Statistical Analysis of ETAS Parameters in the Early Stage of the 2013 Lushan MS7.0 Earthquake Sequence[J].Acta Seismologica Sinica,2013,35(5):661-669.(in Chinese)

[9]蒋长胜,吴忠良,韩立波,等.地震序列早期参数估计和余震概率预测中截止震级MC的影响:以2013年甘肃岷县漳县6.6级地震为例[J].地球物理学报,2013,56(12):4048-4057.

JIANG Chang-sheng,WU Zhong-liang,HAN Li-bo,et al.Effect of Cutoff Magnitude MCof 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[J].Chinese Journal of Geophysics,2013,56(12):4048-4057.(in Chinese)

[10]蒋长胜,韩立波,郭路杰.新疆于田地区2008年以来3次地震序列参数的早期特征[J].地震学报,2014,36(2):165-174.

JIANG Chang-sheng,HAN Li-bo,GUO Lu-jie.Parameter Characteristics in the Early Period of Three Earthquake Sequences in the Yutian,Xinjiang Since 2008[J].Acta Seismologica Sinica,2014,36(2):165-174.(in Chinese)

[11]Wang Q,Jackson D D,Zhuang J C.Missing Links in Earthquake Clustering Models[J].Geophys Res Lett,2010,37:L21307.doi:10.1029/2010GL044858.

[12]Ogata Y.Increased Probability of Large Earthquakes near Aftershock Regions with Relative Quiescence[J].J Geophys Res,2001,106(B5):8729-8744.

[13]Omori F.On Aftershocks of Earthquakes[J].J Coll Sci Imp Univ Tokyo,1894,7:11-200.

[14]Utsu T.A Statistical Study of on the Occurrence of Aftershocks[J].Geophys Mag,1961,30:521-605.

[15]Akaike H.A New Look at the Statistical Model Identification[J].IEEE Trans Automat Control,1974,AC-19:716-723．

[16]Iwata T.Low Detection Capability of Global Earthquakes After the Occurrence of Large Earthquakes:Investigation of the Harvard CMT Catalogue[J].Geophys J Int,2008,174(3):849-856.

[17]Huang Q.Search for Reliable Precursors:A Case Study of the Seismic Quiescence of the 2000 Western Tottori Prefecture Earthquake[J].J Geophys Res,2006,111(B4):B04301.doi:10.1029/2005JB003982.

[18]蒋长胜,吴忠良.玉树MS7.1地震前的中长期加速矩释放(AMR)问题[J].地球物理学报,2011,54(6):1501-1510.

JIANG Chang-sheng,WU Zhong-liang.Intermediate-term Medium-range Accelerating Moment Release(AMR)Priori to the 2010 Yushu MS7.1 Earthquake[J].Chinese J Geophysics,2011,54(6):1501-1510.(in Chinese)

[19]Zhuang J,Harte D,Werner M J,et al.Basic Models of Seismicity:Temporal Models[G]//Community Online Resource for Statistical Seismicity Analysis,2012.doi:10.5078/corssa-79905851.Available at http://www.corssa.org.

[20]Schoenberg F P,Chu A,Veen A.On the Relationship between Lower Magnitude Thresholds and Bias in ETAS Parameter Estimates[J].J Geophys Res,2010,115(B4):B04309.doi:10.1029/2009JB006387.

Analysis of ETAS Model Parameters and Associated Variation Characteristics for 2010 Yushu,Qinghai MS7.1 Earthquake Sequence

YU Na1,JIANG Chang-sheng2,MA Yu-hu1

(1.Earthquake Administration of Qinghai Province,Xining 810000,Qinghai,China;2.Institute of Geophysics,CEA,Beijing 100081,China)

The maximum likelihood method is used and the epidemic-type aftershock sequence (ETAS)model is fitted in the analysis of aftershock sequence characteristics of the April 14,2010 Yushu (MS7.1)earthquake that occurred in Qinghai province.To examine parameter stability,ETAS parameters and their standard errors are estimated with different cutoff magnitudes,MC,and different ending times for the fitting interval.The overall parameters are observed to beα=0.948 2,p=1.059 6,and b=0.817 3 with setting MC=ML1.5. Compared to other M>7.0 earthquake sequences on the continental China,the Yushu MS7.1 earthquake sequence is characterized by a weak triggering capability in generating secondary aftershocks,and a quick decay rate of aftershocks.A systematic approach for using multiple cutoff magnitudes above the magnitude of completeness is adopted in the calculation process to examine stability of the sequence parameters and to analyze the degree to which they are affected.The cutoff magnitude,MC,is observed to have a certain effect on the value ofα,k,and p (with an increase in MCthe value of k decreases,there is a general increasing trend ofα,but less influence on p).ETAS model parameters with different sequence durations are studied in this paper,together with their temporal variations and stability in the early stage after the shock.There is variation in the range of p andαshortly after the earthquake stops for different cutoff magnitudes within the same earthquake,whereas there is a trend of p andαis towards convergence and stability after stabilization of seismic activity.Furthermore,there is evident change in the temporal variation of ETAS parameters 14 days after the mainshock,when the parameters become relatively stable.

Yushu MS7.1 earthquake; aftershock; earthquake sequence; ETAS model

2015-07-24

2016年度震情跟踪定向任务(2016010133);青海省科技厅基础研究项目(2013-Z-710)

余娜(1984-),女,工程师,从事数字地震资料的应用研究。E-mail:yuna2003@126.com。

蒋长胜(1979-),男,研究员,主要从事数字地震学和地震预测研究。E-mail:jiangcs@cea-igp.ac.cn。

P315.8

A

1000-0844(2016)04-0609-07

10.3969/j.issn.1000-0844.2016.04.0609