地球物理学报  2019, Vol. 62 Issue (11): 4300-4312   PDF    
华北地区地震序列参数的分布特征
毕金孟1, 蒋长胜2     
1. 天津市地震局, 天津 300201;
2. 中国地震局地球物理研究所, 北京 100081
摘要:为系统地考察华北地区地震序列参数的分布特征,以期构建适合区域地震活动特征的短期概率预测模型和评估地震危险性,本文利用当前较为前沿的ETAS模型和R-J模型,采用连续滑动、多时段拟合的方式,对华北地区1970年以来的16个地震序列进行了参数拟合,并对参数的整体情况、参数之间的关系、参数与大地热流之间的关系等特征进行了分析研究.研究结果表明,主震发生后的早期阶段,地震序列参数变化相对较为剧烈,误差也较大;两种地震预测模型的序列参数呈现出一种优势分布特征,主要参数的平均值分别为αETAS=1.7404±0.3420,pETAS=0.9769±0.1396,aOML=-1.6638±0.5284,bOML=0.8312±0.1658,pOML=0.9053±0.1527,这与国际上其他区域的研究具有很强的一致性;大多数情况下,ETAS模型参数pETAS高于R-J模型参数pOML,平均偏高0.0716;序列参数αETAS与大地热流值整体上呈现一种负相关关系,激发次级余震的能力与大地热流具有一定的相关性.地震序列参数的分布特征以及与大地热流之间的关系对地震序列类型的判断,及基于地震序列参数构建定量化的短期概率预测模型,均有重要的应用价值.
关键词: 地震序列      ETAS模型      R-J模型      模型参数      大地热流     
Distribution characteristics of earthquake sequence parameters in North China
BI JinMeng1, JIANG ChangSheng2     
1. Tianjin Earthquake Agency, Tianjin 300201, China;
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract: The purpose of this work is to systematically examine the distribution characteristics of earthquake sequence parameters in North China, which is helpful in establishing a short-term probability prediction model and assessing seismic hazard suitable for regional seismicity. Using currently advanced ETAS and R-J models and continuous sliding and multi-period fitting, the parameters for 16 earthquake sequences since 1970 in North China are fitted, and some parameters characteristics are analyzed such as the overall state of the parameters, the relationship between the parameters and the relationship between the parameters and the terrestrial heat flow. Results show intense changes and large errors in the early stage after the occurrence of the main shock. The sequence parameters of the two earthquake forecasting models exhibit a dominant distribution characteristic. The average values of the main parameters are αETAS=1.7404±0.3420, pETAS=0.9769±0.1396, aOML=-1.6638±0.5284, bOML=0.8312±0.1658, pOML=0.9053±0.1527, which are highly consistent with those in other regions of the world. In most cases, the ETAS model parameter pETAS is higher than the R-J model parameter pOML, with an average difference 0.0716. Sequence parameter αETAS and heat flow are of a negative correlation and the capability of stimulating secondary aftershocks have certain correlation with terrestrial heat flow. Understanding the distribution characteristics of earthquake sequence parameters and their relationship with the heat flow would help to better determining the type of earthquake sequence and constructing a quantitative short-term probability forecasting model based on earthquake sequence parameters.
Keywords: Earthquake sequence    Epidemic type aftershock sequence (ETAS) model    Reasenberg-Jones model    Model parameter    Terrestrial heat flow    
0 引言

地震序列参数尤其是参数的早期特征,包含了大量的震源以及周围构造等方面的信息,可对序列类型快速判定(Guo and Ogata, 1997; 蒋海昆等, 2007)、研究震源区构造应力水平和区域大地热流值(Kagan et al., 2010)以及理解地震过程和评估后续地震危险性均具有重要的参考价值.在实际研究和应用中,目前尚存在地震序列参数计算的科学性、区域研究的系统性以及震后早期阶段序列参数的稳定性等系列问题.上述问题的解决,对地震序列参数进一步在地球动力学和断层物理学等研究中提供一定的科学依据.

大震之后的余震活动信息可以很好的在地震序列参数中体现出来,而余震的活动特点是影响震后救灾最重要的因素之一.通过国际“地震可预测性合作研究”(CSEP)计划近10年的发展以及“可操作的地震预测”(Operational Earthquake Forecast,简写为“OEF”)等概念的提出,地震预测结果的效能评价技术逐渐得到完善,对地震的可预测属性的认识也得到了极大提高(Schorlemmer et al., 2018).如Han等(2017)提出的一种结合地球物理观测中的预测信息和地震目录中的预测信息的混合模型.而正在开展的CSEP计划中“向前”预测的余震预测模型包括R-J模型(Reasenberg and Jones, 1989)、ETAS模型(Ogata, 1988)等都是基于序列参数构建的,而基于这些定量化的短期概率预测模型,开展的地震预测及预测策略研究(Guo et al., 2015; 蒋长胜等, 2015; 2017; 2018;毕金孟和蒋长胜, 2017; Nicolis et al., 2017; Ogata, 2017; Taroni et al., 2018),成为震后有效的指导科学救灾的重要依据.

蕴涵大量地球物理、地球动力学以及地质信息的大地热流是地球内部热动力在地表上最直接的显示(Sclater et al., 1980; Furlong and Chapman, 1987; Pollack et al., 1993),反映了岩石圈内部的热平衡状态,此外,大地热流的分布与岩浆活动、构造以及地震活动性等都密切相关(Chapman and Rybach, 1985).中国大陆地区大地热流测点分布极不均匀且存在大部分空区,总体呈现东高、西低,西南高、西北低的整体布局(胡圣标等, 2001; 汪集旸等, 2012),实测的热流值变化范围为23.4~319 mW·m-2,算术平均值为61.5±13.9 mW·m-2,东部整体热流值在40 mW·m-2以上,平均热流值在60 mW·m-2左右,中部较低,在40~60 mW·m-2之间(姜光政等, 2016).Enescu等(2009)研究了美国南加州地区的大地热流值与地震序列参数的分布关系后发现,大地热流与ETAS模型参数的α值呈现一种负相关关系,表明大地热流值越大的地区,地球内部活动性相对偏高,更容易激发出自己的次级余震.此外,一些研究中发现(Mogi, 1962; Kisslinger and Jones, 1991),p值与热流值存在正相关关系,并被解释为高温度对断层区介质应力释的加速作用.这些模型参数与大地热流值之间关系对于该区域余震模型的构建、地震危险性评估以及余震预测策略的研究具有很好的指示意义.

华北地区作为我国经济、政治、文化中心所在地,同时,华北地区的地震活动强烈,该区历史上发生过多次中强地震.自1970年以来,该区域发生6.0级以上地震17余次,包括1975年2月4日辽宁海城的MS7.3、1976年7月28日河北唐山的MS7.8、1998年1月10日河北张北MS6.2等多次大震,造成了重大的人员伤亡和财产损失.在华北地区,尚未系统的开展过有关序列参数总体特征以及与大地热流的相关研究,本文针对上述问题,采用ETAS模型和R-J模型,利用最大似然法对所选取的华北地区4.5级以上16个地震序列参数进行拟合,得到稳定的地震序列参数,并与大地热流进行对比分析.此次利用两种短期预测模型统计其序列参数特征,这对于构建适合该区域“可操作的余震预测”(Operational Aftershock Forecasting,简称OAF)模型,具有较好的科学价值.

1 研究区和所用数据

本次研究选用了人口稠密,经济文化发达而构造活动复杂的华北地区(34°N—43°N,109°E—124°E)为研究区,该区是在太平洋板块、印度板块和欧亚板块的相互作用下形成的仍在强烈活动的块体(李三忠等, 2011).华北地区上分布着张渤地震带、河北平原地震带、山西地震带、郯庐地震带等复杂的地震带,且该区构造活动强烈、断裂和褶皱发育、地震活动频繁,给人民的生命财产造成了极大的损失.受所选模型对地震序列质量的要求,本研究中使用了华北地区1970年以来的记录较好的4.5级以上地震序列.

研究中使用中国地震台网中心提供1970年1月1日到2018年03月31日期间的《全国统一正式编目》地震目录(http://10.5.160.18/uniteWeekCatalog/index.action?bg=true&menu_index=1,2018-04-12).此次使用Gardner-Knopoff方法(Gardner and Knopoff, 1974)对完整性震级Mc=3.0的地震目录进行筛选,参与计算的地震目录为10145条,经过余震删除,此时地震目录2629条.为考察华北地区的地震序列参数的空间分布与大地热流的关系,确保地震序列在空间上覆盖充分,需要选用足够的地震序列,这里设定主震震级满足MS≥4.5即为初步的选取对象,经G-K删除余震后的主震为159个.为选定各地震序列事件,通过纬度-时间图、经度-时间图以及震中分布图三者相结合的方式,即“自然边界法”来选取地震序列展开研究.此外,根据ETAS模型和R-J模型的计算要求,地震序列中需一定量的余震事件并满足统计显著性要求,参照毕金孟和蒋长胜(2017)的选取规则,初步的地震序列选取过程中要求地震事件不低于60个,最终在进行自然边界法和满足完整性震级以上的条件下,地震序列中的地震事件数据不低于30个才进行参数拟合.图 1给出了1981年8月13日内蒙古丰镇MS5.6地震序列的选取示意图,根据上述选取规则,共选取了华北地区的16个地震序列,这些地震序列的主震参数如表 1所示.

图 1 1981年8月13日内蒙古丰镇MS5.6地震序列的选取示意图 (a)纬度-时间图;(b)经度-时间图;(c)震中分布图;(d)震级-序号图.图中黄色五角星代表主震,黑点代表所选地震序列的余震事件,灰点代表未作为序列事件的“背景地震”.图上水平虚线和垂直虚线分别标出了McC0的位置. Fig. 1 Schematic diagrams showing selection of 13 August 1981 Fengzhen, Inner Mongolia, MS5.6 earthquake sequence (a) Latitude-time; (b) Longitude-time; (c) Epicenter distribution; (d) Magnitude-number. The yellow star is the main shock. Black dots are aftershocks in the selected earthquake sequence. Grey dots are "background earthquakes" in the earthquake sequence. Horizontal dashed lines indicate the cutoff magnitude, Mc, while the vertical dashed lines indicate the starting time of model fitting, C0.
表 1 华北地区地震序列的主震参数及模型参数值 Table 1 Main shock parameters and model parameters of earthquake sequences in North China

因每个地震序列的独特性,也为了保证足够的余震事件参与计算,此次并没有设定“固定”的完整性震级Mc.分别使用“震级-序号”法(Huang, 2006; 蒋长胜和吴忠良, 2011; Zhuang et al., 2012)确定所选地震序列的完整性震级Mc.“震级-序号”法将地震事件按发震时间的前后进行排序,然后根据震级和序号的分布情况来确定Mc.此外,在参数拟合过程中,需同时设定两个模型拟合的起始时刻C0McC0的选取具有相关性,一般Mc越大,C0相对越小,通过调节两个参数的关系,找出适合每个地震序列的McC0值,以保证足够多的地震参与到拟合和计算中,每个地震序列McC0值的如表 1所示.

2 研究所用方法 2.1 时间ETAS模型及其参数估计

与以往对产生余震的传统认识方式不同,Ogata(1988)进一步的将自相似思想引入,认为无论是背景地震还是触发地震都可以激发原震的余震,一次主震事件后能产生大量的直接余震以及间接余震(余震的余震),从而构建具有分支点过程特征的“传染型余震序列”ETAS模型.其条件强度函数可表示为

(1)

其中τi=t-ti代表地震事件i的离逝时间,KETAS为归一化的常数,决定了有Mi事件直接触发余震的预期地震数目,αETAS代表地震事件激发次级余震大小的一种能力(Ogata, 1989; 1992),相对于孤立型和主余型地震,震群型地震序列的αETAS更小一些,一般αETAS < 1(Ogata, 2001),pETAS代表地震序列的衰减程度,pETAS值的大小与序列的衰减快慢呈正相关关系.参数cETAS用于调节主震发生后极短时间内余震记录的不完整性,是一个数值极小的正常数,也避免了分母出现0的情况.参数μETAS表示场地本身的地震活动性,即背景地震的发生率,在计算过程中,为更好的确保参数拟合的稳定性,当区域内背景地震发生率相对不高时,一般通过设置μETAS=0来降低拟合的难度.

利用最大似然法(MLEs)对ETAS模型中的参数进行估计,似然函数L的表示为

(2)

将(1)式代入(2)式,对待估参数[μETAS, KETAS, cETAS, αETAS, pETAS]进行最大似然估计.通过反演多元函数二阶偏导数的Hessian矩阵获得相应参数的标准差(Ogata, 1978).在获得5阶Hessian矩阵后,将2倍的各对角线元素的平方根与参数最大似然值的平方根进行乘积,从而获得5个参数的标准差.

2.2 Reseanberg-Jones(R-J)模型及其参数估计

对于简单地震序列的余震活动特征,常用“Omori- Utsu”公式(Omori, 1894; Utsu, 1961)进行描述,很多模型用此模型预测余震的发震时间,表示为

(3)

而描述震级频度关系的G-R定律(Gutenberg and Richter, 1944)可以很好的用来预测余震的发震数目,地震序列的频次N与震级M的关系可表示为

(4)

Reasenberg和Jones(1989)在“Omori-Utsu”公式和G-R定律的基础上发展了余震短期概率预测的Reseanberg-Jones(R-J)模型,来开展余震短期发生率的预测.因而,在地震序列中t时刻主震震级为Mm的、震级不低于M的余震发生率可以表示为

(5)

其中,参数KOMLcOMLpOML可利用最大似然法对大森-宇津公式拟合获得(Ogata, 1983),并可求取相应的标准差.相应的对数似然函数可以表示为

(6)

对上式进行最大似然估计获得参数KOMLcOMLpOML后,同时,利用最大似然法的渐进式,也即“费舍尔信息矩阵”(Fisher information matrix)估算这些参数的标准差δKOMLδcOMLδpOML.

bOML值利用下式最大似然法(Utsu, 1965; Aki, 1965)求得

(7)

其中M的为包括截止震级M0以上地震事件的平均震级,ΔMbin为地震目录的震级滑动宽度,一般为ΔMbin=0.1.

bOML值的标准差,可采用Shi和Bolt(1982)给出的方法计算:

(8)

(8) 式中的n为参与计算的地震样本数量.对于aOML值,需要利用KOMLbOML值进行计算推导出:

(9)

aOML值的标准差可利用误差传递理论计算得出:

(10)

通过以上关系即可推算出R-J模型中的参数KOMLcOMLpOMLaOMLbOML值以及这些参数的标准差δKOMLδcOMLδpOMLδaOMLδbOML.

2.3 序列参数的稳定性分析

利用地震活动模型进行余震序列类型判定和短期概率预测时,稳定的序列参数将更有助于判断序列类型以及进行余震发生率、余震概率预测等方面的工作(蒋长胜等, 2013; 2015;2017;2018).对所选的16个地震序列进行了连续、滑动的拟合,根据序列之间的差异,对每个序列选取了不同拟合时间窗的起始时刻C0,但拟合时间窗结束时刻/序列持续时间t2统一为从主震后0.50天开始,以0.05天步长直至序列结束,进行了多个时段的滑动拟合.以1976年7月28日河北唐山MS7.8地震序列为例,拟合获得的αETAS值、pETAS值、aOML值、bOML值以及pOML值随序列持续时间的变化如图 2所示.在地震序列早期,地震序列参数αETAS值、pETAS值、aOML值、bOML值以及pOML值均变化较为剧烈,相应的标准差的数值范围也较大,显示了较不稳定的拟合状态和序列发生的复杂性,其他地震序列也表现出了相同的这种特性.因此,将序列参数用于震后地震序列类型快速判断、地震预测等研究中需谨慎对待,早期序列参数的剧烈变化也正反映了地震主震发生后震源区应力的快速调整过程.

图 2 1976年07月28日河北唐山MS7.8地震序列模型拟合参数随序列持续时间的变化 (a)参数αETAS值随序列持续时间的变化;(b)参数pETAS值随序列持续时间的变化;(c)参数aOML值随序列持续时间的变化;(d)参数bOML值随序列持续时间的变化;(e)参数pOML值随序列持续时间的变化.序列截止震级Mc=ML4.0,起始时刻固定为C0=0.0132天,序列持续时间t2以0.05天步长,自主震后0.50天增加至72.00天进行滑动拟合.图中灰色区域为相应的参数的标准差的范围 Fig. 2 Model fitting parameters varying with duration time of the sequence of the 28 July 1976 Tangshan MS7.8 earthquake. (a) αETAS; (b) pETAS; (c) aOML; (d) bOML; (e) pOML The cut-off magnitude is Mc=ML4.0 in model fitting, the starting time fixed by C0=0.0132 day and the ending time from 0.50 day to 72.00 days by sliding with a 0.05 day step. The gray area is the standard deviation range of the parameters shown in each sub-graph.
3 华北地区地震序列参数的分布特征

地震序列参数是震后趋势研判和余震危险性分析的重要研究对象,对断层带结构、区域地球动力学以及构造演化研究等均有一定的参考作用.地震序列参数是地震活动特征在一定时间、空间范围内的集中反映,包括地震发生区域内的地球动力学信息,例如地震序列参数的空间分布与区域深部介质环境(蒋海昆等, 2006)以及大地热流的分布存在一定的关系.此外,大震后震源周边构造场的应力调整、主震的破裂特征、主震所在断层的愈合速度等诸多方面也会使反映地震序列的衰减特征、激发余震能力的参数值等出现明显差异.

3.1 华北模型参数的总体分布特征

前人对参数特征的研究很早就已经开展,如Utsu(1961)应用修正的的Omori公式确定了不同地区的51个余震序列的pOML值,得到pOML在0.9~1.8之间,其中以1.1~1.4最多.Utsu等(1995)对1962—1995年全球已发生的200多个计算结果分析,pOML值分布在0.6~2.5之间,均值为1.1.赵志新等(1992)的研究指出,中国大陆的32个序列的pOML值在0.63~1.54之间,均值为0.95.日本13次地震序列bOML值的数值分布范围为0.51~1.33,中值为0.83;南加州10次地震序列bOML值的数值分布范围为0.46~1.00,中值为0.82;希腊10次地震序列bOML值的分布范围为0.56~1.36,中值为0.82.蒋海昆等(2007)给出了中国大陆294次MS≥5.0中强震余震序列参数的ETAS模型拟合结果以及对中国不同地区、不同主震破裂类型、不同主震震级范围、不同类型序列的αETAS值、pETAS值和bETAS值等早期参数特征.Ogata等(2018)基于G-R关系的b值(b=0.9)对三种不同预测未来地震震级大小的震级预测模型进行了探讨和评价.

此次利用ETAS模型和R-J模型对所选的16个地震序列进行研究,在序列参数达到稳定时,序列参数的总体特征如图 3所示,总体呈现一种优势分布.此次,利用箱线图(Boxplot)对各个序列参数进行了统计,由图 3可知,各个参数分布相对较为集中,具有很好的对称性,五个序列参数统计中均未出现异常值,表明华北地区的序列参数具有较好的一致性.序列中主要参数αETAS值、pETAS值、aOML值、bOML值、pOML值的统计特征值,包括最大值、最小值、平均值、中值以及标准差如表 2所示,主要参数的平均值分别为αETAS=1.7404±0.3420,pETAS=0.9769±0.1396,aOML=-1.6638±0.5284,bOML=0.8312±0.1658,pOML=0.9053±0.1527,总体来看,此次研究结果比蒋海昆等(2007)对中国大陆294次MS≥5.0中强震余震序列参数的ETAS模型拟合结果偏大,这可能与余震序列长度的选取有关.此次研究所得到的结果与前人在全球范围内得到的结果具有可比性,如Reasenberg和Jones(1989)研究了加利福尼亚州1933—1987年间62次M≥5.0级以上地震序列,参数平均值为aOML=-1.76±0.07,bOML=0.90 ±0.02,pOML=1.07±0.03;Lolli和Gasperini(2003)利用意大利1981—1996年的20个M≥4.0的地震序列,得到了意大利地区的参数平均值为aOML= -1.66±0.72,bOML=0.96±0.18,pOML=0.93±0.21.

图 3 华北地区1970年以来16次MS≥4.5地震序列的模型参数分布特征 图中利用了“箱线图”对各个参数数据进行了统计,红、蓝五条竖实线依次代表了各个序列参数的最小值、第一四分位数、中位数、第三四分位数与最大值. Fig. 3 Distribution characteristics of model parameters for 16 MS≥4.5 earthquake sequences in North China since 1970 The "Box-plot" is used to make statistics on each parameter data. The five vertical solid lines represent the minimum, first quartile, median, third quartile and maximum values of each sequence parameter, respectively.
表 2 华北地区地震序列参数整体统计特征 Table 2 Statistical characteristics of seismic sequence parameters in North China

通过统计这些模型参数,无论对进一步构建单一的短期概率预测模型,还是国际较为前沿的“混合模型”都具有一定的潜在应用价值.Reasenberg和Jones(1989)根据模型参数中值(a=-1.67,b=0.91,p=1.08,c=0.05)构成了“普通加州模型”的概率预测模型,用于强余震和大地震的概率预测.

3.2 华北地区地震序列参数之间的对比分析

为了探究模型参数之间的相互关系,这里分别考察了pETAS值与pOML值之间的关系,如图 4所示,图中黑色虚线为pETASpOML的等值线.在本文所选用的16个地震序列中,pETAS>pOML的序列为12个,占75%,表明ETAS模型参数pETAS在数值上总体高于R-J模型参数pOML,与Guo和Ogata(1997)利用ETAS模型研究的34个地震的余震序列所得到的结论一致.这可能与ETAS模型考虑了高阶余震活动、pETAS包含了更多的次级余震衰减信息.由图 4还可见,pETAS值与pOML值与主震震级无明显的相关性,也与前人研究相一致(Utsu et al., 1995; Lomnitz, 1966; Papazachos, 1974).此外,其他序列参数与主震震级也不存在明显的相关性.

图 4 ETAS模型参数pETAS与R-J模型参数pOML对比 图中黑色虚线为pETAS=pOML的等值线,颜色的深浅代表主震震级的大小. Fig. 4 Comparison of ETAS model parameter pETAS and R-J model parameter pOML The black dotted line is the contour of pETAS=pOML.Colors represent the magnitude of the main shock.
3.3 华北地区地震序列模型参数与大地热流的关系

地震活动性的背后都有物理机制的作用,大量观测体现的统计特征往往是区域地壳岩石的物理化学条件的反映,例如激发次级余震的能力与热流值呈正相关关系(Enescu et al., 2009; Yang and Ben-Zion, 2009; Zaliapin and Ben-Zion, 2013).Hu等(2000)总结发布了截至1999年的中国大陆热流值结果,其中渤海湾盆地热流值为44~106 mW·m-2,均值69±13.3 mW·m-2;鄂尔多斯地块热流值在35~72 mW·m-2之间,均值60±8.7 mW·m-2;汾渭地堑热流值为63~98 mW·m-2,均值73± 9.2 mW·m-2;郯庐断裂带热流值44~79 mW·m-2,均值69±9.8 mW·m-2;太行山—燕山山体热流值30~79 mW·m-2,均值50±11.3 mW·m-2,具体大地热流的分布情况如图 5a所示.考虑到部分主震震中附近大地热流数据比较稀疏,将震中为圆心、半径200 km的圆形区域内的平均值作为震中的大地热流值.获得的震中附近地区大地热流值如表 1中“HF栏”所示,相应的数值范围为43.2429~ 64.3312 mW·m-2,平均值为56.1307±8.1207 mW·m-2.华北地区地震序列的αETAS数值分布范围为1.2415~2.3568,空间分布如图 5b所示.为考察次级余震激发能力的大地热流值的关系,将αETAS与大地热流值进行了对比和线性回归分析,如图 5c所示.图中结果显示,αETAS与大地热流呈现较弱的负相关趋势,与Enescu等(2009)给出高热流地区地震序列激发次级余震的能力更强的认识较为一致.然而由于本研究未涉及α值更小的震群型地震序列、α值变化范围较小,上述的负相关关系仍需要补充更多的序列震例进一步确认.

图 5 地震序列稳定时序列参数αETAS与大地热流(q)的关系分布图 (a)大地热流值得空间分布图;(b) αETAS值空间分布图;(c) αETAS值与大地热流值的关系图.图(a)颜色的深浅代表大地热流值的大小,大地热流数据来源于姜光政等(2016);图(b)中不同大小的矩形代表主震震级的大小,不同的颜色代表地震序列参数稳定时的模型参数值的大小;图(c)中黑色实线代表了数据的拟合情况,虚线代表了一倍的均方差. Fig. 5 The relationship between the stable sequence parameter αETAS and the terrestrial heat flow (a) Spatial distribution map of heat flow; (b) Spatial distribution map of the αETAS value; (c) Parameter αETAS of the ETAS model versus the heat flow. The shade of the color in (a) represents the size of the heat flow. The data of the terrestrial heat flow are from Jiang et al. (2016). The different sizes of rectangles in (b) represent the size of the main shock, and different colors represent the size of the stability model parameters. The black solid line in (c) represents the overall situation of data fitting, and the dashed line represents the 1-fold standard deviation.
4 讨论和结论

在大地热流值与地震序列参数的相关性讨论中存在两个主要的问题(Kisslinger, 1996):首先,地表热流值能否准确代表震源深度处的温度值是不一定的;其次,余震区地表热流值的观测较少,利用一定区域内的平均值代表主震区的热流值可能存在较大的不确定性.p值分布与热流关系并不显著,这可能是因为地热观测点分布稀疏且极不均匀,所以单个地震序列的衰减特征很难直接与大地热流值建立强相关关系,还需要进一步的积累地震序列和大地热流数据观测值.

为了进一步探究区域性pETAS值、pOML值与大地热流之间的关系,根据不同的区域特点,特将此研究区划分为A(小华北区域)、B(营-海-岫地区)、C(唐山老震区)、D(晋冀蒙交界地区)、E(河北平原地震带南部)五部分展开研究,如图 6a图 6b所示,五部分的平均热流值依次为63.6814 mW·m-2、69.7636 mW·m-2、58.1375 mW·m-2、57.6667 mW·m-2、61.5600 mW·m-2,ETAS模型的pETAS平均值为pETAS(A)= 0.9769±0.0931、pETAS(B)= 0.9822±0.1492、pETAS(C)= 0.8710±0.0759、pETAS(D)= 1.0795±0.0968、pETAS(E)= 0.8769±0.0666,R-J模型的pOML平均值为pOML(A)= 0.8990±0.0886、pOML(B)= 0.9885± 0.0822、pOML(C)= 0.8248±0.1078、pOML(D)= 0.9519±0.0757、pOML(E)= 0.8100±0.0925,大地热流与模型参数的关系如图 6c图 6d所示.研究可知,五个区域的数据关系总体上落在Kisslinger(1991)的统计关系里,即p值与热流值存在正相关关系,高温度对断层区介质应力释放具有加速作用,在高热流地区地震序列衰减速度较快,但受于数据的限制,上述结果还需继续积累更多的震例和序列参数进一步研究.

图 6 地震序列稳定时序列参数p值与大地热流值的关系分布 (a)pETAS值空间分布图;(b) pOML值空间分布图;(c) pETAS值与大地热流值的关系图;(d) pOML值与大地热流值的关系图;图中不同大小的矩形代表主震震级的大小,不同的颜色代表地震序列参数稳定时的模型参数值.图(c,d)中蓝色、浅红色分别代表了A、B、C、D、E五个区域稳定的序列参数p的平均值与大地热流平均值的关系,浅绿色数据来源于Kisslinger和Jones(1991)统计的南加州地区的33个有效序列参数p与大地热流的关系图(此处为R-J模型中的p值,与ETAS模型的p值仅做粗略讨论),黑色实线代表了数据的拟合情况,虚线代表了一倍的均方差. Fig. 6 The relationship between the stable sequence parameter p and the terrestrial heat flow (a) The spatial distribution map of the pETAS value; (b) The spatial distribution map of the pOML value; (c) Parameter pETAS of the ETAS model versus the heat flow; (d) Parameter pOML of the R-J model versus the heat flow. The different sizes of rectangles represent the size of the main shock, and different colors represent the size of the stability model parameters. The blue and light red in (c) and (d) represent the mean value of the stable sequence parameter p versus the mean value of the heat flow in the five regions A, B, C, D and E, respectively. Light green data are derived from Kisslinger′s (1991) diagram of the relationship between the 33 effective sequence parameters p and heat flow in Southern California (here it refers to the p value in the R-J model, which is roughly discussed with the p value of the ETAS model). The black solid line represents the overall situation of data fitting, and the dashed line represents the 1-fold standard deviation.

为系统地考察华北地区ETAS模型与R-J模型参数的统计特征以及与大地热流的关系,本文利用G-K方法和自然边界法筛选出1970年以来的16个地震序列,且每个地震序列完整性震级以上的余震数目满足N≥30,对挑选的16个地震序列进行了参数拟合.在得到稳定的序列参数并进行分析后,与区域内大地热流之间的关系进行了讨论.获得了以下几个方面的认识:

(1) 通过对所选的16个序列在震后0.5天开始,以0.05天为步长利用连续滑动、多时段拟合的方式对地震序列进行了参数拟合,发现序列参数在震后的早期变化较为剧烈,当序列趋于稳定时,各个主震事件的序列参数如表 1所示.因此,将序列参数用于震后地震序列类型快速判断、地震预测等研究中需谨慎对待,早期序列参数的剧烈变化也正反映了地震主震发生后震源区应力的快速调整过程.

(2) 利用ETAS模型和R-J模型对所选的16个地震序列进行研究,在序列参数达到稳定时,序列参数的总体特征如表 2所示,模型各个参数分布较为集中,呈现一种优势分布,且与主震震级大小关系不大.

(3) 通过对16个地震序列模型参数对比发现,ETAS模型参数pETAS大于R-J模型参数pOML的序列为12个,占75%,即大多数情况下地震序列ETAS模型参数pETAS大于R-J模型参数pOML,整体性平均偏高0.0716,这为以后研判该区域地震早期序列参数特征提供了科学依据.

(4) 地震序列参数αETAS与大地热流值整体上呈现较弱的负相关关系,可能表明大地热流值较高的地方激发次级余震的能力相对较强.

本文获得的华北地区序列参数统计特征,对构建适合该区域的地震预测模型以及地震危险性评估,都具有潜在的应用价值.由于华北地块可获得的大地热流的数值分布范围相对集中,难于判断单个序列参数pETASpOML与大地热流值的相关性.但从划分的五个子区域的平均参数值与大地热流值的分布关系看,存在正相关趋势,基本符合Kisslinger和Jones(1991)给出的美国南加州地区统计规律,也即在较高的大地热流地区,余震序列衰减速率可能更快.但受限于热流值的观测数据有限、分布不均匀等客观条件,序列参数与大地热流值的关系尚需要进一步的验证.此外,在震后的早期阶段,序列参数的稳定性实际上涉及到震源区应力调整和恢复,以及断层带愈合问题.由于本论文从统计地震学角度讨论问题,未涉及更多的断层带物理问题,因此对相关问题的讨论还不够深入,在未来的研究中,需要结合应力降、地球动力学的地壳厚度以及断层带结构等因素进一步深入开展相关研究.

致谢  本文属中国地震局2018年度震情跟踪定向工作任务经费项目.国际地震可预测性合作研究(CSEP)计划中国检验中心筹备组对本研究给予了指导.研究中使用了中国地震台网中心全国地震编目系统提供的统一正式目录,研究中使用的ETAS模型拟合程序由日本统计数理研究所庄建仓教授所提供,在此一并表示感谢.
References
Aki K. 1965. Maximum likelihood estimate of b in the formula logN=a-bM and its confidence limits. Bulletin of the Earthquake Research Institute, University of Tokyo, 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 (in Chinese), 32(1): 8-17. DOI:10.6038/pg20170102
Chapman D S, Rybach L. 1985. Heat flow anomalies and their interpretation. Journal of Geodynamics, 4(1-4): 3-37. DOI:10.1016/0264-3707(85)90049-3
Enescu B, Hainzl S, Ben-Zion Y. 2009. Correlations of seismicity patterns in Southern California with surface heat flow data. Bulletin of the Seismological Society of America, 99(6): 3114-3123. DOI:10.1785/0120080038
Furlong K P, Chapman D S. 1987. Thermal state of the lithosphere. Reviews of Geophysics, 25(6): 1255-1264.
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.
Guo Y C, Zhuang J C, Zhou S Y. 2015. An improved space-time ETAS model for inverting the rupture geometry from seismicity triggering. Journal of Geophysical Research:Solid Earth, 120: 3309-3323. DOI:10.1002/2015JB011979
Guo Z Q, Ogata Y. 1997. Statistical relations between the parameters of aftershocks in time, space, and magnitude. Journal of Geophysical Research:Solid Earth, 102(B2): 2857-2873. DOI:10.1029/96JB02946
Gutenberg R, Richter C F. 1944. Frequency of earthquakes in California. Bulletin of the Seismological Society of America, 34: 185-188.
Han P, Hattori K, Zhuang J C, 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
Hu S B, He L J, Wang J Y. 2000. Heat flow in the continental area of China:a new data set. Earth and Planetary Science Letters, 179(2): 407-419. DOI:10.1016/S0012-821X(00)00126-6
Hu S B, He L J, Wang J Y. 2001. Compilation of heat flow data in the China continental area (3rd edition). Chinese Journal of Geophysics (in Chinese), 44(5): 611-626.
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:Solid Earth, 111(B4): B04301. DOI:10.1029/2005JB003982
Jiang C S, Wu Z L. 2011. Intermediate-term medium-range Accelerating Moment Release (AMR) priori to the 2010 Yushu MS7. 1 earthquake. Chinese Journal of Geophysics (in Chinese), 54(6): 1501-1510. DOI:10.3969/j.issn.0001-5733.2011.06.009
Jiang C S, Wu Z L, Han L B, et al. 2013. 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 (in Chinese), 56(12): 4048-4057. DOI:10.6038/cjg20131210
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 (in Chinese), 58(11): 4163-4173. DOI:10.6038/cjg20151123
Jiang C S, Zhuang J C, Wu Z L, et al. 2017. Application and comparison of two short-term probabilistic forecasting models for the 2017 Jiuzhaigou, Sichuan, MS7. 0 earthquake. Chinese Journal of Geophysics (in Chinese), 60(10): 4132-4144. DOI:10.6038/cjg20171038
Jiang C S, Bi J M, Wang F C, et al. 2018. Application of the Omi-R-J method for forecast of early aftershocks to the 2017 Jiuzhaigou, Sichuan, MS7. 0 earthquake. Chinese Journal of Geophysics (in Chinese), 61(5): 2099-2110. DOI:10.6038/cjg2018M0113
Jiang G Z, Gao P, Rao S, et al. 2016. Compilation of heat flow data in the continental area of China (4th edition). Chinese Journal of Geophysics (in Chinese), 59(8): 2892-2910. DOI:10.6038/cjg20160815
Jiang H K, Qu Y J, Li Y L, et al. 2006. Some statistic features of aftershock sequences in Chinese mainland. Chinese Journal of Geophysics (in Chinese), 49(4): 1110-1117.
Jiang H K, Zheng J C, Wu Q, et al. 2007. Earlier statistical features of ETAS model parameters and their seismological meanings. Chinese Journal of Geophysics (in Chinese), 50(6): 1778-1786.
Kagan Y Y, Bird P, Jackson D D. 2010. Earthquake patterns in diverse tectonic zones of the globe. Pure and Applied Geophysics, 167(6-7): 721-741. DOI:10.1007/s00024-010-0075-3
Kisslinger C, Jones L M. 1991. Properties of aftershock sequences in southern California. Journal of Geophysical Research:Solid Earth, 96(B7): 11947-11958. DOI:10.1029/91JB01200
Kisslinger C. 1996. Aftershocks and fault-zone properties. Advances in Geophysics, 38: 1-36. DOI:10.1016/S0065-2687(08)60019-9
Li S Z, Zhang G W, Zhou L H, et al. 2011. The opposite Mesor-Cenozoic intracontinental deformations under the superconvergence:Rifting and extension in the North China Craton and shortening and thrusting in the South China Craton. Earth Science Frontiers (in Chinese), 18(3): 79-107.
Lolli B, Gasperini P. 2003. Aftershocks hazard in Italy Part I:Estimation of time-magnitude distribution model parameters and computation of probabilities of occurrence. Journal of Seismology, 7(2): 235-257. DOI:10.1023/A:1023588007122
Lomnitz C. 1966. Magnitude stability in earthquake sequences. Bulletin of the Seismological Society of America, 56(1): 247-249.
Mogi K. 1962. Study of elastic shocks caused by the fracture of heterogeneous materials and its relations to earthquake phenomena. Bulletin of the Earthquake Research Institute, University of Tokyo, 40: 125-173.
Ogata Y. 1978. The asymptotic behaviour of maximum likelihood estimators for stationary point processes. Annals of the Institute of Statistical Mathematics, 30(2): 243-261. DOI:10.1007/BF02480216
Ogata Y. 1983. Estimation of parameters in the modified Omori formula for aftershock frequencies by the maximum likelihood procedure. Journal of Physics of the Earth, 31: 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. 1992. Detection of precursory relative quiescence before great earthquakes through a statistical model. Journal of Geophysical Research:Solid Earth, 97(B13): 19845-19871. DOI:10.1029/92JB00708
Ogata Y. 2001. Increased probability of large earthquakes near aftershock regions with relative quiescence. Journal of Geophysical Research:Solid Earth, 106(B5): 8729-8744. DOI:10.1029/2000JB900400
Ogata Y. 2017. Statistics of earthquake activity:models and methods for earthquake predictability studies. Annual Review of Earth and Planetary Sciences, 45: 497-527. DOI:10.1146/annurev-earth-063016-015918
Ogata Y, Katsura K, Tsuruoka H, et al. 2018. Exploring magnitude forecasting of the next earthquake. Seismological Research Letters, 89(4): 1298-1304. DOI:10.1785/0220180034
Omori F. 1894. On the aftershocks of earthquakes. Journal of the College of Science, Imperial University of Tokyo, 7: 111-200.
Papazachos B C. 1974. On the time distribution of aftershocks and foreshocks in the area of Greece. Pure and Applied Geophysics, 112(3): 627-631. DOI:10.1007/BF00877298
Pollack H N, Hurter S J, Johnson J R. 1993. Heat flow from the earth's interior:Analysis of the global data set. Reviews of Geophysics, 31(3): 267-280.
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
Schorlemmer D, Werner M J, Marzocchi W, et al. 2018. The Collaboratory for the Study of Earthquake Predictability:Achievements and Priorities. Seismological Research Letters, 89(4): 1305-1313. DOI:10.1785/0220180053
Sclater J G, Jaupart C, Galson D. 1980. The heat flow through oceanic and continental crust and the heat loss of the Earth. Reviews of Geophysics, 18(1): 269-311.
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.
Taroni M, Marzocchi W, Schorlemmer D, et al. 2018. Prospective CSEP evaluation of 1-day, 3-month, and 5-yr earthquake forecasts for Italy. Seismological Research Letter, 89(4): 1251-1261. DOI:10.1785/0220180031
Utsu T. 1961. A statistical study on the occurrence of aftershocks. The 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 for earthquakes. Geophysical Bulletin of Hokkaido University, 13: 99-103.
Utsu T, Ogata Y, Matsu'ura R S. 1995. The centenary of the Omori formula for a decay law of aftershock activity. Journal of Physics of the Earth, 43(1): 1-33.
Wang J Y, Hu S B, Peng Z H, et al. 2012. Estimate of geothermal resources potential for hot dry rock in the continental area of China. Science & Technology Review (in Chinese), 30(32): 25-31.
Yang W Z, Ben-Zion Y. 2009. Observational analysis of correlations between aftershock productivities and regional conditions in the context of a damage rheology model. Geophysical Journal International, 177(2): 481-490. DOI:10.1111/j.1365-246X.2009.04145.x
Zaliapin I, Ben-Zion Y. 2013. Earthquake clusters in southern California II:Classification and relation to physical properties of the crust. Journal of Geophysical Research: Solid Earth, 118(6): 2865-2877. DOI:10.1002/jgrb.50178
Zhao Z X, Matsumura K, Oike K, et al. 1992. The p value of aftershock activity in Chinese mainland. Acta Seismologica Sinica (in Chinese), 14(1): 9-16.
Zhuang J, Harte D, Werner M J, et al. 2012. Basic models of seismicity: Temporal models.Community Online Resource for Statistical Seismicity. Analysis.https://doi.org/10.5078/corssa-79905851. http://www.corssa.org. Accessed 8 Apr 2018.
毕金孟, 蒋长胜. 2017. 晋冀蒙交界地区余震短期发生率的预测效能评估和预测策略研究. 地球物理学进展, 32(1): 8-17. DOI:10.6038/pg20170102
胡圣标, 何丽娟, 汪集旸. 2001. 中国大陆地区大地热流数据汇编(第三版). 地球物理学报, 44(5): 611-626. DOI:10.3321/j.issn:0001-5733.2001.05.005
蒋长胜, 吴忠良. 2011. 玉树MS7. 1地震前的中长期加速矩释放(AMR)问题.地球物理学报, 54(6): 1501-1510. DOI:10.3969/j.issn.0001-5733.2011.06.009
蒋长胜, 吴忠良, 韩立波, 等. 2013. 地震序列早期参数估计和余震概率预测中截止震级Mc的影响:以2013年甘肃岷县漳县6. 6级地震为例.地球物理学报, 56(12): 4048-4057. DOI:10.6038/cjg20131210
蒋长胜, 吴忠良, 尹凤玲, 等. 2015. 余震的序列参数稳定性和余震短期发生率预测效能的连续评估-以2014年云南鲁甸6. 5级地震为例.地球物理学报, 58(11): 4163-4173. DOI:10.6038/cjg20151123
蒋长胜, 庄建仓, 吴忠良, 等. 2017. 两种短期概率预测模型在2017年九寨沟7. 0级地震中的应用和比较研究.地球物理学报, 60(10): 4132-4144. DOI:10.6038/cjg20171038
蒋长胜, 毕金孟, 王福昌, 等. 2018. 利用早期余震预测的Omi-R-J方法对2017年四川九寨沟MS7. 0地震的应用研究.地球物理学报, 61(5): 2099-2110. DOI:10.6038/cjg2018M0113
姜光政, 高堋, 饶松, 等. 2016. 中国大陆地区大地热流数据汇编(第四版). 地球物理学报, 59(8): 2892-2910. DOI:10.6038/cjg20160815
蒋海昆, 曲延军, 李永莉, 等. 2006. 中国大陆中强地震余震序列的部分统计特征. 地球物理学报, 49(4): 1110-1117. DOI:10.3321/j.issn:0001-5733.2006.04.024
蒋海昆, 郑建常, 吴琼, 等. 2007. 传染型余震序列模型震后早期参数特征及其地震学意义. 地球物理学报, 50(6): 1778-1786. DOI:10.3321/j.issn:0001-5733.2007.06.018
李三忠, 张国伟, 周立宏, 等. 2011. 中、新生代超级汇聚背景下的陆内差异变形:华北伸展裂解和华南挤压逆冲. 地学前缘, 18(3): 79-107.
汪集旸, 胡圣标, 庞忠和, 等. 2012. 中国大陆干热岩地热资源潜力评估. 科技导报, 30(32): 25-31. DOI:10.3981/j.issn.1000-7857.2012.32.002
赵志新, 尾池和夫, 松村一男, 等. 1992. 中国大陆性地震的余震活动p值. 地震学报, 14(1): 9-16.