0 引言

地震序列是地震活动的主要表现形式之一，是统计地震学的重要研究对象。不同的构造区域、地震序列类型、震源区构造应力水平和区域大地热流值等均可能表现为地震序列参数的差异(Kagan et al，2010)，且余震序列的衰减特征和激发余震的能力可通过序列的统计参数来表征，因此，获得准确可靠的地震序列参数对快速判定地震序列类型、研究震源区特征和评估后续地震危险性等均具有重要的参考价值(蒋海昆等，2007；Guo et al，1997)。目前国际上对地震序列参数的计算主要采用时间序列的“传染型余震序列”(epidemic type aftershock sequence，简称ETAS)模型(Ogata，1988、1989；Zhuang，2011)。蒋海昆等(2006)系统研究了中国大陆5级以上地震序列的统计特征，并从物理上讨论了不同条件下的序列衰减和余震激发问题。蒋长胜等(2013a、2013b、2013c)对1976年唐山地震序列、2013年岷县漳县M_S6.6地震序列和芦山M_S7.0地震序列利用ETAS模型进行了研究。余娜等(2016a、2016b)采用时间序列的ETAS模型计算了2010年玉树M_S7.1、2016年门源M_S6.4地震序列。

青海地区是西北地震活动最为活跃的地区之一，是地震学家高度关注的地震危险区域，历史上多次强震的发生，积累该区域内更多、更可靠的地震序列参数，有利于震后地震趋势估计和余震危险性研究。

1 研究资料 1.1 时间序列ETAS模型

用于拟合地震序列的ETAS模型(Ogata，1988、1989；Zhuang，2011)，其假设所有的余震均可按照大森-宇津公式(Omori，1894；Utsu，1961)激发自己的余震

$f\left(t \right) = \mathit{\mu + K/}{\left({t + c} \right)^p} $

(1)

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

$\mathit{\lambda }\left(t \right) = \mathit{\mu + K}\sum\limits_{{t_\mathit{i}} < t} {\frac{{{\rm{e}}\mathit{\alpha }\left({{M_i} - {M_0}} \right)}}{{{{\left({t - {t_i} + c} \right)}^p}}}}, {M_i} > {M_0} $

(2)

其中，M_i和t_i分别为第i个事件的震级和发生时间；M₀为参考震级，一般可取截止震级。式(2)中的参数p表示地震序列衰减的快慢，p越大衰减越快，反之则越慢；参数α表示触发次级余震的能力(Ogata，1989)，α越大表示触发次级余震的能力越弱，反之亦然。对于震群型序列，一般情况下α＜1，而当地震序列中无明显的被激发的次级余震时，α一般大于1(Ogata，2001)；参数c为数值极小的正常数，一般用于调节当地震序列在主震发生后极短时间内余震记录不完整，以及避免式(2)中分母为0的情况。为确保参数拟合的稳定性，当区域内背景地震发生率不高时，可通过设置背景地震的发生率μ=0来降低拟合难度。

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

本文选取的“青海地区”的空间范围为32.0°~39.0°N，89.5°~103.0°E。该研究区域位于青藏高原中北部及东北缘，地质构造规模大、活动性强，主要有阿尔金断裂带、祁连山-河西走廊活动断裂带、昆仑山断裂带、巴颜喀拉断裂带和唐古拉断裂带等，以NNW、NW-NWW向断裂带最为发育。地震具有频次高、强度大、分布广的特点(马玉虎等，2006；曾秋生，1999)。2009年以来研究区主震震级M_S≥5.0地震有12个，其中，5.0≤M_S＜6.0的有8个，6.0≤M_S＜7.0的有3个，M_S≥7.0的有1个。受ETAS模型计算对地震序列质量的要求，本文使用了记录较好的5级以上地震序列共9个(图 1)，分别为2009年8月28日大柴旦M_S6.4、2010年3月24日西藏聂荣M_S5.8、2010年4月14日玉树M_S7.1、2011年6月26日囊谦M_S5.2、2013年1月30日杂多M_S5.2、2013年9月20日门源M_S5.1、2015年11月23日祁连M_S5.0、2016年1月21日门源M_S6.4和2016年10月17日杂多M_S6.2地震。

本研究所采用的地震目录为中国地震台网中心提供的《全国统一正式目录》。青海地区地震目录存在着明显的空间不均匀现象，说明不同时期、不同地区地震记录水平存在明显的差异。在实际的地震序列选取中，采用了蒋长胜等(2014)区分地震序列与背景地震的时-空分布界限的“自然边界法”，通过纬度-时间图、经度-时间图及震中分布图相结合的方式选取地震序列。图 2给出了2009年大柴旦M_S6.4地震序列的选取结果。依照上面的选取方法，选取了其余8个地震序列。此外，由于在一些强震发生后的短期内，主震的波形振幅较大，面波等波列持续时间较长，可出现将随后大量发生的震级较小的余震“淹没”而导致余震区甚至更大范围内的地震监测能力显著降低(Akaike，1974)的现象。为确保地震序列的完整性，采用“震级-序号”法(Huang，2006；蒋长胜等，2014；Zhuang et al，2012)确定了每个地震序列完整性震级M_c(图 3)。“震级-序号”法按地震发生时间的先后顺序排序，地震密度较大的区域的连线大致为截止震级M_c的时序变化。

2 计算结果 2.1 地震序列ETAS模型拟合情况

ETAS模型参数一般使用最大似然法进行估计。在拟合时间[S，T]范围内(0＜S＜T)，似然函数L的形式为

$\lg \mathit{L = }\sum\limits_{i:S \le {t_\mathit{i}} < T} {\lg \mathit{\lambda }\left({{t_\mathit{i}}} \right)} - \int_S^T {\mathit{\lambda }\left(t \right){\rm{d}}\mathit{t}} $

(3)

将式(2)带入式(3)，即可对序列参数[μ，K，c，α，p]进行最大似然估计。作为示例，图 4给出了2009年海西M_S6.4地震序列在M_c=M_L1.9情况下的ETAS模型参数结果，获得ETAS模型参数为μ=0.0000，K=0.0280，c=0.0035，α=1.2411和p=1.2650。为进一步考察ETAS模型的拟合效果，常使用“转换时间”域的“残差分析”法(Ogata，1988)。“残差分析”方法可将复杂的地震序列{t_i}转换为齐次泊松过程，并在“转换时间”域{τ}中进行拟合效果分析，相应的转换过程为

$t \mapsto \mathit{\tau = \Lambda }\left(t \right) = \int_0^t {\mathit{\lambda }\left(u \right){\rm{d}}\mathit{u}} $

(4)

这种经过转化的{τ}域的数据被称为“残差点过程”(RPP)。若在残差分析中，转换时间域的RPP在{τ}上的累计频次近似直线，并与理论直线较为接近(Zhuang et al，2012)，则ETAS模型拟合效果较好。以2009年大柴旦M_S6.4地震序列为例，图 4(a)给出了主震之后0~91天和M_c=M_L1.9时的“残差分析”结果，可以看出，序列累计地震数与理论ETAS拟合曲线较为接近，ETAS模型对此地震序列的拟合较为理想；图 4(b)为相应的{τ}域的M-τ分布。研究区9个地震序列参数的估算结果见表 1。

2.2 序列参数稳定性分析

ETAS模型参数与截止震级的选取有一定的关系(蒋长胜等，2013a、2013b)。本文以2009年8月28日大柴旦M_S6.4、2010年4月14日玉树M_S7.1地震为例，研究了不同截止震级对地震序列参数的影响。2009年大柴旦M_S6.4地震选取截止震级M_c=M_L1.8、1.9、2.1，发现地震序列参数α值和p值在M_c=M_L2.1时变化较小，误差也小(图 5)。2010年玉树M_S7.1地震选取截止震级M_c=M_L1.5、1.6、…、2.1，当截止震级M_c从M_L1.5到2.1逐渐增大时，α值和p值随截止震级的变化均有较大变化，参数k值逐渐减小，α值总体上有增加的趋势，p值随截止震级的增大而减小，但是变化幅度不大。因此，在考察地震序列参数的早期特征时，须考虑在何种截止震级条件下进行讨论。

由于板内地震所处的构造运动速率不同，故地震序列持续时间存在明显差异。因此，所选取的地震序列的拟合截止时间不同，获得的ETAS模型参数可能不同，拟合参数的标准差也将存在差异(Iwata，2008)。以2010年4月14日青海玉树7.1级地震为例，设定t_end=[1.0，2.0，3.0，…，42.0]，当M_c=M_L1.5时，ETAS模型参数α值和p值均在震后14天出现突跳变化，标准差幅值较大，尤其是激发次级余震能力的α值变化最为剧烈，p值在震后18天内有一定幅度波动。b值变化相对稳定，但在震后14天内有逐渐增加的变化，随着序列持续时间的增加α值减小。地震早期ETAS模型参数的α值、p值变化幅度较大，其后，趋于相对稳定。

3 青海地区地震序列参数的时变特征

地震序列参数的早期变化比较剧烈，反映了地震序列主震发生后应力快速调整的过程，余震也主要集中发生在地震早期阶段。为了研究地震序列参数在时间上的变化特征，本文将地震序列不同截止条件下的ETAS模型参数稳定时间和震级进行统计(表 2)，可以看出，地震序列以不同的截止震级作为研究对象时，地震序列的稳定时间T_α和T_p是不同的，也就是说，在一个地震序列中，地震序列参数α值和p值的稳定时间T_α和T_p与截止震级没有明显的正相关关系，当M_c=M_L1.8时，T_α=19，T_p=19；当M_c=M_L1.9时，T_α=17，T_p=21；当M_c=M_L2.1时，T_α=23，T_p=19。地震序列参数α值和p值的稳定时间T_α和T_p与主震震级间没有明显的相依关系，主震震级较大的其稳定时间也可能较短，震级较小的地震稳定时间反而较长。比如玉树M_S7.1地震，地震序列参数α值和p值的稳定时间为：当M_c=M_L2.0时，T_α=15，T_p=15；杂多M_S5.2地震，当M_c=M_L2.0时，T_α=18，T_p=13；反映了玉树7.1地震震后应力迅速调整的特征。

研究区不同震级档地震序列参数的最长稳定时间和最短稳定时间分别为：地震序列主震震级M_S5.0~6.0，T_α，min=9，T_α，max=47，T_p，min=13，T_p，max=45；主震震级M_S6.0~7.0，T_α，min=17，T_α，max=57，T_p，min=11，T_p，max=26。

4 讨论与结论

为系统考察青海地区地震序列的ETAS模型，并对地震序列的震后早期特征进行分析，本文运用中国地震台网中心提供的《全国统一正式目录》，对青海地区的地震序列进行了筛选，发现青海地区地震目录存在着明显的空间不均匀现象，不同时期、不同地区地震记录水平存在明显的差异。在实际的地震序列选取中采用“自然边界法”选取了研究区2009年以来的9个地震序列，利用“震级-序号”法确定了每个地震序列的完整性震级M_c，使用最大似然法对每个地震序列进行参数估算，获得了如下认识：

(1) 为检测地震序列参数的稳定性，考察了截止震级和拟合截止时间的选取对ETAS模型参数影响。研究发现，截止震级M_c对α、k和p值有一定的影响。在考察地震序列参数的早期特征时，须考虑在何种截止震级下进行讨论；地震早期ETAS模型参数的α值、p值变化幅度较大，反映了地震主震发生后震源区应力的快速调整过程，其后，地震参数的α值和p值趋于相对稳定。

(2) 青海地区地震序列的ETAS模型参数α值、p值的稳定时间T_α、T_p与截止震级间没有明显的正相关关系，与主震震级间也无明显的相依关系，主震震级较大的地震事件稳定时间可能较短，震级较小的地震事件其稳定时间也有可能较长。

(3) 青海地区不同震级档地震序列参数的最长稳定时间和最短稳定时间分别为：地震序列主震震级M_S5.0~6.0，T_α，min=9，T_α，max=47，T_p，min=13，T_p，max=45；主震震级M_S6.0~7.0，T_α，min=17，T_α，max=57，T_p，min=11，T_p，max=26。

由于研究区域强震数目和ETAS模型计算条件所限，仅对9个地震序列进行了估算，并对其时间特征进行了讨论，且所得出的结论仅反映部分统计特征，今后还需要继续积累更多的震例进行更深入研究。

致谢: 研究中使用了中国地震台网中心“全国地震编目系统”提供的“统一正式目录”，日本统计数理研究所庄建仓教授编制、提供了本文所使用的ETAS模型和余震短期概率预测程序，在此一并表示感谢。