震级-频度关系中 b值的极大似然法估计及其影响因素分析
吴果, 周庆, 冉洪流
中国地震局地质研究所, 活动构造与火山重点实验室, 北京 100029
*通讯作者: 周庆, 男, 研究员, E-mail: zqcsb@163.com

〔作者简介〕 吴果, 男, 1988年生,2018年于中国地震局地质研究所获固体地球物理学博士学位, 主要研究方向为地震活动性与地震危险性分析, 电话: 18810404034, E-mail: wgfirst@foxmail.com

摘要

b值在地震活动性研究和地震危险性分析中起着十分重要的作用, 通常拟合 b值的方法有最小二乘法(Least Squares Method)和极大似然法(Maximum Likelihood Estimation)。最小二乘法简单易行, 得到了广泛的应用。然而很多研究表明该方法存在一定的局限性, 极大似然法在特定条件下可以作为最小二乘法的一种可行的替代或补充方法。前人对极大似然法的研究非常繁杂, 提出了各种各样的方程式, 每个方程式的隐含假设和求解方式各不相同。文中对主要方程式进行了简要的回顾, 并按照是否考虑震级的归档效应、 是否设定有限最大震级、 是否对不同震级档数据取不同的观察时段和是否具有解析解这4个方面, 对这些方程式进行了分类和总结。进而对震级的归档效应、 震级的测量误差、 样本量、 震级跨度、 最小完整震级和前余震共6个可能影响极大似然法估计 b值的因素进行了分析和总结。最后对正确使用这些方程式提出了合理的建议。文中的分析和总结有助于更准确地理解和使用不同的极大似然法估计 b值的方程式, 以供相关研究者参考。

关键词: 震级-频度关系; b; 最小二乘法; 极大似然法; 震级的归档效应; 震级跨度
中图分类号:P315.5 文献标志码:A 文章编号:0253-4967(2019)01-0021-23
THE MAXIMUM LIKELIHOOD ESTIMATION OF b-VALUE IN MAGNITUDE-FREQUENCY RELATION AND ANALYSIS OF ITS INFLUENCING FACTORS
WU Guo, ZHOU Qing, RAN Hong-liu
Key Laboratory of Active Tectonics and Volcano, Institute of Geology,China Earthquake Administration,Beijing 100029, China
Abstract

b-value in the magnitude-frequency(G-R)relationship plays a vital role in seismicity research and seismic hazard analysis, and the most commonly used techniques to simulate it are least square approach and maximum likelihood method. Least square method is simple and easy to apply, therefore widely used in China. However, many researches show that there exist some limits in least square estimation of b-value. Earthquakes with different magnitudes are not equally weighted in this method, and larger events have higher weights, so b-value is vulnerable to the fluctuation of several big earthquakes; meanwhile, least square method needs to divide magnitude intervals artificially. With a small sample size, data points could be not enough if the magnitude interval is too wide, and events in a magnitude interval may be lacking if it is divided to be too narrow. Especially for incremental G-R relationship, it is possible that N( Mi)equals 0 in an interval with large magnitude, so log( N( Mi))loses meaning and has to be ignored, resulting in a low b-value. Therefore, under certain conditions, maximum likelihood method is recommended as an effective substitution or supplementary for least square estimation of b-value. Among numerous previous researches on maximum likelihood estimation of b-value, lots of equations have been provided, based on varied implicit assumptions and different ways of solution. A brief overview is first presented for these equations, and classification and summary are provided based on whether taking account of the effect of binned magnitude, with finite maximum magnitude, using unequal observation periods for different magnitude intervals, and with analytic solution or not. Following this, a total of 6 influential factors are analyzed, such as binning magnitude, measurement errors of magnitude, sample size, magnitude span, minimum completeness magnitude and fore- and aftershocks. At last, reasonable suggestions are provided for using those equations properly. The equations ofAki(1965),Utsu(1965),Page(1968)andKijko and Smit(2012)are based on assumption that magnitudes are continuous random variables, and have no corrections for this, so these equations are not recommended here. For simplicity, the equations ofUtsu(1966)orTinti and Mulargia(1987)can be used, but magnitude span should be greater than 2.5 due to without finite maximum magnitude in the formulas. For researchers having capability to write code and calculate numerically,Weichert(1980)orBender(1983)’s algorithm could be utilized. Especially when it is required to apply data with different observation periods for varied magnitudes, the formula ofWeichert(1980)is recommended. This study contributes to more accurately understand and use different formulas of estimating b-value by maximum likelihood technique, which can be used as reference for peers.

Keyword: magnitude-frequency relation; b-value; maximum likelihood method; least squares method; binned magnitude; magnitude span
0 引言

震级-频度(G-R)关系, 自Gutenberg等(1944)提出至今, 在地震学研究的各个领域有着广泛的应用(Cornell, 1968; 黄玮琼等, 1989; Wyss et al., 2000; Schorlemmer et al., 2005; Helmstetter et al., 2007; Petersen et al., 2008, 2015; 易桂喜等, 2013; 刘静伟等, 2016)。而b值作为G-R关系中最重要的参数, 在地震活动性研究和地震危险性分析中影响巨大, 因此如何合理估计b值一直以来都是地震学研究的热点(Aki, 1965; Utsu, 1965; Page, 1968; Weichert, 1980; 张建中等, 1981; Bender, 1983; Tinti et al., 1987; 秦长源, 2000; Kijko et al., 2012)。

常用的估计b值的方法为最小二乘法(Least Squares Method)(黄玮琼等, 1989; 陈培善等, 2003; 易桂喜等, 2008, 2013; 雷建成等, 2010; 邵延秀等, 2011), 使用极大似然法(Maximum Likelihood Estimation)估计b值的文献相对较少(龙锋等, 2009; 吴果, 2014; 吴果等, 2014; 张盛峰等, 2014)。然而国际上很多研究表明用最小二乘法估计b值存在一定的不足和局限性, 并且会带来一定的偏差, 在特定条件下推荐使用极大似然法作为替代或者补充(Suzuki, 1958; Page, 1968; Weichert, 1980; 张建中等, 1981; Dong et al., 1984; Sandri et al., 2007)。而中国对b值的极大似然法估计仍缺乏相对完整的分析和介绍。

本文简要地回顾了7个最具影响力的采用极大似然法估计b值及其标准差的方程式, 这些方程式的隐含假设和求解方式各不相同。针对是否考虑震级的归档效应(binned magnitude)、 是否设定有限最大震级、 是否对不同震级档的地震数据取不同的观察时段以及是否存在解析解这4个方面, 对这些方程式进行了统一的分类和总结。

b值的估计容易受到多种因素的影响, 国内对此有一些研究先例。例如黄玮琼等(1989)研究了G-R关系的两端 “ 掉头” 现象对于b值拟合的影响; 秦长源(2000)研究了地震的震级误差对b值的影响; 潘华等(2006)分析了样本统计时段和样本处理方法等的影响; 任雪梅等(2011)研究了地震目录的完整性对b值的影响。但是, 目前仍缺乏系统和全面的总结和分析。本研究全面回顾了前人关于极大似然法估计b值的影响因素的相关论述, 总结了震级的归档效应、 震级的测量误差、 样本量、 震级跨度、 最小完整震级和前余震共6个因子的作用方式, 并提出了相应的规避建议。

基于上述分析和总结, 本文还给出了对各种方程式的使用建议, 以期为广大地震工作者深入理解和准确应用各种不同的极大似然估计b值的方程式提供参考。

1 G-R关系
1.1 G-R关系简介

G-R关系, 即震级-频度关系最早由Gutenberg等(1944, 1949)提出。某地区在某段时间内发生震级≥ M的地震数N(M)可以用方程式(1)表示:

logN(M)=a-bM(1)

参数 a表示该地区的地震活动水平, 参数 b描述大、 小地震之间的比例关系。G-R关系在地震研究中具有极其重要的作用, 因为其不仅可以用来描述构造地震, 也可以用来描述诱发地震, 同时其对于不同的时间尺度和较大的震级间隔也同样适用(Kijko et al., 2012)。如果令 α=a·ln10=2.302i6a, β=b·ln10=2.302i6b, 则方程式(1)可以写为

lnN(M)=α-βM(2)

1.2 震级M的概率密度函数

假设某地区在某段时间内发生的地震的震级 MminMMmax, 那么由方程式(2)可以得到该震级范围内总的地震数:

Ntotal=N(Mmin)-N(Mmax)=eα-βMmin-eα-βMmax(3)

进一步构建震级M的概率分布函数:

F(M)=N(Mmin)-N(M)Ntotal=N(Mmin)-eα-βMNtotal(4)

那么, 如果假设震级M为1个连续随机变量, 对 F(M)M的导数可得到其概率密度函数 f(M):

f(M)=dF(M)dM=βe-βMe-βMmin-e-βMmax,  MminMMmax(5)

f(M)也可以写成以下形式:

f(M)=bln(10)10-bM10-bMmin-10-bMmax,  MminMMmax(6)

MmaxMmin时, 一般要求 Mmax-Mmin3(Marzocchi et al., 2003), 方程式(5)可以简化成以下形式:

f(M)=βe-β(M-Mmin)(7)

相应地, 方程式(6)可以写为

f(M)=bln(10)·10-b(M-Mmin)(8)

当震级范围较小时, 这一简化将不再适用(Knopoff, 2000)。

1.3 b值的含义

自方程式(1)被提出至今, b值被广泛地应用于各种地震学研究, 尤其是地震活动性模拟(Ogata et al., 2006; Felzer, 2009)、 地震预测(Kagan et al., 1987; Helmstetter et al., 2007)和地震危险性及风险性评估(Cornell, 1968; Beauval et al., 2004; Petersen et al., 2008, 2015; 潘华等, 2013; 李正芳等, 2014)中。

直接从方程式中理解可知, b值指示大、 小地震之间的比例关系。假设有40i000条M4.0以上地震的记录, 当b =1.0时, M7.0以上地震记录为10条; 当b=0.9时, M7.0以上地震记录则为20条, 为前者的2倍。Bender(1983)通过数组地震危险性分析示例发现, 当地震总数不变, b值发生一定幅度的变化时, 概率地震危险性图中的地震动加速度值变化的比例更大。由此可见, b值的影响巨大。

一些研究表明b值具有清晰的物理含义, 并且b值会随着时间和空间变化。例如 Mogi(1962a, b, 1963)通过实验证实, 对于前震序列和余震序列, 或者2个区域不同或震源深度不同的余震序列, b值都会存在差异。这一差异与材料属性或者震源区的应力分布有关。Papazachos等(1967)的研究表明在希腊附近区域, 随着震源深度的增加, 余震序列的b值会随之减小。Scholz(1968)的研究表明, 在实验室的岩石压力试验中, 微观破裂的发生也遵从类似的规律, 其b值与压应力呈负相关。对于一些特定的地震序列, b值会在主震前后发生变化, 与余震相比, 前震序列的b值更小(Suyehiro et al., 1964; Suyehiro, 1966; Gibowicz, 1973; Robinson, 1979; Smith, 1981)。这一特性与初始的高应力状态随着主震的应力降而降低相一致, 该结论被认为有助于地震预测(Scholz, 1968; 李全林等, 1978)。Cao等(2002)的研究显示, 日本岛弧东北部的b值从1984— 1990年的0.86降低到1991— 1995年的0.73, 造成该现象的原因可能与该区域的地震矩释放速度增加和火山活动更加频繁有关。

还有一些研究表明G-R关系只在一定的震级范围内有效(黄玮琼等, 1989; Pacheco et al., 1992a, b; Scholz, 1997; Triep et al., 1997; Knopoff, 2000)。Schwartz等(1984)提出了特征地震模型, “ 特征地震” 是指在同一震源区或特定断裂上重复发生的大小相差不大的地震。Youngs等(1985)认为对于独立的断层, 特征地震模型比G-R模型更能有效反映其地震活动规律, 故建议对于研究程度较高的具有发生中高震级地震能力的断层采用特征地震模型。类似的关于b值和G-R关系的研究过于繁杂, 这里不再展开介绍。

综上所述, b值不仅指示大、 小地震之间的比例关系, 同时也与地应力的状态密切相关。 b值高, 代表地应力水平较低, 而b值低, 则代表地应力水平较高。而另一些研究则表明G-R关系只适用于一定的震级范围。由此可见, b值含义丰富, 应用广泛, 影响深远, 合理估计b值显得尤为重要。

2 最小二乘法估计b值的局限性及其影响

最小二乘法(Least Squares Method)简单易行, 是回归分析中估计参数常用的一种方法, 被大量用于求取G-R关系中的b值(Pacheco et al., 1992a, b; Karnik et al., 1993; Okal et al., 1995; Scholz, 1997; Triep et al., 1997; Main, 2000; Working Group on California Earthquake Probabilities, 2003)。然而, 很多研究认为其存在一定的不足和局限性, 这可能会为结果带来一定的偏差, 因此建议在一定条件下采用极大似然法作为替代或者补充(Suzuki, 1958; Page, 1968; Weichert, 1980; 张建中等, 1981; Shi et al., 1982; Bender, 1983; Dong et al., 1984; Sandri et al., 2007)。

一般而言, 常见的用来估计b值的最小二乘法是一元线性最小二乘法, 其拟合对象分为2种: 第一种为累计震级-频度关系, 即如方程式(1)或方程式(2)所示, N(M)表示震级≥ M的地震数; 第二种为非累计震级-频度关系, 方程式形式与累计震级-频度关系相同, 但N(M)的含义变为震级M所在震级档( M-ΔM/2M< M+ΔM/2)内的地震数。2种震级-频度关系的a值不同, 但是b值相同(Sandri et al., 2007)。下面在介绍最小二乘法拟合b值的局限性及其影响时, 将针对这2种情况分别予以讨论。

(1)大、 小地震权重不一致: 累计和非累计震级-频度关系, 都需要划分震级档来获取每个震级档的数据点(Mi, log(N(Mi)))。对于普通的最小二乘法, 每个数据点(Mi, log(N(Mi)))之间都是等权重的。然而, 现实中越小的地震个数越多, 越大的地震个数越少, 大、 小地震之间的个数符合指数关系。随着震级档Mi由小到大变化, log(N(Mi))包含的地震数逐渐减少。所以越大的地震在最小二乘法中所占的权重越大, b值很容易受少数几个大地震变化的影响, 而在极大似然法中所有地震都具有相等的权重(Suzuki, 1958)。使用加权最小二乘法则会依赖更多不能验证的假设(Weichert, 1980)。

(2)地震样本量的限制: 如上所述, 最小二乘法需要划分震级档来获取每个震级档的数据点(Mi, log(N(Mi)))。当地震样本量较小时, 如果震级档间隔( ΔM)划分过大, 总的数据点数则会不足。而震级档间隔划分过小, 则每个震级档内的地震数可能会不足或者缺失。尤其当拟合对象为非累计震级-频度关系时, 随着震级不断增大, 极有可能出现某个震级档内地震数N(Mi)=0的情况, 这时其对数形式 log(N(Mi))则无意义, 在最小二乘法中只能被忽略, 这会导致拟合的b值偏小(Weichert, 1980)。极大似然法的原理不同于最小二乘法, 虽然其在拟合b值时样本量同样不能太小(详见后文), 但是不需要人为划分震级档, 可以在一定程度上避免上述问题。

Sandri等(2007)从理论分析和数值拟合2个方面论述了最小二乘法拟合b值带来的影响, 还将最小二乘法应用的对象分为累计和非累计震级-频度关系2种。蒙特卡洛模拟结果显示: 将最小二乘法应用于G-R关系的累计形式时, 地震数的累计过程会产生一种过滤效应, 这会导致对b值的不确定性被低估。如图 1所示, 由于需要对地震数取对数形式, 无论是累计还是非累计震级-频度关系, 最小二乘法拟合的b值都偏小, 其中前者的偏差相对较小, 随着地震数的增加偏差减小。与对地震数取对数形式带来的误差相比, 震级自身的测量误差可以忽略不计。

图 1 用最小二乘法计算1i000个人工拟合地震目录的b值的平均值(虚线)随样本量的变化曲线
拟合对象为累计(a)和非累计(b)震级-频度关系。每个平均值附95%的置信区间, 实线代表b值的真实值(图a、b中真实b值均为1, 改自Sandri et al., 2007)
Fig. 1 Average of b-value calculated through least squares technique in 1i000 synthetic catalogs as a function of the catalog size for the cases of cumulative(a)and binned(b)G-R relationship.

前人对最小二乘法和极大似然法拟合b值的效果进行过对比分析, Weichert(1980)通过对加拿大西部3个地震带的数据进行分析发现, 对于受到较好约束的数据, 与极大似然法相比, 最小二乘法能得到相近的结果。但是, 对于质量较差的数据, 由于受主观绘图和权重取值的影响, 最小二乘法得到的结果则与极大似然法存在较为显著的差异。总体上, 极大似然法的效果较好, 同时最小二乘法不能反映人们预先对最大震级的判断, 因此建议使用极大似然法估计b值。

张建中等(1981)用蒙特卡洛模拟的方法比较了极大似然法、 线性最小二乘法和非线性最小二乘法共3种算法的均方误差、 运算量和相关性。结果显示, 极大似然法估计b值的均方误差最小, 其次为线性最小二乘法, 均方误差最大的是非线性最小二乘法; 三者运算量的大小与均方误差的规律相同; 在与极大似然法结果的相关性上, 线性最小二乘法同样优于非线性最小二乘法; 随着样本量N的不断增大, 3种方法结果的差异逐渐减小。

上述总结和分析重点介绍了利用最小二乘法求取b值的一些局限性。事实上, 最小二乘法和极大似然法互有优劣, b值的最小二乘法估计的最大优点是简单易行, 且结果与极大似然法差异不大。当观察点(xi, yi)之间相互独立, 随机误差 εi~N(0, σ2)时, 最小二乘法和极大似然法给出的结果是相同的。但是在特定条件下, 例如样本量较小时, 极大似然法可以作为最小二乘法的一种有效的替代和补充方法。而国内对极大似然法的介绍相对较少, 因此后文重点探讨基于极大似然法拟合b值的重要方程式的隐含假设和应用条件, 以及各种可能对结果带来不确定性的影响因素。

3 极大似然法估计b值及其标准差

基于极大似然法估计b值并分析其不确定性一直都是地震学者研究的热点, 历史上各研究者提出了众多的方程式。这些方程式基于的假设条件有的明显, 有的不明显, 适用范围也各不相同, 极易混淆。本文先以时间顺序重点介绍前人的代表性工作, 然后再根据不同方程式的隐含假设和求解方式, 统一进行归类和总结分析。

3.1 Aki(1965)和Utsu(1965)

Aki(1965)最早提出用极大似然法估计b值, 并给出了方程式。用极大似然法推导b值的过程如下: 假设震级M是1个连续随机变量, 当MmaxMmin时, 其概率密度函数可以用方程式(7)来描述。则 Mi发生的概率为

f(Miβ)=βe-β(Mi-Mmin)(9)

震级样本( M1, M2, M3, …, MN)的联合概率密度函数为

f(M1, , MNβ)=i=1Nβe-β(Mi-Mmin)(10)

构建似然函数, 令:

L=ln[f(M1, , MNβ)]=i=1N[ln(β)-β(Mi-Mmin)](11)

令似然函数对 β的偏导数为0:

Lβ=i=1N1β-(Mi-Mmin)=0(12)

可以得到 β的极大似然法估计值

β˙=11Ni=1NMi-Mmin=1M̅-Mmin(13)

式中, M̅=1Ni=1NMi为样本均值, 于是得到b值的极大似然估计值:

b˙=log(e)M̅-Mmin(14)

由于只有记录完整的地震数才符合G-R关系, 所以 Mmin通常被定义为最小完整震级(Minimum Completeness Magnitude, 简记为Mc), 有时也被称为临界震级(Threshold Magnitude)。因此, 方程式(14)通常被记为

b˙=log(e)M̅-Mc(15)

式中, Mc为最小完整震级。 b值的极大似然估计是它的渐进无偏估计(Page, 1968; 张建中等, 1981)。Aki(1965)利用中心极限定理推得 β服从高斯分布, 均值为其极大似然估计值 β˙, 标准差为 β˙/N, 从而可以得到b值的标准差为

σ˙b˙=b˙N(16)

Utsu(1965)利用总体(Population)和样本(Sample)的一阶矩(First Moment)相等的方法得到了相同的方程式。同样也假设震级M是连续随机变量, 当MmaxMmin时, 由方程式(7)可以得到震级M总体的一阶矩:

E(M)=MminM·βe-β(M-Mmin)dM=Mmin+1β(17)

而震级样本( M1, M2, M3, …, MN)的样本矩 M̅=1Ni=1NMi, 令 M̅=E(M), 同样可以得到方程式(13)。因此, 方程式(15)有时被称为b值的Aki-Utsu估计值(Kijko et al., 2012)。

3.2 Utsu(1966)

Aki(1965)和Utsu(1965)得到的估计值基于的假设为震级是连续的随机变量, 同时没有设定震级的最大值。然而, 实际中的震级并非连续的随机变量。由于地震震级的精度有限, 对于仪器记录的地震, 通常为0.1个震级单位, 对于历史地震记录, 则≥ 0.5个震级单位。地震震级通常以 “ 0.1” 或者 “ 0.5” 为间隔跳跃式显示, 如 “ M6.0” 、 “ M6.1” 、 “ M6.2” , 相邻的震级被归并到与之最近的 “ 震级bin” (Magnitude bin)中。例如, 震级 Mi实际代表 Mi-ΔM/2M< Mi+ΔM/2范围内的地震, 其中 ΔM为 “ 震级bin” 的宽度。通常在英文文献中用 “ Binned magnitude” 或 “ Grouped in magnitude classes” 来描述这一现象, 本文将其翻译为 “ 震级的归档效应” 。

因此, 假设震级精度为 ΔM, 那么最小完整震级 Mc实际代表的震级范围为 Mc-ΔM/2M< Mc+ΔM/2。Utsu(1966)最早提出可以对方程式(15)进行修改, 用真实的最小完整震级 Mc-ΔM/2代替 Mc, 相应地方程式(15)变为

b˙=log(e)M̅-(Mc-ΔM/2)(18)

虽然方程式(18)提高了估计b值的精度, 但是在实际中, 应用最广的仍然是方程式(15)(Utsu, 1999)。还存在使用方程式(18), 但是却引用Aki(1965)的研究作为参考文献的情况(Gerstenberger et al., 2001; 黄亦磊等, 2016)。方程式(18)的一个重要应用是用Matlab编写的地震活动性分析软件ZMAP(Wiemer et al., 2000; Wiemer, 2001; Woessner et al., 2005)。

3.3 Page(1968)

Aki(1965)和Utsu(1965, 1966)的方程式有一个共同的问题, 即没有设定震级上限 Mmax, 因此他们的方程式无法体现人们对于某地区最大震级的认识。Page(1968)同样假设震级M为连续的随机变量, 但是引入了有限最大震级 Mmax, 推导了此条件下用于求解b值的极大似然估计方程式:

1β=M̅-Mmin-Mmaxexp(-β(Mmax-Mmin))1-exp(-β(Mmax-Mmin))(19)

式中, M̅是震级样本均值, β=b·ln10=2.302i6bMmax-Mmin> 2.0时, 该方程式可以简化为Aki(1965)和Utsu(1965)的形式。Page(1968)还给出了求解b值95%置信区间的方程式:

N1b'2+Mspan22-eb'Mspan-e-b'Mspan-1/2·1b'+Mmin-Mmaxe-b'Mspan1-e-b'Mspan-M̅=±1.96(20)

式中, b'=bln(10), Mspan为震级跨度(Magnitude Span), 定义为经过筛选的最小完整震级以上的地震数据的最大震级减去最小震级, 即 Mspan=Mmax-Mmin

3.4 Weichert(1980)

人类通过仪器记载的地震目录不过只涵盖数十年至一百年, 这与大地震几百年乃至上千年的复发周期相比显得十分短暂。同时, 一个地区的地震监测能力往往随着时间的推移逐步提高, 所以震级越小的地震具有相对完整记录的时间越短。因此, 有限的地震目录资料显得尤为宝贵。为了充分利用地震目录资料, Weichert(1980)将大小不同的震级档对应不同的完整时段, 对地震目录进行分档并予以利用。

图2 不同震级档对应不同观察时段示意图Fig. 2 A schematic illustration of unequal observation periods for different magnitude intervals.

如图 2所示, 假设震级M满足 MminMMmax, 将震级范围划分为J个震级档 M1, M2, M3, , MJ, 震级档Mj实际代表 Mj-ΔMj2M< Mj+ΔMj2内的所有地震。各个震级档的地震记录对应的完整时段的时长分别为 T1, T2, T3, , TJ, 每个震级档在完整时段内记录到的地震数分别为 (k1, k2, k3, , kJ)据此, Weichert推导了新的基于极大似然法估计 β(b·ln10)的方程式:

jTjMjexp(-βMj)jTjexp(-βMj)=kjMjN(21)

式中, N为各个震级档在相应完整时段内记录到的地震数的和, N=j=1Jkj等式右侧即样本均值 M̅式中 β没有解析解, 可以用牛顿迭代法求解。或者采用枚举法, 用间隔为0.001的不同 β进行尝试, 取令方程式(21)两侧之差最小的 β值为结果。 β相应的标准差为

σ˙β˙=1N(jθj)2(jθjMj)2-jθjjθjMj2(22)

式中, θj=Tjexp(-β˙Mj)

方程式(21)基于3个基本假设: 1)将地震目录按震级档予以划分, 震级档的跨度可以各不相等, 也可以等间隔( ΔMj=ΔM, j=1, , J)。这里的震级档 Mj还可以用来处理震级的归档效应(Binned magnitude), 因此该方程式针对的震级是实际中被归档过的随机变量, 而不是原始的连续随机变量; 2)对于不同的震级档, 根据地震记录完整的时间不同, 取不同的观察时段 Tj; 3)该方程式可以设定最大震级, 令 Mmax=Mmin+j=1JΔMj即可。

3.5 Bender(1983)

Bender(1983)同样假设震级为归档后的随机变量, 将震级范围 MminMMmaxΔM为间隔划分为J个震级档, 震级档 Mj实际代表 Mj-ΔM/2M< Mj+ΔM/2内的所有地震。假设一共记录到N个地震 (M1, M2, M3, , MN), 将其按照震级大小分配到J个震级档中, 每个震级档内的地震数记为 (k1, k2, k3, , kJ)。推导得到用于求解β(b·ln10)的极大似然估计值的方程式为

q1-q-JqJ1-qJ=j=1J(j-1)kjN(23)

式中, q=exp(-βΔM), kj为震级档 Mj内的地震数, 有 j=1Jkj=N该方程式也只能通过数值方法求解。

Bender(1983)的方程式基于的假设与上述Weichert(1980)的方程式(21)基本相同, 区别在于: 1)方程式(23)以等间距 ΔM划分震级档, 而方程式(21)可以按照等间距或不等间距划分震级档; 2)方程式(23)中每个震级档都取同样的完整时段, 因此方程式中没有 Tj这一因子, 而方程式(21)针对每个震级档可以取不同的完整时段。基于这2点分析, 方程式(23)似乎是方程式(21)的一种特殊情况, 功能没有后者强大。但是Bender(1983)强调方程式(23)可以合理地考虑最大震级、 震级的归档效应和没有观察值的震级档, 只有当样本量较大时才能体现方程式(21)的应用价值, 而当样本量较小的时候, 使用方程式(23)即可。

3.6 Tinti等(1987)

Tinti等(1987)假设震级为等间隔 (ΔM)取值的归档后的随机变量, 同时不设定可能的最大震级, 推导了新的方程式。 b值的极大似然估计值为

b˙=1ln(10)ΔMln1+ΔMM̅-Mc(24)

相应的b值的标准差为

σ˙b˙=1-pln(10)ΔMNp(25)

式中, p=1+ΔMM̅-Mc。当 ΔM较小时, 方程式(24)和(25)可以在样本量非常小(N≈ 30)的情况下, 较好地估计b值及其标准差, 其另1个优点是存在解析解, 不需要进行数值拟合。

3.7 Kijko等(2012)

Kijko等(2012)认为虽然Weichert(1980)的方程式可以使用不同震级档对应不同完整记录时段的数据, 但是方程式(21)没有解析解, 不便于使用。于是, 他们在方程式(15)(Aki, 1965; Utsu, 1965)的基础上推导了新的方程式, 使用的假设条件也与方程式(15)相同, 即假设震级为连续随机变量, 同时没有设定震级上限。

图3 根据最小完整震级随时间变化划分次级地震目录示意图(改自Kijko et al., 2012)Fig. 3 A schematic illustration of the division of subcatalogs based on variation of minimum completeness magnitude versus time(Adapted after Kijko et al., 2012).

Kijko等(2012)按照地震监测水平的差异, 将总记录时长为T的地震目录分成L个时段, 每个时段的次级地震目录对应时长为 Tl, 有 l=1LTl=T, 如图 3所示。由于地震记录的完整性水平随着时间变化, 每个次级地震目录的最小完整震级 Mcl各不相同, 在最小完整震级以上记录到的地震数分别为 Nl, 有 l=1LNl=N, N为总目录中的地震数。推导得到总目录的 β的极大似然估计值为

β˙=N·N1β˙1+N2β˙2+N3β˙3++NLβ˙L-1(26)

式中, (β˙1, β˙2, β˙3, , β˙L)为各个次级地震目录的 β值, 可以用方程式(13)求解。 β的标准差为

σ˙β˙=β˙N(27)

对比图2和图3可以发现, 方程式(26)和方程式(21)不仅基于的基本假设不同, 而且在处理不同大小的震级对应不同完整时段时的策略也不同。方程式(21)是将地震目录分成多个震级档, 每个震级档可分别取不同时长的地震记录数据。而方程式(26)则根据地震完整记录水平随着时间的变化, 将地震目录分成多个次级地震目录。Kijko等(2012)开发了基于Matlab平台的程序包AUE, 可以方便地使用方程式(26)估计 β值, 同时还可以估算最大可能震级和年发生率等地震活动性参数。

3.8 估计b值标准差的重要方程式

前文在介绍各种基于极大似然法估计b值的方程式时, 已经给出了与之相应的求解b值标准差的方程式, 其中应用最广泛的是方程式(16)(Aki, 1965)。这里再补充介绍2个求解b值标准差的重要方程式。

Utsu(1966)最早在其研究中构建χ 2分布来研究b值, 张建中等(1981)进一步通过构建 χ2分布研究b值的标准差, 得到:

σ˙b˙(ZS)=bN+2(N-1)(N-2), N> 2(28)

并且给出了b值在置信水平为 (1-α)时的区间估计:

χα/22(2N)/2NM̅< b< χ1-α/22(2N)/2NM̅(29)

式中, χα2(2N)为自由度等于 2N概率为 αχ2分布的下侧分位点。需要指出的是, 方程式(28)基于的假设与方程式(16)相同, 即震级是连续随机变量且未设最大震级。

另1个重要的方程式是Shi等(1982)提出的:

σ˙b˙(SB)=2.30·b2i=1N(Mi-M̅)2N(N-1)(30)

ZMAP软件在计算b值的标准差时即使用方程式(30)(Woessner et al., 2005)。

由于绝大部分文献都没有直接提及G-R关系中的a值及其标准差的计算, 故在此作简单介绍。在求得b值的极大似然估计值 b˙及其标准差 σ˙b˙的基础上, 方程式(1)中的参数a可用方程式(31)求取(De Santis et al., 2011):

a˙=log(N)+b˙·Mc(31)

式中, N为样本量, Mc为最小完整震级。在已知震级精度 ΔM的情况下, 通过高斯误差传播理论可以求得相应的a值的标准差为

σ˙a˙=(Mc·σ˙b˙)2+(b˙·ΔM)2(32)

3.9 估计b值的方程式分类

本文一共总结了7个最常见的基于极大似然法估计b值的方程式, 这些方程式的假设和应用条件各不相同, 现对其予以分类和总结。

这些方程式大致有3个隐含假设, 每个假设有相互对立的2种选择: 1)假设震级M为归档后的随机变量(Binned Random Variable, 简称BRV, 详细论述见3.2节), 否则为通常的连续随机变量(Continuous Random Variable); 2)设定有限最大震级(Mmax), 否则Mmax趋向于无穷大, 或者要求MmaxMmin; 3)不同的震级档对应不同的观察时段(记为Ti), 否则每个震级档的观察时段相等, 即Ti=T

此外, 方程式的求解方式也分为存在解析解(Analytic solution, AS)和只能求数值解(Numerical solution, NS)。7组方程式的假设条件和求解方法具体见表1, 若相应参考文献中明确给出了b值或者 β的标准差方程式时, 其也被归纳进表1中。

表1 基于极大似然法估计b值的方程式及其基本假设和求解方式总结 Table1 Summary on the implicit assumptions and the way of solving the equations estimating b-value by maximum likelihood method

在这些方程式中, 应用最广的是假设震级为连续随机变量的方程式(15)。那些假设震级为归档后的随机变量的方程式可能由于求解困难反而应用很少(Marzocchi et al., 2003), 如Weichert(1980)和Bender(1983)提出的方程式都只有数值解。

求解b值标准差的方程式有时会和估计b值的方程式搭配使用, 例如在ZMAP软件中使用方程式(18)(Utsu, 1966)估计b值, 而相应标准差则用Shi等(1982)提出的方程式(30)来估算。

4 利用极大似然法估计b值的影响因素

b值在地震活动性预测和地震危险性分析中有着极其重要的作用, 然而其自身会受到多种因素的影响, 例如样本量、 样本的震级跨度以及余震等。本文总结前人的工作, 现将各种可能影响极大似然法估计b值的因素一一予以介绍。

4.1 震级的归档效应

表1所示, 很多用极大似然法估计b值的方程式都假设震级是连续随机变量, 然而这与实际情况是不相符的。如前文所述(见3.2节), 现实中的震级存在归档效应(Binned magnitude), 相邻的震级被归并到离其最近的 “ 震级bin” 。例如, 震级 Mi实际代表 Mi-ΔM/2M< Mi+ΔM/2范围内的地震, 其中 ΔM为 “ 震级bin” 的宽度。

Marzocchi等(2003)指出, 如果简单地假设震级为连续随机变量会带来2类偏差: 第一类, 为幂律分布的连续随机变量的震级均值 μ与同一样本但震级被归并到各个 “ 震级bin” 的震级均值存在差异; 第二类, 最小完整震级 Mc与真实值之间也存在一定差异。为了便于区分, 将这2类偏差分别记为 θ1θ2

Bender(1983)对第一类偏差进行了详细的论述。其指出归并到 “ 震级bin” 后的震级均值要系统地大于连续随机变量的震级均值。因为地震数在1个 “ 震级bin” 内并不是均匀分布的, Mi-ΔM/2~Mi内的地震数比 Mi~Mi+ΔM/2内的地震数多。当 ΔM=0.1时, 对结果的影响较小, 可以忽略不计。随着 ΔM变大(例如0.6), 影响将会变得十分显著。

对于第二类偏差, Utsu(1966)最早提出可以对方程式(15)进行1个小的修改。由于最小的 “ 震级bin” , 即最小完整震级 Mc所对应的 “ 震级bin” , 实际包含 Mc-ΔM/2 M< Mc+ΔM/2内的所有地震, 于是可以用真实的最小完整震级 Mc-ΔM/2代替 Mc, 相应的方程式(15)变为方程式(18)(见3.2节)。

Bender(1983)指出对于同样的样本, 采用不同的方程式估计b值, 结果将显示出明显的区别。当震级间隔 ΔM=0.1时, 用于处理连续随机变量的方程式(15)(Aki, 1965; Utsu, 1965)和方程式(19)(Page, 1968)能得到与方程式(23)近似的结果。然而, 当震级间隔 ΔM较大时, 例如 ΔM=0.6时, 方程式(15)和方程式(19)得到的结果则明显偏小。

Marzocchi等(2003)将方程式(24)(Tinti et al., 1987)与Aki(1965)的方程式(15)以及Utsu(1966)的方程式(18)进行了系统的对比。如图 4所示, Aki(1965)的方程式估计的b值存在严重的偏差。Utsu(1966)的方程式虽然只做了1个简单的矫正, 却大大提高了精度, 在不同样本量下不存在显著的偏差。方程式(24)与Utsu(1966)的方程式相比, 当真实值b=1.0时, 结果十分近似。当真实值b=2.0时, 方程式(24)的结果略优。

图 4 用3个方程式计算1i000个人工拟合地震目录的b值中位数和95%置信区间随样本量的变化
虚线、 点线和粗实线分别对应方程式(15)(Aki, 1965; Utsu, 1965)、 方程式(18)(Utsu, 1966)和方程式(24)(Tinti et al., 1987)。细实线代表b值的真实值, 图a和b的真实b值分别为1.0和2.0(改自Marzocchi et al., 2003)
Fig. 4 Medians and 95%confidence bands of b-value calculated through 3 kinds of equations from 1i000 synthetic catalogs, as a function of the catalog size.

4.2 震级的测量误差

震级的测量误差不同于上述震级的归档效应。震级的归档是指震级自身作为1个连续随机变量被强行归并到各 “ 震级bin” 中, 类似于 “ 四舍五入” 的过程。而震级的测量误差是指地震台站的在测量震级时的所包含的随机误差 ε, 该误差无法消除。如果假设 “ 真” 震级为 Mreal, 其严格按照方程式(1)分布, 不带随机误差, 那么对于测量得到的 “ 视震级” M(目录震级)有以下关系式:

M=Mreal+ε(33)

式中, ε为震级的测量误差, 一般取高斯分布(Tinti et al., 1985; Rhoades, 1996; Miao et al., 2007; Sandri et al., 2007)。

Marzocchi(2003)通过理论和数值模拟2种方式研究了震级的归档效应和震级的测量误差对于极大似然法估计b值的影响, 结果显示与对震级的归档效应不当处理相比, 例如方程式(15)(Aki, 1965)直接假设震级为连续随机变量, 震级的测量误差带来的影响可以忽略不计。这可能也是大部分作者在拟合估计b值的方程式时, 没有特别关注震级测量误差的原因(Weichert, 1980; Kijko et al., 2012)。Sandri等(2007)在研究地震数取对数形式和震级的测量误差对于最小二乘法估计b值的影响时, 也有类似的结论, 与前者相比, 震级的测量误差带来的影响可以忽略不计。

4.3 样本量

影响极大似然法估计b值的另外一个重要因素是样本量, 需要指出的是这里的样本量是指经过筛选的最小完整震级(Mc)以上的地震数, 而不是直接下载得到的地震目录中的地震数。

样本量早期就被用于估计b值的标准差。Aki(1965)利用中心极限定理推导得出 β服从高斯分布, 均值为其极大似然估计值 β˙, 标准差为β˙/N, N为样本量。需要注意的是, 该方程式不能用于地震数过少的情况(Weichert, 1980)。事实上, Aki(1965)给出的示例表中 N50

Shi等(1982)用方程式(15)(Aki, 1965; Utsu, 1965)估计美国加州中部的b值, 结果显示, 对于该地区, 要得到相对稳定的结果, 筛选样本所需的时空窗中大约需要包含100个地震。Bender(1983)指出当样本量过小时, 无论使用何种方程式, 结果都将变得不可信。当样本量N=100时, b值的标准差 σb0.1b准确估计b值至少需要25个震级样本, 对应的 σb约小于0.25b

Kutliroff(2017)用蒙特卡洛模拟方法对方程式(18)(Utsu, 1966)和方程式(24)(Tinti et al., 1987)进行检验, 发现二者的结果十分相近, 当样本量 N> 25时, 结果较为可信。

4.4 震级跨度

早期的研究者没有针对震级跨度进行专门的研究, 然而一些方程式的推导过程事实上包含对震级跨度的隐含假设。如1.2节所述, MmaxMmin, 震级的概率密度函数可以由方程式(5)简化为没有最大震级项的形式(方程式(7))。若满足这一条件, 则一般要求Mmax-Mmin ≥ 3(Marzocchi et al., 2003)。所以, 当使用基于方程式(7)推导的方程式时, 例如方程式(15)(Aki, 1965; Utsu, 1965)、 方程式(18)(Utsu, 1966)以及方程式(26)(Kijko et al., 2012), 所使用的地震数据要满足震级跨度≥ 3.0这一条件。

另一处对震级跨度有隐含假设的是Page(1968)推导的带有最大震级项的方程式(19), 其指出当Mmax-Mmin> 2时, 该方程式可以简化成Aki(1965)和Utsu(1965)提出的不带最大震级项的形式。

Kutliroff(2017)首先对震级跨度的影响做了系统的分析。蒙特卡洛模拟显示(图5), 对于方程式(18)(Utsu, 1966)和方程式(24)(Tinti et al., 1987), 当样本量的影响可以忽略时(N ≥ 6i837), 要求震级跨度Mspan> 2.5, 结果才较为可信, 误差约在± 0.015倍真实值之间。而方程式(15)(Aki, 1965)的结果则系统地偏高。

图 5 3种方法拟合人工地震目录的b值随震级跨度的变化(改自Kutliroff, 2017)
bAbUbTM分别对应方程式(15)(Aki, 1965; Utsu, 1965)、 方程式(18)(Utsu, 1966)和方程式(24)(Tinti et al., 1987)。每个人工地震目录包含10i000个地震, b值为1; Mmin固定为2.0, Mmax在2.5~6.5之间变动, 且Δ M=0.1; 在震级范围内, 地震数在6i837~10i000之间变动
Fig. 5 Three estimates of the b-value of synthetic catalogs versus the magnitude span Mspan(Adapted after Kutliroff, 2017).

4.5 最小完整震级

本文在描述样本量N时, 特别指出是经过筛选的最小完整震级(Mc)以上的地震数。这是因为如果选择了错误的最小震级, 将会对b值的极大似然估计值产生显著的影响。如图6a所示, 在G-R关系的小震级端, 由于地震记录不全, 将导致该震级段曲线斜率变缓, 中文文献中习惯将这一现象称为 “ 掉头” (黄玮琼等, 1989; 任雪梅等, 2011)。如果错误地选择了过小的最小完整震级, 会导致算得的b值偏小。在Woessner等(2005)的示例中(图6b), 当真实的Mc被人为地从M1.4降低到M0.7时, 相应的b值从真实值1.02降低到0.8。

图 6 a 北加州地震台网(NCSN)地震目录的子目录的G-R关系; b b值(b)、 b值的均值(bave)和标准差(σ b)随小震级端的截断震级Mco 的变化(改自Woessner et al., 2005)
a 最小完整震级Mc由最大曲率法(MAXC)求得; b 箭头指向利用b值稳定性法(MBS)求得的Mc
Fig. 6 Frequency-magnitude distribution of the subset of the Northern California Seismic Network(NCSN)catalogue(a), and b-value, bave and the uncertainties σ b as a function of cutoff magnitude Mco(b)(Adapted after Woessner et al., 2005).

事实上, 根据b值是否能稳定取值还可以用来估计最小完整震级。Cao等(2002)提出了1种估计最小完整震级的方法: 如图6b所示, 将人为选择的小震级端的截断震级(Cutoff magnitude)记为Mco, 真实的最小完整震级记为Mc。当Mco< Mc时, 随着Mco 增大, b值不断增大; 当McoMc时, 随着Mco 增大, b值保持相对稳定; 当McoMc时, 随着Mco 增大, b值再次开始增大。如果使用的McoMc, 则求得的b值会明显偏小。黄亦磊等(2016)对比了多种分析地震目录完整性的方法, 结果显示Cao等(2002)的方法与多数其他方法相比效果较好。

由上述分析可知, 准确选择最小完整震级对于估计b值来说十分重要。而估计最小完整震级的方法多种多样(Woessner et al., 2005; Mignan et al., 2012; 黄亦磊等, 2016), 只利用单一方法, 例如G-R关系的线性程度进行估计, 往往存在缺陷, 准确估计Mc需要结合使用多种方法(任雪梅等, 2011; 吴果等, 2014)。

4.6 前余震

在1.3节论述b值的含义时已经指出, 很多前人的研究表明前余震会影响b值的取值, 这与估计b值的方法无关, 属于客观规律。例如: Mogi(1962a, b, 1963)通过实验证实, 对于前震序列和余震序列, 或者2个区域不同或震源深度不同的余震序列, b值都会存在差异。这一差异与材料属性或者震源区的应力分布有关。Scholz(1968)的研究表明, 在实验室的岩石压力试验中, 微观破裂的发生也遵从类似的规律, 其b值与压应力呈负相关。对于一些特定的地震序列, b值会在主震前后发生变化, 与余震相比, 前震序列的b值更小(Suyehiro et al., 1964; Suyehiro, 1966; Gibowicz, 1973; Robinson, 1979; Smith, 1981)。这一特性与初始的高应力状态随着主震的应力降而降低相一致, 该结论被认为有助于地震预测(Scholz, 1968; 李全林等, 1978)。

考虑到绝大部分的地震活动性研究和地震危险性分析都基于主震目录(Helmstetter et al., 2006, 2007, 2014; Wang et al., 2011; Petersen et al., 2015; Talebi et al., 2017), 本文建议在利用极大似然法估计b值之前, 先用合理的方法删除前余震, 常用的算法有Gardner等(1974)和Reasenberg(1985)等提出的方法。

4.7 影响因素总结

基于上述分析, 现对上述6个方面的影响因素的影响方式进行总结, 并提出合理的规避建议, 详见表2

表2 影响极大似然法估计b值的6种影响因素的作用方式和规避建议总结 Table2 Summary on the effects of 6 influential factors of maximum likelihood estimation of b-value, and suggestions to avoid them
5 结论

本文回顾了震级-频度(G-R)关系中震级M的概率密度函数以及b值的含义, 分析了国内常用的基于最小二乘法估计b值的一些局限性及其影响, 系统地总结了前人关于极大似然法估计b值及其标准差的研究成果, 并对各种方程式的隐含假设、 求解方式以及影响因素等进行了分析和归类, 主要结论如下:

(1)最小二乘法简单易行, 应用广泛, 但是利用其估计b值会受到一定的局限, 例如大、 小地震之间权重不一致, 样本量小时震级档的划分问题等。因此, 在特定条件下, 极大似然法可以作为一种较为有效的替代或者补充。

(2)本文将利用极大似然法估计b值的隐含假设分为3大类: 首先, 假设震级M为连续随机变量, 或考虑震级的归档效应; 其次, 方程式中存在有限最大震级项, 否则则要求最大震级远大于最小震级; 最后, 是否对不同震级档的地震数据取不同的观察时段。同时方程式的求解方式可以分为存在解析解与只有数值解。按照上述隐含假设和求解方式对7个方程式进行了分类, 详细情况见表1

(3)本文详细总结了影响极大似然法估计b值的重要因素, 一共分为6个方面: 震级的归档效应、 震级的测量误差、 样本量、 震级跨度、 最小完整震级和前余震。其中除了震级的测量误差, 其他因素都有可能对结果造成较大影响, 不容忽视。详情见表2

(4)基于本文的分析和讨论, 对7个方程式的推荐使用情况如下: 方程式(15)(Aki, 1965; Utsu, 1965)和方程式(19)(Page, 1968)都是基于假设震级为连续随机变量, 并且没有做任何矫正, 建议不要使用。而方程式(26)(Kijko et al., 2012)对每个次级地震目录的b值的求取是基于方程式(15)的, 故也建议不要采用。 为使计算简单方便, 可以使用方程式(18)(Utsu, 1966)或者方程式(24)(Tinti et al., 1987), 但是要求震级跨度Mspan> 2.5, 因为二者都没有有限最大震级项; 如果有编程进行数值计算的能力, 可以使用方程式(21)(Weichert, 1980)或者方程式(23)(Bender, 1983)。当需要对不同震级档的地震数据取不同的观察时段时, 可以使用方程式(21)。

(5)对前人研究b值的方法进行总结, 可以发现: 除了基础的理论推导, 最重要的研究方法是蒙特卡洛模拟(张建中等, 1981; Bender, 1983; Marzocchi et al., 2003)。蒙特卡洛模拟可以人工生成大量的参数取值不同的地震目录, 然后用多种方法分析各个目录的b值, 再与真实值作对比。这在研究震级的归档效应、 样本量、 震级跨度等因素的影响时能很好地发挥作用。

The authors have declared that no competing interests exist.

参考文献
[1] 陈培善, 白彤霞, 李保昆. 2003. b值和地震复发周期[J]. 地球物理学报, 46(4): 510519.
CHEN Pei-shan, BAI Tong-xia, LI Bao-kun. 2003. b-value and earthquake occurrence period[J]. Chinese Journal of Geophysics, 46(4): 510519(in Chinese). [本文引用:1]
[2] 黄玮琼, 时振梁, 曹学锋. 1989. b值统计中的影响因素及危险性分析中 b值的选取[J]. 地震学报, 11(4): 351361.
HUANG Wei-qiong, SHI Zhen-liang, CAO Xue-feng. 1989. Factors influencing the estimation of b-value and the selection of b-value in hazard analysis[J]. Acta Seismologica Sinica, 11(4): 351361(in Chinese). [本文引用:4]
[3] 黄亦磊, 周仕勇, 庄建仓. 2016. 基于地震目录估计完备震级方法的数值实验[J]. 地球物理学报, 59(4): 13501358.
HUANG Yi-lei, ZHOU Shi-yong, ZHUANG Jian-cang. 2016. Numerical tests on catalog-based method to estimate magnitude completeness[J]. Chinese Journal of Geophysics, 59(4): 13501358(in Chinese). [本文引用:2]
[4] 雷建成, 高孟潭, 吕红山, . 2010. 地震带划分方案对地震活动性参数的影响[J]. 地震学报, 32(4): 457465.
LEI Jian-cheng, GAO Meng-tan, Hong-shan, et al. 2010. Effects of seismic belt division scheme on seismicity parameters[J]. Acta Seismologica Sinica, 32(4): 457465(in Chinese). [本文引用:1]
[5] 李全林, 陈锦标, 于渌, . 1978. b值时空扫描: 监视破坏性地震孕育过程的一种手段[J]. 地球物理学报, 21(2): 101125.
LI Quan-lin, CHEN Jin-biao, YU Lu, et al. 1978. Time and space scanning of b-value: A method for monitoring the development of catastrophic earthquakes[J]. Chinese Journal of Geophysics, 21(2): 101125(in Chinese). [本文引用:2]
[6] 李正芳, 周本刚. 2014. 利用断裂带上的低 b值识别凹凸体方法的探讨: 以龙门山断裂带和鲜水河断裂带为例[J]. 震灾防御技术, 9(2): 213225.
LI Zheng-fang, ZHOU Ben-gang. 2014. Asperity identification based on low b-value: Application to the Longmenshan and Xianshuihe fault zone[J]. Technology for Earthquake Disaster Prevention, 9(2): 213225(in Chinese). [本文引用:1]
[7] 刘静伟, 吕悦军. 2016. 川滇地区 b值空间分布特征及其与震源类型关系的初步探讨[J]. 震灾防御技术, 11(3): 561572.
LIU Jing-wei, Yue-jun. 2016. Spatial distribution of b values and its relationship with the type of focal mechanism in the Sichuan-Yunnan area[J]. Technology for Earthquake Disaster Prevention, 11(3): 561572(in Chinese). [本文引用:1]
[8] 龙锋, 闻学泽, 倪四道. 2009. 区域最小完整性震级时空分布的确定: 以龙门山断裂带为例[J]. 地震, 29(3): 2736.
LONG Feng, WEN Xue-ze, NI Si-dao. 2009. Determination of temporal-spatial distribution of the regional minimum magnitudes of completeness: Application to the Longmenshan fault zone[J]. Earthquake, 29(3): 2736(in Chinese). [本文引用:1]
[9] 潘华, 高孟潭, 谢富仁. 2013. 新版地震区划图地震活动性模型与参数确定[J]. 震灾防御技术, 8(1): 1123.
PAN Hua, GAO Meng-tan, XIE Fu-ren. 2013. The earthquake activity model and seismicity parameters in the new seismic hazard map of China[J]. Technology for Earthquake Disaster Prevention, 8(1): 1123(in Chinese). [本文引用:1]
[10] 潘华, 李金臣. 2006. 地震统计区地震活动性参数 b值及 v4不确定性研究[J]. 震灾防御技术, 1(3): 218224.
PAN Hua, LI Jin-chen. 2006. Study on uncertainties of seismicity parameters b and υ4 in seismic statistical zones[J]. Technology for Earthquake Disaster Prevention, 1(3): 218224(in Chinese). [本文引用:1]
[11] 秦长源. 2000. 地震震级误差对 b值的影响[J]. 地震学报, 22(4): 337344.
QIN Chang-yuan. 2000. The effect of the magnitude uncertainty on the b-value[J]. Acta Seismologica Sinica, 22(4): 337344(in Chinese). [本文引用:2]
[12] 任雪梅, 高孟潭, 冯静. 2011. 地震目录的完整性对 b值计算的影响[J]. 震灾防御技术, 6(3): 257268.
REN Xue-mei, GAO Meng-tan, FENG Jing. 2011. Effect of completeness of earthquake catalogue on calculating b-value[J]. Technology for Earthquake Disaster Prevention, 6(3): 257268(in Chinese). [本文引用:3]
[13] 邵延秀, 袁道阳, 王爱国, . 2011. 西秦岭北缘断裂破裂分段与地震危险性评估[J]. 地震地质, 33(1): 7990. doi: DOI: 103969/j. issn. 0253-4967. 2011. 01. 008.
SHAO Yan-xiu, YUAN Dao-yang, WANG Ai-guo, et al. 2011. The segmentation of rupture and estimate of earthquake risk along the north margin of western Qinling fault zone[J]. Seismology and Geology, 33(1): 7990(in Chinese). [本文引用:1]
[14] 吴果. 2014. 基于多模型的概率地震危险性分析方法实例研究 [D]. 北京: 中国地震局地质研究所.
WU Guo. 2014. Case study of probabilistic seismic hazard analysis method based on multi-model [D]. Institute of Geology, CEA, Beijing(in Chinese). [本文引用:2]
[15] 吴果, 周庆, 冉洪流. 2014. 中亚地震目录震级转换及其完整性分析[J]. 震灾防御技术, 9(3): 368383.
WU Guo, ZHOU Qing, RAN Hong-liu. 2014. Magnitude conversion of earthquake catalog in central Asia and the completeness analysis[J]. Technology for Earthquake Disaster Prevention, 9(3): 368383(in Chinese). [本文引用:1]
[16] 易桂喜, 闻学泽, 苏有锦. 2008. 川滇活动地块东边界强震危险性研究[J]. 地球物理学报, 51(6): 17191725.
YI Gui-xi, WEN Xue-ze, SU You-jin. 2008. Study on the potential strong-earthquake risk for the eastern boundary of the Sichuan-Yunnan active faulted-block, China[J]. Chinese Journal of Geophysics, 51(6): 17191725(in Chinese). [本文引用:1]
[17] 易桂喜, 闻学泽, 辛华, . 2013. 龙门山断裂带南段应力状态与强震危险性研究[J]. 地球物理学报, 56(4): 11121120.
YI Gui-xi, WEN Xue-ze, XIN Hua, et al. 2013. Stress state and major-earthquake risk on the southern segment of the Longmen Shan fault zone[J]. Chinese Journal of Geophysics, 56(4): 11121120(in Chinese). [本文引用:2]
[18] 张建中, 宋良玉. 1981. 地震 b值的估计方法及其标准误差: 应用蒙特卡罗方法估计 b值精度[J]. 地震学报, 3(3): 7887.
ZHANG Jian-zhong, SONG Liang-yu. 1981. On the method of estimating b-value and its stand ard error: The Monte Carlo method of estimating the accuracy of b-value[J]. Acta Seismologica Sinica, 3(3): 7887(in Chinese). [本文引用:7]
[19] 张盛峰, 吴忠良, 房立华. 2014. 双差(DD)定位地震目录能用于地震序列的统计地震学参数计算吗?——云南鲁甸 MS6. 5地震序列 b值的空间分布[J]. 地震地质, 36(4): 12441259. doi: DOI: 103969/j. issn. 0253-4967. 2014. 04. 024.
ZHNAG Sheng-feng, WU Zhong-liang, FANG Li-hua. 2014. Can the DD-relocation earthquake catalogue be used for the statistical parameters of an earthquake sequence?A case study of the spatial distribution of b-value for the aftershocks of the 2014 Ludian MS6. 5 earthquake[J]. Seismology and geology, 36(4): 12441259(in Chinese). [本文引用:1]
[20] Aki K. 1965. Maximum likelihood estimate of b in the formula log N= a- bM and its confidence limits[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 43(2): 237239. [本文引用:12]
[21] Beauval C, Scotti O. 2004. Quantifying sensitivities of PSHA for France to earthquake catalog uncertainties, truncation of ground-motion variability, and magnitude limits[J]. Bulletin of the Seismological Society of America, 94(5): 15791594. [本文引用:1]
[22] Bender B. 1983. Maximum likelihood estimation of b-values for magnitude grouped data[J]. Bulletin of the Seismological Society of America, 73(3): 831851. [本文引用:7]
[23] Cao A, Gao S S. 2002. Temporal variation of seismic b-values beneath northeastern Japan island arc[J]. Geophysical Research Letters, 29(9): 48-1—48-3. [本文引用:1]
[24] Cornell C A. 1968. Engineering seismic risk analysis[J]. Bulletin of the Seismological Society of America, 58(5): 15831606. [本文引用:2]
[25] De Santis A, Cianchini G, Favali P, et al. 2011. The Gutenberg-Richter law and entropy of earthquakes: Two case studies in Central Italy[J]. Bulletin of the Seismological Society of America, 101(3): 13861395. [本文引用:1]
[26] Dong W M, Bao A B, Shah H C. 1984. Use of maximum entropy principle in earthquake recurrence relationships[J]. Bulletin of the Seismological Society of America, 74(2): 725737. [本文引用:2]
[27] Felzer K R. 2009. Simulated aftershock sequences for a M7. 8 earthquake on the southern San Andreas Fault[J]. Seismological Research Letters, 80(1): 2125. [本文引用:1]
[28] Gardner J K, Knopoff L. 1974. Is the sequence of earthquakes in Southern California, with aftershocks removed, Poissonian?[J]. Bulletin of the Seismological Society of America, 64(5): 13631367. [本文引用:1]
[29] Gerstenberger M, Wiemer S, Giardini D. 2001. A systematic test of the hypothesis that the b-value varies with depth in California[J]. Geophysical Research Letters, 28(1): 5760. [本文引用:1]
[30] Gibowicz S J. 1973. Variation of the frequency-magnitude relation during earthquake sequences in New Zealand [J]. Bulletin of the Seismological Society of America, 63(2): 517528. [本文引用:2]
[31] Gutenberg B, Richter C F. 1944. Frequency of earthquakes in California[J]. Bulletin of the Seismological Society of America, 34(4): 185188. [本文引用:1]
[32] Gutenberg B, Richter C F. 1949. Seismicity of the Earth and Associated Phenomena [M]. Princeton University Press: 273. [本文引用:1]
[33] Helmstetter A, Kagan Y Y, Jackson D D. 2006. Comparison of short-term and time-independent earthquake forecast models for southern California[J]. Bulletin of the Seismological Society of America, 96(1): 90106. [本文引用:1]
[34] Helmstetter A, Kagan Y Y, Jackson D D. 2007. High-resolution time-independent grid-based forecast for M≥5 earthquakes in California[J]. Seismological Research Letters, 78(1): 7886. [本文引用:3]
[35] Helmstetter A, Werner M J. 2014. Adaptive smoothing of seismicity in time, space, and magnitude for time-dependent earthquake forecasts for California[J]. Bulletin of the Seismological Society of America, 104(2): 809822. [本文引用:1]
[36] Kagan Y Y, Knopoff L. 1987. Statistical short-term earthquake prediction[J]. Science, 236(4808): 15631567. [本文引用:1]
[37] Karnik V, Klima K. 1993. Magnitude-frequency distribution in the European-Mediterranean earthquake regions[J]. Tectonophysics, 220(1-4): 309323. [本文引用:1]
[38] Kijko A, Smit A. 2012. Extension of the Aki-Utsu b-value estimator for incomplete catalogs[J]. Bulletin of the Seismological Society of America, 102(3): 12831287. [本文引用:7]
[39] Knopoff L. 2000. The magnitude distribution of declustered earthquakes in Southern California[J]. Proceedings of the National Academy of Sciences, 97(22): 1188011884. [本文引用:2]
[40] Kutliroff J R. 2017. Estimating the proportions of large to small earthquakes in seismic regions with a short span of earthquake magnitudes [D]. The University of Memphis, America. [本文引用:1]
[41] Main I. 2000. Apparent breaks in scaling in the earthquake cumulative frequency-magnitude distribution: Fact or artifact?[J]. Bulletin of the Seismological Society of America, 90(1): 8697. [本文引用:1]
[42] Marzocchi W, Sand ri L. 2003. A review and new insights on the estimation of the b-value and its uncertainty[J]. Annals of Geophysics, 46(6): 12711282. [本文引用:5]
[43] Miao Q, Langston C A. 2007. Empirical distance attenuation and the local-magnitude scale for the central United States[J]. Bulletin of the Seismological Society of America, 97(6): 21372151. [本文引用:1]
[44] Mignan A, Woessner J. 2012. Estimating the magnitude of completeness for earthquake catalogs [EB/OL], Community Online Resource for Statistical Seismicity Analysis. doi: DOI:10.5078/corssa-00180805 [本文引用:1]
[45] Mogi K. 1962 a. Study of elastic shocks caused by the fracture of heterogeneous materials and its relations to earthquake phenomena[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 40(1): 125173. [本文引用:1]
[46] Mogi K. 1962 b. Magnitude-frequency relation for elastic shocks accompanying fractures of various materials and some related problems in earthquakes(2nd paper)[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 40(4): 831853. [本文引用:1]
[47] Mogi K. 1963. The fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena(1st paper)[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 41: 595614. [本文引用:1]
[48] Ogata Y, Zhuang J. 2006. Space-time ETAS models and an improved extension[J]. Tectonophysics, 413(1-2): 1323. [本文引用:1]
[49] Okal E A, Kirby S H. 1995. Frequency-moment distribution of deep earthquakes; implications for the seismogenic zone at the bottom of slabs[J]. Physics of the Earth and Planetary Interiors, 92(3-4): 169187. [本文引用:1]
[50] Pacheco J F, Scholz C H, Sykes L R. 1992 a. Changes in frequency-size relationship from small to large earthquakes[J]. Nature, 355(6355): 7173. [本文引用:1]
[51] Pacheco J F, Sykes L R. 1992 b. Seismic moment catalog of large shallow earthquakes, 1900 to 1989[J]. Bulletin of the Seismological Society of America, 82(3): 13061349. [本文引用:1]
[52] Page R. 1968. Aftershocks and microaftershocks of the great Alaska earthquake of 1964[J]. Bulletin of the Seismological Society of America, 58(3): 11311168. [本文引用:6]
[53] Papazachos B, Delibasis N, Liapis N, et al. 1967. Aftershock sequences of some large earthquakes in the region of Greece[J]. Annals of Geophysics, 20(1): 143. [本文引用:1]
[54] Petersen M D, Frankel A D, Harmsen S C, et al. 2008. Documentation for the 2008 update of the United States national seismic hazard maps [R]. United States Geological Survey Open File Report, 2008-1128. [本文引用:2]
[55] Petersen M D, Moschetti M P, Powers P M, et al. 2015. The 2014 United States national seismic hazard model[J]. Earthquake Spectra, 31(S1): 130. [本文引用:3]
[56] Reasenberg P. 1985. Second-order moment of central California seismicity, 1969-1982[J]. Journal of Geophysical Research: Solid Earth, 90(B7): 54795495. [本文引用:1]
[57] Rhoades D A. 1996. Estimation of the Gutenberg-Richter relation allowing for individual earthquake magnitude uncertainties[J]. Tectonophysics, 258(1-4): 7183. [本文引用:1]
[58] Robinson R. 1979. Variation of energy release, rate of occurrence and b-value of earthquakes in the main seismic region, New Zealand [J]. Physics of the Earth & Planetary Interiors, 18(3): 209220. [本文引用:2]
[59] Sand ri L, Marzocchi W. 2007. A technical note on the bias in the estimation of the b-value and its uncertainty through the least squares technique[J]. Annals of Geophysics, 50(3): 329339. [本文引用:5]
[60] Scholz C H. 1968. The frequency-magnitude relation in microfracturing in rock and its relation to earthquakes[J]. Bulletin of the Seismological Society of America, 58: 399415. [本文引用:2]
[61] Scholz C H. 1997. Size distributions for large and small earthquakes[J]. Bulletin of the Seismological Society of America, 87(4): 10741077. [本文引用:2]
[62] Schorlemmer D, Wiemer S, Wyss W. 2005. Variations in earthquake-size distribution across different stress regimes[J]. Nature, 437(7058): 539542. [本文引用:1]
[63] Schwartz D P, Coppersmith K J. 1984. Fault behavior and characteristic earthquakes: Examples from the Wasatch and San Andreas fault zones[J]. Journal of Geophysical Research: Solid Earth, 89(B7): 56815698. [本文引用:1]
[64] Shi Y, Bolt B A. 1982. The stand ard error of the magnitude-frequency b-value[J]. Bulletin of the Seismological Society of America, 72(5): 16771687. [本文引用:1]
[65] Smith W D. 1981. The b-value as an earthquake precursor[J]. Nature, 289(5794): 136139. [本文引用:2]
[66] Suyehiro S. 1966. Difference between aftershocks and foreshocks in the relationship of magnitude to frequency of occurrence for the great Chilean earthquake of 1960[J]. Bulletin of the Seismological Society of America, (1): 185200. [本文引用:2]
[67] Suyehiro S, Asada T, Ohtake M. 1964. Foreshocks and aftershocks accompanying a perceptible earthquake in central Japan: On the peculiar nature of foreshocks[J]. Papers in Meteorology & Geophysics, 15(1): 7188. [本文引用:2]
[68] Suzuki Z. 1958. A statistical study on the occurrence of small earthquakes Ⅲ[J]. Science Reports of the Tohoku University Ser 5, Geophysics, 10(1): 1527. [本文引用:3]
[69] Talebi M, Zare M, Peresan A, et al. 2017. Long-term probabilistic forecast for M≥5. 0 earthquakes in Iran[J]. Pure and Applied Geophysics, 174(4): 15611580. [本文引用:1]
[70] Tinti S, Mulargia F. 1985. Effects of magnitude uncertainties on estimating the parameters in the Gutenberg-Richter frequency-magnitude law[J]. Bulletin of the Seismological Society of America, 75(6): 16811697. [本文引用:1]
[71] Tinti S, Mulargia F. 1987. Confidence intervals of b-values for grouped magnitudes[J]. Bulletin of the Seismological Society of America, 77(6): 21252134. [本文引用:6]
[72] Triep E G, Sykes L R. 1997. Frequency of occurrence of moderate to great earthquakes in intracontinental regions: Implications for changes in stress, earthquake prediction, and hazards assessments[J]. Journal of Geophysical Research: Solid Earth, 102(B5): 99239948. [本文引用:2]
[73] Utsu T. 1965. A method for determining the value of b in formula log N= a- bM showing the magnitude-frequency relation for earthquakes[J]. Geophysical bulletin of Hokkaido University, (13): 99103. [本文引用:7]
[74] Utsu T. 1966. A statistical significance test of the difference in b-value between two earthquake groups[J]. Journal of Physics of the Earth, 14(2): 3740. [本文引用:7]
[75] Utsu T. 1999. Representation and analysis of the earthquake size distribution: A historical review and some new approaches[J]. Pure and Applied Geophysics, 155(2-4): 509535. [本文引用:1]
[76] Wang Q, Jackson D D, Kagan Y Y. 2011. California earthquake forecasts based on smoothed seismicity: Model choices[J]. Bulletin of the Seismological Society of America, 101(3): 14221430. [本文引用:1]
[77] Weichert D H. 1980. Estimation of the earthquake recurrence parameters for unequal observation periods for different magnitudes[J]. Bulletin of the Seismological Society of America, 70(4): 13371346. [本文引用:10]
[78] Wiemer S. 2001. A software package to analyze seismicity: ZMAP[J]. Seismological Research Letters, 72(3): 373382. [本文引用:1]
[79] Wiemer S, Wyss M. 2000. Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the western United States, and Japan[J]. Bulletin of the Seismological Society of America, 90(4): 859869. [本文引用:1]
[80] Woessner J, Wiemer S. 2005. Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty[J]. Bulletin of the Seismological Society of America, 95(2): 684698. [本文引用:3]
[81] Working Group on California Earthquake Probabilities. 2003. Earthquake probabilities in the San Francisco Bay region: 2002—2031 [R]. United States Geological Survey Open File Report, 03-214. [本文引用:1]
[82] Wyss M, Wiemer S. 2000. Change in the probability for earthquakes in southern California due to the Land ers magnitude 7. 3 earthquake[J]. Science, 290(5495): 13341338. [本文引用:1]
[83] Youngs R R, Coppersmith K J. 1985. Implications of fault slip rates and earthquake recurrence models to probabilistic seismic hazard estimates[J]. Bulletin of the Seismological Society of America, 75(4): 939964. [本文引用:1]