2019, Vol. 62
, GAO MengTan, CHEN XueLiang
, WU Jian, ZHOU Hong, LI TieFei, LI ChangLong, WU Qing, BU ChunYao
2017年8月8日21时，四川省阿坝州九寨沟县发生MS7.0地震，地震共造成25人死亡，500多人受伤.据中国地震台网中心测定：震源位置位于(33.22°N，103.83°E)(图 1)，震源深度20 km.另据美国地质调查局(USGS)测定此次地震的地震矩M0为7.228×1018 N·m，震源机制参见网站https://earthquake.usgs.gov/earthquakes/search/.九寨沟地震发生在南北地震带中北段，青藏高原东缘巴颜喀拉块体内部，邻近该块体北边界断裂带(易桂喜等，2017；邓起东等，2010), 围绕巴颜喀拉块体边界断裂带区域已经发生过多次较大的破坏性地震，例如10年玉树MS7.1地震、2008年汶川MS8.0地震、2008年和2014年的两次MS7.3地震等(易桂喜等，2017；邓起东等，2010).
|
图 1 九寨沟MS7.0地震地形图 红色实心点表示县市位置，黄色菱形表示本文中模拟的台站.图中震源处色带表示断层走向. Fig. 1 The terrain and city distribution of Jiuzhaigou Valley MS7.0 earthquake The red solid point means counties and cities. The yellow diamond means simulated stations. The colour strip near the source is fault strike. |
众多学者对地震的震源破裂过程和震源机制、发震构造等进行了广泛的研究.郑绪君等(2017)采用IDS方法反演强震数据来确定九寨沟地震的破裂过程，表明在当前的数据条件下，即使定位结果存在不确定性，通过反演强震数据也能够独立可靠地确定中强地震的破裂过程.杨宜海等(2017)搜集了四川地震台网的波形资料，采用全波形反演九寨沟地震序列震源机制解.房立华等(2018)基于三维速度模型的定位方法, 测定了九寨沟地震的主震位置，估算了主震初始破裂点的深度不浅于14.3 km.徐锡伟等(2017)对九寨沟地震的发震断层属性及青藏高原东南缘现今应变状态进行了分析，青藏高原东缘至东南缘各块体主干边界活动断层现今处于中等偏高的应变积累状态，未来发生高震级地震的危险性不容忽视.单斌等(2017)利用弹性位错理论和分层岩石圈模型，计算了汶川地震和九寨沟地震引起的应力场变化，研究表明汶川地震导致的同震及震后应力变化显著增强了九寨沟地震震源处的应力积累，有效的促进了九寨沟地震的发生.单新建等(2017)基于InSAR技术，反演了九寨沟地震的断层滑动分布，基于三种不同接收断层计算了同震库伦应力变化.
关于九寨沟强地面运动数值模拟的研究工作相对较少.地震发生后快速获取地震动强度的时空分布特征对于震后应急救援意义重大，地震动强度预测对于未发生破坏性地震地区的建筑物抗震设防以及未来城市规划等意义非常重大，是防震减灾中一项很重要的基础工作.本文采用经验格林函数方法对九寨沟地震的强地面运动特征进行数值模拟，评估缺失强震记录的九寨章扎台站的地震动强度.因未得到九寨沟地震的余震，本文初次尝试将汶川地震和定西地震的余震作为格林函数预测九寨沟地震.本文的创新点是在缺乏余震记录的情况下用经验格林函数法对九寨沟地震的地震动特征进行预测，通过将其他大震的余震记录作为格林函数对九寨沟地震进行模拟(本文将采用定西MS6.6地震以及汶川MS8.0地震的部分余震记录作为格林函数)，本研究也为经验格林函数方法在缺乏小震记录地区的使用积累了经验.
1 方法与数据 1.1 经验格林函数方法Hartzell(1978)最早提出了经验格林函数方法，用大地震的前震或余震作为经验格林函数合成大地震.由于小震记录本身己经包含了传播介质的影响，所以用小震记录合成的大地震时程也考虑了传播介质的复杂性，并能克服计算理论格林函数的困难.
Kanamori(1979)提出了用大小地震地震矩的比值取整确定所需小震数目的方法，并针对所用的小震记录与大震不是取自同一地震和同一台站的情况提出了修正公式.Irikura等(1983, 1986)依据对大小地震震源参数的统计研究结果，提出了假定大小地震满足相似律并将大地震中的子源和所利用的小震位错函数都取为斜坡函数，则大震产生的地面运动可由对子源在断层上的二重求和转换为对小震在断层上和时间域内的三重求和.此方法放宽了对小震记录的限制并且转变了人们对“经验格林函数”法的认识，把问题从选择什么样的小震作为格林函数转换为如何合理地对小震记录进行修正以得到大震中的子源贡献.罗奇峰(1989)考虑大小地震应力降的差异，将该方法进一步推广，使该方法在大小地震断层形状不相似的情况下也可使用.Kamae等(1998)将随机方法和经验格林函数方法相联系，提出一种混合方法模拟地震动，即低频长周期部分用三维有限断层方法模拟，高频短周期部分的小震记录首先由Boore的随机方法生成，再用Irikura的经验格林函数方法模拟(用到的小震是由随机方法生成的)高频短周期地震动.
Irikura等(Irikura and Kamae, 1994; Irikura and Miyake, 2011; Irikura et al., 2017)系统提出了用经验格林函数方法预测未来地震动的想法，结合多个震例模拟了大量地震动并总结了用该方法预测未来地震动的一般步骤. Kurahashi等(2010)出于没有获得汶川地震余震记录的客观原因而用日本的余震记录成功模拟了汶川大地震的主震记录.经过众多学者的不断努力，经验格林函数方法已经发展成为一种比较完善和成熟的模拟强地震地面运动的方法，已被广大学者认可.用其他大地震的余震记录作为格林函数预测某地未来的破坏性大地震的地震动也是非常实用和有意义的研究方向，随着方法的不断完善和模拟经验的越来越丰富，必定可以为建筑物抗震设防提供实用可靠的参考依据.
本文所用的经验格林函数法采用Haskell(图 2)的震源模型并结合Haskell(1969)凹凸体的震源模型，合成地震动的核心思想是用小震记录作为格林函数合成大震地震动.将大震震源看作是由一系列子震震源构成，选择一个大小合适的余震或者前震记录作为格林函数，将小震等同于子震，按照一定的破裂方式，将这些经验格林函数叠加得到大震地震动时程.
|
图 2 Haskell震源模型 rij为第(i, j)个子断层到观测点的距离, r0为破裂初始点到观测点的距离; ξij是第(i, j)子断层到初始破裂位置的距离. Fig. 2 Haskell source model The rij is the distance between (i, j) element fault and observed point while r0 is the distance between initial rupture point and observed point. ξij is the distance between (i, j) element fault and initial rupture point. |
该方法合成主震的公式如式(1)—(4)所示(Irikura, 1986; Miyake et al., 2003)：
|
(1) |
|
(2) |
|
(3) |
|
(4) |
式(1)—式(4)中U(t)是合成的大震记录，u(t)是小震的地震记录(格林函数).rij和rs分别是子断层的震源距和小震的震源距.ξij是第(i, j)子断层到初始破裂位置的距离.T是大震的上升时间，上升时间与修正函数的持时有关.F(t)是修正小震的修正函数(Hiroe Miyake, 2003).N和C分别是大小地震断层维度的比值以及大小地震应力降的比值，M0和m0分别是大震和小震的地震矩.VS和Vr分别是S波靠近震源区的波速和断层面上的破裂传播速度.Δσlarge和Δσsmall分别是大震和小震的应力降.小震记录包含了震源的破裂过程、传播介质和场地条件复杂性的影响，所以用小震记录合成的大震记录也考虑了震源、传播介质和场地条件的复杂性，同时可以克服理论格林函数的计算困难(Hartzell, 1978).该方法在日本应用较多，Irikura等经过多次的震例反演证明了该方法在地震动数值模拟方面的实用性，这也是选择格林函数法作为数值模拟工具的原因.
1.2 地震数据本文使用的地震数据全部来自国家强震动观测台网中心，包括九寨沟地震的加速度数据，汶川地震和定西地震部分余震的加速度记录.根据国家地震科学数据共享中心提供的数据，本次地震共收集到46个台站的强震数据，地震动峰值加速度(Peak Ground Acceleration，简称PGA)最大值出现在九寨百河(051JZB)，为185.0 Gal，出现在NS方向分量.EW方向分量和UD方向分量的最大值也是在九寨百河台站，分别为129.5 Gal和124.7 Gal.此次获取的强震台站数据显示，大部分台站的PGA普遍小于10 Gal，地震动分量中超过10 Gal的台站共有9个台站.本文重点选择了PGA超过10 Gal或者比较靠近震源的10个强震台站进行了地震动数值模拟，各台站基本信息表 1所示，台站分布如图 1所示.
|
表 1 地震动模拟所选台站及备选余震基本信息 Table 1 The information about aftershocks and stations of simulation |
中国地震台网中心监测到九寨沟地震的余震接近4500个，但绝大部分余震是在MS3.0以下，MS4.0以上余震只有三个，最大余震为MS4.8.在向国家强震动观测台网中心申请九寨沟强震数据时只获得了九寨沟主震的加速度记录，未申请到九寨沟地震的余震记录，无法用九寨沟地震的余震做格林函数.因此本文初次尝试用汶川MS8.0地震和定西MS6.7地震的余震来做格林函数，即通过使用其他大震的余震作为格林函数来反演九寨沟地震动特征以及预测章扎镇地震动强度特征.由于未能找到能够被所选的10个台站同时记录到的余震信息，故本文中每个台站所选择的余震记录不一样(表 1).大部分台站的余震是MS4.0到MS5.0的小震，但在台站迭部(062DIB)和舟曲(062ZHQ)两个台站分别选的定西地震的MS6.7主震和汶川地震的MS6.4余震作为格林函数，原因是上述两个台站的地震记录只有定西地震和汶川地震的主震，但这至少也保证了余震与目标主震的场地条件的一致性.
2 关键参数的计算过程本文在预测地震动强度时需要设定的参数很多，例如震级、地震矩、断层破裂尺度、初始破裂位置等.本次地震动预测的设定震级为MS7.0，但是不同的研究机构给出的震源深度值却不尽相同.参考张旭等(2017)对九寨沟余震精定位的结果以及汇总现有对九寨沟地震的研究成果信息，本文取九寨沟地震的震源深度为20 km.
地震矩取USGS的研究成果，M0为7.228×1018 N·m.参考龙门山断裂带上已经发生的汶川地震的剪切波速，取值为VS为3.6 km·s-1，破裂速度取Vr为2.9 km·s-1.不同震级的余震地震矩、断层破裂面积等参数根据李宗超(2017)论文中对局部区域地震震源参数的经验关系进行计算(文献的图 6).震源模型采用Haskell模型和凹凸体相结合的混合震源模型，参考Somerville等(1999)，李宗超(2017)等对凹凸体的研究成果，将凹凸体的面积设定为断层破裂面积的22%，地震矩能量设定为主震地震矩的44%.
震源上升时间的计算如式(5)所示：
|
(5) |
S为主震的破裂面积，β为S波的速度，本文计算得到的震源上升时间为4.26 s.
经验格林函数法模拟地震动过程中对结果影响较大的两个参数是大小地震的应力降和地震矩(李宗超，2017).应力降的具体的表现形式为大小地震应力降的比值，用大写C表示，地震矩的表现形式为大小地震矩的比值除以C值再开立方然后取整，用大写N表示，即沿着断层走向和倾向划分子断层的个数(图 2).C值和N值自身没有具体的物理意义，但作为重要的间接参数表征大小地震应力降和地震矩之间的能量关系，对模拟结果影响非常大.经验格林函数法的基本思想是小震合成大震，需要分析大小地震之间的比例关系来确定需要叠加多少次小震来合成大震.
Kurahashi和Irikura(2010)及Xia等(2015)均用经验格林函数法对汶川地震进行了地震动模拟，并且给出了C和N值的计算方式.例如Xia等(2015)给出了详细的计算N值和C值的过程：
首先根据矩震级与地震矩的经验关系(Hanks and Kanamori, 1979):
|
(6) |
计算获得该震级的地震矩.然后计算拐角频率，通常把振幅高频谱渐进趋势(包络线)和低频趋势(零频水平)的交点称为拐角，与之相应的频率称为拐角频率.拐角频率是远场位移谱的一种估计，余震的拐角频率f根据拟合ω2模型获得，再根据公式(7)—(8)：
|
(7) |
|
(8) |
得到断层破裂面的等效半径r和余震的应力降Δσ，再根据公式(3)—(4)计算C值和N值.Susumu Kurahashi在模拟汶川地震时也是基于上述的方法计算的这两个参数.但是事实上拐角频率只是远场位移谱的一个估计，在根据位移的震源谱来确定拐角频率时存在较大的误差(图 3)，振幅高频谱渐进趋势(包络线)和低频趋势(零频水平)的交点不稳定.
|
图 3 拐角频率的取值有一定的误差范围 振幅高频谱渐进趋势(包络线)和低频趋势(零频水平)的交点所示的拐角频率位置并不明显, 不同包络线和零频水平线的交点可以不同，如图中四个不同颜色的点所示(黑色底图摘自Xia等(2015)文章中的位移震源谱，谱的拐角频率是通过拟合ω2获得). Fig. 3 The value of corner frequency has error range The corner frequency position shown at the intersection of the high frequency′s amplitude spectrum asymptotic trend (envelope line) and low frequency trend (zero frequency level) is not obvious (the black base map from Xia et al., 2015). The intersections between different envelope line and zero frequency level could be very different. |
本文对C值和N值的计算方法做了修正，主要是采取对应力降的直接计算方式，绕过了对拐角频率的依赖.地震矩M0和断层破裂面积S参考李宗超(2017)博士论文中通过统计局部区域震例获得的面波震级MS与地震矩M0、破裂面积S间的经验关系，具体计算过程如下所示：
计算过程中继续采用“布龙模型”，该模型将地震断层面等效成一个圆盘模型(Brune, 1970)，把地震看作是圆盘形断层面上剪切应力的突然释放.C值和N值的确定按照下列步骤进行计算(图 4)，各台站进行地震动模拟所用的相关参数计算列于表 2中.
|
图 4 应力降比值C和子断层个数相关的参数N值的计算流程 M(Asperity)为凹凸体的地震矩，M(small)为小震的地震矩，MLarge为全部的地震矩.ΔσAsperity为凹凸体的应力降，Δσsmall小震的应力降. Fig. 4 The calculation process of C and N values M(Asperity) is the seismic moment of Asperity, M(small) is the seismic moment of small earthquake while MLarge is the seismic moment of large earthquake. |
|
|
表 2 余震相关参数的计算汇总 Table 2 Calculation summary of aftershock related parameters |
通过经验格林函数法的模拟得到10个台站的时程的观测值与模拟值的对比情况(图 5, 图 6)，以及10个台站反应谱的观测值与模拟值的对比情况(图 7).为方便直观的对比分析，此处去掉观测值和模拟值时程两端为零或者趋近于零的时段，所有时程均截取40 s的时间段进行对比.
|
图 5 台站051JZB, 台站051JZW, 台站051JZY, 台站051PWM, 台站062DIB时程的观测值与模拟值的比较情况 黑色时程表示观测值，红色时程表示模拟值.为方便比较，观测值截取了20~60 s时间段的时程，去掉两端加速度为0的时程段，模拟值选取的时长为40 s的时程记录. Fig. 5 Comparison between observed values and synthetic values of station 051JZW, station 051JZB, station 051JZY, station 051PWM, station 062DIB Red line means synthetic value while black line means observed values. For convenient comparing, the time interval of the 20~60 s time period was cut off from the observed accelerations, and the time history segments with accelerations of 0 at both ends were removed. The time length of the simulation values was selected about 40 s. |
|
图 6 台站062MXT, 台站062WEX, 台站062ZHQ, 台站062SHW, 台站62ZM2时程的观测值与模拟值的比较情况(图注同图 5) Fig. 6 Comparison between observed values and synthetic values of station 062MXT, station 062WEX, station 062ZHQ, station 062SHW, station 62ZM2 |
|
图 7 台站反应谱的观测值与模拟值的比较 红色表示模拟值加速度反应谱, 蓝色表示观测值加速度反应谱. Fig. 7 Station response spectrum comparison of observed and simulated values Red indicates simulated acceleration response spectrum while blue indicates observed acceleration response spectrum. |
通过对以上10个台站的模拟值与观测值的对比，结果显示所有台站的时程和反应谱整体上拟合的较好.PGA的强度、地震动持时的长度的匹配上整体较好，尤其是反应谱周期在小于1.0 s内模拟结果程度更好，这也验证了经验格林函数在模拟地震动高频成分时效果较好.模拟结果整体上可以体现此次九寨沟地震的地震动三要素的整体特征，即振幅、持时、频谱.但不可否认，有些台站在局部时程段或者局部反应谱周期段拟合不太好，例如台站062SHW的E-W方向分量的反应谱，其在周期0.1~1.0 s期间，反应谱观测值的幅值与模拟值的幅值差距较大.台站062WEX的反应谱在周期大于0.3 s时，观测值反应谱比模拟值反应谱要小很多.某些台站在时域范围内存在非平稳现象，例如062MXT台站.可以推断以上现象出现的主要几个原因：(1)余震受数据数量和质量的限制，台站的余震选取的欠佳，余震与目标主震的传播路径差异较大；(2)地震发生位置在地形较为复杂的九寨沟山区，实际地震波在传播时，可能会发生复杂的反射或折射现象，或者受到类似盆地的放大效应的影响，而在进行地震动模拟时很难将这些复杂的因素完整考虑在内.
但总的来说，时程和反应谱的模拟结果基本可以再现九寨沟地震主要的地震动特征，也可以验证经验格林函数法模拟九寨沟地震的可行性，况且此次模拟是在没有九寨沟地震自身余震的情况下进行的，用其他主震的余震作为格林函数输入，倘若有本次地震合适的余震，模拟结果会更好.随着国内地震记录的日趋丰富，各地区的余震数据会越来越多，未来可以建立国内各区域的小震记录数据库以备选取合适的格林函数，此举可以摆脱缺乏地震记录地区不能使用经验格林函数法的限制，使该方法更好的应用在地震工程中.
4 九寨沟漳扎镇地震动讨论九寨沟地震的震中位于九寨沟景区内，距离震中最近的人口聚集区是漳扎镇，距离震中约11.4 km(图 1).漳扎镇位于九寨沟景区北部，镇上的强震台站为九寨漳扎(103.8°E, 33.3°N)，编号为051JZZ.地震发生后，由于九寨漳扎台站的台基出现垮塌现象且记录到的峰值加速度异常，因此漳扎镇和九寨沟景区内的地震记录缺失.本节的目的是预测该区域的地震动强度.因为未搜集到九寨章扎台站的任何地震记录，故本文用临近的九寨百合(051JZB)台站的余震记录作为格林函数进行漳扎镇地震动的强度预测.台站九寨百合与台站九寨漳扎的震源距较为接近，九寨章扎台站的震源距为11.4 km，九寨百合的震源距为20.3 km，且两个台站所处的地形及场地条件相似.由于地震记录的缺失导致格林函数的选取较为被动，可以想象漳扎镇的地震动强度预测会存在较大的不确定性，很难给出精确的地震动强度值，所以本节的预测结果给出的是一个可能的取值范围而不是固定的值，我们选取了台站051JZB记录到的汶川地震的两个余震记录作为格林函数分别进行地震动强度预测(表 3).
|
|
表 3 格林函数余震信息 Table 3 Information about aftershock as Green function |
进行数值模拟之前，需要对选择的2个余震记录进行修正.通过对汶川地震MS4.0~MS5.5地震的地震事件的地震动特征的统计，在震源10~15 km范围内，PGA的值可能出现的取值范围大约在10~20 Gal之间.考虑到这2个汶川余震的震源距都超过100 km，要用这2个余震作为格林函数进行漳扎镇地震动强度预测时，其强度需要修正到震源距小于15 km内的地震动强度，这就需要给每一条余震记录乘上一个合适的放大倍数或者修正系数.所以将统计的PGA取值范围(10~20 Gal)除以2个余震的实际PGA得到余震需要放大的倍数范围.将修正后的余震记录作为格林函数运用经验格林函数法进行漳扎镇地震动强度预测，模拟的结果列于表 4中，分别给出了2个余震在最小PGA放大倍数和最大PGA放大倍数时得到的两个方向分量的PGA的模拟结果，得到的主震的地震动时程如图 8所示，主震的反应谱如图 9所示.
|
|
表 4 两个余震作为格林函数进行地震动模拟得到的PGA取值范围 Table 4 PGA of two aftershock as Green functions |
|
图 8 MS4.7余震和MS5.0余震合成的九寨沟地震漳扎镇的时程幅值范围 MS4.7余震合成的漳扎镇地区PGA范围为118~200 Gal; MS5.0余震合成的漳扎镇地区PGA范围为66.7~152 Gal. Fig. 8 The time history amplitude range of the JiuZhaigou earthquake in Zhangzha town, which was composed of the MS4.7 and MS5.0, respectively The range of the synthetic PGA in Zhangzha town by MS4.7 is in the range of 118~200 Gal while the MS5.0 is 66.7~152 Gal. |
|
图 9 两个余震在不同修正系数下合成的大地震的反应谱图 (a) EW方向分量的反应谱，MS4.7余震在修正系数分别为1.3(黑色线)和2.2(红色线)时的反应谱，MS5.0余震在修正系数分别为2.0 (蓝色线)和3.5(粉色线)时的反应谱; (b)NS方向分量的反应谱，MS4.7余震在修正系数分别为1.3(黑色线)和2.2(红色线)时的反应谱，MS5.0余震在修正系数分别为2.0(蓝色线)和3.5(粉色线)时的反应谱. Fig. 9 The response spectrum of two aftershock by different correction factors (a) The response spectrum of E-W component. MS4.7 with correction factors 1.3(black line) and 2.2 (red line); MS5.0 with correction factors 2.0 (blue line) and 3.5 (pink line); (b) The response spectrum of N-S component. MS4.7 with correction factors 1.3 (black line) and 2.2 (red line); MS5.0 with correction factors 2.0 (blue line) and 3.5 (pink line). |
漳扎镇地震动强度的预测结果是两个余震在最大、最小修正系数下得到的主震的地震动强度.从图 8的时程图来看.MS4.7级余震作为格林函数得到的主震的PGA的取值范围：EW方向分量为118~200 Gal，NS方向分量为117~199 Gal；MS5.0级余震作为格林函数得到的主震的PGA的取值范围：EW方向分量为65.1~114 Gal，NS方向分量为87.5~153 Gal，综合两个余震的模拟结果，漳扎镇的地震动PGA的E-W方向分量取值范围为65.1~200 Gal, N-S方向分量取值范围为87.5~199 Gal.漳扎镇周围几个台站的PGA(强震动观测简报.国家强震动台网中心. 2017(强震动观测简报，2017)(4)):九寨白河(051JZB)的E-W方向分量为129.5 Gal, N-S方向分量为185.0 Gal；九寨勿角(051JZW)的E-W方向分量为73.8 Gal, N-S方向分量为91.7 Gal；九寨永丰(051JZY)的E-W方向分量为45.8 Gal, N-S方向分量为66.7 Gal.考虑三个台站的地震动强度以及现有的衰减关系，我们推测九寨沟漳扎镇PGA的取值大于九寨白河(051JZB)、九寨勿角(051JZW)、九寨永丰(051JZY)台站的PGA，故九寨章扎台站的PGA取值范围约为180~200 Gal，震中位置的PGA为200 Gal左右.
图 9所示的是加速度反应谱在相同水平分量时，两个余震在不同修正系数下合成的大地震的反应谱.9a为EW方向分量的反应谱，MS4.7余震在修正系数分别为1.3(黑色线)和2.2(红色线)时合成大震的反应谱，MS5.0余震在修正系数分别为2.0(蓝色线)和和3.5(粉色线)时合成大震的反应谱；9b为NS方向分量的反应谱，MS4.7余震在修正系数分别为1.3(黑色线)和2.2(红色线)时合成大震的反应谱，MS5.0余震在修正系数分别为2.0(蓝色线)和和3.5(粉色线)时合成大震的反应谱.漳扎镇地震动的反应谱模拟值取值范围应该在红色和蓝色线之间，综合现有台站的PGA的观测值和模拟值，可以得到在研究区域(103.2°E—105.5°E；32.6°N—34.5°N)范围内PGA的观测值和模拟值的空间分布特征(图 10)，从图中可以较明显的看出，PGA的模拟结果与观测值拟合较好.
|
图 10 九寨沟地震PGA空间分布特征(漳扎镇和震源位置没有观测值PGA，所以在绘图时都取得模拟值的最高值200 Gal) (a) PGA观测值在E-W方向分量的分布图; (b) PGA观测值在N-S方向分量的分布图; (c) PGA模拟值在E-W方向分量的分布图; (d) PGA模拟值在N-S方向分量的分布图. Fig. 10 Jiuzhaigou earthquake′s PGA distribution characteristic (note: there is no observed PGA in Zhangzha town and source location, so we choose the highest synthetic value 200 Gal as PGA value both observed and synthetic values.) (a) Observed PGA distribution in E-W component; (b) Observed PGA distribution in N-S component; (c) Synthetic PGA distribution in E-W component; (d) Synthetic PGA distribution in N-S component. |
受台站数量的限制且分布不均匀特征的影响，震中西侧以及西南方向没有合适的台站记录，PGA的空间分布特征在椭圆长轴方向没有体现出与断层走向一致的特征.
5 结论以及格林函数选择标准的讨论 5.1 结论综合上述九寨沟地震时程图、反应谱图以及九寨沟地震PGA的空间分布特征，本文可以得到如下的结论：
(1) 九寨沟地震动模拟结果基本可以再现主要的地震动特征，10个台站的模拟结果整体上可以反映各台站地震动的强度特征，尤其是地震动高频成份拟合较好.模拟值的地震动峰值加速度、时程数据、反应谱等与观测值拟合较好.可以验证经验格林函数法模拟九寨沟地震的可行性.
(2) 预测结果显示漳扎镇的地震动峰值加速度值约为180~200 Gal.结果也表明在缺少大震的余震记录时，经验格林函数法使用其他大震的余震同样可以再现目标地震的强度特征.本研究也为经验格林函数方法在缺乏小震记录地区的使用积累了经验.
(3) 本文提出的计算子断层划分数目以及应力降比值的计算方式是可行的，输入N值和C值后合成的地震动可以再现该震级主要的地震动特征.
5.2 格林函数选择标准的讨论经验格林函数发展初期，在选择小震时最大的先决条件是主震与小震的震源机制要一致或近似，小震的初始破例位置与主震的破裂位置尽可能的接近，这样保证了大小地震的震源机制、传播路径和场地条件的高度一致，所以对于余震的选取非常苛刻.Hartzell(1978)最初假定如果在同一台站上记录到来自目标断层上足够多的点源所产生的小震记录，那么用此法可模拟出该台站的大震记录.此假定大大限制了该方法的实际应用，因为更多的情况是所模拟的大震断层附近的台站仅有少量记录. Kanamori(1979)提出对个别小震记录进行距离和辐射图案校正的方法，来修正Hartzell经验格林函数法的缺点.
Irikura等(1983, 1986)依据对大小地震震源参数的统计研究结果，提出了假定大小地震满足相似律并将大地震中的子源和所利用的小震的位错函数都取为斜坡函数，把问题从选择什么样的小震作为格林函数转换为如何合理地对小震记录进行修正使其可以等效成为大震中的一个子震源进而合成大地震，即对满足与大地震相似性特征的小震级记录加以修正即可作为格林函数合成大地震.罗奇峰(1989)等考虑大小地震应力降的差异，使该方法在大小地震断层形状不相似的情况下也可使用，可选为格林函数的小震记录更加容易得到.
Kurahashi和Irikura(2010)等在进行汶川地震动模拟时，将小震选取的条件设置为：(1)大小地震的震源机制一致；(2)震源距对传播路径的影响；(3)大小地震的场地条件尽可能的一致.这三条标准涵盖了地震的三个主要环节的影响，即震源、传播路径和场地条件的影响.但他选取的余震是发生在日本的地震，模拟结果同样可以表征大震的主要地震动特征.原因在于选择余震时完全是按照汶川地震的震源以及场地条件选择的.Kueahashi等的工作无意中证明了不同大震的余震记录可作为格林函数合成目标地震的地震动，是一次很有意义的尝试.受此启发，本文在合成九寨沟地震动时用的汶川地震和定西地震的余震记录.
Kamae等(1998)将随机方法和经验格林函数方法联系起来，提出一种混合方法来模拟地震动，即低频长周期部分用三维有限断层方法模拟，高频短周期部分的小震记录首先由Boore的随机方法生成，再用Irikura的经验格林函数方法模拟(用到的小震是由随机方法生成的)高频短周期地震动.Irikura等(Irikura and Kamae, 1994; Irikura and Miyake, 2011; Irikura et al., 2017)系统提出了用经验格林函数方法预测未来地震动的想法(李启程和景立平，2012)，结合多个震例模拟了大量地震动并总结了用该方法预测未来地震动的一般步骤.
地震动数值模拟未来发展的趋势是应用在预测未来破坏性地震的强度特征并据此给出科学合理的抗震设计的参考意见.先前的地震动模拟是对已发生地震的强度的反演，格林函数选取时要严格按照已发生大震的震源机制和震源参数为参考标准，无形中增加了诸多限制条件，但预测未发生的大地震时没有震源机制和参数的条件限制，毕竟主震的震源机制和参数都是未知的，这种发展趋势反而更加放宽了格林函数的选取限制，用随机的小震去预测未知震源机制的主震特征未尝不是一种可能的办法.综上所述本文给出几点未来进行强地震动预测时小震的选择标准及优先级，以供参考：
(1) 余震要选自研究区域内的台站的记录，即余震与目标主震要来自同一个台站，优先级也最高；(2)震级选择时原则上不选择破坏性震级(5.5级以上地震)的地震作为格林函数，因为首先高震级作格林函数缺乏模拟的意义，其次震级越高所包含的震源的特征过于明显，不具备作为格林函数的普适性；(3)震级太小的余震也不适合选择，其信噪比太低，噪声太大；(4)余震记录选择最好要位于目标地震的断层面上，这样有助满足与目标主震相似的传播路径；(5)只有在小震记录缺乏的时候再考虑修正其他地区的小震记录做为格林函数进行地震动强度预测.以上是本文对未来强地震动预测时格林函数选择标准的建议.
致谢 本文中所有的地震数据皆来自国家强震台网中心(由中国地震局工程力学研究所负责)，在此表示感谢.
Brune J N. 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes. Journal of Geophysical Research, 75(26): 4997-5009. DOI:10.1029/JB075i026p04997 |
China Strong Motion Network Centre. 2017. Strong ground motion observation bulletin (in Chinese). 2017(4).
|
Deng Q D, Gao X, Chen G H, et al. 2010. Recent tectonic activity of Bayankala fault-block and the Kunlun-Wenchuan earthquake series of the Tibetan Plateau. Earth Science Frontiers (in Chinese), 17(5): 163-178. |
Fang L H, Wu J P, Su J R, et al. 2018. Relocation of mainshock and aftershock sequence of the MS7.0 Sichuan Jiuzhaigou earthquake. Chinese Science Bulletin, 63(7): 649-662. DOI:10.1360/N972017-01184 |
Hanks T C, Kanamori H. 1979. A moment magnitude scale. Journal of Geophysical Research:Solid Earth, 84(B5): 2348-2350. DOI:10.1029/JB084iB05p02348 |
Hartzell S H. 1978. Earthquake aftershocks as Green's functions. Geophysical Research Letters, 5(1): 1-4. DOI:10.1029/GL005i001p00001 |
Haskell N A. 1969. Elastic displacements in the near-field of a propagating fault. Bulletin of the Seismological Society of America, 59(2): 865-908. |
Irikura K. 1983. Semi-empirical estimation of strong ground motion during large earthquake. Bull Disas Prev Res, 33: 63-104. |
Irikura K. 1986. Prediction of strong acceleration motion using empirical Green's function.//Proceedings of the 7th Japan Earthquake Engineering Symposium. Tokyo: Architectural Institute of Japan, 151-156.
|
Irikura K, Kamae K. 1994. Estimation of strong ground motion in broad-frequency band based on a seismic source scaling model and an empirical Green's function technique. Annals of Geophysics, 37(6): 1721-1743. |
Irikura K, Miyake H. 2011. Recipe for predicting strong ground motion from crustal earthquake scenarios. Pure and Applied Geophysics, 168(1-2): 85-104. DOI:10.1007/s00024-010-0150-9 |
Irikura K, Miyakoshi K, Kamae K, et al. 2017. Applicability of source scaling relations for crustal earthquakes to estimation of the ground motions of the 2016 Kumamoto earthquake. Earth, Planets and Space, 69: 10. DOI:10.1186/s40623-016-0586-y |
Kamae K, Irikura K, Pitarka A. 1998. A technique for simulating strong ground motion using hybrid Green's function. Bulletin of the Seismological Society of America, 88(2): 357-367. |
Kanamori H. 1979. A semi-empirical approach to prediction of long-period ground motions from great earthquake. Bulletin of the Seismological Society of America, 69(6): 1654-1670. |
Kurahashi S, Irikura K. 2010. Characterized source model for simulating strong ground motions during the 2008 Wenchuan earthquake. Bulletin of the Seismological Society of America, 100(5B): 2450-2475. DOI:10.1785/0120090308 |
Li Q C, Jing L P. 2012. A review of ground motion simulation with empirical Green's function. World Earthquake Engineering (in Chinese), 28(4): 89-94. |
Li Z C. 2017. The uncertainty factors research of near-field ground motion numerical simulation[Ph.D.thesis] (in Chinese). Beijing: Institute of Geophysics, China Earthquake Administration.
|
Luo Q F. 1989. Semi empirical synthesis of Near field acceleration[Ph.D.thesis] (in Chinese). Harbin: Institute of Engineering Mechanics, China Earthquake and Ministration.
|
Miyake H, Iwata T, Irikura K. 2003. Source characterization for broadband ground-motion simulation:Kinematic heterogeneous source model and strong motion generation area. Bulletin of the Seismological Society of America, 93(6): 2531-2545. DOI:10.1785/0120020183 |
Shan B, Zheng Y, Liu C L, et al. 2017. Coseismic Coulomb failure stress changes caused by the 2017 M7.0 Jiuzhaigou earthquake, and its relationship with the 2008 Wenchuan earthquake. Science China Earth Sciences, 60(12): 2181-2189. DOI:10.1007/s11430-017-9125-2 |
Shan X J, Qu C Y, Gong W Y, et al. 2017. Coseismic deformation field of the Jiuzhaigou MS7.0 earthquake from Sentinel-1A InSAR data and fault slip inversion. Chinese Journal of Geophysics (in Chinese), 60(12): 4527-4536. DOI:10.6038/cjg20171201 |
Somerville P, Irikura K, Graves R, et al. 1999. Characterizing crustal earthquake slip models for the prediction of strong ground motion. Seismological Research Letters, 70(1): 59-80. DOI:10.1785/gssrl.70.1.59 |
Xia C, Zhao B M, Horike M, et al. 2015. Strong ground motion simulations of the MW7.9 Wenchuan earthquake using the empirical Green's function method. Bulletin of the Seismological Society of America, 105(3): 1383-1397. DOI:10.1785/0120140001 |
Xu X W, Chen G H, Wang Q X, et al. 2017. Discussion on seismogenic structure of Jiuzhaigou earthquake and its implication for current strain state in the southeastern Qinghai-Tibet Plateau. Chinese Journal of Geophysics (in Chinese), 60(10): 4018-4026. DOI:10.6038/cjg20171028 |
Yang Y H, Fan J, Hua Q, et al. 2017. Inversion for the focal mechanisms of the 2017 Jiuzhaigou M7.0 earthquake sequence using near-field full waveforms. Chinese Journal of Geophysics (in Chinese), 60(10): 4098-4104. DOI:10.6038/cjg20171034 |
Yi G X, Long F, Liang M J, et al. 2017. Focal mechanism solutions and seismogenic structure of the 8 August 2017 M7.0 Jiuzhaigou earthquake and its aftershocks, northern Sichuan. Chinese Journal of Geophysics (in Chinese), 60(10): 4083-4097. DOI:10.6038/cjg20171033 |
Zhang X, Feng W P, Xu L S, et al. 2017. The source-process inversion and the intensity estimation of the 2017 MS7.0 Jiuzhaigou earthquake. Chinese Journal of Geophysics (in Chinese), 60(10): 4105-4116. DOI:10.6038/cjg20171035 |
Zheng X J, Zhang Y, Wang R J. 2017. Estimating the rupture process of the 8 August 2017 Jiuzhaigou earthquake by inverting strong-motion data with IDS method. Chinese Journal of Geophysics (in Chinese), 60(11): 4421-4430. DOI:10.6038/cjg20171128 |
邓起东, 高翔, 陈桂华, 等. 2010. 青藏高原昆仑-汶川地震系列与巴颜喀喇断块的最新活动. 地学前缘, 17(5): 163-178. |
房立华, 吴建平, 苏金蓉, 等. 2018. 四川九寨沟MS7.0地震主震及其余震序列精定位. 科学通报, 63(7): 649-662. |
国家强震动台网中心. 2017.强震动观测简报. 2017(4).
|
李启成, 景立平. 2012. 经验格林函数方法模拟地震动研究现状. 世界地震工程, 28(4): 89-94. DOI:10.3969/j.issn.1007-6069.2012.04.014 |
李宗超. 2017.大震近场地震动数值模拟不确定性研究[博士论文].北京: 中国地震局地球物理研究所.
|
罗奇峰. 1989.近场加速度的半经验合成[博士论文].哈尔滨: 中国家地震局工程力学研究所. http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y129975
|
单斌, 郑勇, 刘成利, 等. 2017. 2017年M7.0级九寨沟地震同震库仑应力变化及其与2008年汶川地震的关系. 中国科学:地球科学, 47(11): 1329-1338. |
单新建, 屈春燕, 龚文瑜, 等. 2017. 2017年8月8日四川九寨沟7.0级地震InSAR同震形变场及断层滑动分布反演. 地球物理学报, 60(12): 4527-4536. DOI:10.6038/cjg20171201 |
徐锡伟, 陈桂华, 王启欣, 等. 2017. 九寨沟地震发震断层属性及青藏高原东南缘现今应变状态讨论. 地球物理学报, 60(10): 4018-4026. DOI:10.6038/cjg20171028 |
杨宜海, 范军, 花茜, 等. 2017. 近震全波形反演2017年九寨沟M7.0地震序列震源机制解. 地球物理学报, 60(10): 4098-4104. DOI:10.6038/cjg20171034 |
易桂喜, 龙锋, 梁明剑, 等. 2017. 2017年8月8日九寨沟M7.0地震及余震震源机制解与发震构造分析. 地球物理学报, 60(10): 4083-4097. DOI:10.6038/cjg20171033 |
张旭, 冯万鹏, 许力生, 等. 2017. 2017年九寨沟MS7.0级地震震源过程反演与烈度估计. 地球物理学报, 60(10): 4105-4116. DOI:10.6038/cjg20171035 |
郑绪君, 张勇, 汪荣江. 2017. 采用IDS方法反演强震数据确定2017年8月8日九寨沟地震的破裂过程. 地球物理学报, 60(11): 4421-4430. DOI:10.6038/cjg20171128 |