引言

古登堡和里希特在1941年美国地质学会会刊中提出全球的地震活动服从经验关系lgN=a-bM(式中M表示震级，N表示震级≥M的地震次数，a、b是常数)。一般认为a值代表地区的地震活动总体水平(Rundle，1989)；b值代表地震活动的大小地震数量的比例，是地震活动研究中的重要参量。研究表明，b值具有明确的物理意义，与地壳的介质特性、应力状态和不均匀性有关，能反映所研究区域的地质构造特征(王熠熙等，2015；谢卓娟等，2015)以及地震的震源特征(Schorlemmer等，2005；Gulia等，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   模拟地震目录

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

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

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

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

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

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

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

式中，$ \beta = b{\rm{ln}}(10)$，M₀为计算所用截止的震级，M_u为震级上限。

设定b=1、M₀=2.0、M_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年以来汾渭地震带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  汾渭地震带范围及震中分布

Figure 1.  The distribution of earthquakes in the Fen-Wei seismic tectonic zone

下载: 全尺寸图片 幻灯片

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

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

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

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

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

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

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

Figure 2.  The fitting curve of mean b values calculated by maximum likelihood method under different sample sizes

下载: 全尺寸图片 幻灯片

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

Figure 3.  The fitting curve of mean b values calculated by least squares method under different sample sizes

下载: 全尺寸图片 幻灯片

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

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

表 1  最大似然法模拟结果

Table 1.  Simulation results by the maximum likelihood method

N	b1	Pr-	Pr	Pr+	σ
10	1.221	0.46	0.14	0.40	0.42
20	1.211	0.46	0.25	0.29	0.38
50	1.145	0.44	0.30	0.16	0.16
100	1.061	0.53	0.37	0	0.13
200	1.039	0.52	0.48	0	0.11
400	1.036	0.47	0.53	0	0.09
500	1.025	0.45	0.55	0	0.08
1000	1.003	0.30	0.70	0	0.08

下载: 导出CSV 

| 显示表格

表 2  最小二乘法模拟结果

Table 2.  Simulation results by the least squares method

N	b1	Pr-	Pr	Pr+	σ
10	0.870	0.79	0.08	0.13	0.52
20	0.895	0.76	0.07	0.17	0.48
50	0.918	0.68	0.27	0.05	0.42
100	0.941	0.64	0.28	0.08	0.37
200	0.942	0.46	0.48	0.06	0.32
400	0.955	0.48	0.52	0	0.28
500	0.982	0.45	0.50	0.05	0.16
1000	0.996	0.20	0.80	0	0.12

下载: 导出CSV 

| 显示表格

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

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

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

Figure 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

震级档M	地震个数	年平均发生率
2.0—2.4	5038	119.95
2.5—2.9	2006	47.76
3.0—3.4	734	17.48
3.5—3.9	230	5.61
4.0—4.4	110	2.62
4.5—4.9	91	2.17
5.0—5.4	15	0.36
5.5—5.9	10	0.23
6.0—8.5	3	0.07
注：年平均发生率指震级≥M的年均地震数，代表地震活动水平。

下载: 导出CSV 

| 显示表格

根据表 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

样本量	均匀抽样计算b值
	最小二乘法	最大似然法
10	0.60	0.897
50	0.62	0.923
100	0.63	0.816
200	0.63	0.728
300	0.65	0.767
500	0.69	0.752
700	0.70	0.756
1000	0.72	0.746
2000	0.74	0.751
5000	0.75	0.753

下载: 导出CSV 

| 显示表格

由表 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种方法计算结果的差异越来越小。但由于半对数坐标下不同震级档数据权重的不对等性，国外已极少使用最小二乘法。