据美国地质调查局测定：2016年3月2日20时49分(北京时)左右在印尼苏门答腊岛海域发生7.8级地震(以下简称印尼地震)，震中位于东经94.330°/南纬4.952°, 深度约24 km.印尼地震是印度板块和澳大利亚板块共同作用的结果，震中位于印尼苏门答腊以西海域，沃顿盆地(Wharton Basin)以北的离散变形带之间(图 1).震中区附近300 km内历史地震稀少，但沃顿盆地北部和东部发生过多次有巨大影响的地震，包括2004年12月26日M_W9.0逆冲型地震和2012年4月11日M_W8.6走滑型事件.其中，2012年M_W8.6地震被认为是有地震记录以来全球最大的走滑型地震之一.该事件是一次发生在板块内的大型走滑型事件.其发震断层面相互交错，破裂过程异常复杂(Yue et al., 2012; Meng et al., 2012; Wang et al., 2012; Wei et al., 2013).2016年印尼地震是在2012年M_W8.6地震之后，时隔4年又一次7.5级以上地震，引人注目.

图 1

Fig. 1

图 1 印尼地震基于点源模型的震源机制解和研究区地质构造.选取本研究区内自1976年以来历史地震矩张量资料, 数据来源于哥伦比亚大学Global CMT目录.震源球的颜色表示震源的深度. Fig. 1 The focal mechanism solution of the Indonesia earthquake and tectonic setting in the studied area. The historical centroid moment tensor catalog has been chosen since 1976, which is derived from the Global CMT catalog provided by Columbia University. The color of the beachball indicates its hypocenter depth.

印尼地震发生后，国内多家机构和学者对震源机制解和破裂过程进行快速测定，服务于大震应急.其中，本文作者采用W震相(W phase)波形反演方法，震后约50分钟正式发布此地震震源机制解结果，其中断层面节面1参数为5°/倾角88°/滑动角4°，判定为一次走滑型为主事件.张勇在震后1 h左右给出了破裂过程结果.上述结果已在中国地震信息网上发布(http://www.csi.ac.cn/manage/eqDown/05LargeEQ/201603022049M7.8/zonghe.html)，为地震海啸预警和科学抗震减灾提供重要参考.由于印尼地震震后数小时内余震较少(截止4月22日记录到20次较大余震，震级范围M4.1-5.6)，且震源区地质构造信息尚不清楚，因此，有关印尼地震真实发震断层、破裂方式以及是否存在超剪切波速破裂等是值得深入研究的问题.

考虑到破裂过程反演结果的可靠性依赖于反演前断层面几何参数、破裂速度及破裂持续时间等参量(或范围)设置，因此，本文首先根据W震相测定矩张量方法，获得比较可靠的震源机制解，用来约束初始断层面模型.然后，根据反投影分析获得破裂速度和破裂持续时间较合理范围，为获得可靠细致的有限断层模型提供更好的约束.

本文在前人工作基础上，根据中国和全球地震台网记录的波形数据，采用W震相矩张量测定、反投影分析及有限断层模型反演方法，详细地研究了2016年印尼M_W7.8地震动态破裂时空过程，据此分析讨论印尼地震震源运动学特征.本研究对增进海洋地震震源特性的了解和认识，探讨本区域岩石圈大型转换断层孕震背景等有着重要意义.

本文使用全球和中国地震台网记录的远震波形，采用基于Kanamori和Rivera(2008)提出了基于W震相测定质心矩张量的方法，求解印尼地震的震源机制解.W震相是在S波之前到达的一种较明显的长周期波(100~1000 s)，最早于1992年日本记录的尼加拉瓜海啸地震的位移波形中被辨识出(Kanamori，1993).根据地震射线理论，在选定初至P波开始至15Δ秒(Δ表示震中距，单位为(°))时窗中，W震相可解释为长周期P、PP、SP和S等多个震相的叠加.采用的长周期W震相在测定大震震源参数优势主要有两点：(1)W震相传播速度快，在周期100~1000 s范围，群速度为4.5~9.0 km·s^-1，明显高于面波，因此相比于传统的面波波形反演方法，更适合大地震快速响应；(2)对远震台，W震相为长周期的体波，主要在地幔中传播，浅层大陆地壳结构对其影响不大.因此，此方法逐渐被美国国家地震信息中心、美国太平洋海啸警报中心及中国地震台网中心等国内外多所机构所使用(Hayes et al., 2009；Duputel et al., 2012；赵旭等，2014).本文基于一维速度模型PREM，采用简振正型叠加法计算远震(5°~85°)W震相的理论波形.

我们从中国地震台网中心和IRIS收集中国和全球地震台网记录的震中距5°~ 85°范围内地震波形，数据预处理主要步骤：(1)波形质量控制，舍去信噪比不高和断点等“坏”波形；(2)重采样，将不同采样率的波形数据重采样为1 s；(3)时间域反卷积，将原始的速度记录转化为位移记录；(4)采用4阶巴特沃斯(Butterworth)滤波器带通滤波，滤波频段为150~600 s；(5)选用初至P波开始至15Δ秒时间窗提取W震相.

W震相矩张量反演采用的主要技术思路类似与Dziewonski等(1981)和Ekström等(2005)提出的标准质心矩张量(Centroid moment tensor)反演方法，但两种方法有着四方面不同：(1)选用地震波时窗；(2) W震相矩张量反演使用长周期地震波形；(3)前者采用时间域反卷积方法，能对时间序列每个抽样点进行实时处理，使得W震相可在面波限幅的记录中无失真恢复.从而该方法可充分使用强震近震(震中距5°~30°)波形.而后者主要使用频率域反卷积方法提取位移波形；(4)确定质心水平位置和质心深度的反演算法.

W震相反演未知震源参数向量m可由式(1)表示：

(1)

其中：地震矩张量f=[M_rr, M_pp, M_tt, M_rp, M_rt, M_pt]^t和表征质心时空坐标η_c=[θ_c, ϕ_c, h_c, τ_c]^t，其中θ_c, ϕ_c, h_c, τ_c分别表示质心余纬度、质心经度、质心深度以及质心相对发震时刻的时间偏移，共10个未知参数.t表示转置运算符.

为了确定最佳的质心时空坐标η_c，可通过网格搜索的方法使定义的误差函数(见式2)最小：

(2)

其中向量s_W表示W震相的理论波形，向量d_W为W震相的实际观测.

计算中最终选用了震中距16° ~84°内68个台站共104道波形(图 2d，限于篇幅仅显示部分台站拟合波形).经反演后，最后获得实际观测波形与理论波形拟合最优的矩张量解，如表 1所示.图 2d给出实际观测波形和理论波形对比.大多数台站的实际观测波形与理论波形拟合较好，波形拟合误差均方根值(RMS)仅为0.016 mm.

图 2

Fig. 2

图 2 反演结果和波形拟合情况 (a)反演获得的质心水平位置最佳解.红色三角形表示求得的质心水平位置，黑色十字为美国地质调查局测定的震中位置，不同颜色不同大小圆圈为相对于最小波形拟合误差均方根值的百分比; (b)质心时间偏移τc最佳解(红色三角).横坐标表示质心时间偏移(单位：s)，纵坐标表示相对于最小波形拟合误差均方根值的百分比; (c)质心深度最佳解(红色三角).横坐标表示质心深度(单位：km)，纵坐标表示相对于最小波形拟合误差均方根值的百分比; (d)实际观测波形(黑线)和理论波形(红线)对比(仅显示部分台站).波形附近两个红色圆圈表示选用的时间窗.图上:台网名、台站名、位号、分向、方位角及震中距/°.图左:振幅值/mm.图右:震中(蓝色五角星)和台站(红色圆圈)分布. Fig. 2 Inverted results and waveform comparisons (a) The optimal centroid horizontal location was obtained from grid searching. The triangle denotes the inverted centroid horizontal location, and the cross is the epicentral location determined by the USGS. The different circles show the percentage relative to the RMS of the minimum waveform fit error; (b) The optimal centroid time shift τc was obtained from the grid searching. The horizontal axis indicates the centroid time shift τc. The vertical axis shows the percentage relative to the RMS of the minimum waveform fit error; (c) The optimal centroid depth was obtained from the grid searching. The horizontal axis indicates the centroid depth. The vertical axis shows the percentage relative to the RMS of the minimum waveform fit error; (d) Comparison of observed (black line) with synthetic waveforms (red line) from part of the used stations. The red circles near each waveform denote the selected time window. At each top, they are network name, station name, location, channel, azimuthal angle and epicenter distance respectively. On the left, it shows amplitude(mm). On the right, the blue star denotes the epicenter and the red circle represents the corresponding station.

表 1 Table 1

表 1 印尼地震震源参数对比 Table 1 Comparison of source parameters of the Indonesia earthquake 结果来源 MW 质心水平位置/(°) 深度/km 矩张量/1020Nm 节面Ⅰ 节面Ⅱ Mrr Mtt Mpp Mrt Mrp Mtp 走向/倾角/滑动角 走向/倾角/滑动角 本研究 7.8 4.85S/94.43E 21 0.29 -1.14 0.85 -0.52 -0.19 -5.93 5°/85°/-2° 95°/88°/-175° Global CMT 7.8 4.75S/94.22E 37.5 0.16 -1.11 0.94 0.01 0.29 -5.83 5°/90°/3° 275°/87°/180° 美国地质调查局 (体波测定) 7.8 _ 24 - - - - - - 184°/90°/-5° 275°/85°/-180°

表 1 印尼地震震源参数对比 Table 1 Comparison of source parameters of the Indonesia earthquake

结果表明：印尼地震的总标量地震矩M₀为6.04×10²⁰Nm，矩震级约M_W7.79，最佳双力偶解的节面Ⅰ和节面2(走向/倾角/滑动角)为5°/85°/-2°、95°/88°/-175°.通过网格搜索最佳波形拟合的质心时间偏移τ_c约16 s，见图 2b.假设质心时间偏移τ_c能表征震源破裂持续时间的二分之一(Kanamori and Rivera, 2008)，推断震源破裂大多数能量释放的持续时间约32 s.质心水平位置(东经94.43°/南纬4.85°，图 2a)，位于破裂起始点东北.通过网格搜索10~70 km范围深度，获得波形拟合误差最小的深度约21 km，浅于破裂起始点的深度.

本文采用中国地震台网记录的远震P波波形，基于反投影法，分析印尼地震破裂过程.本文使用Wang等(2016)提出的方法，其主要思路为：在选定某个时窗，通过波形互相关技术，叠加与某一可能的震源位置所对应的波形，以抵消噪音和次声波之影响，进而放大震源处释放的有用信号，然后将叠加的能量投影至相对应的震源位置.对震源区所有可能位置皆作能量后投影后，即可获得此时窗内的能量辐射图像.采用全球一维速度模型IASPEI 1991(Kennett, 1991)计算每个台站记录的P波理论到时.

本文从中国地震台网选用152个台站记录的高质量垂直向P波波形(图 3a)，震中距范围为27°~48°，方位角范围7°~30°.选用的台站主要分布在河北、河南、湖北、四川及云南等地，记录的原始波形采样率皆为100 Hz.使用巴特沃斯滤波器带通滤波，最终使用的波形频段范围为0.5~2.0 Hz.在实际计算中，对各个台站的P波前6 s波形与参考地震台(成都台，台站代码CD2)做相关运算以对齐所有台站的初始波形.破裂初始点选用美国地质调查测定的结果，我们分别采用10 s的时窗长度，2 s的滑动窗对200 s的数据进行计算.

图 3

Fig. 3

图 3 基于中国地震台网数据的反投影分析结果 (a)台站分布；(b)地震破裂迁移及能量叠加图；(c)叠加的高频能量释放值随时间变化曲线图. Fig. 3 The back projection analysis results obtained from seismic waveforms observed by the China Digital Seismic Network (a) The station distribution; (b) The rupture distribution and the cumulative stacked energy during this earthquake; (c) The release value of stacked high-frequency energy varies with time.

通过反投影分析，我们得到了印尼地震辐射高频能量的震源位置分布(图 3b).不同的圆圈大小表示每个时窗最大振幅相对应的震源位置，圆圈的颜色表示不同的破裂时刻.根据中国地震台网计算的结果表明：印尼地震最大能量释放区域位于震中以北.地震发生后，破裂主要向北北东向扩展，至18 s左右，破裂方向未较大改变，一直延伸至最北端.18~40 s期间，破裂在最北端附近沿着南东至北西方向徘徊，但扩展区域面积并不大.从图 3c可以看出，整个破裂总时间时间小于40 s，能量峰值出现在震后20 s左右.通过测定破裂长度与破裂时间大小，估计破裂平均速度约2.0 km·s^-1.

为了更细致地了解印尼地震破裂过程信息，本文联合中国和全球地震台网记录的远场体波和面波波形，基于有限断层模型，重建此地震时空破裂图像.根据波形资料的信噪比和方位角覆盖情况，在震中距30°~90°范围内选取了30个P波和18个SH波位移数据(滤波频段0.0033 Hz至1 Hz).为了更准的对齐数据，人工重新拾起了P波到时.另外，增加了49个长周期的面波资料(周期范围为167~333 s)，以更好约束此地震总破裂方式和标量地震矩大小.

本文采用一种基于小波变换的波形反演方法(姚振兴和纪晨，1997; Ji et al., 2002; 王卫民等，2008; Hao et al., 2013)，利用模拟退火算法求解每个子断面滑动量、滑动角、破裂开始时间和震源时间函数.该方法主要优点有两方面：(1)采用小波变换技术能在时域和频域对地震波形数据分析，能够更有效地提取震源破裂过程信息; (2)联合利用远场体波和面波，强地面运动等多种地震观测资料，可重建更可靠的震源破裂过程物理图像.在实际计算中，我们结合质心水平位置和高频能量移动轨迹，能够确定断层面(走向5°/倾角85°)为实际破裂面.经过多次模拟计算，最终建立一个沿走向108 km和沿倾向38.5 km的断层面.考虑震中附近地形和海水层厚度的影响，设定断层面上端距地表约5.2 km.将断层面划分为126个6.0 km×5.5 km的子断层, 每个子断层的滑动量搜索范围为0~14 m.基于前文获得此地震的震源机制，在此将滑动角范围设定为-30°~30°, 间隔为3°.破裂速度允许在0.5~3.5 km·s^-1范围内变化, 每个子断层的上升时间为0.6~6 s, 间隔为0.6 s.

本文中印尼地震震源区地壳结构采用了CRUST2.0模型, 地壳以下使用PREM模型，基于前面设定反演参数, 我们反演得到了此次印尼7.8级地震震源破裂过程图像(图 4).结果表明：此地震为一次对称的双侧破裂事件，破裂分量主要以走滑为主.地震在深度24 km处开始破裂，最大的滑动量位于初始破裂点北东，沿着断层面走向约30 km处，大小高达11 m左右.反演得到的标量地震矩约为6.19×10²⁰ Nm，换算成矩震级M_W为7.79，该结果与点源模型获得的地震矩大小(6.04×10²⁰ Nm)吻合.整个破裂总持续时间约35 s，其中主要能量在5~25 s内释放(见图 4c)，约占总能量的97%，换算成矩震级M_W约7.5.在地震发生后约2 s内，破裂扩展缓慢.2 s后，破裂以更快的速度，沿北东和南西方向向两侧同时扩展.沿断层走向北东向，破裂在5~25 s内，在深度5~21 km处形成高滑动能量集中带.在断层南端，破裂在15~20 s内也集中释放较多应变能.整个断层面上的滑动分布并不均匀.在震源上方浅部区域，沿走向方向，形成了震中以北约50 km和以南约20 km长的高滑动量集中带.破裂平均速度约2.2 km·s^-1, 约为起始破裂点深度处剪切波传播速度(4.5 km·s^-1)的一半.此次地震破裂最终形成了总长度约90 km的断层.从实际观测波形和理论预测波形拟合情况来看, 如图 5所示, 无论体波和面波，波形拟合情况皆很好, 吻合度很高，说明结果可信度比较高.

图 4

Fig. 4

图 4 印尼地震最佳的同震滑动位移模型 (a)滑动位移在地面上的投影.黄色五角星示意震中位置，红色五角星代表本文获得的质心水平位置, 红色震源球表示4个较大余震(MW4.8~5.7)，其震源机制解资料来源于Global CMT目录；(b)同震滑动位移沿走向和深度的分布.黄色五角星表示破裂初始位置(94.330°E/4.952°S/深度24千米)，不同颜色代表滑动量大小，白色箭头表示上盘相对下盘的滑动方向.黑色等值线表示破裂传播时间；(c)地震矩率随时间变化. Fig. 4 The preferred coseismic slip models of the Indonesia earthquake (a) The surface projection of the slip distribution. The rupture initiation is located at the point (94.330°E/4.952°S/) with a depth of 24 km. The red star indicates the horizontal centroid location obtained by the centroid moment tensor inversion. The red beachballs indicate large aftershocks with their source mechanism solutions derived from the Global CMT catalog; (b) The coseismic slip is distributed along the strike and depth directions, respectively. The yellow star shows the hypocenter. The color indicates the slip amplitude, and white arrows indicate the motion direction of the hanging wall relative to the footwall. Black contour lines indicate the rupture initiation time; (c) The moment rate varies with time.

图 5

Fig. 5

图 5 所用体波和面波的实际观测波形和理论预测波形比较 (a)实际观测体波记录(黑线)和理论波形(灰线)对比.滤波频段为0.0033 Hz至1 Hz.台站名、震中距、方位角和波形记录的最大幅度(单位：μm)分别表示在每列的左边、中下、中上和右边；（b)实际位移面波(4~ 6 mHz，黑线)和理论预测波形(灰线)比较.台站名、震中距、方位角和波形记录的最大幅度(单位：mm)分别表示在每列的左边、中下、中上和右边. Fig. 5 Waveform comparison of teleseismic observations and synthetics predicted by our preferred model (a) Comparison of the teleseismic body-wave displacement observations (black lines) and synthetics (grey lines). These data were band-pass filtered from 0.0033 Hz to 1 Hz. The black lines indicate the data, whereas the grey lines indicate the synthetics. Both the data and synthetics are normalized to the peak amplitude of the data in μm, which is shown in the end of each trace. The numbers above and below the beginning of each trace are the source azimuth and epicentral distance, respectively, (b) Comparison of long-period (4~6 mHz) surface wave observations and synthetics. The peak displacement of the observation in mm is shown above the end of each trace

本文使用中国和全球地震台网记录的波形数据，采用W震相波形矩张量反演、反投影分析及有限断层模型反演方法，研究了2016年苏门答腊岛海域印尼7.8级地震动态破裂时空过程.质心矩张量反演结果表明此地震是一次纯走滑型事件，质心深度约21 km.获得的震源机制解结果同Global CMT给出结果比较，两者基本吻合，但质心深度相差约16 km(见表 1).有限断层模型反演结果表明印尼地震破裂时在深度5~21 km处形成高滑动量集中区，质心深度约18 km.该结果与本文矩张量反演的深度21 km更为吻合.另外，震后3天，4个较大余震(M_W4.8~5.7)皆发生在主震上方深度为10~14 km浅部区域，也佐证本文结果的可信度.

质心时间偏移τ_c用于研究甄别具有特殊震源性质的“异常”地震(如海啸型地震)，也是估算震源破裂持续时间大小有效手段之一(Duputel et al., 2013).在近似条件下，震源破裂持续时间可用质心时间偏移τ_c的2倍来表征(Kanamori and Rivera, 2008).本文通过网格搜索，获得最佳波形拟合的质心时间偏移τ_c约16 s，推断震源破裂大多数能量释放的持续时间约32 s.该结果同反投影分析和有限断层模型反演破裂持续时间结果吻合.根据Duputel等(2012)提出的质心时间偏移τ_c和地震矩M₀统计关系：τ_c=1.2×10^-8×M₀^1/3，其中M₀单位dyne·cm，由本文反演的标量地震矩估算的质心时间偏移(约22 s)，是本文获得质心时间偏移τ_c1.4倍.相比于相同规模的地震，研究结果表明印尼地震破裂持续时间较短.

破裂过程反演结果表明：此地震为一次对称的双侧破裂事件，北北东─南南西向的断层节面(走向5°/倾角85°)为其发震断层面.标量地震矩约为6.19×10²⁰ Nm，换算成矩震级为7.79，最大的滑动量约11 m, 位于初始破裂点北东，沿着断层走向约30 km处.破裂总持续时间约35 s，破裂在5~25 s内释放大多数能量，约占总能量的97%.最终形成了总长度90 km左右断层.该结果同国外学者发表的结果比较一致(Lay et al., 2016).相对于其他位于大型转换型断层附近的海洋型大震(M_W≥7.5)，例如，2013年11月17日M_W7.8 Scotia Ridge地震破裂尺度长达200 km(Ye et al., 2014)，以及1906年美国旧金山M7.7~7.9地震破裂尺度300~470 km(Wald et al., 1993; Thatcher et al., 1997)，印尼地震破裂尺度显得比较小，分别仅为前两个地震破裂尺度的45%和25%.说明印尼地震破裂形成的高滑动能量带相对集中.

通过反投影分析和有限断层模型反演，表明印尼地震破裂平均速度约2.0~2.2 km·s^-1，约为起始破裂点深度处剪切波传播速度(4.5 km·s^-1)的一半.与同样位于沃顿盆地以北，距本地震约700 km的2012年M_W8.6地震相比，两者破裂速度相差较大.已有的研究结果表明：2012年M_W8.6地震破裂过程中可能存在超剪切波速破裂现象，多个断层面上的平均破裂速度高达5 km·s^-1，超过了本地的剪切波速(Wang et al., 2012).另外，两个地震破裂模式也不同.2012年M_W8.6地震破裂过程异常复杂，破裂在震源区附近多个相互交错的断层上扩展，被认为是一次发生板块内共轭断层上的大型走滑型事件(Yue et al., 2012; Meng et al., 2012; Wang et al., 2012; Wei et al., 2013).相比于前者，2016年M_W7.8地震破裂过程并不太复杂，为一次对称的双侧破裂事件.破裂在前18 s主要在北北东─南南西向的断层上扩展.反投影分析结果表明：破裂约18 s后，其方向似乎发生了一些变化.此时，破裂由开始北北东向扩展，逐步转变为朝南东至北西方向徘徊，但扩展区域面积并不大，推测2016年M_w7.8地震并未明显在第二个交错的断层发生较大破裂.

本文通过地震波形分析，深入地研究了印尼地震动态破裂过程，探究其运动学特征，对增进本区域岩石圈海洋型地震的震源特性的了解和认识，具有重要参考意义.

谨此祝贺姚振兴先生从事地球物理教学科研工作60周年.文中数据处理主要使用地震信号处理通用软件SAC2000.图件制作采用GMT软件包.文中高质量的波形数据来源于中国地震台网中心地震台网部和IRIS.王墩教授在讨论中给予宝贵建议，在此一并感谢.最后感谢两位评审专家提出的建议.