• 首页关于本刊投稿须知期刊订阅编委会期刊合作查询检索English
应急演练中基于G-R关系的主余型地震余震震级模拟
应急演练中基于G-R关系的主余型地震余震震级模拟
王青平*, 王辉山, 肖健, 林岩钊, 周施文, 张树君
(福建省地震局,福州 350003)
 [收稿日期]: 2017-04-21
摘要

在现有的日常地震演练过程中,与应急救援密切相关的地震余震信息产品较为缺乏,直接影响发震构造的判断以及影响场修正等关键环节。本文从计算机系统提供的均匀分布随机数出发,运用反函数法模拟生成余震序列,并进行系统检验,证实该方法产生的余震序列满足G-R频次关系。模拟生成的余震震级数据既有助于增强地震应急救援演练的现实性,也有助于丰富地震应急宣传产品,提升地震部门的履职能力。



引言

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

周蕙兰等(1982)利用主震所释放能量占全序列所释放总能量的比例RE对地震序列进行划分:当RE≥99.99%为孤立型(IET,Isolated Earthquake Type),其主要特点是几乎没有前震,也几乎没有余震;当90%≤RE<99.99%为主余型(MAT,Mainshock-Aftershock Type),其主要特点是主震震级突出,主震和最大前震、最大余震的震级相差显著。当RE<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)。基本方法如下:

对于任意初始整数I0,随机数序列由如下递推公式确定:

$ \left\{ \begin{align} &{{I}_{i}}=\left(a\cdot {{I}_{i-1}}+c \right)\rm{mod}\ \mathit{M} \\ &{{\mathit{x}}_{\mathit{i}}}={{I}_{i}}/M \\ \end{align} \right.\ \ i=1, 2, \cdots $ (1)

其中M是模数(modulus),M>0;初始值I0也称种子(Seed),0≤I0Ma是乘子(multiplier),0<aMc是增量(additive constant),0≤cM;mod是取模运算。xi为均匀伪随机数,0≤xi<1。通常用表达式LCG(MacI0)表示上述随机数生成器。考虑到随机数的产生效率,模数M通常取2的整数次方(在32位计算机上,M可取231),取模操作可以用位移操作来代替,提高运算速度。

通过(1)式可以得到[0,1]区间内的均匀随机数,通过下列的线性变换可得到[AB]区间内的均匀随机数:

$ \xi_i =A+\left(B-A \right){{x}_{i}} $
1.3 任意概率分布随机数

对于[AB]区间内的均匀分布随机数x,其归一化条件可以表示为:

$ \int_{A}^{B}{\rm{d}\mathit{x=}\rm{1}} $ (2)

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

$ \int_{A}^{B}{\mathit{f}\left(x \right)\rm{d}\mathit{x=}\rm{1}} $ (3)

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

$ \int_{F\left(A \right)}^{F\left(B \right)}{{\rm{d}}\left(F\left(x \right) \right)=\int_{F\left(A \right)}^{F\left(B \right)}{{\rm{d}}y=1}} $ (4)

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

(1)求f(x)的不定积分F(x);

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

(3)求解F(x)=y,得到x=F-1(y)。

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

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

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

$ f\left(M \right)=N={{10}^{a-bM}} $ (5)

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

$ f\left(M \right)=\int{f\left(M \right){\rm{d}}\mathit{M}}=\int{{{10}^{a-bM}}{\rm{d}}\mathit{M}=\frac{{{10}^{a-bM}}}{-b\cdot \rm{ln}\ \rm{10}}} $ (6)

其反函数为:

$ M=\frac{a-{\rm{lg}}\left[ -b\cdot F\left(M \right)\cdot \rm{ln}\ \rm{10} \right]}{b} $ (7)
2.2 模拟步骤

模拟一个主震为M0级的地震的余震序列(最大余震记为Mmax,最小可监测的余震记为Mmin)的步骤为:

(1)使用计算机程序语言产生一个[0,1]的均匀随机数x

(2)通过线性变换转换将x为[F(Mmin),F(Mmax)]的随机数F

$ \begin{align} &F=F\left({{M}_{\rm{min}}} \right)+\left(F\left({{M}_{\rm{max}}} \right)-F\left({{M}_{\rm{min}}} \right) \right)\cdot x \\ &\ \ \ =\frac{{{10}^{a-b{{M}_{\rm{min}}}}}+\left({{10}^{a-b{{M}_{\rm{max}}}}}-{{10}^{a-b{{M}_{\rm{min}}}}} \right)\cdot x}{-b\cdot \rm{ln}\ \rm{10}} \\ \end{align} $

(3)本次生成的余震震级为:

$ M=\frac{{\rm{lg}}\left[ {{10}^{-b{{M}_{\rm{min}}}}}+\left({{10}^{-b{{M}_{\rm{max}}}}}-{{10}^{-b{{M}_{\rm{min}}}}} \right)\cdot x \right]}{b} $ (8)

(4)重复第(1)步,产生下一次余震震级,直到所有余震模拟结束。

3 模拟结果及分析

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

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

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


图 1 随机数相关性检验 Fig. 1 Correlation test of random numbers

图 2 随机数频次分布图 Fig. 2 Frequency distribution of random numbers
3.2 模拟参数

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

3.3 余震序列模拟

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


图 3 余震震级分布图 Fig. 3 Magnitude distribution of aftershocks

图 4 余震概率密度函数和累积分布函数 Fig. 4 PDF and CDF of aftershocks

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

$ \rm{lg}\ \mathit{N}=\rm{5}\rm{.5583}-\rm{0}\rm{.7289}\mathit{M} $ (9)

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

表 1 模拟的余震频次及拟合频次对比表 Table 1 Comparison of simulated with fitted frequency of aftershocks

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


图 5 拟合的余震频次分布图 Fig. 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
3.5 随机数种子的影响

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

表 3 不同随机数种子的拟合效果 Table 3 The fitting results of different seeds
4 结语

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

参考文献
傅征祥, 吕晓健, 邵辉成, 等, 2008. 中国大陆及其分区余震序列b值的统计特征分析[J]. 地震, 28(3): 1-7.
蒋海昆, 曲延军, 李永莉, 等, 2006. 中国大陆中强地震余震序列的部分统计特征[J]. 地球物理学报, 49(4): 1110-1117.
刘善琪, 李永兵, 田会全, 等, 2013. 影响b值计算误差的Monte Carlo实验研究[J]. 地震, 33(4): 135-144.
马宏生, 刘杰, 张国民, 等, 2005. 基于G-R关系的应变积累释放模型研究中国大陆强震的分区活动[J]. 地震学报, 27(4): 355-366.
任雪梅, 高孟潭, 刘爱文, 等, 2009. 1900年以来我国西南地区强余震统计特征[J]. 震灾防御技术, 4(2): 200-208.
王菊, 王和明, 徐海龙, 等, 2015. 任意分布随机数的FPGA实现[J]. 火力与指挥控制, 40(4): 173-175, 180.
王伟锞, 李志强, 李晓丽, 2011. 利用余震法快速判定宏观震中的研究[J]. 震灾防御技术, 6(1): 36-48.
吴开统, 焦远碧, 吕培苓, 等, 1990. 地震序列概论. 北京: 北京大学出版社.
中华人民共和国国家质量监督检验检疫总局, 中国国家标准化管理委员会, 2005. GB/T 18207.2-2005防震减灾术语第2部分:专业术语. 北京: 中国标准出版社.
周蕙兰, 房桂荣, 章爱娣, 等, 1982. 余震序列的持续时间[J]. 地震学报, 4(1): 45-54.
Bath M., 1965. Lateral in homogeneities of the upper mantle[J]. Tectonophysics, 2(6): 483-514. DOI:10.1016/0040-1951(65)90003-X
Eob B. C., Soo S. J., 1992. Monte-Carlo simulation of earthquake sequence in the time and magnitude space[J]. The Journal of Engineering Geology, 2(2): 147-154.
Gutenberg B., Richter C. F., 1944. Frequency of earthquakes in California[J]. Bulletin of the Seismological Society of America, 34(4): 185-188.


Magnitude Simulation of Aftershocks Based on G-R Relationfor MAT Earthquake in Earthquake Emergency Response Exercise
Wang Qingping*, Wang Huishan, Xiao Jian, Lin Yanzhao, Zhou Shiwen, Zhang Shujun
(Fujian Earthquake Agency, Fuzhou 350003, China)
Abstract

Insufficient aftershock information, which is one of the key factors to emergency rescue, directly affects the essential steps such as seismogenic structure determination and influence field correction in the daily exercise against earthquake. In this paper we simulate and generate aftershock sequence by using inverse function method on uniformly distributed random numbers generated by computer system, Then we carry out system examination with the simulation results, and finally prove that the aftershock sequence generated by this method satisfies the magnitude distribution from G-R relation. The simulated aftershock sequence not only helps to enhance the reality of emergency rescue exercise against, but also contributes to enrich promotional products of earthquake emergency response and improves the executing capacity of earthquake departments.



主办单位:中国地震台网中心
版权所有:《震灾防御技术》编辑部
地址:北京西城区三里河南横街5号,   邮编:100045
邮箱:zzfy2006@126.com   电话:59959251;59959139;59959462
访问人数:1384039
应急演练中基于G-R关系的主余型地震余震震级模拟
王青平*, 王辉山, 肖健, 林岩钊, 周施文, 张树君
《震灾防御技术》, DOI:10.11899/zzfy20180117