〔作者简介〕 吴果, 男, 1988年生,2018年于中国地震局地质研究所获固体地球物理学博士学位, 主要研究方向为地震活动性与地震危险性分析, 电话: 18810404034, E-mail: wgfirst@foxmail.com。
b值在地震活动性研究和地震危险性分析中起着十分重要的作用, 通常拟合 b值的方法有最小二乘法(Least Squares Method)和极大似然法(Maximum Likelihood Estimation)。最小二乘法简单易行, 得到了广泛的应用。然而很多研究表明该方法存在一定的局限性, 极大似然法在特定条件下可以作为最小二乘法的一种可行的替代或补充方法。前人对极大似然法的研究非常繁杂, 提出了各种各样的方程式, 每个方程式的隐含假设和求解方式各不相同。文中对主要方程式进行了简要的回顾, 并按照是否考虑震级的归档效应、 是否设定有限最大震级、 是否对不同震级档数据取不同的观察时段和是否具有解析解这4个方面, 对这些方程式进行了分类和总结。进而对震级的归档效应、 震级的测量误差、 样本量、 震级跨度、 最小完整震级和前余震共6个可能影响极大似然法估计 b值的因素进行了分析和总结。最后对正确使用这些方程式提出了合理的建议。文中的分析和总结有助于更准确地理解和使用不同的极大似然法估计 b值的方程式, 以供相关研究者参考。
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 of
震级-频度(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值的方程式提供参考。
G-R关系, 即震级-频度关系最早由Gutenberg等(1944, 1949)提出。某地区在某段时间内发生震级≥ M的地震数N(M)可以用方程式(1)表示:
参数
假设某地区在某段时间内发生的地震的震级
进一步构建震级M的概率分布函数:
那么, 如果假设震级M为1个连续随机变量, 对
当Mmax ≫Mmin时, 一般要求
相应地, 方程式(6)可以写为
当震级范围较小时, 这一简化将不再适用(Knopoff, 2000)。
自方程式(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值显得尤为重要。
最小二乘法(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所在震级档(
(1)大、 小地震权重不一致: 累计和非累计震级-频度关系, 都需要划分震级档来获取每个震级档的数据点(Mi,
(2)地震样本量的限制: 如上所述, 最小二乘法需要划分震级档来获取每个震级档的数据点(Mi,
Sandri等(2007)从理论分析和数值拟合2个方面论述了最小二乘法拟合b值带来的影响, 还将最小二乘法应用的对象分为累计和非累计震级-频度关系2种。蒙特卡洛模拟结果显示: 将最小二乘法应用于G-R关系的累计形式时, 地震数的累计过程会产生一种过滤效应, 这会导致对b值的不确定性被低估。如图 1所示, 由于需要对地震数取对数形式, 无论是累计还是非累计震级-频度关系, 最小二乘法拟合的b值都偏小, 其中前者的偏差相对较小, 随着地震数的增加偏差减小。与对地震数取对数形式带来的误差相比, 震级自身的测量误差可以忽略不计。
前人对最小二乘法和极大似然法拟合b值的效果进行过对比分析, Weichert(1980)通过对加拿大西部3个地震带的数据进行分析发现, 对于受到较好约束的数据, 与极大似然法相比, 最小二乘法能得到相近的结果。但是, 对于质量较差的数据, 由于受主观绘图和权重取值的影响, 最小二乘法得到的结果则与极大似然法存在较为显著的差异。总体上, 极大似然法的效果较好, 同时最小二乘法不能反映人们预先对最大震级的判断, 因此建议使用极大似然法估计b值。
张建中等(1981)用蒙特卡洛模拟的方法比较了极大似然法、 线性最小二乘法和非线性最小二乘法共3种算法的均方误差、 运算量和相关性。结果显示, 极大似然法估计b值的均方误差最小, 其次为线性最小二乘法, 均方误差最大的是非线性最小二乘法; 三者运算量的大小与均方误差的规律相同; 在与极大似然法结果的相关性上, 线性最小二乘法同样优于非线性最小二乘法; 随着样本量N的不断增大, 3种方法结果的差异逐渐减小。
上述总结和分析重点介绍了利用最小二乘法求取b值的一些局限性。事实上, 最小二乘法和极大似然法互有优劣, b值的最小二乘法估计的最大优点是简单易行, 且结果与极大似然法差异不大。当观察点(xi, yi)之间相互独立, 随机误差
基于极大似然法估计b值并分析其不确定性一直都是地震学者研究的热点, 历史上各研究者提出了众多的方程式。这些方程式基于的假设条件有的明显, 有的不明显, 适用范围也各不相同, 极易混淆。本文先以时间顺序重点介绍前人的代表性工作, 然后再根据不同方程式的隐含假设和求解方式, 统一进行归类和总结分析。
Aki(1965)最早提出用极大似然法估计b值, 并给出了方程式。用极大似然法推导b值的过程如下: 假设震级M是1个连续随机变量, 当Mmax ≫Mmin时, 其概率密度函数可以用方程式(7)来描述。则
震级样本(
构建似然函数, 令:
令似然函数对
可以得到
式中,
由于只有记录完整的地震数才符合G-R关系, 所以
式中, Mc为最小完整震级。 b值的极大似然估计是它的渐进无偏估计(Page, 1968; 张建中等, 1981)。Aki(1965)利用中心极限定理推得
Utsu(1965)利用总体(Population)和样本(Sample)的一阶矩(First Moment)相等的方法得到了相同的方程式。同样也假设震级M是连续随机变量, 当Mmax ≫Mmin时, 由方程式(7)可以得到震级M总体的一阶矩:
而震级样本(
Aki(1965)和Utsu(1965)得到的估计值基于的假设为震级是连续的随机变量, 同时没有设定震级的最大值。然而, 实际中的震级并非连续的随机变量。由于地震震级的精度有限, 对于仪器记录的地震, 通常为0.1个震级单位, 对于历史地震记录, 则≥ 0.5个震级单位。地震震级通常以 “ 0.1” 或者 “ 0.5” 为间隔跳跃式显示, 如 “ M6.0” 、 “ M6.1” 、 “ M6.2” , 相邻的震级被归并到与之最近的 “ 震级bin” (Magnitude bin)中。例如, 震级
因此, 假设震级精度为
虽然方程式(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)。
Aki(1965)和Utsu(1965, 1966)的方程式有一个共同的问题, 即没有设定震级上限
式中,
式中,
人类通过仪器记载的地震目录不过只涵盖数十年至一百年, 这与大地震几百年乃至上千年的复发周期相比显得十分短暂。同时, 一个地区的地震监测能力往往随着时间的推移逐步提高, 所以震级越小的地震具有相对完整记录的时间越短。因此, 有限的地震目录资料显得尤为宝贵。为了充分利用地震目录资料, Weichert(1980)将大小不同的震级档对应不同的完整时段, 对地震目录进行分档并予以利用。
如图 2所示, 假设震级M满足
式中, N为各个震级档在相应完整时段内记录到的地震数的和,
式中,
方程式(21)基于3个基本假设: 1)将地震目录按震级档予以划分, 震级档的跨度可以各不相等, 也可以等间隔(
Bender(1983)同样假设震级为归档后的随机变量, 将震级范围
式中,
Bender(1983)的方程式基于的假设与上述Weichert(1980)的方程式(21)基本相同, 区别在于: 1)方程式(23)以等间距
Tinti等(1987)假设震级为等间隔
相应的b值的标准差为
式中,
Kijko等(2012)认为虽然Weichert(1980)的方程式可以使用不同震级档对应不同完整记录时段的数据, 但是方程式(21)没有解析解, 不便于使用。于是, 他们在方程式(15)(Aki, 1965; Utsu, 1965)的基础上推导了新的方程式, 使用的假设条件也与方程式(15)相同, 即假设震级为连续随机变量, 同时没有设定震级上限。
Kijko等(2012)按照地震监测水平的差异, 将总记录时长为T的地震目录分成L个时段, 每个时段的次级地震目录对应时长为
式中,
对比图2和图3可以发现, 方程式(26)和方程式(21)不仅基于的基本假设不同, 而且在处理不同大小的震级对应不同完整时段时的策略也不同。方程式(21)是将地震目录分成多个震级档, 每个震级档可分别取不同时长的地震记录数据。而方程式(26)则根据地震完整记录水平随着时间的变化, 将地震目录分成多个次级地震目录。Kijko等(2012)开发了基于Matlab平台的程序包AUE, 可以方便地使用方程式(26)估计
前文在介绍各种基于极大似然法估计b值的方程式时, 已经给出了与之相应的求解b值标准差的方程式, 其中应用最广泛的是方程式(16)(Aki, 1965)。这里再补充介绍2个求解b值标准差的重要方程式。
Utsu(1966)最早在其研究中构建χ 2分布来研究b值, 张建中等(1981)进一步通过构建
并且给出了b值在置信水平为
式中,
另1个重要的方程式是Shi等(1982)提出的:
ZMAP软件在计算b值的标准差时即使用方程式(30)(Woessner et al., 2005)。
由于绝大部分文献都没有直接提及G-R关系中的a值及其标准差的计算, 故在此作简单介绍。在求得b值的极大似然估计值
式中, N为样本量,
本文一共总结了7个最常见的基于极大似然法估计b值的方程式, 这些方程式的假设和应用条件各不相同, 现对其予以分类和总结。
这些方程式大致有3个隐含假设, 每个假设有相互对立的2种选择: 1)假设震级M为归档后的随机变量(Binned Random Variable, 简称BRV, 详细论述见3.2节), 否则为通常的连续随机变量(Continuous Random Variable); 2)设定有限最大震级(Mmax), 否则Mmax趋向于无穷大, 或者要求Mmax ≫Mmin; 3)不同的震级档对应不同的观察时段(记为Ti), 否则每个震级档的观察时段相等, 即Ti=T。
此外, 方程式的求解方式也分为存在解析解(Analytic solution, AS)和只能求数值解(Numerical solution, NS)。7组方程式的假设条件和求解方法具体见表1, 若相应参考文献中明确给出了b值或者
在这些方程式中, 应用最广的是假设震级为连续随机变量的方程式(15)。那些假设震级为归档后的随机变量的方程式可能由于求解困难反而应用很少(Marzocchi et al., 2003), 如Weichert(1980)和Bender(1983)提出的方程式都只有数值解。
求解b值标准差的方程式有时会和估计b值的方程式搭配使用, 例如在ZMAP软件中使用方程式(18)(Utsu, 1966)估计b值, 而相应标准差则用Shi等(1982)提出的方程式(30)来估算。
b值在地震活动性预测和地震危险性分析中有着极其重要的作用, 然而其自身会受到多种因素的影响, 例如样本量、 样本的震级跨度以及余震等。本文总结前人的工作, 现将各种可能影响极大似然法估计b值的因素一一予以介绍。
如表1所示, 很多用极大似然法估计b值的方程式都假设震级是连续随机变量, 然而这与实际情况是不相符的。如前文所述(见3.2节), 现实中的震级存在归档效应(Binned magnitude), 相邻的震级被归并到离其最近的 “ 震级bin” 。例如, 震级
Marzocchi等(2003)指出, 如果简单地假设震级为连续随机变量会带来2类偏差: 第一类, 为幂律分布的连续随机变量的震级均值
Bender(1983)对第一类偏差进行了详细的论述。其指出归并到 “ 震级bin” 后的震级均值要系统地大于连续随机变量的震级均值。因为地震数在1个 “ 震级bin” 内并不是均匀分布的,
对于第二类偏差, Utsu(1966)最早提出可以对方程式(15)进行1个小的修改。由于最小的 “ 震级bin” , 即最小完整震级
Bender(1983)指出对于同样的样本, 采用不同的方程式估计b值, 结果将显示出明显的区别。当震级间隔
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)的结果略优。
震级的测量误差不同于上述震级的归档效应。震级的归档是指震级自身作为1个连续随机变量被强行归并到各 “ 震级bin” 中, 类似于 “ 四舍五入” 的过程。而震级的测量误差是指地震台站的在测量震级时的所包含的随机误差
式中,
Marzocchi(2003)通过理论和数值模拟2种方式研究了震级的归档效应和震级的测量误差对于极大似然法估计b值的影响, 结果显示与对震级的归档效应不当处理相比, 例如方程式(15)(Aki, 1965)直接假设震级为连续随机变量, 震级的测量误差带来的影响可以忽略不计。这可能也是大部分作者在拟合估计b值的方程式时, 没有特别关注震级测量误差的原因(Weichert, 1980; Kijko et al., 2012)。Sandri等(2007)在研究地震数取对数形式和震级的测量误差对于最小二乘法估计b值的影响时, 也有类似的结论, 与前者相比, 震级的测量误差带来的影响可以忽略不计。
影响极大似然法估计b值的另外一个重要因素是样本量, 需要指出的是这里的样本量是指经过筛选的最小完整震级(Mc)以上的地震数, 而不是直接下载得到的地震目录中的地震数。
样本量早期就被用于估计b值的标准差。Aki(1965)利用中心极限定理推导得出
Shi等(1982)用方程式(15)(Aki, 1965; Utsu, 1965)估计美国加州中部的b值, 结果显示, 对于该地区, 要得到相对稳定的结果, 筛选样本所需的时空窗中大约需要包含100个地震。Bender(1983)指出当样本量过小时, 无论使用何种方程式, 结果都将变得不可信。当样本量N=100时, b值的标准差
Kutliroff(2017)用蒙特卡洛模拟方法对方程式(18)(Utsu, 1966)和方程式(24)(Tinti et al., 1987)进行检验, 发现二者的结果十分相近, 当样本量
早期的研究者没有针对震级跨度进行专门的研究, 然而一些方程式的推导过程事实上包含对震级跨度的隐含假设。如1.2节所述, Mmax ≫Mmin, 震级的概率密度函数可以由方程式(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)的结果则系统地偏高。
本文在描述样本量N时, 特别指出是经过筛选的最小完整震级(Mc)以上的地震数。这是因为如果选择了错误的最小震级, 将会对b值的极大似然估计值产生显著的影响。如图6a所示, 在G-R关系的小震级端, 由于地震记录不全, 将导致该震级段曲线斜率变缓, 中文文献中习惯将这一现象称为 “ 掉头” (黄玮琼等, 1989; 任雪梅等, 2011)。如果错误地选择了过小的最小完整震级, 会导致算得的b值偏小。在Woessner等(2005)的示例中(图6b), 当真实的Mc被人为地从M1.4降低到M0.7时, 相应的b值从真实值1.02降低到0.8。
事实上, 根据b值是否能稳定取值还可以用来估计最小完整震级。Cao等(2002)提出了1种估计最小完整震级的方法: 如图6b所示, 将人为选择的小震级端的截断震级(Cutoff magnitude)记为Mco, 真实的最小完整震级记为Mc。当Mco< Mc时, 随着Mco 增大, b值不断增大; 当Mco ≥ Mc时, 随着Mco 增大, b值保持相对稳定; 当Mco ≫ Mc时, 随着Mco 增大, b值再次开始增大。如果使用的Mco ≪ Mc, 则求得的b值会明显偏小。黄亦磊等(2016)对比了多种分析地震目录完整性的方法, 结果显示Cao等(2002)的方法与多数其他方法相比效果较好。
由上述分析可知, 准确选择最小完整震级对于估计b值来说十分重要。而估计最小完整震级的方法多种多样(Woessner et al., 2005; Mignan et al., 2012; 黄亦磊等, 2016), 只利用单一方法, 例如G-R关系的线性程度进行估计, 往往存在缺陷, 准确估计Mc需要结合使用多种方法(任雪梅等, 2011; 吴果等, 2014)。
在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)等提出的方法。
本文回顾了震级-频度(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] |
|
[2] |
|
[3] |
|
[4] |
|
[5] |
|
[6] |
|
[7] |
|
[8] |
|
[9] |
|
[10] |
|
[11] |
|
[12] |
|
[13] |
|
[14] |
|
[15] |
|
[16] |
|
[17] |
|
[18] |
|
[19] |
|
[20] |
|
[21] |
|
[22] |
|
[23] |
|
[24] |
|
[25] |
|
[26] |
|
[27] |
|
[28] |
|
[29] |
|
[30] |
|
[31] |
|
[32] |
|
[33] |
|
[34] |
|
[35] |
|
[36] |
|
[37] |
|
[38] |
|
[39] |
|
[40] |
|
[41] |
|
[42] |
|
[43] |
|
[44] |
|
[45] |
|
[46] |
|
[47] |
|
[48] |
|
[49] |
|
[50] |
|
[51] |
|
[52] |
|
[53] |
|
[54] |
|
[55] |
|
[56] |
|
[57] |
|
[58] |
|
[59] |
|
[60] |
|
[61] |
|
[62] |
|
[63] |
|
[64] |
|
[65] |
|
[66] |
|
[67] |
|
[68] |
|
[69] |
|
[70] |
|
[71] |
|
[72] |
|
[73] |
|
[74] |
|
[75] |
|
[76] |
|
[77] |
|
[78] |
|
[79] |
|
[80] |
|
[81] |
|
[82] |
|
[83] |
|