古登堡-里希特定律中的b值统计样本量研究1

李世杰 吕悦军 刘静伟



古登堡-里希特定律中的值统计样本量研究1

李世杰 吕悦军 刘静伟

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

值是研究地震活动的重要指标，其广泛应用于地震危险性分析和地震预测研究之中，与实际资料的完整性、样本量的大小、计算方法等因素有着重要的关系。常见的值计算方法有最小二乘法和最大似然法，样本量的大小对这2种方法影响很大。本文利用蒙特卡罗模拟地震目录和汾渭地震带实际目录作为样本，从中抽取不同大小的样本量进行计算，研究不同样本量下这2种方法计算得到的值与设定值或真实值之间的差别。结果表明，最小二乘法需要的最低样本量为1000，最大似然法为200；当样本量达不到要求时，计算出的值是不可靠的；由于对样本量的要求不同，前者适用于计算区域的整体值，而后者在研究某区域值在时间轴上的变化方面更有优势。本研究为确定2种值计算方法对样本量的最低要求提供了参考依据。

值 最小二乘法 最大似然法 样本量大小

引言

古登堡和里希特在1941年美国地质学会会刊中提出全球的地震活动服从经验关系lg=-（式中表示震级，表示震级≥的地震次数，、是常数）。一般认为值代表地区的地震活动总体水平（Rundle，1989）；值代表地震活动的大小地震数量的比例，是地震活动研究中的重要参量。研究表明，值具有明确的物理意义，与地壳的介质特性、应力状态和不均匀性有关，能反映所研究区域的地质构造特征（王熠熙等，2015；谢卓娟等，2015）以及地震的震源特征（Schorlemmer等，2005；Gulia等，2010；刘静伟等，2016；张广伟，2016）。因此，值广泛应用于地震危险性分析和地震预测研究之中。在地震危险性分析中，值和地震年发生率共同用于确定地震活动的水平（胡聿贤，1999），其取值对地震危险性分析结果的影响较大（鄢家全等，1996；黄玮琼等，1998；谢卓娟等，2013）；而在地震预测研究中，值作为基本的地震活动性参数，成为地震预测的常用指标参数（韩渭宾，2003；沈建文等，2007）。

通常根据实际地震资料统计得到值。目前，常用的值估算方法是最大似然估计法（Aki，1965）和最小二乘法等，这些方法不仅对实际资料的完整性和精度有一定要求，同时也需要足够的地震资料。在实际工作中，对于历史地震资料短缺或地震活动水平低的地区，在计算时常常将现代小震资料与历史地震资料联合使用，以弥补地震资料样本量的不足（黄玮琼等，1989；鄢家全等，1996；胡聿贤，1999；潘华等，2006）。统计值对地震资料样本量需求的定量研究，国内尚无专门的研究报道。在国外，Nava等（2017）利用蒙特卡罗模拟地震目录进行抽样估计，研究最大似然法计算值时对样本量的需求，得出计算值时样本量和精度之间的相互关联关系，但没有应用实际地震目录进行分析。国内仅有少部分学者的研究涉及相关内容，如韩晓明等（2016）研究河套地震带值时空变化特征的文章中，探讨了地震前后值的变化规律，分别使用最小二乘法和最大似然法对值时间和空间进行了扫描计算，在最小二乘法计算中，设定每次计算的窗长内包含的样本数目不少于100，而最大似然法的扫描窗内包含的样本数目不少于20；刘方斌等（2017）在鲁西南聊考断裂带地震危险性评价与活动性分布的研究中，利用最大似然法和最小二乘法估算值并进行了对比研究，但是没有给出样本量的具体数目，且没有用分震级段的方法进行值估算。

本文采用Utsu（1965）提出的最大似然法和最小二乘法，利用模拟地震目录和实际地震目录，定量分析最小二乘法和最大似然法计算值时分别对地震资料样本量的需求。

1 资料概况

1.1 模拟地震目录

蒙特卡罗法是以抽样和随机数的产生为基础的随机性方法，也称为随机抽样法、计算机随机模拟法等。蒙特卡罗方法的基本原理是通过数字模拟试验，得到所要求解的出现某种事件的概率作为问题的近似解。其基本思想是：为了求解数学、物理、工程技术以及管理等方面的问题，首先建立概率模型或随机过程，使用相应的参数，得到某些问题（如概率分布或数学期望等问题）的解；然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，并用算术平均值作为所求解的近似值。因此，只要尽可能完整地表示地震震源模型，即获得某一区域内地震事件的发生时间和空间上的分布规律，那么就可以直接利用蒙特卡罗方法产生合成地震序列。

蒙特卡罗方法的基本依据是随机模拟次数足够多的情况下，事件发生的频率可以反映事件发生的概率。当模拟次数有限时，计算结果和真实值之间必然存在误差，这种误差随着模拟次数的增加而减少。

蒙特卡罗模拟中随机变量的简单子样（1，2，3，……，X）是独立分布的，即每次模拟中事件发生的次数与其它任意1次中事件发生的次数无关，那么随机变量就是服从泊松分布的，当随机变量的期望值趋向于无穷时，泊松分布趋近于正态分布。

本文利用蒙特卡罗方法模拟地震目录（张建中，1974a，1974b；任雪梅等，2011），研究不同样本量下计算方法对值的影响。

首先，构造概率分布模型：

其中，（）为震级的概率密度函数，表示为：

其次，将公式（2）带入公式（1），得到关于模拟震级的函数：

设定=1、0=2.0、u=8.0，利用MATLAB随机数程序模拟生成不同样本量的地震目录，样本数分别为5000、4500、4000、3500、3000、2500、2000、1750、1500、1250、1000、750、500、250、100、80、60、40、20、10，每个目录模拟10000组。由于实际大地震非常稀缺，所以震级上限的限制对计算结果没有显著影响。

1.2 汾渭地震带实际地震目录

本研究使用的实际地震资料来自中国地震台网中心的地震数据库。研究表明地震资料的完整性对地震活动分析非常重要（Wiemer等，2000；Rotondi等，2002）。经统计分析（焦远碧等，1990；黄伟琼等，1994a，1994b；周公威等，2007；谢卓娟等，2012a；徐伟进等，2014），1970年以来汾渭地震带2.0以上的地震资料基本完整，故本研究选取1970年1月—2010年12月，≥2.0的8184次地震资料进行研究分析，其中2.0—2.9地震7044次、3.0—3.9地震964次、4.0—4.9地震148次、5.0—5.9地震25次、6.0—6.9地震3次，所用的地震目录均已删除前震、余震。汾渭地震带4.7级以上地震的震中分布见图1。

图1 汾渭地震带范围及震中分布

2 模拟地震目录的结果分析

利用最小二乘法和最大似然法对模拟的地震目录进行定量研究分析，比较2种方法在不同的样本量下计算得到的平均值以及样本量的大小对计算结果的影响。

文中使用的最大似然法为Utsu（1965）提出的公式：

采用上述最大似然法和最小二乘法计算，分别得到平均值随样本量的变化，如图2、3所示。

图2 不同样本量下最大似然法计算得到的平均b值拟合图

图3 不同样本量下最小二乘法计算得到的平均b值拟合图

由图2、3可知，当样本量大于300时，最大似然法计算的平均值能取得符合预期的数值；样本量大于1000时，最小二乘法计算的平均值趋于稳定，向理论值收敛。

表1 最大似然法模拟结果

表2 最小二乘法模拟结果

当设定=1时，不同样本下得到的值分布直方图见图4。由图可知，当样本量≥200时，最大似然法计算的值很好地向设定值收敛；样本量≥500时，最小二乘法计算的值有较好的收敛效果。从估值概率来看，取得的值一般都偏低。

根据上述分析，2种方法虽然都能得到符合预期的数值，但最小二乘法计算的平均值通常比最大似然法的小，而标准差比最大似然法的大；样本量小于1000时，利用最小二乘法得到的平均值精确度较低；样本量大于1000时，2种计算方法得出的结果差别不大。进一步分析认为，虽然2种方法都需要一定的样本量，且样本量越大、得到的结果越准确，但最大似然法对样本量的要求要比最小二乘法低，样本量大于200时，计算得到的平均值与设定值一致性较好；样本量大于1000时，基本等于设定值。

3 实际地震目录的计算结果分析

根据汾渭地震带地震资料的完整性和可靠性研究，自1970年以来，台站记录得到的≥2.0地震目录基本完整，故本研究以汾渭地震带1970年1月—2010年12月的地震记录为计算样本，研究实际地震目录样本量对值计算的影响。表3给出该地震带不同震级档的地震数目，并计算出不同震级档的年平均发生率（黄玮琼等，1989；潘华等，2006；吴兆营等2005；谢卓娟等，2012b）。

图4 不同样本量的b值分布直方图

表3 汾渭地震带地震分档统计和年平均发生率

续表

震级档M地震个数年平均发生率 5.0—5.4150.36 5.5—5.9100.23 6.0—8.530.07

注：年平均发生率指震级≥的年均地震数，代表地震活动水平。

根据表3中统计的各震级档的地震个数和年平均发生率，利用最小二乘法和最大似然法计算得到的值分别为0.75和0.75335。

为研究实际地震目录下值计算对样本量的要求，采用任雪梅（2011）的抽样原则，对汾渭地震带的地震目录进行均匀抽样分析，每次抽样完将样本放回进行下一次抽样，得到样本量分别为10、50、100、200、300、500、700、1000、5000的地震目录，每个目录重复抽样10000次。同样，采用2种方法分别计算不同样本量下的平均值（表4），并分析不同样本量下的值变化。

表4 汾渭地震带b值拟合情况（1500—2010年）

由表4可以看出，2种计算方法得到值随着样本量的增加逐渐接近真实值；当样本量大于200时，最大似然法能够得出相对稳定可靠的值；当样本量大于500时，最小二乘法才能得到相对稳定可靠的值。因此，对实际的地震目录，最小二乘法对样本量的要求也比较高，且受地震目录完整性的影响较大。

4 结论

通过上述研究和对比分析，得到以下认识：

（1）利用最大似然法计算值时，样本量至少要在200以上；对于最小二乘法，样本量要求不少于1000；一般来说，计算的值都低于设定值。

（2）最小二乘法利用震级-频度的线性关系来拟合计算值，样本量不足会影响其线性关系。因此，最小二乘法受样本量影响较大，样本量小于1000时，值的计算值与设定值相差较大，数值也不稳定。

（3）最大似然法方便快捷，受样本量影响小，计算出的值相对稳定，但误差估计值偏大，其利用平均震级计算，只与地震个数有关，受数据质量的影响较小；最小二乘法受样本量影响较大，在样本量充足的情况下计算的值比较准确，而在样本量不足时值波动较大，误差也随着样本量的减少而增加。在实际应用中，以地震目录充足为前提，研究不同区域的值优先选用最小二乘法；而研究某区域值时间上的变化时，采用最大似然法估算的值相对稳定，更能体现值长期变化的趋势。

（4）从误差和计算量来看，最大似然法比最小二乘法要小，但随着样本量的增加，2种方法计算结果的差异越来越小。但由于半对数坐标下不同震级档数据权重的不对等性，国外已极少使用最小二乘法。

韩渭宾，2003．值在地震预测中的三类应用及其物理基础与须注意的问题．四川地震，（1）：1—5．

韩晓明，张文韬，王树波等，2016．河套地震带的值时空变化特征分析．中国地震，32（3）：522—532．

胡聿贤，1999．地震安全性评价技术教程．北京：地震出版社．

黄玮琼，时振梁，曹学锋，1989．值统计中的影响因素及危险性分析中值的选取．地震学报，11（4）：351—361．

黄玮琼，李文香，曹学锋，1994a．中国大陆地震资料完整性研究之一——以华北地区为例．地震学报，16（3）：273—280．

黄玮琼，李文香，曹学锋，1994b．中国大陆地震资料完整性研究之二——分区地震资料基本完整的起始年分布图象．地震学报，16（4）：423—432．

黄玮琼，李文香，1998．地震区划中值统计时空范围的确定．地震学报，20（5）：2—6．

焦远碧，吴开统，杨满栋，1990．我国地震台网监测能力及台网观测条件质量评定．中国地震，6（4）：3—9．

刘方斌，曲均浩，田兆阳等，2017．鲁西南聊考断裂带地震危险性评价与活动性分布．大地测量与地球动力学，37（8）：802—807．

刘静伟，吕悦军，2016．川滇地区值空间分布特征及其与震源类型关系的初步探讨．震灾防御技术，11（3）：561—572．

潘华，李金臣，2006．地震统计区地震活动性参数值及4不确定性研究．震灾防御技术，1（3）：218—224．

任雪梅，2011．地震区划中值统计的若干问题研究．北京：中国地震局地球物理研究所．

任雪梅，高孟潭，冯静，2011．地震目录的完整性对值计算的影响．震灾防御技术，6（3）：257—268．

沈建文，余湛，邱瑛，2007．地震安评中地震活动性的统计区域与值．国际地震动态，（3）：1—6．

王熠熙，张辉，刘双庆等，2015．河北平原地震带值时空变化特征．地震工程学报，37（01）：188—195.

吴兆营，薄景山，刘志平等，2005．东北地震区值和地震年平均发生率的统计分析．东北地震研究，21（3）：27—32，72．

谢卓娟，吕悦军，张力方，等，2012a．基于现代地震资料确定汾渭地震带分区及其地震活动性参数．地球物理学进展，27（3）：894—902．

谢卓娟，吕悦军，彭艳菊等，2012b．东北地震区小震资料完整性分析及其对地震活动性参数的影响研究．中国地震，28（3）：256—265．

谢卓娟，吕悦军，兰景岩等，2013．值和4的统计分析及其不确定性对地震危险性分析结果的影响研究．地震研究，36（1）：86—92．

谢卓娟，李山有，吕悦军，2015．滇西南地区主要活动断裂的值空间分布特征．地球科学（中国地质大学学报），40（10）：1755—1766.

徐伟进，高孟潭，2014．中国大陆及周缘地震目录完整性统计分析．地球物理学报，57（9）：2802—2812.

鄢家全，韩炜，高孟潭，1996．地震活动性参数的不确定性及其对区划结果的影响．中国地震，12（S1）：71—77．

张广伟，2016．云南地区地震的重新定位及值研究．中国地震，32（1）：54—62．

张建中，1974a．蒙特卡洛方法(Ⅰ)．数学的实践与认识，（01）：28—40.

张建中，1974b．蒙特卡洛方法(Ⅱ)．数学的实践与认识，（02）：43—56.

周公威，张伯明，吴忠良等，2007．中国数字地震台网（CDSN）的近期发展．地球物理学进展，22（4）：1130—1134．

Maximum likelihood estimate ofin the formula log=-and its confidence limits. Bulletin of the Earthquake Research Institute, 43(2): 237—239.

Gulia L., Wiemer S., 2010. The influence of tectonic regimes on the earthquake size distribution: A case study for Italy. Geophysical Research Letters, 37(10): L10305.

Nava F. A., Márquez-Ramírez V. H., Zúñiga F. R., et al., 2017. Gutenberg-Richter-value maximum likelihood estimation and sample size. Journal of Seismology, 21(1): 127—135.

Rundle J. B., 1989. Derivation of the complete Gutenberg-Richter magnitude-frequency relation using the principle of scale invariance. Journal of Geophysical Research: Solid Earth, 94(B9): 12337—12342.

Rotondi R., Garavaglia E., 2002. Statistical Analysis of the Completeness of a Seismic Catalogue. Natural Hazards, 25(3):245—258.

Schorlemmer D., Wiemer S., Wyss M., 2005. Variations in earthquake-size distribution across different stress regimes. Nature, 437(7058): 538—542.

Utsu T., 1965. A method for determining the value ofin a formula log=-showing the magnitude-frequency relation for earthquakes. Geophysical Bulletin of Hokkaido University, (13): 99—103.

Wiemer S., Wyss M., 2000. Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, the Western United States, and Japan. Bulletin of the Seismological Society of America, 90(4): 859—869.

The Study of Sample Size on-value Statistics in the Gutenberg-Richter's Law

Li Shijie, Lü Yuejun and Liu Jingwei

(Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China)

The-value is an important indicator for evaluating the level of seismicity, and is widely used in seismic hazard analysis and earthquake prediction research. The-value is basically effected by some factors, such as the actual data integrity, the quantity of earthquake samples, calculation methods and so on. The most widely used two methods for calculatingvalues include the least square method (LSM) and the maximum likelihood method (MLM), and the results from both methods can remarkably vary with different the quantity of earthquake samples. In this paper, based on the real catalog of the Fen-Wei seismic zone and the Monte-Carlo simulated earthquake catalog, we use different number of samples to calculate thevalues and try to find the differences between the calculated results and the set values/true values. It turns out that the threshold value of samples is 1000 and 200 for LSM and MLM, respectively. When the sample size does not meet the requirements, the calculated-value is not reliable. It suggets that LSM is suitable for calculating the-value of the whole region, while MLM is more advantageous in studying the-value of a region on the time axis. This study provides a reference for the threshold values of the samples for two methods, which is of great significance to the seismic hazard analysis and earthquake prediction.

-value; The least square method; The maximum likelihood method; Sample size

李世杰，吕悦军，刘静伟，2018．古登堡-里希特定律中的值统计样本量研究．震灾防御技术，13（3）：636—645．

10.11899/zzfy20180315

中国地震局地壳应力研究所中央级公益性科研院所基本科研业务专项（ZDJ2016-03）

2018-02-28

李世杰，男，生于1991年。硕士研究生。主要从事工程地震研究。E-mail：729883299@qq.com

吕悦军，男，生于1966年。研究员。研究领域：工程地震。E-mail：luyj1@263.net