2017, Vol. 60
2. Canada Center for Mapping and Earth Observation, Natural Resources Canada, Ottawa, K1A0E4
2. Canada Center for Mapping and Earth Observation, Natural Resources Canada, Ottawa, K1A0E4
根据中国地震台网测定，北京时间2017年8月8日21时19分于我国四川省九寨沟县发生了MS7.0级地震(下文称之为“2017年九寨沟MS7.0地震”)，震中位于33.20°N，103.82°E，震源深度20 km (图 1).截至8月13日20时，此次地震共造成25人死亡，525人受伤.
|
图 1 地震区域构造及震后早期余震分布 五角星为不同机构测定的主震震中位置(参见表 1)，其中红色五角星来自中国地震台网中心(CENC)，绿色五角星来自美国地质调查局(USGS)，紫色五角星来自全球矩心矩张量解(GCMT).沙滩球表示不同机构给定的主震震源机制解(参见表 2).红色圆圈为震后约13 h内ML≥0.5的余震(数据来自中国地震台网中心).蓝色虚线框展示了本研究使用的同震InSAR资料覆盖区域.白色实线表示区域内活断层, 塔藏断裂(TZF)，岷江断裂(MJF)和虎牙断裂(HYF). Fig. 1 Tectonic settings and aftershocks distribution Stars are the epicenter locations of the mainshock from different institutions, with the red one being from CENC, the green from USGS and the purple from GCMT. Beach balls show the focal mechanism solutions of the mainshock from different institutions (see Table. 2 for details). Red circles are the aftershocks (ML≥0.5) within the first 13 hours (from CENC). Blue dashed frame shows the spatial coverage of the InSAR data used in this study. White solid curves are active faults, Tazang fault (TZF), Minjiang fault (MJF) and Huya fault (HYF). |
|
表 1 主震震源位置 Table 1 Hypocenter locations of the mainshock |
|
|
表 2 主震震源机制 Table 2 Focal mechanism solutions of the mainshock |
除中国地震台网外，国内外其他地震研究机构也发布了震源位置与震源机制解信息(表 1和表 2；图 1).中国地震局地球物理研究所利用双差定位方法(Waldhauser and Ellsworth, 2000)迅速确定了早期余震的位置(房立华等，2017)，还有一些研究机构发布了基于远场体波的震源破裂过程初步反演结果(罗翊菁和赵里，2017；王卫民等，2017；张勇等，2017)，这些结果之间仍有明显的差异存在.随后，一些国内外学者利用InSAR同震位移反演得到了断层面上的静态位错分布(冯万鹏，私人交流；单新建等，2017)，这些结果与地震学反演结果明显不同.
为了确认2017年九寨沟MS7.0地震的发震断层面的几何属性、断层面位错的静态分布以及震源破裂的时空特征，从而认识发震断层的复杂性，我们首先利用双差定位方法(Waldhauser and Ellsworth, 2000)对主震后约13个小时内的余震进行了定位，然后利用远场体波反演、近场同震InSAR资料反演以及二者的联合反演技术确定了震源的破裂过程，最后利用基于有限动态破裂模型的烈度估计方法(许力生等，2016)估计了此次地震的烈度分布.
2 余震双差定位为了确认发震断层的几何属性，我们利用双差定位技术对中国地震台网中心提供的震后约13 h内的震相数据进行了处理，重新确定了293次余震的相对位置.从图 2可以清楚地看到，重新定位前后震源位置发生了明显变化，无论在地表分布或是垂直剖面的分布，定位后的余震更显集中，成簇特征更加凸现.
|
图 2 双差定位前(红色)后(蓝色)的余震分布 (a)双差定位前后余震震中分布，其中插图展示了用于双差定位的台站(蓝色三角形)分布; (b)双差定位前后余震在AA′剖面上的投影; (c)双差定位前后余震在BB′剖面上的投影.红色五角星均表示主震. Fig. 2 Distribution of the aftershocks before (red) and after (blue) the DD-relocation (a) Aftershock epicenters. The inset shows the stations (blue triangles) used in doing the DD-relocation; (b) Aftershocks projected on the profile AA′; (c) Aftershocks projected on the profile BB′. Note: the stars refer to as the mainshock. |
假设发震断层为平面断层且余震发生在主震断层面及其附近，通过线性拟合我们分别确定了断层面的走向和倾角.如图 2所示，最佳拟合走向约为142°，最佳拟合倾角为81°.拟合走向与远场震源机制反演结果相差约10°(见表 2)，但拟合倾角与远场震源机制反演结果非常接近(见表 2).
基于拟合结果，我们在走向方向和倾向方向分别建立两个垂直剖面AA′和BB′，并将重定位的余震投影其上.从图 2可以看出，大部分余震的震源深度分布在5~15 km之间.从图 2b可以看出，余震沿发震断层展布约35 km，其中震中东南15 km，震中西北20 km；西北段余震较浅，东南段余震较深.从图 2c可以看出，定位后余震分布更为集中，且条状分布更加明显.这种条状特征似乎表明发震断层向西南倾斜，且深部较陡浅部较缓.
3 主震破裂过程 3.1 数据与初始模型构建远场地震资料来自从IRIS数据中心下载的震中距在30°~90°范围内的全球分布的台站.根据记录的信噪水平和空间分布的均匀性，我们挑选了其中36个台站(如图 3a所示)的垂直向记录，并利用0.02~0.2 Hz的带通滤波对其进行了滤波以压制非同震长周期和短周期信号干扰.参与反演的波形长度为40 s，其中P波初动前后分别为10 s和30 s.计算格林函数时，我们采用了ak135全球一维速度模型(Kennett et al., 1995)和反透射系数法(Wang, 1999)，并将其保持在与观测资料相同的频带范围内.
|
图 3 地震台站分布及反演参数的确定 (a)地震台站(蓝色三角形)与主震震中(红色沙滩球)分布; (b)相对失配度随震源深度的变化.红色圆圈对应于最佳震源深度.(c)方差降随最大破裂速度与子断层上升时间的变化.红色圆圈对应于最佳破裂速度与上升时间; (d)方差降随InSAR资料相对于地震资料的相对权重值的变化.红色圆圈对应于最佳权重值. Fig. 3 Distribution of the teleseismic stations and determination of the inversion parameters (a) Distribution of the teleseismic stations (blue triangles) and the epicenter location of mainshock (red beach ball); (b) Variation of the relative misfit with the source-depths. The red circle is used to emphasize the preferred depth; (c) Variation of the variance reduction with the maximal rupture velocity and rise time. The red circle is used to highlight the preferred velocity and time; (d) Variation of the variance reduction with the weight of the InSAR data relative to the seismic data. The red circle refers to as the preferred value. |
同震InSAR数据来自欧洲空间局(ESA).我们收集了Sentinel-1A沿128轨道(T128) 扫描干涉模式的两景雷达图像(表 3).其中只有覆盖震区的第3个子条带中3个Burst被提取出来并进行了干涉处理(Feng et al., 2016, 2017a, 2017b).形变区内最大视线向(Light of Sight, LOS)形变达~-22 cm.近断层两侧出现部分连续失相干，不能直接判定最大跨断层视线向位错量.干涉图像经用四叉树算法(Simons et al., 2002)降采样为2376个点(表 3).同样基于ak135全球一维速度模型(Kennett et al., 1995)，利用Fortran程序EDGRN/EDCMP(Wang et al., 2003)计算了在InSAR观测资料同名点的静态位移格林函数.
|
|
表 3 InSAR数据信息 Table 3 Information on the InSAR data |
初始破裂模型的选择对于反演至关重要.众所周知，近场波形对初始破裂时间和位置以及断层面参数都很敏感，近场静态观测对断层面参数很敏感但与初始时间无关，而远场波形对起始时刻和起始位置并不是很敏感.本研究只涉及远场波形资料和近场InSAR资料，所以，我们首先采用了中国地震台网中心测定的震源位置和USGS确定的W-phase矩张量反演提供的断层面参数，让断层面穿过33.20°N/103.82°E，保持断层面走向153°以及断层面倾角84°.断层面上边界位于地表，沿断层面向下拓展30 km作为下边界；东南边界距起始破裂点12 km，西北边界距起始破裂点20 km.对于子断层的尺度，我们进行多次尝试，最终选定为2 km×2 km.由于震源深度的不确定性素来较大，事实上多家机构给出的震源深度也不同(表 1)，所以我们采用波形反演方法(Zhang et al., 2012)对震源深度进行搜索(如图 3b所示)，获取的最佳深度为11 km.
所幸的是，利用上述初始破裂点的位置和断层面参数可以顺利地进行反演，并且正如下文所述，反演结果不但较好地解释了远场波形资料，而且很好地解释了近场InSAR资料，这表明如此确定的初始模型暂且无可挑剔.
3.2 反演方法基于滑动角可变(Zhang et al., 2012)的震源破裂过程时间域反演方法(Chen and Xu, 2000; Xu et al., 2002)，张旭(2016)发展了基于地震波形资料和静态大地测量资料的有限断层破裂时空过程联合反演方法.该反演方法不仅可以单独反演地震波形和大地测量资料，而且可以联合反演两类数据.该方法不需要预先给定子断层震源时间函数形状，而是通过共轭梯度法(Ward and Barrientos, 1986)迭代反演子断层震源时间函数，从而避免了先验假定给反演结果带来的影响，同时该反演方法允许子断层的滑动方向在总体的滑动角附近(±45°)发生瞬时变化.
另外，为了稳定反演结果以使其具有可接受的物理意义，我们还引入了时间和空间光滑约束(Yagi et al., 2004; Zhang et al., 2012)以及标量地震矩最小约束(Hartzell and Heaton, 1983; Antolik and Dreger, 2003; Zhang et al., 2012).时间光滑约束用于抑制子断层震源时间函数的过高频率，空间光滑约束用于消除相邻子断层间位错的不连续性，而标量地震矩最小约束则用于压制较弱的过低频滑动.
3.3 远场体波反演利用地震波形资料反演震源破裂过程时，最大破裂速度和子断层上升时间需事先给定，然而破裂速度与子断层上升时间之间存在“折衷(trade-off)”效应，最终影响反演结果(e.g., Lay et al., 2010; Gusman et al., 2015).为了弱化这种效应，我们采用网格搜索策略，尝试不同的破裂速度和上升时间，评价模型对观测资料的拟合度，选择拟合度最好的参数对.如图 3c所示，选取3.0 km·s-1的最大破裂速度和6 s的子断层上升时间，观测资料的拟合度最佳.
图 4展示了远场体波的反演结果.结果表明，此次地震破裂持续时间~15 s，释放总标量地震矩~6.68×1018 Nm，相当于矩震级MW6.5，最大滑动量达~1.2 m.同震滑动区主要位于震源附近，破裂区几乎呈圆盘形，走向方向约15 km，倾向方向约18 km，且西北方向略浅，东南方向略深，这一趋势性特征与余震的分布特征一致(图 2b).滑动方式以走滑为主，但具有一定的正滑分量(总体滑动角约为-14.7°).图 4c展示了反演结果对观测资料的拟合情况，观测波形与合成波形的平均相关系数达0.82.
|
图 4 基于地震资料的反演结果及其对观测资料的解释 (a)地震矩率函数；(b)位错分布.蓝色实线分别表示滑动量为0.24 m，0.48 m，0.72 m和0.96 m的等值线.红色五角星表示起始破裂点；(c)观测波形与合成波形的比较.每幅子图左侧从上至下依次为台站名，震中距和方位角，右侧从上至下依次为相关系数，震相名和台站分量. Fig. 4 The inverted results of the mainshock based on the teleseismic data and the explanation to the observed data (a) Moment rate function; (b) Static slip distribution. The blue solid curves show the slip contours of 0.24 m, 0.48 m, 0.72 m and 0.96 m, respectively, and the red star refers to as the initial point; (c) Comparison of the observed and synthetic waveforms. On the left side of each subplot are station code, epicentral distance and station-azimuth, respectively, and on the right side are correlation coefficient, phase name and channel, respectively. |
从图 1可以看出，我们获取的Sentinel-1的InSAR观测对震区的覆盖非常理想.因此，我们对InSAR同震地表位移进行了单独反演，以了解断层面静态位错的空间分布.
反演结果(如图 5所示)表明，此次地震破裂释放的标量地震矩约6.67×1018 Nm，相当于矩震级MW6.5，和地震资料反演结果几乎一致.不同的是，断层面上出现两个凹凸体，其中一个位于起始破裂点南东较浅处，水平方向约12 km，垂直方向约10 km，最大滑动量约1.1 m，且该凹凸体破裂以走滑机制为主(总体滑动角约为-3.9°)；另一个位于起始破裂点北西方向相对较深处，水平方向约17 km，垂直方向约10 km，最大滑动量约1.2 m，该凹凸体破裂具有明显的正滑分量(总体滑动角约为-38.7°).两个凹凸体总体的滑动角约为-15.7°, 以走滑为主, 兼具少量正滑分量.
|
图 5 基于InSAR资料的反演结果及其对观测资料的解释 (a)反演得到的位错分布.蓝色实线分别表示滑动量为0.24 m，0.48 m，0.72 m和0.96 m的等值线.红色五角星表示起始破裂点; (b)观测位移; (c)合成位移; (d)合成位移与观测位移间的残差; (b)—(d)中红色五角星表示主震震中. Fig. 5 The inverted result of the mainshock based on the InSAR data and the explanation to the observed data (a) Static slip distribution. The blue solid curves show the slip contours of 0.24 m, 0.48 m, 0.72 m and 0.96 m, respectively, and the red star refers to as the initial point; (b) Observed Line-of-Sight (LOS) displacements; (c) Model-predicted LOS displacements; (d) Residuals between the predicted and observed LOS displacements. The red stars in subplots (b)—(d) refer to as the epicenter locations of mainshock. |
图 5b—5d展示了基于反演结果的合成资料与观测资料的比照情况.可以看出，反演得到的模型可以很好地预测观测到的地表位移特征，计算得到的均方根误差约为0.0168 m，合成资料和观测资料的相关度达~0.96.
3.5 远场地震资料和近场InSAR资料联合反演为了融合地震波形资料的时间信息和近场同震InSAR资料的空间信息对震源时空破裂过程进行约束，我们对远场波形资料和近场InSAR位移资料进行了联合反演.
为了确定联合反演中地震资料和InSAR资料的相对权重，我们同样采用网格搜索策略，尝试不同的权重值，评价模型对资料的解释程度.如图 3d所示，当相对权重值为1时，两种资料得到最佳拟合.因此，最终选取的InSAR资料相对于地震资料的权重为1.
联合反演结果(如图 6所示)表明，此次地震破裂持续时间约15 s，释放总标量地震矩约6.61×1018 Nm，相当于矩震级MW6.5，最大滑动约1.0 m.如图 6b所示，此次地震的破裂区可分为两个主要的凹凸体，但不同于单独反演InSAR资料的结果.一个较大，距起始破裂点较近，从起始破裂点北西较深的位置延展到起始破裂点南东较浅的位置，可能破裂到了地表.此凹凸体主要以走滑为主(总体滑动角约为-14.0°).另一个凹凸体较小，距起始破裂点较远，位于起始破裂点北西约15 km且6 km深处.此凹凸体的正滑分量非常明显(总体滑动角约为-52.8°).根据瞬时破裂图像6c，此次地震起始于震源位置，首先向北西和浅部破裂，然后向深部和南东方向扩展，最后终止于北西较远的凹凸体和南东较浅位置.两个凹凸体总体的滑动角约为-19.5°, 以走滑为主, 兼具少量正滑分量.
|
图 6 基于地震资料和InSAR资料的联合反演结果 (a)地震矩率函数; (b)位错分布.蓝色实线分别表示滑动量为0.24 m，0.48 m，0.72 m和0.96 m的等值线; (c)滑动速率快照; (b)—(c)中红色五角星表示起始破裂点. Fig. 6 The jointly inverted results of the mainshock based on the joint inversion of the teleseismic data and the InSAR data (a) Moment rate function; (b) Static slip distribution. The blue solid curves show slip contours of 0.24 m, 0.48 m, 0.72 m and 0.96 m, respectively; (c) Snapshots of the slip rate. The red stars in subplots (b)—(c) refer to as the initial points |
为了验证反演得到的动态破裂模型，我们利用反演结果计算了相应的远场地震波形和近场InSAR同震位移.如图 7所示，合成的地震波形与观测波形的平均相关系数达到0.81，合成的近场InSAR位移与观测位移的相关度达约0.95，且均方根误差仅为0.0193 m.
|
图 7 联合反演结果对观测资料的解释 (a)观测波形与合成波形的比较; (b) InSAR观测位移; (c) InSAR合成位移; (d) InSAR合成位移与观测位移间的残差. Fig. 7 The explanation of the jointly inverted results to the observed data (a) Comparison of the observed waveforms and the synthetic ones; (b) Observed InSAR LOS displacements; (c) Model-predicted InSAR LOS displacements; (d) Residuals between the predicted InSAR LOS displacements and the observed ones. |
烈度是对地震复杂地面运动的简单描述，更是对震害的直接反映(Wald et al., 1999a, 1999b), 震后烈度的快速获取对于地震应急救援十分重要(许力生等，2016).因此，利用基于有限动态源的烈度估计方法(许力生等，2016)，基于反演得到的有限动态源模型，我们对此次地震的烈度分布进行了估计.考虑到近场InSAR位移并不是严格的同震位移，由此反演得到的位错模型也不是严格的同震位错模型，所以，基于单纯的地震资料反演结果和联合反演结果分别计算了两种情况下的烈度分布.
图 8a展示了基于地震资料反演结果计算的烈度分布.可以看出，烈度椭圆长轴方向与主震断层空间展布方向一致，也与余震空间分布方向一致，最大烈度约为Ⅸ度，震中南东方向和北西方向均有分布，但区域很小.Ⅵ-Ⅷ度区大体呈椭圆形，但震中北西侧的面积略大.
|
图 8 基于有限动态源模型估计的烈度分布 (a)基于地震波反演结果的烈度; (b)基于联合反演结果的烈度.红色五角星表示震中位置，白色实线表示活断层，黑色实线表示烈度为Ⅸ度的等震线. Fig. 8 Distribution of the intensity estimated based on the finite-dynamic source models (a) Based on the inverted results of the seismic data; (b) Based on the jointly inverted results of the seismic data and the InSAR data. The red stars refer to as the epicenter locations of the mainshock. The white solid curves are active faults, and the Ⅸ-grade areas are fenced with the black isoseismal contours. |
图 8b展示了基于联合反演结果计算的烈度分布.可以看出，总体特征与基于地震资料反演结果计算的烈度非常相似.但是，细节上却有明显不同.首先，Ⅸ区面积增加，但其他区域明显减小；第二，在东北和西南方向上各等级的烈度区都有不同程度的收缩.这些细节上的变化都是由于破裂区的空间分布的改变以及破裂方向的改变所致.
与野外考察的烈度相比(中国地震局，2017)，估计烈度的总体特征与其基本相同，但无论是基于地震资料反演结果计算的烈度或是基于联合反演结果计算的烈度，其Ⅸ度区明显较小.另外，估计烈度区的非椭圆特征较为明显.
5 讨论 5.1 关于主震震源位置根据中国地震台网测定，北京时间2017年8月8日21时19分于我国四川省九寨沟县发生了MS7.0级地震，震中位于33.20°N，103.82°E，震源深度20 km.美国地震调查局(USGS)也发布了震源位置参数，哈佛大学发布了GCMT的矩心位置(表 1，图 1).毋庸置疑，中国地震台网利用地方和区域台站确定的震源位置可靠性最高.但是，震源深度素来不确定性较大，为此，我们在波形反演过程中考虑震源深度的不确定性，并利用波形反演将震源深度重新确定为11 km.事实上，这个深度和中国地震台网测定的震中位置不但可以较好地解释远场地震波形，还可以很好地解释近场InSAR地表位移，这说明选择这个点作为地震起始破裂点是恰当的.
5.2 关于发震断层参数根据多家研究机构发布的震源机制解(表 2)和余震的空间分布特征，不难判定发震断层走向为南东(150°~153°)，倾角约80°(78°~84°).根据震后13 h的余震精定位结果，发震断层面的走向为142°，倾角为81°.考虑到震源机制解为点源等效力对应的结果，而基于余震空间分布的拟合结果以余震发生在主震断层面及其附近的假设为前提，同时考虑到利用走向为153°而倾角为84°的平面断层对远场地震波和近场InSAR位移进行了很好的解释，我们不妨认定发震断层面为走向153°而倾角84°的等效平面.
5.3 关于断层面上的位错分布从上几节可以看出，建立相同的初始平面模型，单独利用地震资料和InSAR资料反演得到的位错分布具有明显的不同.根据地震资料反演的结果，位错主要分布在震源周围，呈现单独的圆盘形.根据InSAR资料反演的结果，位错可以分为两个区域，一个在起始破裂点左侧(西北方向)，另一个在右上侧(东南方向).目前，普遍认为基于单一资料的反演结果是片面的，而多种资料的联合反演结果更具全面性(Zhou et al., 2004)，更容易被人接受.这里地震资料和InSAR资料联合反演的结果表明，地震的破裂区可分为两个主要的凹凸体，但不同于InSAR资料的反演结果.由于联合反演的结果不但较好地(相关系数达0.81) 解释了远场波形资料，还很好地(相关度达0.95) 解释了近场InSAR资料，所以我们认为联合反演得到的位错分布为2017年九寨沟MS7.0级地震的位错分布.
5.4 关于震源时间函数和标量地震矩无论是地震资料的反演或是地震资料和InSAR资料的联合反演均表明，地震过程持续约15 s，可以分成两次事件，第一次事件持续时间不到10 s, 但释放了大部分标量地震矩，第二次事件持续5s左右，释放标量地震矩较少.整个过程释放标量地震矩6.61×1018 Nm左右，相当于矩震级MW6.5.
5.5 关于地震烈度从前面的反演结果看，单独利用地震资料和单独利用InSAR资料的反演结果截然不同.同时考虑到近场InSAR位移并不是严格的同震位移，由此反演得到的位错模型也不是严格的同震位错模型，所以，基于地震资料反演结果和联合反演结果分别计算了两种情况下的烈度分布.这里估计的烈度分布总体特征与野外考察获取的烈度基本一致，但极震区(Ⅸ度)明显不同.这种不同一方面与我们使用的简单平面断层有关，另一方面与野外考察获取的信息不全有关.需要说明的是，截止目前野外考察人员还未能进入极震区，这里的烈度还有待进一步的科学考察.
5.6 有待深入研究和考证的问题尽管，借助于远场地震资料、近场InSAR资料和震后13 h内的地方和区域台网的震相资料关于主震的震源深度、发震断层参数、断层面上的位错分布、震源的总体破裂历史、地震释放的标量地震矩以及造成的地震烈度得到了上述认识，但仍有一些问题需要进一步的研究和考证.
第一，地震断层是否为一个简单的平面断层？地震区或主要的余震分布区属于无人区，截止目前没有任何关于地表破裂的资料或信息.余震震中在空间上粗略地呈现“线形”排列.但是，如果对余震横断面(BB′)的细节进行分析，可以发现余震并非发生在一个平面上.如果对余震序列进行时窗滑移分析，发震断层很可能由三段组成，而且不在一条直线上.如图 9a所示，震源深度大于10 km的余震都集中在一起，但是小于10 km的余震清楚地分为三支；如图 9b所示，震后1.5 h至5 h时间窗内的余震可以清晰地分为三段.另外，我们没能利用单一机制的格林函数从远场地震波形中提取到可接受的视震源时间函数，且滑动角可变的平面破裂过程模型也只能将观测波形解释到81%的程度(平均相关系数为81%).这些现象表明震源机制具有时空变化，发震断层很可能不是一个平面.
|
图 9 余震空间分布特征 (a)震后13 h的余震在断层横断面BB′上的投影(参看图 2).似乎可以分成三个平面; (b)震后1.5 h至5 h时段内的余震震中分布.似乎可以分成3段.红色五角星为主震. Fig. 9 Spatial feature of the aftershocks (a) Relocated aftershocks within the first ~13 hours projected on the profile BB′ (see Fig. 2 for details). There appear three fault planes being isolated; (b) Epicenter locations of the relocated aftershocks within the first 1.5-to-5.0 hours. There appear three segments being linked. Red stars refer to as the mainshock. |
第二，位错的时空分布特征究竟如何？可以看出, 基于远场地震资料、近场InSAR资料以及两种资料的反演结果，地震的时空破裂过程明显不同.这种差别决不是用近场资料和远场资料的特点差异所能解释.由此引发的第一个疑问是近场InSAR资料揭示的位移是否为同震位移.如表 2所示，InSAR资料的时间窗为2017年7月30日至2017年8月11日.暂且不论震前8天，但震后应力应变会发生急剧调整.震后3天主震断层面发生了多大变化、震区地表变形发生了多大变化均不得而知.另外，在本研究中我们使用一个平面模型，然而，正如上一段的讨论，发震断层很可能有多个平面组成.因此，本研究得到的时空破裂过程仅为等效的平面断层模型.关于震源破裂过程最终模型还需要基于更多资料的更详细的研究.
6 结论基于远场地震资料、近场InSAR资料以及震后数小时的余震震相数据，关于2017年九寨沟MS7.0级地震可以得到如下认识：(1) 主震发生于33.20°N/103.82°E，但震源深度很可能为11 km左右.(2) 发震断层走向~150°，倾角~80°，但很可能是多段多层破裂.(3) 破裂区可分为两个主要的凹凸体.一个较大，距起始破裂点较近，从起始破裂点北西较深的位置延展到起始破裂点南东较浅的位置，似乎破裂到了地表.此凹凸体主要以走滑为主.另一个凹凸体较小，距起始破裂点较远，位于起始破裂点北西15 km且6 km深处.此凹凸体的正滑分量非常明显.但这一认识仅基于单一的平面断层的联合反演.(4) 地震过程持续约15 s，可以分成两次事件，第一次事件持续时间不到10 s, 但释放了大部分标量地震矩，第二次事件持续5 s左右.整个过程释放标量地震矩6.61×1018 Nm左右，相当于MW6.5.(5) 基于地震资料的反演结果以及地震资料和InSAR资料的联合反演结果估计的烈度分布，总体特征与野外考察得到的烈度特征相似，但极震区的烈度明显不同，有待进一步考证.最后，种种迹象表明发震过程很可能更为复杂，还需要后续基于更多资料的详细研究.
致谢感谢欧洲空间局(ESA)无偿共享Sentinel-1雷达数据.远场地震波资料来自IRIS数据中心，余震震相数据由李艳娥副研究员提供.
Antolik M, Dreger D S.
2003. Rupture Process of the 26 January 2001 MW7.6 Bhuj, India, Earthquake from Teleseismic Broadband Data. Bull. Seismol. Soc. Am, 93(3): 1235-1248.
DOI:10.1785/0120020142 |
|
Chen Y T, Xu L S.
2000. A time-domain inversion technique for the tempo-spatial distribution of slip on a finite fault plane with applications to recent large earthquakes in the Tibetan Plateau. Geophys. J. Int, 143(2): 407-416.
DOI:10.1046/j.1365-246X.2000.01263.x |
|
China Earthquake Administration. 2017. http://www.cea.gov.cn/publish/dizhenj/464/476/20170812212207031213560/index.html[2017-09-08].
|
|
Fang L H, Wang W L, Yang T. 2017. http://www.cea-igp.ac.cn/tpxw/275883.html[2017-09-08].
|
|
Feng W, Omari K, Samsonov S V. 2016. An automated insar processing system:potentials and challenges. 2016 IEEE International Geoscience & Remote Sensing Symposium, Beijing, China:IEEE, 3209-3210, doi:10.1109/IGARSS.2016.7729830.
|
|
Feng W, Tian Y, Zhang Y, et al.
2017a. A Slip Gap of the 2016 MW6.6 Muji, Xinjiang, China, Earthquake Inferred from Sentinel-1 TOPS Interferometry. Seismol. Res. Lett, 88(4): 1054-1064.
DOI:10.1785/0220170019 |
|
Feng W, Samsonov S, Tian Y, et al.
2017b. Surface deformation associated with the 2015 MW8.3 Illapel earthquake revealed by satellite-based geodetic observations and its implications for the seismic cycle. Earth Planet. Sci. Lett, 460: 222-233.
DOI:10.1016/j.epsl.2016.11.018 |
|
Gusman A R, Murotani S, Satake K, et al.
2015. Fault slip distribution of the 2014 Iquique, Chile, earthquake estimated from ocean-wide tsunami waveforms and GPS data. Geophys. Res. Lett, 42(4): 1053-1060.
DOI:10.1002/2014GL062604 |
|
Hartzell S H, Heaton T H.
1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake. Bull. Seismol. Soc. Am, 73(6A): 1553-1583.
|
|
Kennett B L N, Engdahl E R, Buland R.
1995. Constraints on seismic velocities in the Earth from traveltimes. Geophys. J. Int, 122(1): 108-124.
DOI:10.1111/gji.1995.122.issue-1 |
|
Lay T, Ammon C J, Hutko A R, et al.
2010. Effects of Kinematic Constraints on Teleseismic Finite-Source Rupture Inversions:Great Peruvian Earthquakes of 23 June 2001 and 15 August 2007. Bull. Seismol. Soc. Am, 100(3): 969-994.
DOI:10.1785/0120090274 |
|
Luo Y J, Zhao L. 2017. http://tec.earth.sinica.edu.tw/new_web/upload/news/EQfile/2017-08-08_M6.5.pdf[2017-09-08].
|
|
Shan X J, Qu C Y, Gong W Y, et al. 2017. http://www.eq-igl.ac.cn/upload/files/2017/8/15151711397.pdf[2017-09-08].
|
|
Simons M, Fialko Y, Rivera L.
2002. Coseismic Deformation from the 1999 MW7.1 Hector Mine, California, Earthquake as Inferred from InSAR and GPS Observations. Bull. Seismol. Soc. Am, 92(4): 1390-1402.
DOI:10.1785/0120000933 |
|
Wald D J, Quitoriano V, Heaton T H, et al.
1999a. Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California. Earthquake Spectra, 15(3): 557-564.
DOI:10.1193/1.1586058 |
|
Wald D J, Quitoriano V, Dengler L A, et al.
1999b. Utilization of the Internet for rapid community intensity maps. Seismol. Res. Lett, 70(6): 680-697.
DOI:10.1785/gssrl.70.6.680 |
|
Waldhauser F, Ellsworth W L.
2000. A double-difference earthquake location algorithm:method and application to the northern hayward fault, California. Bull. Seismol. Soc. Am, 90(6): 1353-1368.
DOI:10.1785/0120000006 |
|
Wang R J.
1999. A simple orthonormalization method for stable and efficient computation of Green's functions. Bull. Seismol. Soc. Am, 89(3): 733-741.
|
|
Wang R J, Martín F L, Roth F.
2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust-FORTRAN programs EDGRN/EDCMP. Computers & Geosciences, 29(2): 195-207.
|
|
Wang W M, He J K, Hao J L, et al. 2017. http://www.itpcas.ac.cn/xwzx/zhxw/201708/t20170809_4840737.html[2017-09-08].
|
|
Ward S N, Barrientos S E.
1986. An inversion for slip distribution and fault shape from geodetic observations of the 1983, Borah Peak, Idaho, Earthquake. J. Geophys. Res, 91(B5): 4909-4919.
DOI:10.1029/JB091iB05p04909 |
|
Xu L S, Chen Y T, Teng T L, et al.
2002. Temporal-Spatial Rupture Process of the 1999 Chi-Chi Earthquake from IRIS and GEOSCOPE Long-Period Waveform Data Using Aftershocks as Empirical Green's Functions. Bull. Seismol. Soc. Am, 92(8): 3210-3228.
DOI:10.1785/0120010173 |
|
Xu L S, Zhang X, Wei Q, et al.
2016. A method for estimating the earthquake intensity caused by a finite-dynamic source. Chinese J. Geophys, 59(10): 3684-3695.
DOI:10.6038/cjg20161015 |
|
Yagi Y, Mikumo T, Pacheco J, et al.
2004. Source Rupture Process of the Tecomán, Colima, Mexico Earthquake of 22 January 2003, Determined by Joint Inversion of Teleseismic Body-Wave and Near-Source Data. Bull. Seismol. Soc. Am, 94(5): 1795-1807.
DOI:10.1785/012003095 |
|
Zhang X. 2016. Study on new methods for analysis of the complexity of source rupture process based on apparent source time functions[Ph. D. thesis] (in Chinese). Beijing:Institute of Geophysics, China Earthquake Administration.
|
|
Zhang Y, Feng W, Chen Y T, et al.
2012. The 2009 L'Aquila MW6.3 earthquake:a new technique to locate the hypocentre in the joint inversion of earthquake rupture process. Geophys. J. Int, 191(3): 1417-1426.
|
|
Zhang Y, Xu L S, Chen Y T. 2017. http://www.cea-igp.ac.cn/tpxw/275883.html[2017-09-08].
|
|
Zhou S, Irikura K, Chen X.
2004. Analysis of the reliability and resolution of the earthquake source history inferred from waveforms, taking Chi-Chi earthquake as an example. Geophys. J. Int, 157(3): 1217-1232.
DOI:10.1111/gji.2004.157.issue-3 |
|
房立华, 王未来, 杨婷. 2017. http://www.cea-igp.ac.cn/tpxw/275883.html[2017-09-08].
|
|
罗翊菁, 赵里. 2017. http://tec.earth.sinica.edu.tw/new_web/upload/news/EQfile/2017-08-08_M6.5.pdf[2017-09-08].
|
|
单新建, 屈春燕, 龚文瑜等. 2017. http://www.eq-igl.ac.cn/upload/files/2017/8/15151711397.pdf[2017-09-08].
|
|
王卫民, 何建坤, 郝金来等. 2017. http://www.itpcas.ac.cn/xwzx/zhxw/201708/t20170809_4840737.html[2017-09-08].
|
|
许力生, 张旭, 魏强, 等.
2016. 一种基于有限动态源的烈度估计方法. 地球物理学报, 59(10): 3684–3695.
DOI:10.6038/cjg20161015 |
|
张旭. 2016. 基于视震源时间函数的震源过程复杂性分析新方法研究[博士论文]. 北京: 中国地震局地球物理研究所.
|
|
张勇, 许力生, 陈运泰. 2017. http://www.cea-igp.ac.cn/tpxw/275883.html[2017-09-08].
|
|
中国地震局. 2017. http://www.cea.gov.cn/publish/dizhenj/464/476/20170812212207031213560/index.html[2017-09-08].
|
|