2017, Vol. 60
概率地震危险性估计是编制地震区划图和场点地震危险性分析的重要工作内容.自概率地震危险性分析方法提出以来(Cornell，1968; McGuire，1976)，绝大数的研究工作都假设地震发生在时间上符合泊松过程.近年来，随着人们对地震构造活动习性研究的深入，时间依从的地震活动性模型也应用到地震危险性分析中来(Petersen et al., 2007).尽管如此，地震发生在时间上符合泊松分布仍然是概率地震危险性分析的基本理论假设.因此，在计算地震发生率前需要删除地震目录中的余震.然而余震时空丛集究竟会对概率地震危险性分析产生多大的影响，这是本文提出的基本科学问题.
事实上，关于余震对概率地震危险性分析结果的影响，已有一些学者作了研究和探讨.Öncel和Alptekin (1999)在North Anatolian Fault Zone分别采用删除余震和未删除余震的地震目录计算了该地区的地震活动性参数，比较了两个参数的差别，定性描述了余震可能对地震危险性结果的影响.Liu等(2002)、杨家亮等(2006)分别对河北地区和邢台地区余震对地震危险性的影响做了研究，认为余震对地震危险性的影响是较大的；Beauval等(2006)采用数值模拟的方法对法国南部余震对概率地震危险性结果的影响做了计算，认为余震对概率地震危险性计算结果没有显著影响.也有学者讨论了余震对地震危险性分析的意义，认为现在有些地区地震是数百年前地震的余震，因此地震危险性分析中可能会高估这些地区的地震危险性程度(Stein and Liu, 2009).
先前关于余震对地震危险性影响的研究或仅限于单个场点或仅作定性的讨论，研究结果不具普适性.本文中我们不改变概率地震危险性分析的理论框架，只考虑余震参与地震发生率计算对概率地震危险性计算的影响.为了避免研究结果的特殊性，我们选取中国东北、华北和川滇三个地震活动性差异较大的地区作为研究区，采用基于Epidemic Type Aftershocks-Sequence(ETAS)(Ogata, 1988, 1998)地震活动性模型的蒙特卡罗模拟方法来研究余震对概率地震危险性分析结果的影响.
2 研究思路本文的基本思路是，使用实际观测到的地震事件估计ETAS地震活动性模型参数，然后根据估计的参数使用蒙特卡罗模式方法模拟地震目录，将模拟得到的目录称为Cluster目录，然后采用余震删除方法将地震目录中的余震删除(本文采用基于ETAS模型的余震删除方法(Zhuang et al., 2002))，将删除余震后的目录称为Poisson目录.使用这两套地震目录进行地震危险性计算，得到两套地震危险性计算结果(PGA)，比较二者差异，研究余震对概率地震危险性计算结果的影响.我们使用的方法和Beauval等(2006)的方法相似.具体思路见流程图 1.
|
图 1 研究思路流程图 Fig. 1 Flow chart of quantifying the impact of aftershocks on PSHA |
文中我们使用了目前国际上研究较多的传染型余震序列(ETAS)模型来进行地震活动性参数的估计和地震目录的模拟.ETAS模型由Ogata(1988, 1998)提出，用于研究地震的时空演化特征.该模型将地震事件看成是由主震(背景)及其触发的余震构成.在ETAS模型中主震触发余震，余震也会产生自己的余震.地震的时空分布被定义为：
|
(1) |
其中κ(M)=Aexp[α(M－Mc)]，
|
κ(M)为震级为M的地震产生余震能力的函数；g(t)为余震发生率随时间衰减的概率密度函数；f(x, y|M)为余震的空间分布概率密度函数.模型中的参数(ν, A, α, c, p, D, q, γ)可采用最大似然法估计，关于ETAS模型参数估计详见Ogata(1998)的文章，这里不再赘述.
模拟地震目录所采用的算法为Ogata(1998)提出的算法，我们将这一算法做了修改，对每个地震及其产生的余震进行了标记，方便选择任何想要的地震事件及其所产生的余震序列.根据Ogata(1998)的算法步骤，我们将模拟地震目录的具体算法重写如下：
(1) 设整数a=b=c=0，i=1.
(2) 设sa=0，生成一个均匀分布的随机数Ub ∈(0, 1]；使Λc=ν0，ua=－log(Ub)/Λc.
(3) 若ua>T，停止模拟.否则使ti=ua，让J=0，转到第8步.
(4) 使b=b+1并生成一个均匀分布的随机数Ub ∈(0, 1]，使a=a+1并让ua=－log(Ub)/Λc.
(5) 使sa=sa－1+ua，若sa>T，停止运行；否则，使b=b+1并生成随机数Ub.
(6) 若Ub>λA(sa)/Λc，使c=c+1并使Λc=λA(sa)，转到第4步.
(7) 使ti=sa，b=b+1并生成随机数Ub，选择一个最小的J使
(8) 若J=0，生成一个点(xi, yi) ∈ R，点的生成是根据细化算法(thinning)(Lewis and Shedler 1979; Ogata 1981) 产生非均匀Poisson分布的点，将在下文中详细介绍.
(9) 若J≠0(J>0)，使b=b+1并生成随机数Ub和Vb，使xi=xJ+r(Ub, J)cos(2πVb)，yi=yJ+r(Ub, J)sin(2πVb)，r(·, ·)将在下文中介绍.
(10) 若(xi, yi)不在研究区域R中，则转到第4步.
(11) 生成一个震级Mi，使c=c+1并使Λc=λA(ti+)；使i=i+1，转到第4步.
上述中：
|
其中
从J的取值能够看出λA(ti)在进行求和时并不包括ti，但是λA(ti+)则是包括了ti项.
上述第9步中r(Ub, J)为点(xi, yi)与点(xJ, yJ)的距离.根据本文中所采用的模型可写为：
|
下面将详细介绍一下上述第8步中的细化算法.设在一个区域内点的发生率函数为λ(x, y)，文中λ(x, y)为根据研究区的实际地震目录，采用Frankel(1995)的方法进行空间光滑得到.设λ*=max{λ(x, y)}，则生成由n个点组成的二维非平稳Poisson点过程步骤如下：
(1) 以λ*为发生率，在研究区域R内生成由n*个点组成的二维平稳Poisson点过程.
(2) 删除不在研究区域内的点，将剩余的点标记为(X1*, Y1*), (X2*, Y2*), …, (Xm*, Ym*)，并使X1* < X2*, …, Xm*.
(3) 设i=1, k=0，生成一个均匀分布的随机数Ui∈(0, 1]，若Ui≤λ(Xi*, Yi*)/λ*，使k=k+1，Xk=Xi*, Yk=Yi*.
(4) 使i=i+1，若i≤n，转到第3步.
(5) 返回点(X1, Y1), (X2, Y2), …, (Xn, Yn).
根据上述的ETAS模型估计了东北、华北、川滇三个研究区域的地震活动性参数(见表 1)，其中最后一列的b值是根据G-R关系模型(Gutenberg and Richter, 1944)估计得到的.根据表 1中的参数，使用上述模拟算法，生成了三个研究区域的模拟地震目录.为了使计算结果稳定，往往模拟多套地震目录进行计算，此外也可以使模拟地震目录的时间尺度足够长.为了节省计算时间同时兼顾计算机的计算能力，本文采用二者相结合的方法，每次模拟10000年的地震目录，模拟100次，从而使三个地震区模拟的地震目录时间尺度为1000000年.图 2为模拟得到的三个研究区域地震目录的空间分布.
|
表 1 华北、东北及川滇地区地震活动性参数 Table 1 Seismicity parameters of North China, Northeast China, and Sichuan-Yunnan regions |
|
图 2 三个研究区模拟地震目录空间分布图 Fig. 2 Spatial distributions of simulated earthquake catalogues of three study areas |
本文采用Frankel(1995)提出的方法进行概率地震危险性计算.即将研究区划分成0.1°×0.1°的网格，统计每个网格中的地震频度，然后对每个网格中地震频度进行空间光滑，得到空间光滑的地震活动性模型，采用该模型进行地震危险性计算.
对于某一场点，计算场点处地震动参数值u超过给定地震动参数值0u的年发生率为(Reiter, 1990)：
|
(2) |
其中m0为起算震级，本文中取4.0级；ni(m0)为第i个潜源区m≥m0的地震年发生率，pi(m)为震级的概率密度函数，pi(r)为场点到潜源距离的概率密度函数；P(u>u0|m, r)为在距离场点r处，震级为m的地震产生的地震动值u超过给定值u0的概率.
对点源(2) 式可简写为：
|
(3) |
其中ri为空间格点到计算场点的距离.
计算条件概率P(u>u0|m, r)需要用到地震动预测方程(衰减关系)，本文中采用的是俞言祥等(2013)建立的地震动预测方程.已知地震动预测方程的方差，则P(u>u0|m, r)可根据高斯误差函数计算.
对于每个场点，我们事先给定一系列地震动参数值，采用(3) 式计算该场点处地震动超过这一系列给定值的概率值，形成一条超越概率曲线(如图 3中红线所示)，若要求取某一超越概率下的地震动值(如图 3中绿线所示)，则可根据超越概率曲线上的数据点进行插值求取.
|
图 3 超越概率曲线示意图 Fig. 3 Hazard curves of annual probability of exceedance versus peak acceleration |
采用上述地震活动性模型和地震危险性分析方法，我们计算了Possion模型和Cluster模型的PGA值.我们将根据Possion地震目录计算得到的PGA值称为标准值，根据Cluster的地震目录计算得到的PGA值称为测量值.则绝对差值为：
|
(4) |
相对差值为：
|
(5) |
图 4为东北地震区分别采用未删除和删除余震的地震目录计算得到的50年超越概率10%的PGA分布(a，b), 及两个结果之间的绝对差值(c)和相对差值(d).从图 4a中可以看出，相对于图 4b PGA>50 Gal的区域有所扩大，但变化并不明显.图 4d为东北地震区余震对PGA影响的空间分布，可以看出，余震对PGA的影响主要在3%~7%之间.余震影响很显著的场点大都分布于研究区域的边界地区，这是由于边界效应引起的.此外我们还发现在PGA较高的地区余震对其影响较PGA较低的地区要小.从图 4c可以看出在东北地震区不删除余震与删除余震计算得到的PGA绝对差值最大不超过5.0 Gal，这说明在东北地震区余震对地震危险性的影响是微小的，不会对PGA值产生实质的影响.
|
图 4 东北地震区(a)采用未删除余震的模拟地震目录计算得到的50年超越概率10%的PGA分布；(b)采用删除余震后的地震目录计算得到的50年超越概率10%的PGA分布；(c)二者的绝对差值；(d)二者的相对差值 Fig. 4 Spatial distributions of PGAs with 10% probability of exceedance in 50 years calculated using Poisson (a) and clustering (b) earthquake catalogues, respectively for Northeast China. (c) The absolute difference between the two sets PGA in (a) and (b). (d) The relative difference between the two sets PGA. |
图 5为东北地震区不同超越概率水平下余震对PGA影响的分布直方图，可以看出余震对地震危险性影响的概率密度分布呈现单峰，且影响集中分布在某一个区间值内，这说明了余震影响的分布来自某个单一概率分布族，因此在同一个区域内我们可以用均值及方差来描述余震对地震危险性影响的程度.
|
图 5 东北地震区余震对不同超越概率水平下PGA影响的分布直方图 Fig. 5 Probability density distribution of the effect of aftershocks on PGA with different probability levels of exceedance in 50 years in Northeast China |
图 6为余震影响与地震50年超越概率水平的关系曲线，可以看出随着地震50年超越概率水平的增大余震对地震危险性的影响呈变大趋势.图中灰色区域为余震影响值90%置信区间，上、下两条虚线分别为95%和5%的边界值，绿色实线为均值曲线.随着50年超越概率值的增大，余震影响的90%置信区间范围也变大，这从余震影响值的概率密度分布也可以看出(图 5).
|
图 6 东北地震区余震影响与50年超越概率的关系曲线 Fig. 6 Relation curves between the effect of aftershocks and the probability of exceedance in 50 years in Northeast China |
从图 6中还可以看出在东北地震区余震对50年超越概率10%PGA的影响均值约为6%，影响的90%置信区间为[4.6%~7.8%].余震对其他概率水平下PGA的影响见表 2.
|
|
表 2 三个研究区不同超越概率水平下余震对PGA的影响(%) Table 2 Effect of aftershocks on PGA with different probability levels of exceedance in three study regions (%) |
图 7为东北地震区选取的地震危险性较强的点(图 4a中A点)与地震危险性较弱的点(图 4a中B点)余震影响与50年超越概率的关系曲线.从图中可以看出随着50年超越概率水平的提高，余震对地震危险性的影响增大.图中兰色曲线绝大部分在紫色曲线之上，这也再次说明了余震对地震危险性弱的地区的影响要大于地震危险性强的地区.
|
图 7 东北地震区A点与B点余震影响与50年超越概率关系 Fig. 7 Relationship between the aftershock effect and the probability of exceedance in 50 years for point A and point B in Northeast China |
图 8为图 4a中A、B两点考虑余震影响和不考虑余震影响的超越概率曲线，可以看出随着峰值加速度的增大，超越概率的差异是逐渐变小的，这说明了随着超越概率水平的降低，余震影响是逐渐变小的.另外我们发现在地震危险性较强的场点(A点)其两个超越概率的差异大于地震危险性较弱的场点，这意味着在地震活动性较强的地区余震对地震危险性影响的绝对差值是较大的，只是由于其本身的标准值较大，从而导致其相对差值较小.
|
图 8 图 4a中A、B两点超越概率曲线 Fig. 8 Curves of annual probability of exceedance versus peak acceleration for points A and B in Fig. 4a |
图 9为A、B两点分别考虑和不考虑余震影响时各震级档、距离档的地震对地震危险性的贡献.可以看出在A点地震危险性主要受近场地震，特别是近场较高震级地震的影响，远场地震对该点地震危险性的贡献是极小的，并且相对于不考虑余震影响(图 9b)，我们发现考虑余震影响后(图 9a)各震级档、距离档地震对该点地震危险性的贡献变化是很小的，图中很难明显看出；而在B点，近场各震级档地震对该点地震危险性的贡献分布相对于A点是较为均匀的，并且B点地震危险性受远场地震的影响相对来说是大于A点的.比较图 9c和9d可以看出B点在考虑余震影响后近场小震对该点地震危险性的影响是明显增大的.这也解释了为何余震影响在地震危险性较弱的地区偏大的原因.
|
图 9 图 4a中A、B点分别考虑余震影响和未考虑余震影响时各震级档、距离档地震对这两点地震危险性的贡献 Fig. 9 Contributions of magnitude and distance to hazard estimations based on Poisson and cluster catalogues at sites A and B in Fig. 4a |
上述计算了东北地震区余震对地震危险性结果的影响，并对余震影响特征取得了一些认识.为了更多地了解余震对我国其他地区的影响，同时也为了验证上述取得的认识的普遍性，我们又对华北地震区和川滇地震区余震对地震危险性的影响做了研究.
从华北地震区和川滇地震区的计算结果可以看出，余震对地震危险性的影响具有和东北地震区相同的特征：如在地震危险性强的地区余震对其影响要小于地震危险性弱的地区；随着50年超越概率水平的提高，余震影响也逐渐增大，置信区间范围也变宽；研究区内余震影响的概率密度分布也都呈现单峰，且具有明显的优势分布值；在地震危险性较强的场点，考虑余震影响后，各震级档、距离档的地震对该点地震危险性的贡献变化不明显，而在地震危险性较弱的场点，考虑余震影响后近场小震对该点地震危险性的贡献明显增大.因此我们认为上述余震对地震危险性的影响特征是普适的(至少在中国地区是这样的).
从表 2中可以看出华北地震区余震对其地震危险性的影响要强于东北地震区，显然这是由于在华北地震区地震触发余震的能力要远强于东北地震区.川滇地震区余震对其地震危险性的影响与东北地震区大体相当而低于华北地震区.华北地震区余震影响大于川滇地震区说明了在地震触发余震能力相当的情况下，地震活动性越强的地区地震危险性受余震的影响越小.此外，从表 1中第8列D值可以看出在三个研究区中华北地震区余震分布最为集中，这也是导致余震影响偏大的原因之一.
6 结论与讨论本文采用蒙特卡罗模拟方法，系统研究了余震对地震危险性结果的影响.结果表明余震对50年超越概率10%水平下峰值加速度的影响均值为6%左右，最大为10%左右.超越概率水平越高，余震对峰值加速度的影响越大.余震对地震危险性较强地区的影响要小于地震危险性较弱的地区，这是由于在地震危险性较弱的地区地震危险性主要受近场小震的影响，在考虑余震后，小震频数显著增加，因此对场点的地震危险性影响相对较大.纵观三个研究区考虑余震和未考虑余震计算得到的峰值加速度的绝对差值，我们发现两个峰值加速度的绝对差值最大约为15 Gal，这说明在考虑余震后地震危险性分析结果不会有显著变化.事实上，在地震危险性分析中，由于模型或参数的不确定性对最终计算的结果的影响也会大于10%(Beauval et al., 2006)，因此我们认为在一般的地震区划或地震危险性分析中可考虑不用删除余震.
文中我们采用基于ETAS模型的余震删除方法去除了地震目录中的余震，有研究表明该方法会比其他方法删除更多的余震(Wang et al., 2010).可以预见若采用地震危险性分析中常用的Gardner和Knopoff(1974)方法进行余震删除，那么计算的余震对概率地震危险性结果的影响将变小，因此文中得到的余震对概率地震危险性计算结果的影响是一个相对保守的值.
| Beauval C, Scotti O, Bonilla F. 2006. The role of seismicity models in probabilistic seismic hazard estimation: comparison of a zoning and a smoothing approach. Geophysical Journal International, 165(2): 584-595. DOI:10.1111/gji.2006.165.issue-2 | |
| Cornell C A. 1968. Engineering seismic risk analysis. Bulletin of the Seismological Society of America, 58: 1583-1606. | |
| Frankel A. 1995. Mapping seismic hazard in the central and eastern United States. Seismological Research Letters, 66(4): 8-21. DOI:10.1785/gssrl.66.4.8 | |
| Gardner J K, Knopoff L. 1974. Is the sequence of earthquakes in southern California, with aftershocks removed, Poissonian?. Bulletin of the Seismological Society of America, 64(5): 1363-1367. | |
| Gutenberg B, Richter C F. 1944. Frequency of earthquakes in California. Bulletin of the Seismological Society of America, 34(4): 185-188. | |
| Lewis P A W, Shedler G S. 1979. Simulation of nonhomogeneous Poisson processes by thinning. Naval Research Logistics Quarterly, 26(3): 403-413. DOI:10.1002/(ISSN)1931-9193 | |
| Liu Z H, Fu Z X, Jin X S, et al. 2002. Effect of aftershocks on assessment of seismic hazard for Hebei seismic zone, China. Earthquake Research in China, 16(2): 160-166. | |
| McGuire R K. 1976. Fortran computer program for seismic risk analysis. US Geological Survey Open-File Report: 76-67. | |
| Ogata Y. 1981. On Lewis' simulation method for point processes. IEEE Transactions on Information Theory, 27(1): 23-31. DOI:10.1109/TIT.1981.1056305 | |
| Ogata Y. 1988. Statistical models for earthquake occurrences and residual analysis for point processes. Journal of the American Statistical Association, 83(401): 9-27. DOI:10.1080/01621459.1988.10478560 | |
| Ogata Y. 1998. Space-time point-process models for earthquake occurrences. Annals of the Institute of Statistical Mathematics, 50(2): 379-402. DOI:10.1023/A:1003403601725 | |
| Öncel A O, Alptekin Ö. 1999. Effect of aftershocks on earthquake hazard estimation: An example from the North Anatolian Fault Zone. Natural Hazards, 19(1): 1-11. DOI:10.1023/A:1008139802609 | |
| Petersen M D, Cao T Q, Campbell K W, et al. 2007. Time-independent and time-dependent seismic hazard assessment for the State of California: Uniform California Earthquake Rupture Forecast Model 1.0. Seismological Research Letters, 78(1): 99-109. DOI:10.1785/gssrl.78.1.99 | |
| Reiter L. 1990. Earthquake Hazard Analysis: Issues and Insights. New York: Columbia University Press. | |
| Stein S, Liu M. 2009. Long aftershock sequences within continents and implications for earthquake hazard assessment. Nature, 462(7269): 87-89. DOI:10.1038/nature08502 | |
| Wang Q, Jackson D D, Zhuang J C. 2010. Are spontaneous earthquakes stationary in California?. Journal of Geophysical Research: Solid Earth, B8(B08310). DOI:10.1029/2009JB007031 | |
| Yang J L, Ye M, Li H Y. 2006. The affection of aftershock activity in seismic risk evaluation-taking Xingtai aftershock area as example. North China Earthquake Sciences, 24(2): 46-49. | |
| Yu Y X, Li Y S, Xiao L. 2013. Development of ground motion attenuation relations for the new seismic hazard map of China. Technology for Earthquake Disaster Prevention, 8(1): 24-33. | |
| Zhuang J C, Ogata Y, Vere-Jones D. 2002. Stochastic declustering of space-time earthquake occurrences. Journal of the American Statistical Association, 97(458): 369-380. DOI:10.1198/016214502760046925 | |
| 杨家亮, 叶明, 李红印. 2006. 地震危险性评价中余震活动的影响——以邢台余震区为例. 华北地震科学, 24(2): 46–49. | |
| 俞言祥, 李山有, 肖亮. 2013. 为新区划图编制所建立的地震动衰减关系. 震灾防御技术, 8(1): 24–33. DOI:10.11899/zzfy20130103 | |