基于Gutenberg-Richter定律快速估算最大余震震级：以2017年九寨沟MS 7.0地震为例

解孟雨 孟令媛 申文豪 史保平

1）中国科学院大学地球科学学院，北京市玉泉路19号（甲） 100049

2）中国地震台网中心，北京 100045

3）中国地震局地壳应力研究所，北京 100085

0 引言

据中国地震台网中心测定，2017年8月8日21时19分四川省阿坝州九寨沟县发生MS7.0地震（以下简称九寨沟地震），震中 33.22°N、103.83°E，震源深度 20km，矩震级MW6.5。地震震中位于岷江断裂、塔藏断裂和虎牙断裂附近，推测其发震构造为塔藏断裂南侧分支或虎牙断裂北段。此次地震是继2008年汶川MS8.0地震和2013年芦山MS7.0地震后在青藏高原东缘发生的又一次破坏性地震，据统计，截至8月13日20时，地震共造成25人死亡，525人受伤，6人失联。九寨沟县17个乡镇的1680户房屋受损，52个风景名胜中的30个受到影响，造成较为严重的人员经济损失。

中强地震发生之后，常伴随着一系列震级小于主震震级的地震，即余震。主震发生后为了避免较大余震产生的2次危害，有必要快速准确地预测主震震中周边余震的时空分布以及可能发生的最大余震。对于最大余震的震级，Båth（1965）指出，主震震级与最大余震震级之差ΔM是不受主震震级影响的常数（即贝特定律）

其中，Mmain为主震震级；MRLA为最大余震震级；且ΔM的平均值约为1.2。研究该定律对于理解地震破裂过程和余震发生规模有着重要意义（Shcherbakov et al，2005；Shearer，2012；Rodríguez-Pérez et al，2016）。众多研究者对 ΔM的大小、影响因素及物理机制等进行了研究，结果显示，贝特定律的存在表明地震触发机制中存在一定程度的自相似过程（Console et al，2003；Shearer，2012；Žalohar，2014）；贝特定律并不严格适合所有余震序列，常数 ΔM仅在平均意义上接近 1.2（Kisslinger et al，1991；Console et al，2003；Shcherbakov et al，2005、2013；Shearer，2012；Tahir et al，2012），其值为 0～3（Kisslinger et al，1991；Helmstetter et al，2003；Shcherbakov et al，2004、2013；蒋海昆，2010），且可能受震源机制（苏有锦等，2008；Tahir et al，2012；Rodríguez-Pérez et al，2016）、断层之间的相互作用情况（Žalohar，2014）、震源深度（Båth，1965；苏有锦等，2014）、地震序列类型（苏有锦等，2014）、区域构造特征（Rodríguez-Pérez et al，2016）等因素的影响。为了更好地估计最大余震震级，我国学者依据 Gutenberg-Richter定律（以下简称G-R定律）首先提出了b值截距法（吴开统等，1984；毛春长，1989；张智等，1989；国家地震局科技监测司，1990）。该方法是将震级-频度图上 G-R定律的拟合曲线与震级轴交点的震级值作为最大余震震级的估计值。之后，研究者一方面不断验证该方法；另一方面，利用该方法估计最大余震震级，同时也计算了在震后不同时间内估计结果的变化情况（刘正荣，1995；钱晓东等，2008；苏有锦等，2008、2014；付虹等，2008）。研究结果显示，虽然使用该方法得到的最大余震震级存在一定的偏差，但考虑到误差的影响，该方法仍可以很好地估计最大余震震级。在国外，Shcherbakov等（2004）依据G-R定律提出了改进贝特定律，并用于对美国加州包括 Landers7.3级地震在内的10个中强地震余震序列的研究。他们提出了推定最大余震震级MILA的概念，并给出了MILA的计算公式，其实质与b值截距法相同，不同的仅是相应的描述方式。与b值截距法一样，对于不同余震序列，计算得到的推定最大余震震级MILA与真实最大余震震级MRLA也基本一致（Shcherbakov et al，2004、2013）。鉴于推定最大余震震级MILA的含义更为明确，故本文采用推定最大余震震级MILA进行方法说明。本文利用推定最大余震震级MILA的计算公式，针对九寨沟地震余震序列，探讨了在不同时间尺度下提取最大余震震级估算值的可行方案，并与实际测量值进行比较。依据对九寨沟地震余震序列的计算和结果分析，本文给出了快速评估最大余震震级的方法。

1 数据及完备性分析

本文选取九寨沟地震发震后至10月5日3时18分（主震发生后约57天）的余震数据用于建模和参数确定，该目录来源于中国地震台网中心，其中包含里氏震级ML和面波震级MS两类数据。若仅考虑里氏震级目录，则共有余震事件5864个，ML≥3.0事件共92个，其中，3.0≤ML＜4.0地震 78次，4.0≤ML＜5.0地震 13次，5.0≤ML＜6.0地震 1次。特别需要强调的是，在所选取的地震目录中，最大震级的余震为8月9日10时17分（距主震发生时刻约0.54天）发生的ML5.2（或MS4.8）地震。图1为九寨沟地震余震的空间分布情况。由图1可见，九寨沟地震震中位于汶川MS8.0地震、芦山MS7.0地震震中的北侧，且九寨沟地震余震主要集中于主震附近，近似条带状分布。在使用地震目录的数据时，首先应确定目录的完备震级Mc，即在给定的时间、空间范围内可全部测定地震的最小震级。如果完备震级Mc的估算值太高，则会导致地震数据使用不充足，不能获得完整的地震信息；而太低的完备震级Mc值又会引入记录不完整的数据，导致错误的模型参数估计，进而得出错误的结论（Mignan et al，2012；蒋长胜等，2013）。因此，对于九寨沟地震的余震序列，我们需要选择和确定快速且准确的完备震级Mc估算方法。完备震级的估算方法主要有最大曲率、GFT、EMR、MBS等方法（Hamdache et al，2017），其中，利用最大曲率（以下简称 MAXC）方法进行计算，能较为适合快速评估的需要。首先，该方法在实际计算时将非累积频率-震级分布中取得最大值的震级档作为完备震级Mc，因而它不需要任何的参数拟合，是一种快速直接且对数据强健的方法（Mignan et al，2012）。此外，Mignan等（2012）还指出，MAXC方法也可在目录数据较少的情况下得出合理的完备震级，适合在主震发生后短时间内估算完备震级Mc的大小。不过MAXC方法会低估地震目录的完备震级，因此，为更接近实际真值，可以添加相应的修正系数，即Mc=Mc（MAXC）＋0.2（Woessner et al，2005）。

图1 九寨沟地震余震分布与周边近期发生的强震

因此，本文采用修正的MAXC方法对九寨沟地震的余震目录进行完备震级的估算。对于估算完备震级的余震目录，我们选择主震后2h（约为0.083天）内的余震序列目录。之所以采用这样的选择，一方面是为了能够在主震发生后很短的时间尺度下快速得到完备震级，另一方面也是由于采用余震序列早期的目录才可得到较为准确的完备震级（Enescu et al，2002；米琦等，2015）。经过对九寨沟地震2类目录（里氏震级ML、面波震级MS目录）的计算得到，对于里氏震级ML的余震目录，其完备震级Mc约为2.4；而震级类型为面波震级MS的余震目录，其完备震级Mc=1.6。2类目录完备震级估算值的不确定度均约为0.007，可忽略不计。

2 最大余震震级的估算

无论是全球尺度的地震还是区域尺度的地震，除去少数个例，基本都满足Gutenberg-Richter定律（简称 G-R定律）（Marzocchi et al，2003；Shcherbakov et al，2004），即在给定区域和时间段内，其震级大于等于M的地震总个数N（≥M），满足

其中，M≥Mc；a、b为正常数，a给出了选定区域和时间段内0级以上地震个数的对数，它依赖于选定区域和观测数据时间尺度的大小（Hamdache et al，2017），而b值则刻画了区域内不同震级地震的相对分布情况（Hamdache et al，2017），其值为 0.6～1.4（Marzocchi et al，2003；Shcherbakov et al，2004；Wang et al，2015；Hamdache et al，2017）。b值的大小与区域应力状态和断层的破裂过程有关，可用来指示区域地震危险性水平（吴开统等，1984；Hamdache et al，2017）。

图2 九寨沟地震发生后1天内余震目录的震级-频度图

对于余震序列而言，G-R定律依旧成立。在这种情况下，N（≥M）表示余震序列中满足震级大于等于M（M≥Mc）的地震的个数，而且对于余震序列，相应b的取值范围也与区域或全球的统计结果一致（Shcherbakov et al，2004）。图2为九寨沟地震发生后1天内余震目录的震级-频度图。由图2可见，对于实际余震目录而言，N=1的最大震级档（M=5.2）即为主震之后最大余震震级MRLA，即N（≥MRLA）=1。据此，Shcherbakov等（2004）定义推定最大余震震级MILA为用G-R定律描述的余震震级-频度分布满足N（≥M）=1时的地震震级档M（图 2），即

于是，贝特定律可改写为

式（4）即为Shcherbakov等（2004）提出的改进贝特定律。Shcherbakov等（2004）利用改进贝特定律分析了加州1987～2003年间发生的10个震级大于5.5级地震的最大余震震级，结果显示，由该定律得到的推定最大余震震级MILA或震级差（ΔM′）与实际值一致。因此，本文中我们使用推定最大余震震级MILA（式（3））估计实际最大余震震级MRLA。

为了计算推定最大余震震级MILA，首先要估算出b、a值。对于b值的估计，可采用最小二乘法和最大似然法。前人的研究指出，利用最小二乘法计算的b值存在系统性偏差（Sandri et al，2007；Hamdache et al，2017）。因此，本文选择 Marzocchi等（2003）描述的最大似然法估计b值。具体而言，b值及不确定度δb（Uncertainty）的计算公式为

而相应的不确定度δa为（Taylor，1997）（完备震级Mc的不确定度约为 0.007；相对于b值的不确定度可忽略不计）

于是，推定最大余震震级MILA的不确定度δMILA可表示为（Taylor，1997）

3 结果分析与讨论

针对九寨沟地震的余震序列，我们分别对里氏震级ML目录和面波震级MS目录计算了不同时间尺度下推定最大余震震级MILA。图3、图4、表1和表2为相应的计算结果。由图3（a）、3（c）、4及表1、2可见，随着选取余震目录时间范围（即从主震发生后到截止时刻经过的时间）的不断扩大，推定最大余震震级MILA逐渐增大并稳定于固定值。对于里氏震级ML目录而言，该值为MILA=5.3±0.1；而对于面波震级MS，则为MILA=5.0±0.1，这与实际测得的最大余震震级基本一致（MRLA=ML5.2或MRLA=MS4.8）。这也说明利用推定最大余震震级MILA可对最大余震震级进行有效估算，并且不受地震目录中震级类型的限制，即可获得与实际观测较一致的数值。而MILA的不确定度随时间的流逝不断减少，这是由于随时间的增加，地震目录不断增加，有效信息不断增多。这也提醒我们，如果能加强台网对余震的识别能力，增加余震数据，就可以给出更为准确的推定最大余震震级。

图3 不同时间尺度下由九寨沟地震序列得到的推定最大余震震级MILA随时间变化的误差棒图

对于较小的截止时间，即较短的余震目录时间，图3（b）、3（d）显示，在0.6天（14.4h）之前推定最大余震震级MILA并不稳定，并且整体上不确定性较大（余震个数较少，信息不足）。不过，随着截止时间的增加，MILA的值趋于稳定。截止时间为1、2天的MILA值分别为5.0±0.3（ML）、4.6±0.3（MS）和 5.2±0.3（ML）、4.9±0.2（MS），这些值与实际测得的最大余震震级MRLA（ML5.2和MS4.8）基本一致。因此，在短时间内，仍可以利用推定最大余震震级MILA估算最大余震震级MRLA，而且也依然不受震级类型的影响。此外，在实际最大震级余震发生前（约为0.53天前），由于地震目录中余震个数较少，利用0.5天的数据得到的推定最大余震震级ML（4.9±0.4）和MS（4.5±0.3）均小于实际最大余震震级MRLA（MRLA=ML5.2或MRLA=MS4.8），不过它们却可以作为MRLA值的合理下界，这也说明利用推定最大余震震级确实可以有效地预测最大余震震级。据此，在主震发生后利用推定最大余震震级，我们可以在1天或2天时得到较为合理的最大余震里氏震级或面波震级的估计值；另一方面，在更短的时间内（0.5天或12h），对于现有台网识别得到的余震数据，利用推定最大余震震级可以给出MRLA合适的下界。随着台网识别能力的提升，一方面可以增加各个震级余震的个数N，降低不同参数的不确定度，进而得到更为准确的MRLA估计值；另一方面也会降低完备震级Mc，在主震后短时间内就可有大量可使用的余震数据，实现在更短时间内得出MRLA的有效估计值。

图4 九寨沟地震发生后不同时间段内余震目录的震级-频度图

表1 不同时间尺度下计算得到的九寨沟地震推定最大余震里氏震级MILA（ML）

表2 不同时间尺度下计算得到的九寨沟地震推定最大余震面波震级MILA（MS）

4 结论

本文以2017年九寨沟MS7.0（或ML7.1）地震震后57天的余震目录（包含ML、MS目录）为资料，通过对余震目录完备震级、a值和b值的计算，并利用推定最大余震震级，在不同时间尺度上估算了九寨沟地震的最大余震震级。结果显示，推定最大余震震级可有效地预测最大余震震级，并给出合理的估计值。具体而言，在长时间尺度上推定最大余震震级可以给出与实际观测最大余震震级一致的估计值；而在短时间尺度上，1、2天的结果也可作为最大余震震级的合适估计，且0.5天（12h）的结果就可作为最大余震震级的合理下界。

此外，利用推定最大余震震级估计实际最大余震震级的方法，不受震级类型的限制，对于MS、ML目录均可得到较好的结果。因此，我们建议可采用以下步骤进行最大余震震级的快速评估：①根据震后2h的目录，利用修正的MAXC方法计算完备震级；②利用推定最大余震震级分别在主震后12h、1天、2天时估算最大余震震级，将其分别作为最大余震震级的下界和估计值。另外，计算结果也显示，如果可以提高余震的识别能力，利用推定最大余震震级则可在短时间内给出更准确的最大余震震级的估计值。