北京时间4月14日7时49分,青海玉树发生Ms7.1级地震.该次地震对灾区造成巨大的经济损失和人员伤亡.玉树地震灾区建筑物多为土木结构,震后玉树县城土木结构房屋几乎全部倒塌;玉树州所在地结古镇是受灾最为严重的区域,位于结古镇西南部的扎西大通村几乎被夷为平地.截止4月25日下午17时,玉树地震已造成2220人遇难,70人失踪,9145人受伤.
大量的震害资料证实,地震动是造成建(构)筑物破坏的主要动力,也是地震时造成一系列地质灾害(例如砂土液化、地面破裂和塌陷、山体崩塌、滑坡等)的主要诱因.在建(构)筑物抗震设防,尤其是大型建(构)筑物的非线性分析中,地震动时程是必不可少的数据.对于缺乏强震记录的国家或地区,往往需要运用多学科的综合信息预测地震的地震动时程[1].玉树地震中,震中方圆350km以内区域没有架设强震动台网,缺乏近场强震动记录.预测该次地震的近断层加速度时程和加速度场,可为灾区的恢复重建提供抗震设防的依据.
地震动是由3个物理过程(地震破裂过程,波在地壳介质中的传播过程,场地反应)组成的一种复杂系统的产物,地震动模拟与预测均是围绕这3个物理过程的建模开展的.目前,地震动模拟与预测的方法主要有3种:确定性方法,随机性方法和混合方法[1, 2].确定性方法基于Aki和Richard[3]的表示定理,将地震动表达为震源时间函数和格林函数的卷积,用于长周期地震动的模拟.常用的方法有三维有限差分方法,离散波数法和有限元方法.随机性方法基于高斯带限白噪声的随机振动理论[4, 5],或者基于小震的经验格林函数方法[6, 7],用于短周期(高频)地震动的模拟.前者对于远源和小震,一般使用Boore的随机点源模型模拟地震动[4, 5];对于近源或近断层,地震动的空间分布强烈地受到发震断层方位(位置、产状、上界埋深)、断层破裂面上滑动分布的不均匀性和破裂过程的影响,需要使用有限断层震源模型[8~10]来考虑震源破裂过程对近源或近断层地震动的影响.经验格林函数方法将目标地震区记录的、适当大小的、与主震震源机制相同的、信噪比高的小震记录作为经验格林函数模拟地震动[7].该方法的优势是小震记录中包含了波在地壳中传播的路径效应和场地效应;不足之处是适合于地震动模拟的小震记录不易获取,尤其是在地震记录稀疏或缺乏的地区.混合方法利用上述两种方法各自的优势,将地震动分成长、短周期两部分分别进行模拟,长周期地震动用确定性方法模拟,短周期(高频)地震动用随机方法模拟,模拟结果分别经低、高通滤波器滤波之后,在时域中叠加,得到宽频带的地震动[11~13].
本文基于有限断层震源、且使用动力学拐角频率的地震动随机模拟方法预测该次地震近断层191个节点的加速度时程,然后取每个结点的加速度峰值绘制该次地震的近断层加速度场.
2 方法 2.1 有限断层震源建模根据Irikura[14],将震源模型参数分为两类:全局震源参数和局部震源参数.前者表征震源区的宏观特征,例如发震断层的空间方位(位置、产状、上界埋深)、断层类型、破裂尺度(长度、宽度、面积)、断层面上的平均滑动和平均破裂速度等;后者表示断层面上的不均匀性或粗糙度,主要用滑动分布来表征.这两类参数的估计方法如下.
全局震源参数:(1)根据地震地质和地震活动性调查以及地球物理勘探等资料确定活断层的空间方位和滑动类型;(2)根据活断层地震危险性评价的结果,获得活断层设定地震的矩震级,根据发震断层破裂尺度及断层面上平均滑动与矩震级之间的地震定标律估计相应的参数值[15, 16];(3)根据震后几天内余震空间分布的特征估计发震断层的空间方位、断层产状以及断层破裂尺度.
局部震源参数:(1)凹凸体(asperity)模型,根据凹凸体参数与矩震级或相应断层参数之间的经验关系式确定凹凸体模型的参数---凹凸体的位置、大小及其之上的平均滑动量[16, 17];(2)k平方滑动模型,根据空间拐角波数与矩震级之间的经验关系式确定相应的空间拐角波数值[16].当波数大于空间拐角波数时,断层面上滑动分布的空间波数谱以k-2衰减.用k平方滑动模型生成随机滑动分布;(3)凹凸体模型与k平方滑动模型叠加在一起生成断层面上的混合滑动分布[2, 10, 16].
2.2 近断层地震动模拟近断层强地震动模拟中,有限断层震源模型的断层面被分成N个大小相等的矩形子断层,每个子断层即为一个点源,亦称子源.破裂过程以一定的破裂速度(一般取0.8倍的剪切波速)从破裂起始点开始呈辐射状向外传播,传播到每个子源的中心时该子源即被触发.每个子源引起的地震动由Boore[4, 5]的点源模型计算.所有子源在观测点引起的地震动在时域中以适当的延迟时间叠加,可获得观测点的地震动时程a(t),有
(1) |
其中NL和NW分别是沿着断层走向和下倾方向的子断层数,NL ×NW=N为子源总数,Δtij包括破裂传播到第ij个子源引起的时间滞后和从第ij个子源至场点间由于传播距离的不同引起的时间滞后.aij(t)是第ij个子源引起的观测点地震动.
每个子源的地震矩由下式计算:
(2) |
其中M0是地震矩,单位是dyne·cm(1dyne=10-5N),Dij是第ij个子断层上的平均滑动,单位是cm.
合成地震动时,震源谱中的拐角频率如果使用Brune[18]的静力学拐角频率,则对于子断层的大小有严格的要求.而且,为了保证地震矩的守恒,每个子源需要触发多次.
Beresnev和Atkinson[19]确定的子断层大小为
(3) |
如果改变子断层的大小,往往导致断层辐射的总辐射能不守恒.而且合成大震时要求子震的震级在5.0~6.5的范围之内,限制了小震地震动的模拟.
为了避免对每个子源的多次触发以及对子断层及矩震级大小的限制,保证地震矩和辐射能的守恒,Motazedian和Atkinson[9]提出动力学拐角频率的概念:
(4) |
式中foij(t)是第ij个子断层的动力学拐角频率,t是第ij个子源被触发的时刻,NR(t)是在时刻t已破裂子断层的累积数,β是震源附近的剪切波速度,单位为km/s,Δσ是应力降,单位为bar(1bar=105Pa),Moave=Mo/N是子断层的平均地震矩.
第ij个子断层的加速度谱Aij(f)是
(5) |
其中
(6) |
Moij是第ij个子断层的地震矩,f是频率,fc是静力学拐角频率.
(7) |
玉树地震发生在巴颜喀拉活动地块南界甘孜-玉树-凤火山断裂带的玉树部分,该断裂带的性质为左旋走滑型.国内外研究机构确定玉树地震的断层类型均为左旋走滑断层,节面Ⅰ参数为该次地震断层的产状和滑动方向(表 1).但各机构确定的震中位置、地震矩、震源深度、断层倾角、滑动角有所不同.中国地震局将震中位置更新为经度96.6°E,纬度33.2°N,将震源深度更新为14km.本文取哈佛大学确定的节面Ⅰ为该次地震的震源机制,即走向120°,倾角90°,滑动角-13°.震源深度取中国地震局确定的14km.矩震级取Mw6.9.
根据6.5 < Mw≤7.0时走滑断层宏观震源参数的半经验关系式[15, 16]
破裂面积:
(8) |
破裂长度:
(9) |
破裂宽度:
(10) |
平均滑动:
(11) |
求得Mw6.9级地震的宏观震源参数并将破裂长度和宽度取整数得:
(12) |
根据6.5 < Mw≤7.0时走滑地震断层面上滑动分布参数的经验关系式[15, 16]
所有凹凸体的面积:
(13) |
最大凹凸体的面积:
(14) |
第二凹凸体的面积:
(15) |
最大凹凸体上的平均滑动:
(16) |
第二凹凸体上的平均滑动:
(17) |
求得凹凸体模型参数如表 2所示.
将凹凸体模型与k平方模型相结合,获得如图 1a所示的玉树地震的滑动分布模型.图 1(b,c)分别是王卫民等[20]和陈运泰研究小组[21]基于全球台网数据反演的玉树地震的滑动分布模型,二者相差较大.
波在地壳中的传播模型包括几何衰减函数和滞弹性衰减函数,近断层地震动模拟中,几何衰减函数一般采用1/R.根据华卫等[22],龙门山地区品质因子与频率的关系式为:
(18) |
玉树震区与龙门山地区均属于巴颜喀拉地块,有相似的滞弹性衰减关系,故本文也使用该式.
地震动模拟中,路径持时采用Atkinson和Boore[23]对于北美东部提出的经验关系式:
(19) |
其中R为观测点到断层破裂面的最近距离.
场地反应包括地表地层的放大作用和衰减作用两部分.玉树灾区属于高海拔地区,地形以山地为主,平均海拔4493 m.玉树地震近断层岩石类型主要为三叠纪砂岩和板岩互层,由于缺乏该区域场地的剪切波速资料,地震动预测中放大模型采用Boore和Joyner[24]的坚硬岩石场地模型;表征地表衰减作用的kappa参数取0.005.
3.3 加速度场预测基于上述分析,玉树地震近断层加速度时程预测中输入的基本参数如表 3.图 2是用于玉树地震近断层加速度场预测的节点位置和编号,在该区域内,共设计了191个节点,并在发震断层周围加密了节点,以更好地体现断层附近加速度峰值的变化.在此基础上,基于上述地震动模拟方法预测了玉树地震近断层191个节点的加速度时程.图 3是分别沿着图 2中A-A'和B-B'剖面的加速度时程.图 4是取各节点的加速度峰值绘制的玉树地震的近断层加速度场.从图中可见:(1)在直立的断层两盘,近断层加速度场对称分布,且在垂直断层方向,加速度峰值衰减很快;(2)沿着断层走向有2个加速度峰值超过800cm/s2的区域,第一个位于堪扣的东南方向附近,最大加速度峰值达1111.65cm/s2(节点112);第二个位于玉树县结古镇和禅古寺之间,最大加速度峰值为936.94cm/s2(节点65).这2个区域均分布在断层线上及其附近,而且与断层面上2个凹凸体投影到地表的区域一致,说明断层面上凹凸体投影到地表的区域附近,是加速度峰值最大的区域,也是建(构)筑物震害最严重的区域.
图 5是基于王卫民等[20]确定的震源机制参数(表 1中中国科学院青藏高原研究所)和反演的滑动分布模型(图 1b)预测的玉树地震近断层加速度场.预测191个节点的加速度时程时,输入参数除断层方位(走向116°,倾角72.8°)、断层沿走向和下倾方向的尺度(72.8km×23.8km)、子断层大小(5.2km× 3.4km)、滑动分布(图 1b)、断层左上角位置外,其他输入参数均与表 3中的相同.从图中可见,沿着发震断层也有2个加速度峰值较大的区域,第一个区域也位于堪扣的东南方向附近,最大加速度峰值达1786.76cm/s2(节点112),高于本文滑动分布模型预测的最大加速度峰值.该区域与本文建立的滑动分布模型预测的最大加速度峰值区域位置一致,但前者范围更大一些,说明王卫民等反演的滑动分布模型中,最大凹凸体位置与本文建立的滑动分布模型一致,但面积比本文的更大.第二个区域位于禅古寺东南方向,加速度峰值在600cm/s2以上的区域很小,说明王卫民等反演的滑动分布模型中,第二凹凸体位置与本文建立的滑动分布模型不一致,而且面积比本文小的多.
需要说明的是,基于本文建立的震源模型和王卫民等[20]反演的震源模型分别预测玉树地震的近断层加速度场时,二者在震源机制参数方面的差异主要是发震断层产状(走向和倾角)的不同,其他参数均取一致.本文建立的震源模型中断层产状采用哈佛大学的结果,即断层走向和倾角分别是120°和90°;王卫民等反演的震源模型中,断层走向和倾角分别是116°和72.8°.二者断层走向相差4°,结果是加速度峰值最大区域的位置发生4°的偏转;二者断层倾角的差异,使得预测的近断层加速度场中断层附近相同峰值加速度等值线的范围后者大于本文.也就是说,随着断层倾角的变缓,断层破裂过程对近断层地震动的影响范围将增大.
在近断层产生地震动的3个物理过程中,震源破裂过程及其滑动分布对地震动的强度和空间变化影响最大,震源模型对于近断层地震动的模拟与预测至关重要.由于震源破裂过程的复杂性,在近断层强地震动模拟与预测中如何建立实用的震源模型一直是世界性的难题,也是地震学家研究的一个重点,进展相当艰难.本文建立的滑动分布模型,是根据大量浅源地震的滑动分布特征提取出来的地震的平均特征,能够反映地震的共同特点[1, 10, 16].基于全球台网长周期地震动波形反演的震源模型,从理论上来说,能够真实地描述震源破裂过程及断层破裂面上滑动分布的大致轮廓;由于反演过程中影响因素较多,对于相同的地震,不同的研究人员即使采用相同的全球台网数据,反演结果也不尽相同,甚至相差甚远.这也是图 1(b,c)相差较大的原因.由于玉树地震灾区缺乏强震动记录,无法验证哪一个震源模型更接近于真实.玉树州、县所在地及其附近是人口密集区,也是这次地震震害最严重的区域,地震烈度达到Ⅸ度.以本文建立的震源模型预测的玉树地震近断层加速度场中第二个加速度峰值高值区域恰好位于该人口密集区,说明该模型是合理的.
由于场地资料的缺乏,在近断层加速度时程预测中将玉树地区的场地均假设为坚硬岩石场地,没有考虑地形和土层场地的放大效应,某些部位的加速度峰值可能偏低.尤其是在发震断层的西南方向,有一个较大的第四纪盆地---玉树南盆地,沉积物厚度不详(图 2).在该盆地以坚硬岩石场地预测的加速度峰值均超过了100cm/s2,但小于300cm/s2;如果考虑土层的放大作用,预测的加速度峰值会更高.根据Beresnev等[26],土层发生非线性反应的阈值为100cm/s2,那么该盆地内的土层会发生非线性反应,其结果是大幅降低土层在高频的放大效应.而在发震断层的东北方向,即四道班周围有一个较小的第四纪盆地(图 2),在该盆地内以坚硬岩石场地预测的加速度峰值较小(小于50cm/s2),本文预测的峰值加速度值可能严重偏低.再者,河道附近沉积的第四纪沉积物也会对地震动放大,预测的玉树地区河道附近的加速度峰值也应该偏低.
长期以来,一直认为:地震时,当破裂传播方向和断层面上的滑动方向均指向一个场地,该场地将发生破裂方向性效应.也就是沿着断层破裂传播方向,地震动振幅增加,逆破裂传播方向,地震动振幅减小.玉树地震是沿着发震断层从西北到东南破裂的单侧破裂,按照上述条件,从破裂开始点开始,沿着断层破裂传播方向断层沿线附近的场地均应该发生破裂方向性效应.从西北到东南方向,地震动峰值应该逐渐增大.但预测结果并非如此(图 3中B-B'剖面).从图 4,5中可见,加速度峰值主要受断层破裂面上凹凸体分布的控制,凹凸体投影到地表的区域附近,也是加速度峰值最大的区域.对于均匀破裂的断层而言,从破裂开始点沿着断层破裂方向会出现地震动振幅逐渐增加的破裂方向性效应,但由于断层破裂过程的复杂性,均匀破裂是很难出现的.但从破裂开始点开始,沿着断层破裂方向局部出现破裂方向性效应是必然的,例如图 4中2个加速度峰值较大区域,沿着破裂方向,从边缘到中心加速度峰值是逐渐增加的.图 3中B-B'剖面中节点136到113是图 4中加速度峰值最大区域中的破裂方向性效应.
5 结论地震动场是震源破裂过程、地震波在地壳介质中的传播过程以及场地反应3个物理过程组成的一种复杂系统的产物,能够客观地反映地震引起的地面振动的强弱.其确定方法有2种:一是在地震活动性较大的区域布设台站密度较大的强震台网,地震时,根据强震动观测数据确定近断层地震动场,但需要投入大量的人力物力;二是根据近场地震学理论模拟与预测近断层地震动场.由于强震动观测数据的严重不足,一次地震过后,强震动观测数据和地震动模拟结果相结合估计近断层地震动场是目前最佳的选择.
玉树地震中,震中方圆350km之内无强震观测台站,只能根据上述第二种方法预测近断层的地震动.本文基于有限断层、且使用动力学拐角频率的地震动随机模拟方法预测了该次地震近断层的加速度场,并得出如下主要结论:
(1)对于走滑型地震,在垂直断层方向,加速度峰值分别沿断层两侧迅速衰减;对于垂直或者近垂直的断层,近断层加速度场在断层两侧对称分布.
(2)近断层加速度场主要受震源破裂过程和断层面上的滑动分布的影响.断层面上凹凸体投影到地表的区域附近,加速度峰值最大,也是震害最严重的区域.
(3)对于走滑型地震,从破裂开始点开始,沿着断层破裂传播方向断层沿线附近的场地并非均会发生破裂方向性效应,发生破裂方向性效应的场地与凹凸体在断层面上的位置有关.
由于场地资料的缺乏,本文在近断层加速度时程预测中将玉树地区的场地均假设为坚硬岩石场地,没有考虑地形和土层场地的放大效应,某些部位的加速度峰值可能偏低.今后的工作中,需要进一步研究该地区地形和土层场地的放大效应以及土层的非线性反应.
致谢感谢中国科学院青藏高原研究所王卫民博士提供基于全球台网数据反演的玉树地震震源模型数据.感谢同行专家对本文的匿名评审和提出的有益问题.
[1] | 王海云, 谢礼立. 近断层地震动模拟现状. 地球科学进展 , 2008, 23(10): 1043–1049. Wang H Y, Xie L L. A review on near-fault ground motion simulation. Advance in Earth Science (in Chinese) , 2008, 23(10): 1043-1049. |
[2] | 王海云, 谢礼立. 近断层强地震动场预测. 地球物理学报 , 2009, 52(3): 703–711. Wang H Y, Xie L L. Prediction of near-fault strong ground motion field. Chinese J. Geophys. (in Chinese) , 2009, 52(3): 703-711. |
[3] | Aki K, Richards P G. Quantitative Seismology:Theory and Methods. New York:W H Freeman & Co Ltd., 1980 |
[4] | Boore D M. Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra. Bull. Seism. Soc. Am. , 1983, 73: 1865-1884. |
[5] | Boore D M. Simulation of ground motion using the stochastic method. Pure Appl. Geophys. , 2003, 160: 635-676. DOI:10.1007/PL00012553 |
[6] | Hartzell S. Earthquake aftershocks as Green's functions. Geophys. Res. Lett. , 1978, 5: 1-4. DOI:10.1029/GL005i001p00001 |
[7] | Irikura K. Semi-empirical estimation of strong ground motions during large earthquakes. Bull. Disas. Prev. Res. Inst. , 1983, 33: 63-104. |
[8] | Beresnev I, Atkinson G. FINSIM:a FORTRAN program for simulating stochastic acceleration time histories from finite faults. Seism. Res. Lett. , 1998, 69: 27-32. DOI:10.1785/gssrl.69.1.27 |
[9] | Motazedian D, Atkinson G M. Stochastic finite-fault modeling based on a dynamic corner frequency. Bull. Seism. Soc. Am. , 2005, 95(3): 995-1010. DOI:10.1785/0120030207 |
[10] | 王海云, 谢礼立, 陶夏新. 近断层强地震动预测中的有限断层震源模型. 地球科学-中国地质大学学报 , 2008, 33(6): 843–851. Wang H Y, Xie L L, Tao X X. Finite fault source model for predicting near-fault strong ground motion. Earth Science-Journal of China University of Geosciences (in Chinese) , 2008, 33(6): 843-851. DOI:10.3799/dqkx.2008.101 |
[11] | Kamae K, Irikura K, Pitarka A. A Technique for simulating strong ground motion using hybrid Green's function. Bull. Seism. Soc. Am. , 1998, 88(2): 357-367. |
[12] | Hartzell S, Harmsen S, Frankel A, Larsen S. Calculation of broadband time histories of ground motion:comparison of methods and validation using strong-ground motion from the 1994 Northridge earthquake. Bull. Seism. Soc. Am. , 1999, 89(6): 1484-1504. |
[13] | Pitarka A, Somerville P, Fukushima Y, et al. Simulation of near-fault strong-ground motion using hybrid Green's functions. Bull. Seism. Soc. Am. , 2000, 90: 566-586. DOI:10.1785/0119990108 |
[14] | Irikura K. Prediction of strong ground motions from future earthquakes caused by active faults--Case of the Osaka Basin. 12th World Conference on Earthquake Engineering (12WCEE), CDROM, No.2687. Auckland, New Zealand During 30 January-4 February 2000 |
[15] | Wang H Y, Tao X X. Relationships between moment magnitude and fault parameters:theoretical and semi-empirical relationships. Earthquake Engineering and Engineering Vibration , 2003, 2(2): 201-211. DOI:10.1007/s11803-003-0004-x |
[16] | 王海云.近场强地震动预测的有限断层震源模型[博士论文].哈尔滨:中国地震局工程力学研究所, 2004 Wang H Y. Finite fault source model for predicting near field strong ground motion[Ph.D.thesis](in Chinese). Harbin:Institute of Engineering Mechanics, China Earthquake Administration, 2004 http://www.oalib.com/references/17436412 |
[17] | 王海云, 陶夏新. 近场强地震动预测中浅源地震的Asperity模型特征. 哈尔滨工业大学学报 , 2005, 37(11): 1533–1539. Wang H Y, Tao X X. Characterizing a shallow earthquake asperity model for predicting near field strong ground motion. Journal of Harbin Institute of Technology (in Chinese) , 2005, 37(11): 1533-1539. |
[18] | Brune J N. Tectonic stress and the spectra of seismic shear waves from earthquakes. J. Geophys. Res. , 1970, 75: 4997-5009. DOI:10.1029/JB075i026p04997 |
[19] | Beresnev I, Atkinson G. Subevent structure of large earthquakes-a ground motion perspective. Geophys. Res. Lett. , 2000, 28: 53-56. |
[20] | 王卫民, 赵连锋, 姚振兴. 2010年4月14日玉树7.1级地震震源破裂过程初步结果.http://www.csi.ac.cn/manage/html/4028861611c5c2ba0111c5c558b00001/_content/10_04/14/1271250817270 Wang W M, Zhao L F, Yao Z X. A preliminary study on source rupture process of the Yushu earthquake (Ms7.1) on the 14 April 2010.http://www.csi.ac.cn/manage/html/4028861611c5c2ba0111c5c558b00001/_content/10_04/14/1271250817270 |
[21] | 陈运泰院士研究小组. 2010年4月14日青海玉树地震震源机制和破裂过程快报(第四版).http://www.csi.ac.cn/manage/html/4028861611c5c2ba0111c5c558b00001/_content/10_04/17/1271488288058 Chen Y T's Study Group. A flash report on both source mechanism and source rupture process of the Yushu earthquake, Qinghai province on the 14 April 2010(fourth edition). http://www.csi.ac.cn/manage/html/4028861611c5c2ba0111c5c558b00001/_content/10_04/17/1271488288058 |
[22] | 华卫, 陈章立, 郑斯华. 2008年汶川8.0级地震序列震源参数分段特征的研究. 地球物理学报 , 2009, 52(2): 365–371. Hua W, Chen Z L, Zheng S H. A study on segmentation characteristics of aftershock source parameters of Wenchuan M8.0 earthquake in 2008. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 365-371. |
[23] | Atkinson G M, Boore D M. Evaluation of models for earthquake source spectra in Eastern North America. Bull. Seismol. Soc. Am. , 1998, 88: 917-934. |
[24] | Boore D M, Joyner W B. Site amplifications for generic rock sites. Bull. Seism. Soc. Am. , 1997, 87(2): 327-341. |
[25] | Saragoni G R, Hart G C. Simulation of artificial earthquakes. Earthq. Eng. Struct. Dyn. , 1974, 2: 249-267. |
[26] | Bresnev I A, Wen K L. The accuracy of soil response estimates using soil-to-rock spectral rations. Bull. Seism. Soc. Am. , 1996, 86(2): 519-523. |