• 首页关于本刊投稿须知期刊订阅编委会期刊合作查询检索English
古登堡-里希特定律中的b值统计样本量研究
古登堡-里希特定律中的b值统计样本量研究
李世杰*, 吕悦军*, 刘静伟
(中国地震局地壳应力研究所,北京 100085)
 [收稿日期]: 2018-02-28
摘要

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



引言

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

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

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

1 资料概况
1.1 模拟地震目录

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

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

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

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

首先,构造概率分布模型:

$ F(M) = \int_{{M_0}}^M {f(M){\rm{d}}M} $ (1)

其中,f(M)为震级的概率密度函数,表示为:

$ f(M) = \frac{{\beta {{\rm{e}}^{ - \beta M}}}}{{{{\rm{e}}^{ - \beta {M_0}}} - {{\rm{e}}^{ - \beta {M_{\rm{u}}}}}}} $ (2)

其次,将公式(2)带入公式(1),得到关于模拟震级M的函数:

$ F(M)({{\rm{e}}^{ - \beta {M_0}}} - {{\rm{e}}^{ - \beta {M_{\rm{u}}}}}) = {{\rm{e}}^{ - \beta {M_0}}} - {{\rm{e}}^{ - \beta M}} $ (3)
$ {{\rm{e}}^{ - \beta M}} = {{\rm{e}}^{ - \beta {M_0}}} - F(M)({{\rm{e}}^{ - \beta {M_0}}} - {{\rm{e}}^{ - \beta {M_{\rm{u}}}}}) $ (4)
$ M = - \frac{1}{\beta }{\rm{ln(}}{{\rm{e}}^{ - \beta {M_0}}} - F(M)({{\rm{e}}^{ - \beta {M_0}}} - {{\rm{e}}^{ - \beta {M_{\rm{u}}}}})) $ (5)

式中,$ \beta = b{\rm{ln}}(10)$M0为计算所用截止的震级,Mu为震级上限。

设定b=1、M0=2.0、Mu=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等,2000Rotondi等,2002)。经统计分析(焦远碧等,1990;黄伟琼等, 1994a, 1994b;周公威等,2007谢卓娟等,2012a徐伟进等,2014),1970年以来汾渭地震带M 2.0以上的地震资料基本完整,故本研究选取1970年1月—2010年12月,M≥2.0的8184次地震资料进行研究分析,其中M 2.0—2.9地震7044次、M 3.0—3.9地震964次、M 4.0—4.9地震148次、M 5.0—5.9地震25次、M 6.0—6.9地震3次,所用的地震目录均已删除前震、余震。汾渭地震带4.7级以上地震的震中分布见图 1


图 1 汾渭地震带范围及震中分布 Fig. 1 The distribution of earthquakes in the Fen-Wei seismic tectonic zone
2 模拟地震目录的结果分析

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

文中使用的最大似然法为Utsu(1965)提出的公式:

$ b = \frac{{\lg {\rm{e}}}}{{\bar M - {\rm{(}}M{\rm{0}} - \Delta M/2)}} $ (6)

其中,$ \bar M$是平均震级;$\Delta M $=0.1,用来对实际地震记录进行校正。当实际地震的震级记录精度为$\Delta M $、计算所用截止的震级为M0时,其实际代表($M{\rm{0}} - \Delta M/2 $)≤M0 < ($ M{\rm{0}} + \Delta M/2$)的地震。为了使拟合的地震目录尽量贴近实际情况,文中模拟的地震目录震级是连续变化的,将震级归并到各自$\Delta M $中,再用最大似然法拟合b值。

最小二乘法使用累积震级-频度关系($\lg N = a - bM $)计算b值,采用0.5作为震级间隔,即震级分档为0.5级。

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


图 2 不同样本量下最大似然法计算得到的平均b值拟合图 Fig. 2 The fitting curve of mean b values calculated by maximum likelihood method under different sample sizes

图 3 不同样本量下最小二乘法计算得到的平均b值拟合图 Fig. 3 The fitting curve of mean b values calculated by least squares method under different sample sizes

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

采用Nava等(2017)提出的方法,分析不同样本量下的b值平均值和设定值之间的差异。设定样本组数Nr=10000,模拟不同容量的样本N来计算b值,并统计得出合理设定值的概率。表 12分别是2种方法在b=1时,不同样本量下的b值计算值、标准差及其正确估值概率,其中,b1表示设定b=1条件下得到的b值平均计算值,精度$\Delta b $=0.025;Pr-表示过低估计的概率(Pr- < 0.975);Pr表示准确估计的概率(0.975≤Pr≤1.025);Pr+表示过高估计的概率(Pr+ > 1.025);$ \sigma $表示标准差。

表 1 最大似然法模拟结果 Table 1 Simulation results by the maximum likelihood method
表 2 最小二乘法模拟结果 Table 2 Simulation results by the least squares method

表 12表明,样本量低于200、b=1.00、精度$\Delta b $=0.025时,得到正确估算值的概率较低(小于0.5),因此要求样本量至少大于200。

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


图 4 不同样本量的b值分布直方图 Fig. 4 The histogram of the b value calculated under different sample sizes

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

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

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

表 3 汾渭地震带地震分档统计和年平均发生率 Table 3 The annual average incidence of different magnitude-class in Fen-Wei seismic zone

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

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

表 4 汾渭地震带b值拟合情况(1500—2010年) Table 4 The b value fitting of earthquakes from 1500—2010 in the Fen-Wei seismic zone

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

4 结论

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

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

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

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

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

参考文献
韩渭宾, 2003. b值在地震预测中的三类应用及其物理基础与须注意的问题[J]. 四川地震, (1): 1-5. DOI:10.3969/j.issn.1001-8115.2003.01.001
韩晓明, 张文韬, 王树波, 等, 2016. 河套地震带的b值时空变化特征分析[J]. 中国地震, 32(3): 522-532. DOI:10.3969/j.issn.1001-4683.2016.03.009
胡聿贤, 1999. 地震安全性评价技术教程[M]. 北京: 地震出版社.
黄玮琼, 时振梁, 曹学锋, 1989. b值统计中的影响因素及危险性分析中b值的选取[J]. 地震学报, 11(4): 351-361.
黄玮琼, 李文香, 曹学锋, 1994a. 中国大陆地震资料完整性研究之一——以华北地区为例[J]. 地震学报, 16(3): 273-280.
黄玮琼, 李文香, 曹学锋, 1994b. 中国大陆地震资料完整性研究之二——分区地震资料基本完整的起始年分布图象[J]. 地震学报, 16(4): 423-432.
黄玮琼, 李文香, 1998. 地震区划中b值统计时空范围的确定[J]. 地震学报, 20(5): 2-6.
焦远碧, 吴开统, 杨满栋, 1990. 我国地震台网监测能力及台网观测条件质量评定[J]. 中国地震, 6(4): 3-9.
刘方斌, 曲均浩, 田兆阳, 等, 2017. 鲁西南聊考断裂带地震危险性评价与活动性分布[J]. 大地测量与地球动力学, 37(8): 802-807.
刘静伟, 吕悦军, 2016. 川滇地区b值空间分布特征及其与震源类型关系的初步探讨[J]. 震灾防御技术, 11(3): 561-572.
潘华, 李金臣, 2006. 地震统计区地震活动性参数b值及ν4不确定性研究[J]. 震灾防御技术, 1(3): 218-224. DOI:10.3969/j.issn.1673-5722.2006.03.006
任雪梅, 2011.地震区划中b值统计的若干问题研究.北京: 中国地震局地球物理研究所.
任雪梅, 高孟潭, 冯静, 2011. 地震目录的完整性对b值计算的影响[J]. 震灾防御技术, 6(3): 257-268. DOI:10.3969/j.issn.1673-5722.2011.03.005
沈建文, 余湛, 邱瑛, 2007. 地震安评中地震活动性的统计区域与b值[J]. 国际地震动态, (3): 1-6. DOI:10.3969/j.issn.0253-4975.2007.03.001
王熠熙, 张辉, 刘双庆, 等, 2015. 河北平原地震带b值时空变化特征[J]. 地震工程学报, 37(01): 188-195. DOI:10.3969/j.issn.1000-0844.2015.01.0188
吴兆营, 薄景山, 刘志平, 等, 2005. 东北地震区b值和地震年平均发生率的统计分析[J]. 东北地震研究, 21(3): 27-32, 72. DOI:10.3969/j.issn.1674-8565.2005.03.002
谢卓娟, 吕悦军, 张力方, 等, 2012a. 基于现代地震资料确定汾渭地震带分区及其地震活动性参数[J]. 地球物理学进展, 27(3): 894-902.
谢卓娟, 吕悦军, 彭艳菊, 等, 2012b. 东北地震区小震资料完整性分析及其对地震活动性参数的影响研究[J]. 中国地震, 28(3): 256-265.
谢卓娟, 吕悦军, 兰景岩, 等, 2013. b值和V4的统计分析及其不确定性对地震危险性分析结果的影响研究[J]. 地震研究, 36(1): 86-92. DOI:10.3969/j.issn.1000-0666.2013.01.013
谢卓娟, 李山有, 吕悦军, 2015. 滇西南地区主要活动断裂的b值空间分布特征[J]. 地球科学(中国地质大学学报), 40(10): 1755-1766.
徐伟进, 高孟潭, 2014. 中国大陆及周缘地震目录完整性统计分析[J]. 地球物理学报, 57(9): 2802-2812.
鄢家全, 韩炜, 高孟潭, 1996. 地震活动性参数的不确定性及其对区划结果的影响[J]. 中国地震, 12(S1): 71-77.
张广伟, 2016. 云南地区地震的重新定位及b值研究[J]. 中国地震, 32(1): 54-62. DOI:10.3969/j.issn.1001-4683.2016.01.005
张建中, 1974a. 蒙特卡洛方法(Ⅰ)[J]. 数学的实践与认识, (01): 28-40.
张建中, 1974b. 蒙特卡洛方法(Ⅱ)[J]. 数学的实践与认识, (02): 43-56.
周公威, 张伯明, 吴忠良, 等, 2007. 中国数字地震台网(CDSN)的近期发展[J]. 地球物理学进展, 22(4): 1130-1134. DOI:10.3969/j.issn.1004-2903.2007.04.018
Maximum likelihood estimate of b in the formula logN=a-bM 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[J]. Geophysical Research Letters, 37(10): L10305.
Nava F. A., Márquez-Ramírez V. H., Zúñiga F. R., et al, 2017. Gutenberg-Richter b-value maximum likelihood estimation and sample size[J]. Journal of Seismology, 21(1): 127-135. DOI:10.1007/s10950-016-9589-1
Rundle J. B., 1989. Derivation of the complete Gutenberg-Richter magnitude-frequency relation using the principle of scale invariance[J]. Journal of Geophysical Research:Solid Earth, 94(B9): 12337-12342. DOI:10.1029/JB094iB09p12337
Rotondi R., Garavaglia E., 2002. Statistical Analysis of the Completeness of a Seismic Catalogue[J]. Natural Hazards, 25(3): 245-258. DOI:10.1023/A:1014855822358
Schorlemmer D., Wiemer S., Wyss M., 2005. Variations in earthquake-size distribution across different stress regimes[J]. Nature, 437(7058): 538-542.
Utsu T., 1965. A method for determining the value of b in a formula logN=a-bM showing the magnitude-frequency relation for earthquakes[J]. 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[J]. Bulletin of the Seismological Society of America, 90(4): 859-869. DOI:10.1785/0119990114


The Study of Sample Size on b-value Statistics in the Gutenberg-Richter's Law
Li Shijie*, Lü Yuejun*, Liu Jingwei
(Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China)
Abstract

The b-value is an important indicator for evaluating the level of seismicity, and is widely used in seismic hazard analysis and earthquake prediction research. The b-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 calculating b values 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 the b values 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 b-value is not reliable. It suggets that LSM is suitable for calculating the b-value of the whole region, while MLM is more advantageous in studying the b-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.



主办单位:中国地震台网中心
版权所有:《震灾防御技术》编辑部
地址:北京西城区三里河南横街5号,   邮编:100045
邮箱:zzfy2006@126.com   电话:59959251;59959462;59959408
访问人数:1765834
古登堡-里希特定律中的b值统计样本量研究
李世杰*, 吕悦军*, 刘静伟
《震灾防御技术》, DOI:10.11899/zzfy20180315