引言

当前,主流的地震区划技术方法是划分潜在震源区的地震危险性概率分析方法(胡聿贤,1999;高孟潭,2002;2006)。该方法最关键的技术环节是根据地震地质构造和历史地震活动性,划分潜在震源区和确定地震活动性参数,要求研究区域内地震地质构造清晰(胡聿贤,1999)。但中强地震活动在空间上具有一定的弥漫性与成丛性的特点,采用划分潜在震源区的地震危险性概率分析方法,评价中强地震的危险性势必会带来较大的不确定性。即使在地震活动性较强的华北地区亦是如此。且在确定潜在震源区的空间分布函数时还要综合考虑众多的判断因子,经验性判断较多,往往会给地震危险性评价结果带来无法预知的人为误差。

由于对美国中东部地震潜源的了解很有限,1995年Frankel以地震目录为唯一输入数据,先将研究区域网格化(0.1°×0.1°),并在每个网格内统计地震个数(地震活动率)(Frankel,1995)。然后将累计地震活动率用圆形高斯光滑函数处理,从而以点源模型评价地震危险性。如此处理,使得潜源就被简化为小方格单元,其最大特点是在地震危险性分析中摒弃了地震构造单元,而使用了点源。圆形空间光滑模型符合普遍认同的原地复发原理,强调历史地震重演,即未来中强震发生在历史地震附近或现代小震密集区。但并不是所有地震发生都遵从该假设,譬如有些小震活动密集区仅处在持续平稳的释放应力阶段,未来并不一定发生破坏性中强地震。中强地震活动的空间分布往往与构造断裂有一定的空间相关性,这就需要在建立地震空间光滑模型时,考虑相关的地震构造因素。于是,Lapajne等(2003)对该方法进行了改进,采用断层导向性的椭圆光滑方法,使评价地震活动性参数显得更加合理。

根据我国地震活动的特征和抗震设防的工作重点,在《中国地震动参数区划图(2001)》的地震区带划分方案中,华北地震区的地震强度和频度仅次于“青藏高原地震区”,位居全国第二。加之华北地区位于我国人口稠密、大城市集中、政治和经济、文化、交通都很发达的地区,地震灾害的威胁极为严重,因此很有必要针对该地区的地震危险性评估进行深入的探索性研究。此外,华北地区作为我国历史地震记录完整性较好,并且现代地震台网的监测能力也较强的地区,为本研究提供了可靠的资料基础。故笔者选择华北地震区作为研究对象,以地震目录为主要输入数据,并划分地震构造背景区,建立断层导向性的椭圆光滑的地震活动模型,计算各网格点的地震发生率。选取合适的地震动衰减关系,利用网格源的地震危险性概率评价方法,尝试给出该地区地震区划结果。

1 总体技术思路

首先对研究区域划分网格单元(20km×20km),并根据地震空间分布和地震构造空间格局等信息建立地震构造背景区,确定各背景区内的地震活动性参数,如:主要断层性质、断层方位角及其权重、b值、震级上下限等。接着以地震构造背景区模型为基础,利用该区内主导断裂的性质、断层方位角及震级-破裂长度关系等,建立断层导向性的椭圆光滑模型,对区域网格内直接统计得到的地震发生率进行处理,得到光滑后的网格点地震年平均发生率。最后,对比选择符合中强地震的地震动特征的衰减关系,计算该区域地震危险性结果(PGA)。

2 地震活动性参数光滑方法及危险性计算
2.1 断层导向性的椭圆空间光滑

根据已有的研究(Lapajne等,2003;王健,2001;张力方等,2008;杨勇等,2008),小震活动密集区和中强震原地复发都受到其相关断裂的影响,即便是弥散地震活动性也与地震构造具有一定的关联性,地震的优势分布区往往沿断层展布的方向。虽然目前尚无法准确掌握两者的对应关系,但可以用一种简化模型——椭圆光滑模型,近似反应这种特征。利用回归分析得到的震级与破裂尺度之间的经验关系,统计得到地震构造单元内影响地震活动的主要断层类型及其走向等参数,确立椭圆光滑模型的相关参数。对网格内直接统计得到的地震发生率 ni(m0)采用椭圆空间光滑函数进行光滑处理,光滑增量值标记为ni(m0)(Lapajne等,2003)。该函数在一定程度上考虑了地震构造对发震的影响,又较符合地震发生的随机性,同时通过平滑处理抵消了地震定位误差。

(1)
n~i(m0)= ni(m0)•e-1/2δTijVTRVδij
j
∑e-1/2δTijVTRVδij
j
(2)
R
|
|
1
σ2
0
|
|
01
τ2
(3)
V
cosα -sinα
sinα  cosα


式中,δij为第i个网格到第j个网格的距离;T为转置运算符;σ为椭圆长半轴,τ为椭圆短半轴,相当于圆形光滑模型的相关距离c,椭圆的长轴方向σ为相关距离,短轴方向τ为相关距离,并且定义断层与正东方向的逆时针夹角为断层的方位角α

根据上述地震构造量化模型和华北地区震级-破裂长度关系:lg(L)=abmu(龙峰等,2006),椭圆模型的长半轴σ和短半轴τ与破裂长度成比例关系,即:σkLk≥1,τωL,ω<1。这一方法明确了光滑中相关距离、平滑半径选取的物理意义,并将断层因素引进光滑函数中,使光滑后的地震空间分布更趋合理,物理意义也更加明确。根据每个网格节点所处的地震构造区,选取相应的构造模型参数(主要断层类型、断层权重、断层走向αστ等几何参数),依据式(1)—式(3)在之前进行的圆形空间光滑的基础上,进一步对各震级档的年发生率进行处理,最后得到起算震级的年发生率。

利用相对较为完整的输入地震目录,在网格内统计大于起始震级m0的地震个数ni(m0),并采用该网格所处的构造背景区的椭圆光滑模型,对其进行空间光滑处理得到ni(m0),进而可以由其统计时段,得到大于起始震级m0的地震发生率V(m0)。由于在地震危险性计算中一般只考虑4级以上地震的影响,故将由不同时段、不同起始震级m0的输入目录光滑后,得到的发生率V(m0)标准化为V4结果,再根据双截断G-R震级分布关系,得到各震级档的发生率。

2.2 地震震级分布

在本文的网格源地震危险性概率评价方法中,两个基本假设为:地震事件服从泊松分布和震级指数分布,即地震事件在空间和时间上随机独立发生。为此,作者对所涉及的地震目录删除前、余震,再进行统计分析得到地震活动性参数(bv等),其中地震频度(或发生率)是概率方法的重要数据。最常用的关系是G-R复发公式,它是地震频度的对数与震级规模的线性关系。通常所说的单截断的G-R关系,即限定震级下限为m0,则大于m0的复发关系为:

n(m)=n(m0)•10b(mm0)     mm0
(4)

根据Cornell等(1969)提出的双截断的地震复发关系,则m0mmu的复发关系为:

(5)
n(m)=n(m0)• 10b(mm0)-10b(mum0)
1-10b(mum0)
   m0mmu

式中,n(m0)=10αb m0

本文采用双截断的复发关系,分别利用直接统计得到的m0m0iΔm……~mu各震级档的发生率,求解起算震级的年发生率。

2.3 地震危险性计算

根据分段泊松分布模型和全概率公式,计算场点处地震动参数值u超过给定地震动参数值u0的年发生率为(Reiter,1990):

(6)
λ(uu0)= ni(m0)
i
 mum0  P(uu0|m,r)pi(m)pi(r)drdm

式中,m0为起算震级;ni(m0)为第i个潜源区mm0的地震年发生率;pi(m)为震级的概率密度函数,即pi(m)= bln10•10b(mm0)
1-10b(mum0)
pi(r)为场点到潜源之间距离的概率密度函数;P(uu0|m,r)为在距离场点处r,震级为m的地震产生的地震动值u超过给定值u0的概率,P(uu0|m,r)中的u值是根据衰减关系计算得到的。

在研究中作者将研究区划分成20km×20km的网格,统计了每个网格中的地震发生率,因此上式可简写为:

(7)
λ(uu0)= ni(m0)
i
 mum0 P(uu0|m,ri)pi(m)dm

式中,ri为空间格点到计算场点的距离。

基于以上点源式的概率地震危险性评价方法,可得到网格节点的地震动参数。

3 实例分析
3.1 建立地震构造背景区模型及参数确定

由于地震活动在空间上具有差异性,地震活动性参数(如地震发生率、b值、震级上限等)在整个研究区域内并不是相同的。应当根据地震活动特征的不同,划分不同的地震构造背景区,统计各区内的地震活动性参数。本文将根据区域构造应力场、活断层分布、地震空间分布特征综合考虑划分的地震构造背景区,然后根据高斯空间光滑法将地震事件分配到空间格点中,以计算地震年发生率,最后计算地震危险性(PGA)。

众多的区域构造研究(沈建文等,1990;闻学泽,2001;龙锋等,2006)认为,我国华北地区以走滑断层占优势,并且根据华北地区的构造应力场信息、震源机制解以及GPS观测数据(徐菊生等,1999;唐方头等,2005)认为,华北地区现今构造应力场模拟的基本特征是优势压应力方向为北东东,同时受北东-南西向挤压和南东-北西向拉张作用。在如此格局下,同时参考《中国地震动参数区划图(2001)》采用的地震区带划分方案,作者给出了9个地震构造单元,亦可称为地震构造背景区,如图1所示。

图 1 研究区域的地震构造单元划分 Fig. 1Division of seismo-tectonic units in study area
3.2 地震构造单元内参数统计

各个统计区内存在明显的构造差异,这种差异很可能导致不同的地震活动水平,表现为地震活动性参数(M0Muvb)的不同,体现在地震空间分布的不均一性。

关于震级-破裂长度关系,国内外很多学者(邓起东等,1992;Wells等,1994;沈建文等,1990;闻学泽,2001;龙锋等,2006)已做过大量研究。根据翟文杰(1995)对华北地区断层性质与各级地震关系的统计,不同的地震断层,其发震规模和强度是不同的,较大的地震大部分发生在走滑断层上。龙峰等(2006)系统收集和整理了1965年以来华北地区发生的34次4.0—8.0级地震,共判定出34次地震的可靠破裂尺度参数。基于这些参数建立起来的华北地区以走滑为主的地震活断层的震级-震源破裂长度经验关系式为:Ms=3.818+1.895lg(L),以此作为区域内断层的震级-破裂长度的回归关系。关于椭圆模型中kω值的确定,根据其统计数据发现,破裂宽度比破裂长度的中值为0.36,在华北地震区椭圆模型中ω值取0.36、k取1时,比较合理。依据此回归关系式进行断层导向的椭圆光滑模型的建立。

如何估计每个地震构造单元内主导断层走向和其对应的权重,综合考虑该地区活动构造的空间展布和不同方向断层的发震能力。具体确定原则与潜在震源区方向性函数的取值原则类似,这里不再赘述。各地震构造背景区中具体参数如表1中所示。

表1 地震构造单元内地震活动性参数 Table 1 Seismicity parameter in each seismo-tectonic unit
3.3 地震发生率

不同的地震目录反应出不同的构造信息,尤其是在本文的方法中地震目录是重要的输入数据之一。根据不同的地震目录本文建立了以下两种计算模型,而两种地震目录都是经过删除前余震处理过的。

模型一(M1):1484—2006年Ms≥4.7级的破坏性地震事件,其中包括历史文献记载的震级Ms≥4¾的破坏性地震事件和地震台网资料中的震级Ms≥4.7的破坏性地震事件。

模型二(M2):1970—2006年Ms≥3.0级的现代仪器记录的现代中小震资料。假设1970年以来的Ms≥3.0级中小地震资料能够反映未来中等强度破坏性地震的发生区域及其年平均发生率。计算研究区域50年超越概率10%地震动峰值加速度结果,并确定两种模型的权系数,加权叠加后得到该地区综合地震区划结果。

图2为未经空间光滑处理的Ms≥4.0级地震的年发生率V4(图2a)和经空间光滑处理的结果V4(图2b)的比较。从图中可以看出,直接统计得到的地震发生率只存在于曾经有地震发生的格网内,而在未发生过历史地震的格网发生率为零,显然这种直接统计结果不符合地震原地复发的实际情况。当经过空间光滑处理后,在保持地震个数不变的情况下,将地震发生率按高斯分布特征分配到周围的格网内,即使地震发生率的峰值和低谷被平均到每个格网单元内,其光滑后的地震发生率分布和该地区的地震构造格局非常一致。

图 2a 由计算模型一格网内直接统计得到的V4结果分布 Fig. 2aDistribution of M≥4 earthquake occurrence result from directly counted in grid cells
图 2b 由计算模型一经过光滑处理得到的V4结果分布 Fig. 2bDistribution of M≥4 earthquake occurrence result from spatially smoothed process
3.4 衰减关系选取

本文主要关注小于7级地震带来的潜在危险性,为此对比分析了国内外根据小于7级地震记录得到的地震动衰减关系。目前,国内外研究较多的地震动衰减关系主要针对大于6级强震动衰减模型,这样往往会低估6级以下中强地震的影响,对此高玉峰等(2000)、王海江(2002)、李小军等(2005)等进行了中小地震动衰减关系的研究,取得了不错的研究结果,为本文的衰减关系的选取提供了一定的参考。

由于本文所选用的地震目录都是小于7.0级的地震,故在加速度、反应谱、持时的回归中不再考虑地震动的近场饱和,即回归的模型选为:lgYc1c2Mc3lgRεlgY 。

王海江(2002)收集了美国西部加州地区中小地震的近场记录(4.0≤M≤6.5,震源距R<70km),采用上述模型将水平方向上的两个分量视为独立的观测值;距离为震中距,Y为峰值加速度(cm/s2);回归结果为:水平向c0=0.4678,c1=0.4709,c3=-0.9807, εlgY=0.29。金星等(2004;2008)利用福建省地震监测台网拥有大量中小地震的速度记录,将速度记录进行仿真处理得到加速度记录,进而利用不同的模型对数据进行统计分析,得到了福建地区中小地震峰值加速度和峰值速度的衰减规律,并用2007年8月29日永春M4.6级地震对该规律进行了验证,说明该关系的合理性。

根据输入地震目录要与地震动衰减关系相匹配的原则,历史地震危险性评价模型(M1)采用王海江(2002)的地震动参数衰减关系,而现代中小地震危险性评价模型(M2)采用金星等(2004;2008)的地震动参数衰减关系。

3.5 地震危险性结果

基于上述网格源的概率地震危险性评价方法,本文得到了9000多个网格节点的水平向峰值加速度结果,并进行了空间插值生成加速度等值线,同时依据表2的分档原则,编制了研究区域的峰值加速度区划图,如图3a、b、c所示,它们分别为模型M1、模型M2和两种模型加权的结果。结果显示,峰值加速度的分布特征与该地区地震活动性分布非常一致,不同模型得到的结果反映了不同地震带来的危险性。在图3a、b、c中的地震空间分布均呈现出华北地区典型的“π”字形地震高危险区分布特征,对应的高峰值加速度主要集中在大同、忻州、太原和临汾盆地内部,其集中了本区大部分0.2g档以上的危险区。这与该地区总体的地震构造相符。

表2 地震动峰值加速度分档值 Table 2 Peak acceleration bracket

由于不同时段、不同震级的地震目录在完整性上的差异和本身所蕴含地震信息的不同,因此得到的区划结果也各有不同,两者是相辅相成的。历史地震模型(M1)结果主要反应了中强地震的影响,其优点是地震记录时间较长,基本能涵盖强震的复发周期;缺点是在某些地区地震目录漏记较严重,如在东部海域地区,造成其危险性评估值偏低。而现代中小地震模型(M2)结果则代表了中短期内的地震危险性,由于其地震目录的可靠性和完整性较高,使存在历史地震漏记地区的地震危险性得到了补充性提高,尤其是在江苏东部和山东半岛北部海域出现两个0.2g区,这是在M1结果中所没有的;但不足的是仅40多年的数字台网记录时间相对于强震的复发周期还太短,但随着观测时间的增加会有所改善。0.1g、0.15g、0.2g和0.3g区明显比M1的分布更广,分别是M1相应面积的2.1、2.6、2.5和3.1倍,对以往低估的中等地震危险区有了显著提高,体现了利用现代中小地震资料进行危险性的优势所在。

如在图3c中,加权结果体现了该地区综合的地震危险性,在进行综合危险性评估时,则根据每个模型中地震资料选取时间段的不同和资料可靠性的差异,对两种模型进行了适当的加权平均,即综合模型MT=0.6M1+0.4M2。关于综合模型中的权重分配,是一个比较复杂的问题,需要建立在对发震构造的活动特征、活动强度、历史地震和古地震等进行全面了解和分析的基础上,目前的取值还有待商榷,这也是今后需进一步研究的问题。图4为北京、太原和南京三地的概率曲线,上述三地的地震危险性分布是:太原最高;北京次之;南京最低。

(a)M1;

(b)M2

(c)加权综合结果

图 3 两种计算模型的结果及两者加权叠加的地表峰值加速度分布(50年超越概率10%) Fig. 3Peak acceleration (gal) with 10% probability of exceedance in 50 years
图 4 太原、北京、南京三地的概率曲线 Fig. 4The hazard curve in Taiyuan, Beijing and Nanjing area

以上结果均为按照20km×20km网格划分得到的,为验证网格大小对地震危险性结果的影响,还采用了10km、15km、30km、50km的网格划分方案,发现当网格大小在≤20km时较为适宜,结果差别很小且结果较稳定。图5a、b分别为首都圈地区当网格大小10km和20km的平均场地的峰值加速度结果。从两图的比较中可发现,当计算网格加密时得到的危险性分布图像在保持总体格局不变的情况下,加强了细节信息的显示,图像更加精细,但划分的网格越小计算量会成指数激增。当网格大小为30km和50km时,相对于网格20km的结果变化较大,可能是由于网格大小过大,造成光滑失效。故一般情况下,网格的大小以不大于地震最大定位误差的三分之一为宜。

(a)网格大小为10km;

(b)网格大小为20km

图 5 首都圈地区根据不同网格划分得到的50年超越概率10%的地表峰值加速度分布结果 Fig. 5Peak acceleration (gal) with 10% probability of exceedance in 50 years in the capital area
4 结论与讨论

考虑地震构造影响的空间光滑地震活动性模型,不仅充分体现了地震活动的空间非均一性,而且也可使构造信息体现在地震危险性结果中。尤其是当评估弥散性中强地震所带来的地震危险性时,避免了划分潜在震源区带来的较大不确定性。对网格点地震发生率的空间光滑过程中,在地震数目不变的情况下,可以有效地抵消地震定位误差。在构造导向性的椭圆光滑处理过程中,用简单易行的方式考虑了地震构造对未来地震危险性的影响作用。

考虑到不同时段、不同震级的地震目录在完整性上的差异和本身所蕴含地震信息的不同,本文建立了以历史强震和现代仪器地震记录为输入数据的计算模型,以此得到了代表不同意义的地震动区划,两者是相辅相成的关系。同时加权的地震动区划结果代表了综合地震危险水平,抵消了系统的不确定性。关于综合模型中的权重分配,需要建立在对发震构造的活动特征、活动强度、历史地震和古地震等多方面的基础上,才能进行全面了解和分析。总之,本文的方法只是对目前地震区划工作中新方法的探索性尝试,还有很多方面需要改进和细化。

参考文献
[1]邓起东,于贵华,叶文华,1992.地震地表破裂参数与震级关系的研究◆见:国家地震局地质研究所编,活动断裂研究(2). 北京:地震出版社,247—264[本文引用:1次]
[2]高孟潭,卢寿德,2006.关于下一代地震区划图编制原则与关键技术的初步探讨. 震灾防御技术,1(1):1—6[本文引用:1次]
[3]高孟潭,2002.中国地震动参数区划图(GB 18306-2001)主要特色与地震区划新动向◆见:中国地震局地球物理研究所等,新世纪地震工程与防震减灾. 北京:地震出版社,173—178[本文引用:1次]
[4]高玉峰,谢康和,2000.中强地震区地震烈度和峰值加速度的衰减规律. 浙江大学学报(自然科学版),34(4):404—408[本文引用:1次]
[5]胡聿贤,1999.地震安全性评价技术教程. 北京:地震出版社[本文引用:2次]
[6]金星,马强,李山有,2004.利用数字强震仪记录实时仿真地动速度. 地震学报工程与工程震动,24(1):49—84[本文引用:2次]
[7]金星,康兰池,2008.福建地区中小地震地震动峰值衰减规律研究. 地震学报,30(3):279—291[本文引用:2次]
[8]李小军,阎秀杰,潘华,2005.中小震近场地震动估计中地震动衰减关系的适用性分析. 地震工程与工程振动,25(1):1—7[本文引用:1次]
[9]龙锋,闻学泽,徐锡伟,2006.华北地区地震活断层的震级—破裂长度、破裂面积的经验关系. 地震地质,28(4):511—535[本文引用:4次]
[10]沈建文,邱瑛,赵志贺,1990.震级-破裂长度关系与断层破裂模型. 地球物理学报,33(2):242—248[本文引用:2次]
[11]唐方头,张培震,邓志辉等,2005.华北地区主要断裂带的现今运动特征. 煤田地质与勘探,33(1):4—6[本文引用:1次]
[12]王健,2001.地震活动性图象处理的网格点密集值计算方法. 地震学报,23(3):262—267[本文引用:1次]
[13]王海江,2000.中小地震地震动衰减关系的研究. 北京:中国地震局地球物理研究所[本文引用:3次]
[14]闻学泽,2001.活动断裂的可变破裂尺度地震行为与级联破裂模式的应用. 地震学报,23(4):380—390[本文引用:2次]
[15]徐菊生,袁金荣,高士钧,1999.利用GPS观测结果研究华北地区现今构造应力场. 地壳形变与地震,19(2):81—89[本文引用:1次]
[16]杨勇,史保平,孙亮,2008.基于华北区域地震活动性分布的地震危险性评价模型. 地震学报,30(2):198—208[本文引用:1次]
[17]翟文杰,1995.地震活动性参数的确定. 东北地震研究,11(2):14—22[本文引用:1次]
[18]张力方,吕悦军,彭艳菊,2008.用空间光滑方法评估弱地震活动区的地震活动性参数. 震害防御技术,3(1):27—36[本文引用:1次]
[19]Cornell C.A. and Vanmarcke E.H., 1969.The major influences onseismic risk. In: Proc. of the Fourth World Conference on Earthquake Engineering, January 1969, Santiago, Chile, A-1, 69—93※ [本文引用:1次]
[20]Frankel, 1995.Mapping seismic hazard in the central and eastern United States.Seismological Research Letters,66 (4): 8—21[本文引用:1次]
[21]Lapajne J., Šket Motnikar, Zupančič P., 2003.Probabilistic seismic hazard assessment methodology for distributed seismicity.BSSA,93 (6): 2502—2515[本文引用:3次]
[22]Reiter L., 1990.Earthquake Hazard Analysis:Issues and Insights.Columbia Press, New York,P254[本文引用:1次]
[23]Wells D.L. and Coppersmith K. J., 1994.New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement.BSSA,84: 974—1002[本文引用:1次]