全球MW≥8.0浅源地震的前震序列研究

薛 艳 解孟雨 刘 杰 庄建仓

1)中国地震台网中心， 北京 100045 2)日本东京数理统计研究所， 东京190—8562

0 引言

近年来全球发生了多次8级以上大地震， 造成了严重的人员伤亡和巨大的财产损失。由于有仪器记录的时间短， 因此有完整小地震记录以来的8级地震震例较少。人们对地震序列的研究主要集中于5级以上地震(Wellsetal.， 1994； 蒋海昆等， 2007； Chenetal.， 2013， 2015)。有学者对全球大地震的地震序列进行了一定的研究(苏有锦等， 2008， 2014； 吕晓健等， 2010)， 但尚未有涉及8级大地震前震序列特征的系统研究。

2011年3月11日日本本州以东海域9.1级地震震前2d发生了7.3级前震序列， 最大前震距主震约43km， 整个前震序列分布在主震NE方向100km范围内。2012年4月11日苏门答腊8.6级地震震前90d发生了7.2级前震序列， 最大前震与主震震中相距约20km， 前震序列分布在距主震震中66km的范围内。2014年4月1日智利伊基克8.2级地震震前15d发生了6.7级前震序列， 最大前震与主震震中相距42km， 前震序列分布在主震周围75km的范围内。这3次大地震的前震活动引起了广泛关注(Meredithetal.， 2011； Xueetal.， 2012； Brodskyetal.， 2014)。

在所有的地震短临前兆中， 前震是学术界公认的预报强震最有效的指标之一(Jonesetal.， 1979)。陈颙等(2015)总结了前震序列的4个特点： 前震序列是一种高频率的地震活动， 其活动频率比背景地震的活动频率高； 前震序列在主震之前(几h—几d)发生； 前震序列发生的位置与主震相同； 主震的震级比前震序列中任何地震的震级都大。由此可见， 前震的确定与其和主震的时间差和距离差的约束有关。Chen等(2015)系统研究了加利福尼亚地区5级以上地震的前震序列特征， 发现48%的5级以上地震在震前30d内， 距离主震5km范围内至少发生了1次前震。

此外， 陈顒(1978)还发现前震序列的震源机制解类似， 而余震序列的震源机制解差异大， 并据此提出利用震源参数的一致性判别一个序列是前震序列还是正常的主震-余震序列的方法。倪四道等(2010)研究了2010年4月14日青海玉树7.1级地震前约2h在距离主震震中仅2km处发生的4.7级地震序列的时空分布及波形特征， 认为该4.7级序列符合陈顒(1978)的理论， 是典型的前震序列。

Scholz(1968)在岩石破裂实验中发现， 微破裂b值与压应力呈负相关。针对一些特定震例的研究表明，b值会在主震前后发生变化， 与余震相比， 前震序列的b值偏低(Smith， 1981； Xueetal.， 2012)。这一特性与初始的高应力状态随着主震的应力降而降低的理论一致， 该结论常被用于识别前震(Scholz， 1968； 李全林等， 1978； 林邦慧等， 1994)。Gulia等(2019)研究认为， 前震b值显著低于区域背景b值， 并提出了实时判别前震的 “交通灯”准则。本文将系统研究1976年以来全球MW≥8.0浅源地震前震序列的统计特征， 应用 “传染型余震序列(epidemic-type aftershock sequence， ETAS)”模型计算前震和余震序列参数， 并讨论参数的差异性。研究中使用了美国国家地震信息中心(NEIC)提供的全球地震目录(1)http://earthquake.usgs.gov。， 震源机制解数据来源于美国哈佛大学(2)http://www.globalcmt.org/CMTsearch.html。。

1 地震目录最小完备性震级的讨论

美国国家地震信息中心(NEIC)提供了1973年以来全球地震目录。该目录早期使用了多个震级标度， 但是近20a以来震级标度主要为mb和MW， 其中4级及以下地震主要为mb震级， 5级以上地震主要为MW震级。本文对该目录中的所有地震统一用震级M表示。

采用定量评估的 “最大曲率MAXC”方法(Wiemeretal.， 2002)分析全球地震目录的最小完备性震级Mc随时间的变化。最大曲率方法(maximum curvature method)是将震级-频度曲线一阶导数的最大值对应的震级作为Mc。结果显示(图 1)： 1973—1994年， 最小完备性震级为4.6～5.0级； 1995—2012年， 最小完备性震级为4.3～4.7级； 2013—2016年， 最小完备性震级为4.2～4.6级。

图 1 基于MAXC方法计算的全球最小完备性震级随时间的变化曲线Fig. 1 Curve of global catalogue minimum completeness magnitude with time calculated by “MAXC” method.

2 前震序列的活动特征

1976—2017年全球共发生独立的MW≥8.0浅源地震29次， 其中逆冲型23次。为了确定前震， 本文先做如下约束： 1)前震序列发生在余震区内， 本文取主震后10d中余震的密集活动区作为余震区； 2)主震前30d， 余震区内的地震活动水平明显高于背景状态， 频度达到背景值的2倍以上； 3)最大前震震级M≥5.0。

依据上述约定， 共挑选出8个前-主-余型地震序列(蒋海昆等， 2007)， 其主震发震断层滑动角均位于69°～110°范围内， 即全部为逆冲型地震， 占8级以上逆冲型地震总数(23次)的34.8%。前-主-余型地震主要分布在环太平洋地震带的西北边界带、 西南边界与东南边界带(图 2)， 这主要受全球8级以上地震空间分布的影响。

图 2 全球1976—2017年MW≥8.0地震的震源机制解Fig. 2 Global shallow focus mechanism solution with MW≥8.0.红色震源球为有前震的地震

之后统计了 8 次大地震前震活动开始的时间、 最大前震与主震的震级差及震中距等参数， 结果显示(表1)， 最大前震与主震的震级差为1.1～2.8级， 震中距为10～53km。87.5%的时间差(最大前震与主震发生的时间差)为2h～15d， 仅2006年11月15日千岛群岛8.3级地震震前45d发生最大前震。

表1 8次大地震的主震和最大前震的基本参数、 最大前震与主震的时间差、 震级差等相关统计Table1 Basic parameters of the 8 great earthquakes and their largest foreshocks， the difference of occurrence time and magnitude between the main shock and the largest foreshock

图 3 是8次前-主-余型地震序列的震中分布， 其中余震序列的持续时间为主震后1个月。图 4 是地震序列沿图 3 所示AB剖面的震中迁移图， 即将地震投影到AB线段上， 绘制每个地震投影点与A点的距离D随时间t的变化情况。可以看出， 前震的空间分布集中， 且集中分布于主震震中附近。图 5 是8次大地震的前震、 主震和余震震源机制解分布图。 可见， 前震序列震源机制解与主震震源机制解显示的断层错动类型及2个节面的走向、 倾角都较为一致， 而余震的震源机制解比较复杂， 不仅破裂类型不一致， 且节面的走向、 倾角都有很大差异， 这与陈颙(1978)提出的前震特征相符。

图 3 8次大地震的前震和余震序列震中分布Fig. 3 Epicenter distribution of foreshock and aftershock sequences of 8 great earthquakes.a 1976-01-14， 克马德克群岛海域， M8.2； b 1985-03-03， 智利中部近海， M8.0； c 1986-05-07， 阿留申群岛， M8.0； d2006-11-15， 千岛群岛， M8.3； e 2007-04-01， 所罗门群岛， M8.1； f 2011-03-11， 日本本州近海， M9.1； g 2013-02-06， 圣克鲁斯群岛， M8.0； h 2014-04-01， 智利北部近海， M8.2。红色圆形、 蓝色圆形和黄色五角星分别表示前震、 余震和主震震中； 图中虚线框表示表2中背景b值的计算范围

图 4 8次大地震的前震和余震序列沿AB剖面(图 3所示)的震中迁移D-t图Fig. 4 Epicenter migration map in AB direction.a 1976-01-14， 克马德克群岛海域， M8.2； b 1985-03-03， 智利中部近海， M8.0； c 1986-05-07， 阿留申群岛， M8.0； d2006-11-15， 千岛群岛， M8.3； e 2007-04-01， 所罗门群岛， M8.1； f 2011-03-11， 日本本州近海， M9.1； 2013-02-06， 圣克鲁斯群岛， M8.0； h 2014-04-01， 智利北部近海， M8.2。横坐标为时间t， 以第1次前震发生的时间为零点， 后续地震距离第1次前震的时间间隔为横坐标值， 单位为d； 图中红色空心圆表示主震

图 5 8次大地震的前震和余震的震源机制解示意图Fig. 5 Sketch map of focal mechanism solution for earthquake sequence.a 1976-01-14， 克马德克群岛海域， M8.2； b1985-03-03， 智利中部近海， M8.0； c1986-05-07， 阿留申群岛， M8.0； d2006-11-15， 千岛群岛， M8.3； e2007-04-01， 所罗门群岛， M8.1； f 2011-03-11， 日本本州近海， M9.1； g 2013-02-06， 圣克鲁斯群岛， M8.0； h 2014-04-01， 智利北部近海， M8.2。红色震源球为前震和主震， 蓝色震源球为余震

为了考察前震序列在临近主震时是否出现加速活动， 统计了主震前半年(180d)余震区内M≥Mc(Mc为最小完备震级)地震累积频度随时间的变化。取余震区及其附近主震前3a左右的4级以上地震， 使用G-R关系确定最小完备震级。结果显示， 1976年克马德克群岛海域8.2级地震余震区附近的最小完备震级为4.9级， 1985年智利中部近海8.0级地震附近的最小完备震级为4.7级， 1986年阿留申群岛8.0级地震附近的最小完备震级为4.6级， 2007年4月1日所罗门群岛8.1级地震附近的最小完备震级为4.2级。另外4次地震， 即2006年千岛群岛8.3级、 2011年日本本州以东海域9.1级、 2013年圣克鲁斯群岛8.0级和2014年智利北部近海8.2级， 最小完备震级分别为4.2级、 4.5级、 4.4级和4.0级。图6a是8次大地震前半年余震区内M≥Mc地震累积频度随时间的变化曲线， 图6b是累积频度归一化的结果。可以看出， 主震发生前半年， 余震区内仅有少量地震活动， 而到主震前约1个月地震活动显著增强， 累计频度出现加速增长的特点。具体而言， 在8次大地震中， 1986年8.0级、 2006年8.3级和2007年8.1级地震的前震序列分别自主震前35d、 45d和43d开始频度显著增加， 其余5次均在主震前15d内开始加速活动。值得注意的是， 1986年和2007年2个前震序列在主震前1d内、 2006年前震在主震前6d又再次出现加速活动(图 6)。

图 6 8次大地震主震前180d余震区内最小完备震级以上(M≥Mc)地震的累积频度曲线Fig. 6 Cumulative frequency curve of earthquakes with M≥Mc in aftershock area 180 days before 8 main earthquakes.横坐标以主震发生时刻为0点， 负值表示主震前

3 应用ETAS模型计算地震序列参数

本文使用ETAS模型分别计算前震和余震序列参数。ETAS模型假设所有余震均可以按照大森-宇津公式(Utsu， 1961)激发自身的余震， 且震级的分布是独立的。假定主震的发生时间为初始零时刻， 在其后某观测时段[0，T]内地震序列{(ti，Mi)；i=1， 2， …，N}的强度函数可表示为(Ogata， 1988)

(1)

其中，t为主震发生后的离逝时间；M0为计算所用序列的截止震级；Mi、ti分别为第i个事件的震级与距主震的时间长度；μ为背景地震发生率；p表示序列衰减快慢；c为主震后余震频次达到峰值的时间长度；K表示余震的活跃程度；α表示触发次级余震的能力(Ogata， 1989， 1992)。

ETAS模型参数可使用最大似然法进行估计。在拟合时段[S，T]内， 似然函数L可表述为

(2)

将式(1)带入式(2)， 即可对参数[μ，K，c，α，p]进行最大似然估计。

为了考察ETAS模型的拟合效果， 一般使用 “残差分析”方法(Daleyetal.， 2003)将地震序列转换为在 “转换时间”(τ)域的分布， 并考察实际地震序列与理论值的拟合情况。将条件强度函数λ(t)采用时间序列{ti}进行如下转换：

(3)

由此将{ti}转化为服从单位速率的稳态泊松分布(3)Zhuang J， Harte D， Werner M J， et al.， 2012， Basic models of seismicity: Temporal models， Community Online Resource for Statistical Analysis， doi: 10.5078/corssa-79905851。{τi}。如果ETAS模型对数据的拟合较好， 则在转换时间域{τi}的地震累积曲线表现为线性， 接近标准稳态泊松过程的理论直线。

此外， 本文利用最大似然法(Aki， 1965； Utsu， 1966)计算b值：

(4)

(5)

式中，n为参与计算的地震数目。

本文采用震级-序号法确定前震和余震的最小完备震级(Huang， 2006)。图 7 是4次大地震的前震和余震的震级-序号图。可以看出， 在主震发生后短时间内， 由于主震波形振幅大、 面波持续时间长， 大量小震级余震被 “淹没”， 使得地震监测能力降低， 余震最小完备震级上升(Iwata， 2008)。本文通过调节ETAS模型拟合的起始时间C0以确保有足够多的地震参与计算， 即主震后—C0时刻期间的数据不参与计算， 取C0时刻后最小完备震级以上的地震计算序列参数(蒋长胜等， 2013)。由图 7 可见， 2006年千岛群岛8.3级地震的前震序列最小完备震级为4.2级(图7a)， 主震后0.139d余震的最小完备震级为4.4级(图7b)； 2011年3月11日日本本州近海9.1级地震前震序列的最小完备震级为4.5级(图7c)， 余震序列在主震后0.529d的最小完备震级是4.5级(图7d)； 2013年2月6日圣克鲁斯群岛8.0级地震前震序列的最小完备震级为4.4级(图7e)， 余震序列在主震发生后0.164d的最小完备震级为4.3级(图7f)； 2014年4月1日智利北部海域8.2级地震前震序列的最小完备震级为4.0级(图7g)， 余震序列在距离主震后0.046d的最小完备震级为4.4级(图7h)。

图 7 4次大地震的前震和余震序列的震级-序号图Fig. 7 The magnitude-number chart of foreshocks and aftershocks of 4 great earthquakes.a和b为2006年11月15日千岛群岛8.3级地震的前震和余震； c和d为2011年3月11日日本本州近海9.1级地震的前震和余震； e和f为2013年2月6日圣克鲁斯群岛8.0级地震前震和余震； g和h为2014年4月1日智利北部近海8.2级地震的前震和余震。 图中2条垂直虚线分别表示最小完备震级Mc和达到最小完备震级的时刻C0， C0为距离主震发生的时间， 单位为d

样本量N和最低计算震级Mj影响地震序列参数的估算结果(Shietal.， 1982； Bender， 1983； 蒋长胜等， 2013； 吴果等， 2019)。为了确保参数计算的可靠性， 本文约定当前震序列满足计算样本量N≥30且最低计算震级Mj≥Mc(Mc为最小完备震级)时， 才同时计算前震和余震序列参数α、p和b值。在8个前震-主震型地震序列中， 有4个可以对比计算其前震和余震序列参数， 其中余震序列的持续时间为主震后1个月(表2)。针对2006年千岛群岛8.3级和2014年智利8.2级地震的前震序列， 我们分别选了2个最低计算震级， 结果显示， 不同的最低计算震级可以得到不同的α、p和b值， 但由此造成的差异与计算误差大体相当， 即ETAS模型的计算结果比较稳定。为了比较， 对同一个地震的前震和余震取相同或相近的最低计算震级。所得结果表明， 前震和余震序列的α值和p值没有规律性差异， 而b值规律性明显， 前震b值显著低于余震(图 8， 表2)。

为了进一步识别前震， 我们计算余震区及附近较长时间的区域背景b值(确定计算时间主要考虑完备震级且避免期间发生大地震， 计算区域为图 3 虚线框所示区域)。经对比发现， 4次前震序列的b值均低于区域背景b值， 而3次的余震b值高于区域背景b值， 1次的余震b值与背景b值相当(表2)。计算4次前震序列b值与区域背景b值的差， 再除以区域背景b值， 结果依次为0.243、 0.163、 0.138和0.103， 即前震序列的b值低于区域背景b值的10%～24%。前震b值与区域背景b值之差依次是区域背景b值标准差的4.2倍、 6.1倍、 7.1倍和2.2倍， 具有显著性。

为了讨论前震b值的稳定性， 对样本量丰富的2个前震序列， 即2006年千岛群岛8.3级和2014年智利北部8.2级地震的前震序列进行累积滑动计算， 最低计算震级分别为4.2级和4.0级。结果显示(图 9)， 在前震发生的最初阶段b值较低， 之后逐渐增大， 当计算样本量N≥70后，b值基本稳定。

图 8 4次大地震的前震和余震b值(a)、 α值(b)和p值(c)计算结果对比Fig. 8 Parameters b-value(a)， α-value(b)and p-value(c)of foreshocks and aftershocks for 4 great earthquakes.图件横坐标为地震编号， 与表2 序号一致

图 9 2006年千岛群岛8.3级(a)和2014年智利北部近海8.2级(b)地震前震序列的b值随计算样本量的变化Fig. 9 Change of b-value with the calculated sample size of foreshock sequences of the 2006-11-15 Kuril Islands M8.3 and the 2014-04-01 off the coast of northern Chile M8.2 earthquake.图中直线和虚线表示背景b值及其误差

4 结论与讨论

综上所述， 得到以下几点结论。

(1)1976—2017年， 全球共发生29次8级以上浅源地震， 其中8次有前震， 占总数的27.6%， 略低于1900—2011年全球8.5级以上地震的前震比例(Xueetal.， 2012)， 明显高于中国大陆地区中强以上地震的前震比例(陈颙等， 2015)。

(2)8次前震序列的主震发震断层的滑动角为69°～110°， 即全部为逆冲型， 占8级以上逆冲型地震总数(23次)的34.8%。

(3)前震序列空间分布集中， 且集中分布在主震震中附近。最大前震与主震的震中距为10～53km， 震级差为1.1～2.8级， 优势发生时间差为2h～15d， 仅2006年11月15日千岛群岛8.3级地震震前45d发生最大前震。

(4)主震前半年， 余震区内有少量的地震活动， 而到主震前约1个月， 特别是主震前15d内， 地震活动显著增强， 累计频度加速增长， 即前震序列具有高频度活动的特点。具体而言， 8次大地震中有5次在主震前15d内出现加速活动； 3次的加速活动开始时间超过1个月(35～45d)， 但在主震前1d和6d又再次出现频次显著增多的现象。

(5)前震的震源机制解一致性好， 且与主震一致， 而余震的震源机制解比较复杂。

(6)应用ETAS模型计算前震和余震序列的α、p和b值， 结果显示， 对于反映激发次级余震能力的α值和序列衰减快慢的p值， 前震和余震没有共性差异； 而b值则共性特征明显， 前震b值显著低于余震。与区域背景b值相比， 前震b值明显偏低； 而余震b值高于背景b值或与其相当。

(7)为了讨论前震b值的稳定性， 计算了2个资料丰富的前震序列b值随样本量的变化。结果显示， 在前震序列的开始阶段b值较低， 之后逐渐增大， 当计算样本量N≥70后，b值基本稳定。

虽然前震是短临地震预报最有效的方法之一， 但前震的识别至今仍较为困难。特别是前震与前兆震群(宋俊高等， 1989)更为相似， 两者均具有空间分布集中、 高频次活动和震源机制解一致的特点， 如乳山震群(刘方斌等， 2018)和盖州震群(王亮等， 2015)。但前兆震群与余震均具有高b值的特点(前兆震群b≥0.65)(宋俊高等， 1989)， 而前震b值则较低(Guliaetal.， 2019)。

b值不仅反映了大、 小地震的比例关系， 更有明确的物理意义。b值的大小主要取决于震源区应力状态和介质性质， 并与应力呈反比(Scholz， 1968； Wyss， 1973； Urbancicetal.， 1992)， 可用来衡量区域相对应力水平(Wyssetal.， 2000； 易桂喜等， 2013)。本文对比研究了4次大地震的前震和余震序列的b值， 考虑了最低计算震级和样本量的影响， 并引入区域背景b值， 将前震b值与其进行对比， 结果显示， 前震b值低于区域背景b值的10%～24%， 前震b值与区域背景b值之差是区域背景b值标准差的2.2～7.1倍， 具有显著性， 这与Gulia等(2019)的研究结果一致。

此外， 本文研究的震例均为板间地震。有研究显示(苏有锦等， 2014)， 对于浅源地震， 板内和板间地震序列具有相似的时、 空、 强统计特征。 因此， 本文的研究结果对板内大地震具有参考意义。