据中国地震台网中心测定，北京时间2017年8月8日21时19分，我国四川省阿坝州九寨沟县发生里氏7.0级地震，震中位于103.82°E，33.20°N，震源深度约20 km.另据美国地质调查局(USGS，2017)利用远场波形给出的震源机制解结果，该地震矩震级为M_W6.5，震源深度为13.5 km，初步推断为一次走滑型地震.USGS结果还显示断层节面Ⅰ走向为246°，倾角为57°，滑动角为-173°；断层节面Ⅱ走向为153°，倾角为84°，滑动角为-33°，根据中国地震局地球物理研究所(2017)余震精定位结果中余震空间展布情况，可以初步判定破裂面为节面Ⅰ.此次地震发生于青藏高原东缘、巴彦喀拉地块东北边界，该区域断裂分布较多，地震活动极为强烈，此次地震震中100 km范围内，历史上就记载有9次6.5级以上地震，最大为1879年武都南8级地震.结合区域地质构造图能够发现，本次地震震中位于岷江断裂、塔藏断裂和虎牙断裂附近，其中岷江断裂、塔藏断裂和虎牙断裂均为全新世断裂，初步推测发震构造为塔藏断裂南侧分支和虎牙断裂北段(中国地震局地质研究所，2017).震区位于高海拔山区，目前还未发现明显的地表破裂带，这为确定发震断层带来一定影响.九寨沟7.0级地震给当地群众的生命财产造成了巨大损失，截至8月13日20时，地震造成25人死亡，525人受伤，6人失联，19768户176492人受灾，73671间房屋不同程度受损，涉及阿坝州九寨沟、松潘、若尔盖、红原等4个县53个乡镇309个村.震区余震不断，截至2017年8月14日8时00分共记录到余震总数为3704个，其中4.0~4.9级地震3个，3.0~3.9级地震27个，目前最大余震4.8级.

地震发生后，主震对余震的影响如何是普遍关心的问题.根据经典的余震静态触发理论，King等(1994)给出了计算主震后区域静态库仑应力变化(ΔCFS)的具体方法，得到的结果可以用来解释余震的分布，这一理论和基于该理论的Coulomb系列软件已经广泛应用于静态库仑应力变化的计算和余震危险性分析(万永革等，2009；单斌等，2012；解朝娣等，2010).此次7.0级地震造成了一定的建筑物破坏及人员伤亡，并引发了崩塌、滑坡等地质灾害，但是相对于芦山、鲁甸等震级相近的地震，此次地震的地震动强度相对偏弱，震害相对较轻.强地面运动是造成建筑物损坏的直接原因，研究主震产生的强地面运动水平则需要分析震源相关物理参数，因此准确评估强地面运动不仅对了解地震的震源机制和力学过程有重要的理论价值，而且在工程抗震和地震灾害的定量化评估、灾后救援及政府决策中有着重要的实际应用价值.

无论是静态库仑应力变化计算还是强地面运动模拟，震源参数尤其是断层滑动模型是必须要确定的参数.地震发生后，陈运泰课题组(中国地震局地球物理研究所，2017)根据远场地震波形快速反演出此次地震震源破裂模型，其结果显示断层面上的空间滑动量分布较为集中，主要滑动都分布在震源附近，有少量破裂传递至地表，但其规模很小，这可能也是地表难以发现破裂带的原因.季灵运等(2017)基于Sentinel-1 SAR影像获取了此次地震的同震形变场，反演获得同震滑动分布，计算了同震位错对余震分布和周边断层的静态库仑应力变化，并对发震构造进行了分析讨论.单新建等(2017)同样基于Sentinel-1 SAR数据反演了断层滑动分布，并基于三种不同接收断层计算了同震库仑应力变化.张旭等(2017)基于九寨沟地震同震InSAR资料和远场地震波反演了震源破裂过程，并利用基于有限动态源的烈度估计方法，对此次地震的烈度分布进行了估计.王阅兵等(2018)计算了震中350 km范围内10个GNSS台站的水平向同震位移量，并反演得到震源滑动分布.上述研究大多基于单一数据源，目前关于九寨沟7.0级地震研究中使用不同种类数据联合反演震源破裂的研究相对较少，尤其是缺乏近场三维同震位移数据与InSAR数据的结合.事实上，近场同震位移对地壳的介质结构不敏感，能对断层大小和滑动分布提供更强约束(金明培等，2014).

本研究利用基线校正方法获得九寨沟7.0级地震震中100 km范围内9个强震台站的同震位移，基于Sentinel-1卫星干涉SAR影像对获取InSAR同震形变场，结合GPS形变数据，利用多种同震形变资料联合反演得到断层滑动模型，基于该模型本文计算了此次四川九寨沟7.0级地震产生的静态库仑应力变化，考察其对余震的触发作用；另外，本研究将滑动模型作为输入参数，基于随机振动有限断层模型(Motazedian and Atkinson, 2005)，模拟九寨沟7.0级地震近断层强地面运动，得到峰值加速度(PGA)、峰值速度(PGV)等强地震动参数，并将其转化为烈度分布图，从物理角度分析震源参数对强地面运动水平可能产生的影响.

九寨沟7.0级地震震中附近强震台分布相对密集，震中100 km范围内共有9个强震台站.需要说明的是距离震中11.4 km的漳扎强震台(51JZZ)在东西向获取了约1.92g的地震动峰值加速度记录，这一记录超越了2013年芦山7.0级地震在宝兴地办强震台(51BXD)东西向产生的1.026g的峰值加速度记录，创造了目前我国在自由地表获得的最高地震峰值加速度记录.从数学上讲，对加速度时程记录经过两次积分可以得到位移时程记录，但在实际操作中往往会产生严重的漂移现象，本研究使用Wang等(2011)开发的SMBLOC程序，对9个强震台站记录进行基线校正，消除漂移现象，获取每个台站的三分量同震位移.图 1给出了51JZB强震台站原始加速度经过基线校正后的加速度、速度和位移时程曲线示意图.

图 1

Fig. 1

图 1 经过基线校正后的51JZB强震记录加速度、速度和位移时程曲线 Fig. 1 Acceleration, velocity and displacement recorded at station 51JZB after baseline correction

经基线校正后获得的强震台同震位移结果如表 1所示，结果显示经过基线校正后的51JZZ台站水平向同震位移高达171.6 mm，据了解该台站建于山腰斜坡处，历史上曾多次产生强震记录远高于衰减关系的异常记录，这与该台站所处地形、场地和构造环境有关.郑绪君等(2017)、王文才等(2017)、陆新征等(2017)在其各自研究中均未使用51JZZ的台站记录，鉴于此，本文中未将其值作为输入参与断层滑动模型反演.图 2比较了水平方向本研究获取的强震同震形变(黑色)和相同区域内王阅兵等(2018)计算的GPS同震形变结果(红色).

图 2

Fig. 2

图 2 基于强震记录(黑色)和GPS(红色)解算得到的水平向同震位移结果 Fig. 2 Horizontal coseismic displacements derived from strong motion data (black) and GPS observation (red)

表 1 Table 1

表 1 基于基线校正的强震同震位移结果 Table 1 Coseismic displacements after baseline correction Lon(°E) Lat(°N) EW(mm) NS(mm) UD(mm) 台站名 103.8 33.2 -63.0 159.6 55.2 51JZZ 104.1 33.3 -1.3 8.8 -1.0 51JZB 104.2 33.0 -16.7 3.4 -4.1 51JZW 104.2 33.2 -17.2 -3.7 -5.8 51JZY 103.6 32.7 -6.9 -16.6 -4.9 51SPC 104.5 33.6 -3.1 -1.9 -4.1 62SHW 104.3 33.8 1.5 -0.4 -1.4 62ZHQ 104.6 32.9 -0.3 1.2 0.1 62WEX 104.5 32.6 4.4 3.3 4.6 51PWM

表 1 基于基线校正的强震同震位移结果 Table 1 Coseismic displacements after baseline correction

图 2显示，强震台与GPS台站分布相对离散，共点观测的台站仅有两对，这为比较同一位置两种不同观测方法所获取的同震形变带来了困难.从宏观上看，两种观测结果均显示离震中越近，形变量越大.在两组共点观测点上，形变幅值和方向具有较好一致性，例如强震台51JZY和GPS站SCJZ共点观测，两者测得水平向同震位移量分别为18.4 mm和10.3 mm，均为各自观测方式测得的最大量.震中偏东方向水平位移场表现为向内聚敛，51JZY、51JZW和62SHW等几个台站最为明显，而在震中北侧和南侧的51JZB、51SPC、SCSP和GSZQ则表现为向外拉张，两种同震位移场均反映出此次地震错动以左旋走滑为主要方式.

本文采用欧空局Sentinel-1双星卫星编队的SAR数据，分别采用升降轨数据，升轨数据的干涉对是2017-07-30—2017-08-11，轨道号为Track 128，降轨干涉对是2017-08-06—2017-08-12，轨道号为Track 62，计算此次地震同震形变场，详细参数见表 2.

表 2 Table 2

表 2 2017年九寨沟7.0级地震干涉图详细信息 Table 2 Data on satellite images of the 2017 Jiuzhaigou MS7.0 earthquake 卫星类型 轨道号 飞行方向 成像模式 主影像 从影像 空间基线 入射角 Sentinel-1 T128 升轨 IW 20170730(S1A) 20170811(S1A) -35 m 44° Sentinel-1 T62 降轨 IW 20170806(S1A) 20170812(S1B) 92 m 42°

表 2 2017年九寨沟7.0级地震干涉图详细信息 Table 2 Data on satellite images of the 2017 Jiuzhaigou MS7.0 earthquake

其中升轨数据采用了两次Sentinel-1A成像数据，而降轨数据采用了Sentinel-1A/B干涉对，因此时间基线缩短到6天.本文采用ISCE软件对Sentinel-1进行干涉处理.外部DEM数据采用了选用30 m分辨率的SRTM-DEM数据，经过地理编码，采用传统的Goldstein滤波和平面拟合方式去除大气和轨道误差，得到九寨沟地震的同震形变场.

图 3展示了九寨沟7.0级地震的InSAR同震形变场，其中图 3a为升轨数据，抬升区域最大的LOS向形变量约10 cm，主要分布在断裂的末梢，最大的LOS向沉降量为22 cm，形变区主要分布于断层的西南侧，在东北侧形变量较小，但是也可以看出其形状分布特征呈现蝴蝶状四象限分布，从该干涉图中也可以初步推断出其为走滑断裂.图 3b是降轨数据，其中形变范围最大的区域与升轨数据中的沉降范围分布基本一致，最大的抬升量达到了10 cm，在断层的南端末梢位置存在一定的沉降信号，最大的沉降量约5 cm.升降轨数据的形变特征显示出不一致性，这主要和升降轨的成像几何特征以及地表形变方向有关.

图 3

Fig. 3

图 3 (a) 九寨沟7.0级地震升轨同震形变图；(b)九寨沟7.0级地震降轨同震形变图 Fig. 3 (a) Coseismic deformation of ascending orbit of the Jiuzhaigou MS7.0 earthquake; (b) Coseismic deformation of descending orbit of the Jiuzhaigou MS7.0 earthquake

反演前采用四叉树算法对InSAR形变场进行降采样，经过处理后升降轨数据量分别降至1630、459个.将降采样处理后的InSAR形变数据与强震同震形变、GPS同震形变(王阅兵等，2018)一起进行联合反演.本文使用德国地学中心SDM软件(Wang et al., 2008)进行震源破裂滑动反演研究，并采用分层介质模型计算格林函数，模型参数取自CRUST 2.0，介质参数如表 3所示.

表 3 Table 3

表 3 本文使用的分层介质模型参数 Table 3 Parameters of the layered medium model used in this study 层数 厚度 (km) VP (km·s-1) VS (km·s-1) 密度 (kg·m-3) 1 0.5 2.5 1.2 2100 2 20.0 6.1 3.5 2750 3 20.0 6.3 3.6 2800 4 9.5 7.2 4.0 3100

表 3 本文使用的分层介质模型参数 Table 3 Parameters of the layered medium model used in this study

根据余震精定位结果(中国地震局地球物理研究所，2017)，此次地震破裂断层相对单一，本研究将断层设定为1段，断层几何模型参考USGS(2017)震源机制解结果.为获得比较精细的滑动分布结果，本研究将断层划分为2 km×2 km的子断层，断层两端点位置设定为(33.412°N，103.730°E)、(33.088°N，103.903°E)，断层倾角设置在70°~90°范围内，得到的断层面长宽为40 km×32 km.子断层滑动角设定在-20°~20°之间，最大滑动量设定为1.0 m.为确定合理的断层位错分布，需选取合适的平滑因子α，本文采用位错模型的粗糙度和数据吻合程度之间的折中曲线来确定，如图 4a所示，最终取光滑因子α=0.25.联合反演的另一个关键问题是参与反演的各数据类别间的权重分配.本研究将强震数据(8个台站×3分量共24个)和GPS数据(7个台站×2分量共14个，无垂向分量)权重定义为1，选用了数据拟合误差作为量度来对强震、GPS数据与InSAR数据(升轨1630个，降轨459个，共2089个)的相对权重进行选择，基本思路就是选用能使得数据拟合残差最小的权重方案.根据拟合误差相对权重的变化结果(图 4b所示)，最终选取单个InSAR相对单个GPS的权重为0.02.按此权重计算，2089个InSAR数据相当于约42个GPS或强震数据，而反演中强震和GPS数据总量为38个，也就意味着反演中InSAR数据与强震和GPS数据所起作用大体相等.此外，GPS和强震数据在垂直方向上的解算精度均不及水平向，因此垂直向上的同震位移数据在反演时应该被赋予较小的权重，在随后实施的滑动模型反演中本研究取垂直向位移的权重值为水平向的1/4(金明培等, 2014, 2017；金明培和汪荣江，2013).

图 4

Fig. 4

图 4 (a) 拟合残差与粗糙度折中曲线图; (b) InSAR相对强震和GPS数据权重方案 Fig. 4 (a) Trade-off curve line between residuals and roughness; (b) Weighting scheme of InSAR with respect to strong motion and GPS data

图 5a和5b分别给出了单独基于InSAR数据和基于强震形变数据、InSAR和GPS数据联合反演得到的九寨沟7.0级地震震源滑动模型.在单独基于InSAR数据反演结果中释放的地震矩总量约为8.06×10¹⁸N·m，换算成矩震级为M_W6.57.而联合反演得到的结果显示断层面最大滑动量为0.74 m，平均滑动量为0.18 m，释放地震矩约为7.60×10¹⁸N·m，换算成矩震级为M_W6.52.两个模型均反映出滑动主要集中于20 km深度范围以上，整体上断层以走滑运动为主，尤其是在最大滑动区断层近乎纯走滑，在地表亦有少量破裂，但规模很小.两种滑动模型也存在差异之处：单独基于InSAR数据的反演结果在20 km深度以下出现了相当规模的滑动量，其释放的地震矩约为1.20×10¹⁸N·m，占释放地震矩总量的14.89%，相当于一次M_W6.0地震的地震矩释放量.本研究分别基于上述两种滑动模型模拟了升轨InSAR同震形变场并与观测值进行了对比，如图 6所示.结果显示，基于两种模型得到的升轨同震形变场与观测值在分布形态和幅值范围上均体现出基本一致，说明两种反演模型的结果是合理的.从形变残差图(图 6b、6d)可以发现，两种结果均不存在系统性残差，残差较大的地方集中于断层两侧位置，基于联合反演模型模拟得到的形变场残差整体上要小于基于InSAR反演模型的模拟结果，整体形变残差的均方根误差分别为2.03 cm和2.55 cm，前者和后者与观测值的相关系数分别为0.96和0.93，说明加入由强震和GPS解算的同震形变后，反演得到的断层滑动模型更好地约束了断层面上滑动分布.在下文中，本研究即以联合反演得到的滑动模型为对象进行分析和应用.

图 5

Fig. 5

图 5 (a) 基于InSAR同震形变场反演获得的九寨沟7.0级地震断层滑动模型；(b)基于强震数据、InSAR和GPS同震形变场联合反演获得的九寨沟7.0级地震断层滑动模型 Fig. 5 Slip model of the Jiuzhaigou MS7.0 earthquake inverted from the InSAR coseismic deformation field; (b) Slip model of the Jiuzhaigou MS7.0 earthquake from joint inversion of strong motion, InSAR and GPS coseismic deformation field

图 6

Fig. 6

图 6 (a)(b)分别为基于InSAR反演滑动模型模拟得到的升轨InSAR形变场及与观测值比较的残差；(c)(d)分别为基于强震数据、InSAR和GPS数据联合反演滑动模型模拟得到的升轨InSAR形变场及与观测值比较的残差 Fig. 6 (a) and (b) Simulated ascending InSAR deformation field based on slip model inverted from InSAR data and its residuals to the observed deformation field. (c) and (d) Simulated ascending InSAR deformation field based on slip model inverted from strong motion, InSAR and GPS data and its residuals to the observed deformation field

在许多研究中都将空间域内的断层滑动模型分布通过二维傅里叶变换转化到波数域内进行研究(申文豪等，2013).Andrews(1980)给出了在波数域中静态滑动位移分布D(k_x, k_y)和静态应力变化Δτ(k_x, k_y)的线性关系：

(1)

k_x和k_y分别为波数域的横轴、纵轴坐标，K(k_x, k_y)为与k_x和k_y有关的系数表达式(Ripperger and Mai, 2004).图 7给出了分别沿走向和沿断层倾角方向滑动模型对应的静态应力变化，正值表示增加，负值表示降低(即意味着能量释放).走向方向应力增加最大值为1.01 MPa，应力降低最大值为2.36 MPa，倾角方向应力增加最大值为0.58 MPa，应力降低最大值为0.99 MPa.对比两个方向上的静态应力变化结果可以发现，走向方向应力变化占据主要地位，在断层主要滑动区出现了较为集中的应力降低区域，说明该区域是能量主要释放区.

图 7

Fig. 7

图 7 (a) 断层面上沿走向方向的静态应力变化；(b)断层面上沿倾角方向的静态应力变化 Fig. 7 (a) Strike-parallel component of static stress change; (b) Up-dip component of static stress change

平均静态应力降是最重要的震源参数之一，它反映了地震动态破裂的平均效应，并且与能量释放及地震辐射相关，在强地面运动模拟中平均静态应力降的取值对地震波中高频成分影响显著(Atkinson and Wald, 2007).考虑到此次九寨沟地震以走滑运动为主，本研究以走滑方向应力变化来计算平均静态应力降.本研究以有效滑动区域内所有子断层应力降值的算术平均求取平均静态应力降，以超过最大滑移量20%的断层面积为有效滑动区域(Causse et al., 2009)，计算得到的平均静态应力降为1.07 MPa.Mohammadioun和Serva(2001)依据全球范围内相关震例统计得到一般走滑型中强地震产生的平均应力降值约为7.7 MPa左右，显然九寨沟7.0级地震平均应力降明显低于相同类型地震的平均水平.

断层滑动模型的另一个重要特征是其分布不均匀性，且同样影响高频地震波辐射(Causse et al., 2010).研究中通常以滑动模型在波数域内的拐角波数值来表示滑动分布不均匀性强弱(Somerville et al., 1999; Mai and Beroza, 2002).对一个给定的滑动模型，拐角波数越小就意味着滑动分布越均匀，而拐角波数越大就意味着滑动分布越粗糙，不均匀性越强.计算时同样需要将滑动分布转化到波数域得到其波数谱(D_obs)，再将其与理论的k^-2波数谱(D_mod)做差异化对比，以差异化参数E取最小值时对应的波数值作为拐角波数(Causs et al., 2010).E的表达式如下(Causs et al., 2010)：

(2)

理论的k^-2波数谱的表达式为：

(3)

式中k为k_x或k_y，k_c为待求拐角波数，k_N为Nyquist波数(Causs et al., 2010).沿走向和沿倾角两个方向的拐角波数计算结果如图 8所示.计算结果显示，沿走向和倾角两个方向的波数谱与理论k^-2波数谱具有一致的衰减趋势，说明联合反演得到的滑动模型符合天然地震断层破裂的运动学特征.沿走向和倾角两个方向拐角波数值分别为0.99×10^-4和1.10×10^-4，这一组数值明显低于Causs等(2010)给出的震级-拐角波数经验关系的预测值.Causse等(2010)的研究发现PGA随拐角波数值增加而增加，即PGA∝k_c^3/2，因此九寨沟地震低水平的拐角波数值可能也是本次地震强地震动水平不高的原因之一.

图 8

Fig. 8

图 8 断层滑动模型拐角波数计算示意图 (a)走向方向；(b)倾角方向.图中黑色实线为理论波数谱，黑色圆圈为滑动模型波数谱.黑色虚线对应的波数即为拐角波数值. Fig. 8 Calculation of the corner wavenumber for fault slip model (a) Strike direction; (b) Dip direction. The solid line indicates the theoretical k-2 model; the circles indicate the slip spectrum and the dashed lines denote the values of corner wavenumber.

最大剪切应力破裂理论认为，当脆性物质中某点的最大剪切应力达到临界值时，破裂即沿最大剪切应力平面产生.基于此可以定义描述临界破裂的库仑破裂应力CFS为(Harris, 1998)：

(4)

其中，|τ|为主震产生的应力张量在断层面上剪切应力的大小，S为内聚应力，μ为断层摩擦系数，σ_n为法向压应力，P_n为孔隙张应力.假设S和μ不随时间变化，库仑应力变化ΔCFS可以表示为(Harris, 1998)：

(5)

Δτ和Δσ_n分别为断层面上剪切应力和正应力变化，μ′称为等效摩擦系数或视摩擦系数，取值范围为0.2~0.8(盛书中等，2015).如果接收断层的几何参数和滑动方向等参数已知，可以计算指定断层面(faults of specified orientation)上的库仑应力变化，如果接收断层面的相关参数未知，则可以计算在主震和区域构造应力共同影响下的最优破裂面(optimally orientated faults)上的库仑应力变化.由于主震后余震的几何参数和滑动方向很难确定，因此通常计算上述库仑应力变化在最优破裂面投影来分析对余震的触发作用(缪淼和朱守彪，2012).本研究基于联合反演得到的滑动模型，利用Coulomb 3.3软件(Lin and Stein, 2004)计算九寨沟地震主震产生的库仑应力变化，计算时采用最优破裂面进行投影.最优破裂面的取向还与区域构造应力相关，本文使用世界应力地图2016版(World Stress Map 2016, http://dataservices.gfz-potsdam.de/wsm/showshort.php?id=escidoc:1680890)提供的区域应力场方位角及倾角数据进行计算(表 4)，计算中地壳的剪切模量取3.3×10¹⁰Pa，泊松比取0.25，等效摩擦系数取0.4(赵立波等，2016；程佳等，2016；缪淼和朱守彪，2012).

表 4 Table 4

表 4 区域构造应力场主要参数 Table 4 Main parameters of the regional tectonic stress field 方位角(°) 倾角(°) 应力值(MPa) 最大主应力 243 38 10 中间主应力 82 50 3 最小主应力 340 9 0

表 4 区域构造应力场主要参数 Table 4 Main parameters of the regional tectonic stress field

图 9显示了九寨沟7.0级地震产生的不同深度上静态库仑应力变化在最优破裂面上的投影结果，图 9a—9d的计算深度分别为5 km、10 km、15 km和20 km，图 9e为前述4种深度静态库仑应力变化的叠加最大值，图 9f为图 9e叠加余震分布示意图，红色圆圈即为余震活动(余震统计时间为主震后至8月15日，深度范围为0~25 km，震级M≥1.0).由图 9f可以看到，余震空间分布呈条带状，优势分布方向为NW，这与计算的库仑应力增加区域空间展布高度符合.统计结果显示，共有83.6%的余震位于库仑应力增加的区域.研究显示当ΔCFS>0.01 MPa即可触发余震(King et al., 1994)，本研究统计得到此次被主震所触发的余震占总数的77.9%，可见主震对后续余震具有显著的触发作用.

图 9

Fig. 9

图 9 九寨沟7.0级地震静态库仑应力变化在最优破裂面投影及余震分布 (a) 5 km深度静态库仑应力变化计算值；(b) 10 km深度静态库仑应力变化计算值；(c) 15 km深度静态库仑应力变化计算值；(d) 20 km深度静态库仑应力变化计算值；(e) 5~20 km深度静态库仑应力变化的叠加最大值；(f)静态库仑应力变化与余震分布. Fig. 9 Map showing projection of Coulomb stress change on preferred rupture plane and aftershock distribution of the Jiuzhaigou MS7.0 earthquake (a) Coulomb stress change at 5 km depth; (b) Coulomb stress change at 10 km depth; (c) Coulomb stress change at 15 km depth; (d) Coulomb stress change at 20 km depth; (e) The stacked maximum Coulomb stress change at 5~20 km depths; (d) Coulomb stress change with aftershock distribution.

本研究使用Motazedian和Atkinson(2005)提出的应用动力学拐角频率的随机振动有限断层模型进行区域强地面运动数值模拟，该模型已经应用于国内外多个震例并取得很好效果(申文豪等，2013；Atkinson and Boore, 2006).模型输入参数分为全局参数和局部参数，全局参数主要用来描述发震断层的宏观几何特征和平均破裂过程，包括断层破裂面积、长度、宽度、埋深、走向、倾角、平均静态应力降、品质因子Q模型、几何衰减模型等.局部震源参数包括断层滑移分布等相关参数.断层的长度、走向、倾角、埋深、位置、滑动模型、平均静态应力降等几何参数可根据本文第1节中结果确定，Q模型参数的确定均来自乔慧珍等(2006)，取Q=334.4f^0.581，另外本研究取平均破裂速度为剪切波速的0.8倍.表 5总结了模型计算参数.

表 5 Table 5

表 5 随机振动有限断层模型模拟九寨沟7.0级地震强地面运动参数 Table 5 Parameters of the stochastic finite fault model for simulating strong ground motion of the Jiuzhaigou MS7.0 earthquake 参数名 参数值 数据来源 震源 走向: 246°倾角: 74° USGS 矩震级 6.52 联合反演结果 平均应力降 1.07 MPa 联合反演结果 断层长宽 40 km×32 km 联合反演结果 破裂速度 0.8β Atkinson and Boore (2006) 脉冲百分比 50% Atkinson and Boore (2006) 剪切波速β 3.7 km·s-1 Q值模型 334.4f 0.581 乔慧珍等(2006) 几何衰减模型 1/R R < 75 km 75-1 75 km≤R < 125 km 75-1(125/R)0.5 R>125 km Motazedian and Atkinson(2005)

表 5 随机振动有限断层模型模拟九寨沟7.0级地震强地面运动参数 Table 5 Parameters of the stochastic finite fault model for simulating strong ground motion of the Jiuzhaigou MS7.0 earthquake

根据野外调查烈度结果，本研究将研究区域划定在32°N—34.5°N，102.5°E—105°E，网格精度为0.05°.本研究模拟计算了区域内每个网格点的峰值加速度(PGA)、峰值速度(PGV)，并将PGA和PGV通过经验关系(丁宝荣等，2017)分别转化为烈度(MMI)，最后取两种烈度的加权平均为最终烈度值，并与野外调查结果进行比较.为更好展示强地面运动参数空间衰减趋势，本研究分别取PGA、PGV以10为底的对数lg(PGA)，lg(PGV)进行绘图，如图 10a、10b所示.图 10c为模拟烈度和中国地震局公布的实际调查烈度的叠加图.

图 10

Fig. 10

图 10 基于联合反演滑动模型模拟得到的九寨沟7.0级地震强地面运动 (a)(b)分别为lg(PGA)，lg(PGV)分布示意图；(c)中红色和白色实线分别为模拟烈度和实际调查烈度分布示意图. Fig. 10 Strong motion simulation for Jiuzhaigou MS7.0 earthquake by the joint inverted slip model (a) and (b) lg(PGA) and lg(PGV) distribution. (c) Simulated and observed MMI (Modified Mercalli Intensity) distribution respectively.

从图 10中可以看到，本文模拟结果在烈度分布的范围以及等级方面与公布结果符合度较高，说明相关模型参数结果是可信的.野外调查烈度图由于其采样点的有限性，各个烈度对应的分布范围大都是经验性地圈划为圆形或者椭圆形.本研究模拟结果中，Ⅸ度区范围要明显小于实际调查，而Ⅵ度区范围要大于实际调查结果.本文在进行模拟时，震源位于断层左上方位置，地震是一个朝东南方向拓展的单侧破裂，因此模拟得到的烈度分布结果还体现出一定的方向性特征：Ⅵ、Ⅶ度区的烈度分布有向东南方向拓展的趋势，这也说明了随机有限断层模型可以有效反映近场地震动特征.

九寨沟7.0级地震发生在东昆仑断裂带与岷江断裂带的交汇处，其周边活动断裂构造十分复杂，相关断裂的地表形迹很不明确，包括1879年武都南8级地震在内的几次历史地震的发震构造及其相互关系也依然不明朗.此次地震震后的应急考察没有发现明显的地表破裂带，而本文联合反演结果也显示，破裂量在地表处规模不大.InSAR数据的测量误差主要来自于卫星轨道误差、DEM误差和大气效应等，本研究干涉处理中已经对各种误差进行了一定的削弱和去除，另外利用四叉树采样方案也在一定程度上抑制了数据误差，突出了有效信息，尽管如此，测量误差的影响仍然存在，使得震源参数估值产生偏差、精度降低.与季灵运等(2017)结果类似，本研究单独使用InSAR数据得到的反演结果在20 km深度以下出现了相当规模的滑动量，而增加强震和GPS数据约束后，反演结果更加合理，模拟形变场残差也更小.本研究使用的SAR辅影像为震后3~4天内获取，包含的震后余滑或黏弹性松弛效应引起的震后形变量较小，说明InSAR数据观测误差影响了断层深部滑动量的估算，直接导致震级高估，而在联合反演中，由于对不同数据采用了合理的权重配置方案，使得InSAR测量误差的影响相对减弱，从而得到更加合理的反演结果.因此，测量误差的识别及合理补偿显得尤为重要.

借鉴USGS震源机制解结果，本文设定滑动角在-20°~20°范围变化，结果显示最大滑动区断层近乎纯走滑，而在其周边区域有逆冲或正断层分量，这与其他反演结果一致.基线校正法最重要的内容是确定事件前漂移、事件中漂移和尾波段漂移的两个特征时间点，而目前这两个特征时间点的确定一般基于单次事件的经验性准则，往往会导致出现个别异常测点，从而给反演带来一定误差.例如，如果将一些受局部非弹性场地效应影响的场点形变值参与反演，极可能造成断层浅部滑动量增加，造成对震级结果的过高估计(金明培和汪荣江，2013).因此，在反演之前要尽量剔除明显异常的观测值.另外，本文反演过程中采取的平滑性条件，也会滤除滑动模型局部的一些细节信息，这都可能会影响最终结果的分辨率.

在基于反演滑动模型计算区域库仑应力变化时，不同的有效摩擦参数取值对于最终结果有重要影响(万永革等，2000).研究表明，有效摩擦参数与断层类型和区域构造背景都有关系(Parsons et al., 1999)，目前研究普遍取值为0.4，而未细致考虑断层类型和地质背景可能给最终结果造成误差.另外，近期的研究显示除主震产生的同震库仑应力变化外，震后余滑对余震也有重要影响(Perfettini and Avouac, 2007；Hsu et al., 2006; Bourouis and Bernard, 2007)，因此在研究余震触发机制时仅考虑同震静态库仑应力变化的影响而忽略震后库仑应力变化显然是不全面的.张旭等(2017)基于有限动态源的烈度估计方法对此次地震的烈度分布进行了估计，其方法基础是烈度衰减经验关系，而本研究使用的随机振动有限断层模型基础是Hanks和McGuire(1981)提出的地震动随机合成理论，两种方法在物理机制上有根本区别.本研究模拟的烈度分布结果与野外调查烈度整体特征相似，张旭等(2017)模拟结果中出现两处不连续的Ⅸ度区且分别位于断层两端，这可能与其反演模型中存在两处凹凸体有关.随机振动模型并未计算地震波在地壳浅层速度结构中的格林函数，因此不会产生体波和面波.另外，该模型仅考虑S波频谱的影响，得到的是水平方向的地震动参数，无法得到垂直方向上的地震动参数，这些都是该模型存在的不足之处(申文豪等，2013).

本研究利用基线校正方法获得九寨沟7.0级地震震中100 km范围内9个强震台站的同震位移，并基于Sentinel-1卫星干涉SAR影像对获取了InSAR同震形变场.结合GPS形变场(王阅兵等，2018)，本研究采用分层介质模型进行震源滑动模型联合反演，结果显示此次地震整体以走滑运动为主，断层面最大滑动量为0.74 m，平均滑动量为0.18 m，释放地震矩约为7.60×10¹⁸N·m，换算成矩震级为M_W6.52，地震在地表有少量破裂，但其规模很小.通过模拟形变场并与观测值对比发现，联合反演结果优于单独基于InSAR形变场的反演结果.依据波数域中静态滑动位移分布和静态应力变化的线性关系，本研究计算了沿走向和倾角方向静态应力变化，结果显示走向方向应力变化占据主要地位，该方向应力增加最大值为1.01 MPa，应力降低最大值为2.36 MPa，倾角方向应力增加最大值为0.58 MPa，应力降低最大值为0.99 MPa，平均静态应力降为1.07 MPa.本研究计算得到反演滑动模型沿走向和沿倾角两个方向拐角波数值分别为0.99×10^-4和1.10×10^-4，且在波数域的傅里叶谱值符合k^-2衰减规律.本文利用反演震源滑动模型计算了主震产生的同震静态库仑应力变化，结果显示共有83.6%的余震位于库仑应力增加的区域，被主震所触发的余震占总数的77.9%，主震对后续余震具有显著触发作用.本研究基于随机振动有限断层模型计算了九寨沟地震区域PGA、PGV以及烈度值并与野外调查结果进行比较，结果显示模拟结果在分布范围以及等级方面与调查烈度结果符合度较高，模拟结果能够很好地反映断层破裂的方向性效应等特征.本研究计算结果显示无论是平均静态应力降的值还是拐角波数的值均低于同类型地震的平均水平，这可能是造成本次地震强地震动水平相对不高的原因.