引言

地震序列是在某一时间段内连续发生在同一震源体内的一组按次序排列的地震。主震是地震序列中最大的地震，前震是序列中主震前所有地震的统称，余震是序列中主震后所有地震的统称（中华人民共和国国家质量监督检验检疫总局等，2005）。

周蕙兰等（1982）利用主震所释放能量占全序列所释放总能量的比例R_E对地震序列进行划分：当R_E≥99.99%为孤立型（IET，Isolated Earthquake Type），其主要特点是几乎没有前震，也几乎没有余震；当90%≤R_E＜99.99%为主余型（MAT，Mainshock-Aftershock Type），其主要特点是主震震级突出，主震和最大前震、最大余震的震级相差显著。当R_E＜90%为震群型或双震型地震（MMT，Multiple MainshockType），其主要特点是没有突出的主震，主要能量是通过多次震级相近的地震释放出来的。蒋海昆等（2006）统计1970年以来中国大陆记录相对完整的294次5.0级以上地震序列，孤立型、主余型以及多震型的比例分别占23%、59%和18%。

在余震强度分布特征方面，Gutenberg等（1944）在研究世界地震活动时，根据全球各大地震区6级以上地震数目的统计发现地震的震级与频度存在以下关系：lgN=a-bM，它反映了地震序列中大小地震的比例关系及其变化，亦称G-R关系式，式中a值和b值分别反映地震活动水平和强度分布特征（吴开统等，1990）；Bath（1965）提出，主震震级与最大余震的震级之差平均为1.2，史称“巴特定律”；马宏生等（2005）根据G-R关系以及能量-震级经验关系，推导出一个相对独立的应变积累速率公式；蒋海昆等（2006）对1970年以来记录相对完整的294次5.0级以上地震序列资料进行研究，余震序列1年内最大余震震级与主震震级正相关。傅征祥等（2008）对1966—1999年中国大陆浅源主震型序列余震的b值进行分析，88个震例样本的b值平均值为0.73。

在余震序列模拟方面，Eob等（1992）基于稳态泊松和马尔科夫过程的统计模型，以余震释放的总能量作为控制余震序列的长度，利用蒙特卡罗（Monte-Carlo）方法对地震序列的时间间隔和震级大小进行协同模拟。整个过程较为复杂，同时并未对模拟的序列进行回归分析验证模拟系统。刘善琪等（2013）在假定地震的发生与时间的关系遵循稳态泊松模型，采用人工生成b值为常数的地震目录，并用蒙特卡罗方法对b值进行统计分析，研究结果显示影响b值计算误差的主要因素为样本数量、震级间隔以及震级误差。

对余震序列的研究是理解地震过程的重要途径。国内外学者对余震的研究一般是基于完整的地震序列来总结余震的类型和主余震之间的统计特征。任雪梅等（2009）根据主震震级与余震震级、时间间隔与震中距之间的经验关系和统计规律进行探索，结果表明对某些主震后的强余震具有一定的参考价值。王伟锞等（2011）提出余震法来快速判定宏观震中，结果表明利用余震法可在震后6小时根据地震破裂性质推断宏观震中，其准确度、时效性都能为震后快速评估提供较为可靠的依据。

依托国家地震社会服务工程项目的建设，福建省地震局的地震应急服务水平得到较大的提升，而在现有的日常演练过程中缺少与应急救援密切相关的地震余震信息产品。本文以实际需求为导向，基于常用的计算机高级程序设计语言提供的连续均匀分布的独立伪随机数，运用反函数法模拟生成满足G-R频次关系的余震震级分布。根据主震震级自动产生的余震序列，为日常地震演练提供可靠的余震震级模拟数据，以丰富宣传产品种类，提升服务水平质量。

1.   随机数的基本理论

1.1   随机数

真正意义上的随机数（Random number）在某次产生过程中是按照实验过程中表现的分布概率随机产生的，其结果是不可预测的。计算机随机函数产生的“随机数”并不随机，是伪随机数（Pseudo-random number），产生的伪随机数序列并不是相互独立的，每个数均依赖于前一个数，但只要这种相关性弱到一定程度就可以认为是相互独立的，产生的序列是“随机”的。好的随机数产生器至少要满足以下3个要求：

（1）周期性大：随机数序列并非是永不重复的，经过一个很长的周期后，序列形成一个循环重复；

（2）随机性好：即随机数序列中所有随机数之间没有大的关联，或自相关系数很小；

（3）算法简单：为了得到好的统计结果，通常需要较大的样本空间，随机数产生器的速度也是一个非常重要的考虑因素。

1.2   均匀随机数

均匀随机数就是在某个区间内任何一点出现的概率是相等的，没有规律可循，它是产生其它概率分布的随机数的基础和关键（王菊等，2015）。目前使用最多、最广的均匀随机数生成器是基于Lehmer在1951年首先提出的线性同余随机数生成器（Linear Congruence random number generator）。基本方法如下：

对于任意初始整数I₀，随机数序列由如下递推公式确定：

其中M是模数（modulus），M＞0；初始值I₀也称种子（Seed），0≤I₀＜M；a是乘子（multiplier），0＜a＜M；c是增量（additive constant），0≤c＜M；mod是取模运算。x_i为均匀伪随机数，0≤x_i＜1。通常用表达式LCG（M，a，c，I₀）表示上述随机数生成器。考虑到随机数的产生效率，模数M通常取2的整数次方（在32位计算机上，M可取2³¹），取模操作可以用位移操作来代替，提高运算速度。

通过（1）式可以得到[0，1]区间内的均匀随机数，通过下列的线性变换可得到[A，B]区间内的均匀随机数：

1.3   任意概率分布随机数

对于[A，B]区间内的均匀分布随机数x，其归一化条件可以表示为：

任意分布f(x)的随机数x的归一化条件可以表示为：

假定有函数y=F(x)，其导数为f(x)，即有d(F(x))=f(x)dx，那么上述归一化条件可以写成：

与均匀分布随机数的归一化条件（2）式对比，显然y是[F(A)，F(B)]内均匀分布的随机数。所以生成[A，B]区间内任意分布f(x)的随机数做法如下：

（1）求f(x)的不定积分F(x)；

（2）生成[F(A)，F(B)]内的均匀分布随机数y；

（3）求解F(x)=y，得到x=F^-1(y)。

此x就是[A，B]区间内满足f(x)分布的随机数。由于涉及求解F(x)的反函数，被称为反函数法。但该方法也存在局限性：只有在f(x)可积并且累积概率函数F(x)的反函数可解析时高效使用。

2.   基于G-R关系的余震模型建立

2.1   模型建立

根据古登堡和里克特给出的余震震级与频度关系：lgN=a-bM，或：

对上述未归一化的概率密度函数进行积分可得累积分布函数：

其反函数为：

2.2   模拟步骤

模拟一个主震为M₀级的地震的余震序列（最大余震记为M_max，最小可监测的余震记为M_min）的步骤为：

（1）使用计算机程序语言产生一个[0，1]的均匀随机数x；

（2）通过线性变换转换将x为[F(M_min)，F(M_max)]的随机数F：

（3）本次生成的余震震级为：

（4）重复第（1）步，产生下一次余震震级，直到所有余震模拟结束。

3.   模拟结果及分析

为了不失一般性，本文以MATLAB软件为平台，模拟7.0级地震的余震序列，为保证实验结果的重现性，使用2017作为随机数的种子，即程序运行前均使用rand（'seed'，2017）命令进行初始化。

3.1   均匀随机数产生及相关性检验

利用rand函数，生成10万个均匀随机数，用连续的两个随机数作为点（x，y）的坐标，形成的99999个点，如图 1所示，未出现条带结构和规则的网格结构，证明随机数之间的相关性较弱。同时图 2给出的频次柱状图，表明随机数在各个区间内的分布是均匀的，10万个随机数均匀地分布在20个小区间内。

图 1  随机数相关性检验

Figure 1.  Correlation test of random numbers

下载: 全尺寸图片 幻灯片

图 2  随机数频次分布图

Figure 2.  Frequency distribution of random numbers

下载: 全尺寸图片 幻灯片

3.2   模拟参数

根据巴特（Bath）定律：主震震级与最大余震的震级之差平均为1.2。假定主震震级为M₀，模拟的余震最大震级M_max=M₀-1.2=5.8。而利用福建省现有的地震台网，可监测网内1.0-8.0级地震，模拟的余震最小震级M_min可取1。同时参考傅征祥等（2008）对88个震例的样本统计，得到b的平均值0.73。

3.3   余震序列模拟

模拟的余震震级分布如图 3所示。绝大多数余震在3级以下，少数余震震级可达5级以上。图 4是模拟的余震的概率密度函数和累积分布函数图，显然1.2级以下地震约占48.96%，2级以下地震约占86.97%，4级以下地震占99.58%。

图 3  余震震级分布图

Figure 3.  Magnitude distribution of aftershocks

下载: 全尺寸图片 幻灯片

图 4  余震概率密度函数和累积分布函数

Figure 4.  PDF and CDF of aftershocks

下载: 全尺寸图片 幻灯片

对模拟的余震序列的震级与频次的对数进行线性拟合（频次与震级进行指数拟合），拟合的b值为0.7289，a值为5.5583。即由模拟得到的余震序列的震级与频次拟合的G-R关系为：

这与模拟设置的b=0.73是一致的，并且相关系数高达0.9997，系统是协调的。表 1为余震序列的频次以及根据（9）式计算的拟合频次，拟合的误差控制均在合理范围内。

表 1  模拟的余震频次及拟合频次对比表

Table 1.  Comparison of simulated with fitted frequency of aftershocks

震级M	1.2	1.6	2.0	2.4	2.8	3.2	3.6	4.0	4.4	4.8	5.2	5.6
模拟频次	48964	25291	12710	6405	3303	1615	859	428	229	105	58	33
拟合频次	48261	24662	12602	6440	3291	1682	859	439	224	115	59	30
相对误差	1.4%	2.5%	0.8%	0.5%	0.4%	4.1%	0	2.6%	2.2%	9.5%	1.7%	9.1%

下载: 导出CSV 

| 显示表格

图 5是拟合的余震频次分布。

图 5  拟合的余震频次分布图

Figure 5.  Frequency distribution of fitted aftershocks

下载: 全尺寸图片 幻灯片

3.4   b值的影响

b值分别选取0.5555、0.6180、0.7123和0.8234生成余震序列，并进行拟合，拟合后的b值和相关系数如表 2所示，初始b值与拟合的b值相差很小，说明使用本文的随机数模拟的余震序列震级频次满足G-R关系式。

表 2  不同b值的拟合效果

Table 2.  The fitting results of different b value

初始b值	0.5555	0.6180	0.7123	0.8234
拟合b值	0.5560	0.6209	0.7084	0.8174
b值相对误差	0.09%	0.47%	0.55%	0.73%
相关系数	0.9999	0.9999	0.9994	0.9991

下载: 导出CSV 

| 显示表格

3.5   随机数种子的影响

上述实验都是以2017作为随机数的种子，下面以12322、350003、1234567、19491001以及20080808作为随机数的初始种子，得到的实验结果如表 3所示，拟合的b值都在0.73附近，相关性很强，进一步说明模拟余震序列是可靠的。

表 3  不同随机数种子的拟合效果

Table 3.  The fitting results of different seeds

随机数种子	12322	350003	1234567	19491001	20080808
拟合b值	0.7214	0.7162	0.7364	0.7379	0.7278
b值相对误差	1.18%	1.89%	0.88%	1.08%	0.30%
相关系数	0.9991	0.9990	0.9994	0.9995	0.9999

下载: 导出CSV 

| 显示表格

4.   结语

大震后通常伴随着频繁的余震活动，余震活动时空分布特征是反映地震断层几何学和力学性质的良好指标，一直是中外地震学家密切关注的领域之一。本文提出的使用反函数法生成余震序列的方法可行、高效，符合G-R震级频次关系，对模拟的地震序列回归得到的b值与初始b值一致，并且不受随机数种子等影响，系统稳定性强。模拟生成的余震震级数据既有助于增强地震应急救援演练的现实性，也有助于丰富地震应急宣传产品，提升地震部门的履职能力。