2016年4月14日以来，日本九州熊本县接连发生多次强烈地震，4月14日12时26分 (UTC)，发生了M_W6.1地震；4月15日16时25分，再次发生M_W7.1地震.日本气象厅 (Japan Meteorological Agency，简称JMA) 已确定4月15日发生的地震为熊本地震序列的“主震”，而在4月14日及之后发生的地震为“前震”.此次地震序列前震和主震的最大烈度 (日本称震度) 均达到7度，这是日本气象厅烈度等级规定的最高等级标准，也是日本九州地方自1923年有观测记录以来首次出现的烈度为7度的地震.根据JMA和日本防灾科学技术研究所 (National Research Institute for Earth Science and Disaster Prevention，简称NIED) 的宽频带地震观测网F-net (Full Range Seismograph Network of Japan；Fukuyama et al., 1996, 1998) 的测定结果 (www.fnet.bosai.go.jp[2016-04-18])，主震震源位于130.7630°E, 32.7545°N，深度12.45 km，发震断层的走向、倾角和滑动角分别是N226°E、84°和-142°，这表明这是一次由走滑断层引起的浅层地壳地震.而日本产业技术综合研究所 (National Institute of Advanced Industrial Science and Technology，简称AIST) 的日本活断层资料库 (www.gbank.gsj.jp/activefault) 和日本地震调查研究推进本部的地震评价报告 (www.jishin.go.jp[2016-05-13]) 指出的发震构造较为复杂，断层应分为南北两段，即震源机制在破裂过程中可能发生了变化.参考震后24 h的余震，其分布特点和震源机制也与此印证.

本文利用NIED强震台网K-NET和Kik-net (Kinoshita，1998；Aoi et al., 2004) 的近场观测记录针对熊本主震进行了震源破裂过程的反演.首先根据F-net的矩张量解，构建了单一震源机制的有限断层模型 (单段式断层模型)，并进行了初步反演；而后参考日本活断层资料库和日本国土地理院 (Geospatial Information Authority of Japan，简称GSI) 提供的合成孔径雷达干涉 (Synthetic Aperture Radar Interferometry，简称InSAR)、全球导航卫星系统 (Global Navigation Satellite System，简称GNSS) 所测量地表形变的反演结果 (www.gsi.go.jp[2016-05-13])，构建了震源机制不同的两个平面断层组成的分段式断层模型 (二段式断层模型)，并再次进行了震源破裂过程的反演.

本文使用了多时间窗函数线性波形反演方法 (Sekiguchi et al., 2000)，这一方法基于Hartzell和Heaton (1983)与Olson和Apsel (1982)的反演方法发展而来：依据震源表示定理 (Aki and Richards, 2002)，空间中某点的运动时程 (位移、速度或加速度) 可表示为，断层面上一系列离散子断层在多个时间窗内破裂效应的加权叠加.这样，反演所需的线性方程可表示为

(1)

其中

式中m(if, is, itm) 表示第if个子断层的第is个滑动方向上第itm个时间窗的系数；为位移的时间导数，即地面运动的速度；Ġ_{kn, l}(x, t-τ; ξ_if) 为格林函数的时间导数.r(=1, …, nt)、n(=1, …, ncmp)、if(=1, …, nf)、is(=1, …, ns) 和itm(=1, …, ntm) 分别代表数据点、数据分量、子断层、滑动方向和时间窗函数的序数.Δt_trig=R/V_r+Δtw(itm-1)，V_r为破裂传播速度 (即首个时间窗函数的激发速度)，R为子断层与破裂起始点的距离，Δtw为时间窗的间隔.每个台站对应的数据和格林函数都以记录三分量中的最大振幅进行了归一化处理.同时，为使反演问题稳定，在反演方程中添加了光滑约束条件，得到的反演方程如下：

(2)

其中，S是光滑约束矩阵，λ为光滑约束强度.λ通过Akaike (1980)提出的赤池-贝叶斯信息准则 (Akaike′s Bayesian Information Criterion，简称ABIC) 选择合适的取值 (Yabuki and Mastsu′ura，1992；Sekiguchi et al., 2000, 2002)，即取ABIC值最小时对应的参数，认为此时在最大程度上从数据和方程中提取了信息.除此之外，限定每个子断层的滑动方向在初始滑动角±45°的范围内，并使用Lawson和Hanson (1974)的非负最小二乘法 (Nonnegative Least Squares，简称NNLS) 求解反演方程.

本文使用了18个近场台站的三分量加速度波形记录，自S波到达前1 s开始截取，时长为25 s.图 1和表 1分别为所采用的18个强震台的分布图和台站信息，主震震中距约为4~64 km，具有较好的方位角覆盖范围.这些台站中有13个K-NET台站和5个KiK-net台站 (www.kyoshin.bosai.go.jp[2016-04-18])，其中KiK-net台站均使用地下钻孔中地震仪的数据记录.台站原始数据需经过预处理，包括滤波、数值积分和重采样，滤波使用Chebyshev一型带通滤波器，滤波范围为0.05~1.0 Hz，重采样的采样率为5 Hz.

图 1

Fig. 1

图 1 主震震中 (五角星) 位置、强震台站及余震分布 图中三角形标示了台站，两个虚线矩形标示了反演使用的二段式断层模型在地面的投影.五角星为主震震中，圆形为主震后24 h (4/15 16:25—4/16 16:25) 内发生的余震震中. Fig. 1 Location of the strong-motion stations (triangle) used in this study and the mainshock (star) with its aftershocks (circle) happened in 24-hour

表 1(Table 1)

表 1 本文中所使用的强震台站 Table 1 Strong-motion stations used in this study 台站 纬度 (°N) 经度 (°E) 深度 (m) 台网 地点 FKO013 33.2252 130.5607 0 K-NET Yame KMM002 33.0185 130.6846 0 K-NET Yamaga KMM003 32.9336 130.5477 0 K-NET Tamana KMM005 32.8761 130.8774 0 K-NET Ohdu KMM006 32.7934 130.7772 0 K-NET Kumamoto KMM009 32.6858 130.9856 0 K-NET Yabe KMM010 32.6136 130.4874 0 K-NET Misumi KMM011 32.6167 130.8652 0 K-NET Tomochi KMM012 32.5078 130.6024 0 K-NET Yatsushiro MYZ001 32.7050 131.3089 0 K-NET Takachiho MYZ020 32.4550 131.1470 0 K-NET Shiiba NGS008 32.9231 130.1843 0 K-NET Kanagai OIT015 32.9722 131.3970 0 K-NET Takeda KMMH02 33.1220 131.0629 123.55 KiK-net Oguni KMMH03 32.9984 130.8301 200 KiK-net Kikuchi KMMH06 32.8114 131.1010 111 KiK-net Hakusui KMMH14 32.6345 130.7521 110 KiK-net Toyono KMMH16 32.7967 130.8199 252 KiK-net Mashiki

表 1 本文中所使用的强震台站 Table 1 Strong-motion stations used in this study

为计算理论格林函数，本文使用Nakamura等 (2008)提出的日本地下P波和S波三维速度结构.基于每个台站不同的地理位置和地质条件，从地下三维速度结构提取出一维水平层状速度结构用于计算，使用的方法是离散波数法 (Bouchon，1981) 和反射/透射系数矩阵法 (Kennett and Kerry, 1979).

为了考查发震断层破裂过程的总体特征，并验证断层模型的几何参数，根据F-net给出的矩张量解，并参考震后24 h的余震分布范围，建立了单一震源机制的有限断层模型 (单段式断层模型)，以此进行了初步反演.断层模型的长度和宽度分别为40 km和18 km，走向N226°E，倾角84°.断层面上破裂起始点位于130.7630°E, 32.7545°N，深度12.45 km处.整个断层面划分成180个大小为2 km×2 km的正方形子断层，每个子断层滑动方向的范围为-142°±45°，滑动时间函数为10个时间窗函数的叠加，时间窗为光滑斜坡函数，上升时间为0.8 s，间隔为0.4 s，这样每个子断层的滑动持续时间为4.4 s.

依据AIST的日本活断层资料库和日本地震调查研究推进本部的地震评价报告，引发熊本地震的是九州地区的Futagawa-Hinagu断层带中Futagawa区段和Mifune区段，两段的走向相差约30°，倾角相差约20°.矩张量解代表了震源机制的总体特征，基于先前使用的单段式断层模型，我们参考GSI的InSAR和GNSS地表形变反演结果改进了断层模型.由于引发地震的断层为两个区段，震源机制在破裂过程中发生变化，因此断层模型采用两个平面断层组合的形式 (二段式断层模型)，总长度依旧为40 km，其中北段长度26 km，南段长度14 km，两段宽度均为18 km，划分为2 km×2 km的子断层，具体参数如表 2所示，其余参数与单段式断层模型一致.

表 2(Table 2)

表 2 二段式断层模型参数 Table 2 Parameters of the two-segment fault model 长度 (km) 宽度 (km) 走向 倾角 滑动角 断层北段 26 18 N235°E 60° -151° 断层南段 10 18 N205°E 72° 176°

表 2 二段式断层模型参数 Table 2 Parameters of the two-segment fault model

依据ABIC准则，选取的破裂传播速度 (即首个时间窗的激发速度) 为2.30 km·s^-1(约为震源深度处S波速度的65%)，光滑约束强度λ为0.004，此时的ABIC值为最小.

基于单一震源机制的有限断层模型，我们进行了初步反演，图 2表示波形反演得到的滑动分布.断层面上出现两个凹凸体，其中一个面积较大，位于破裂起始点东北方向约14 km处，最大滑动量接近10 m；另一凹凸体位于破裂起始点之上接近断层顶端的位置.整个断层在破裂过程中释放地震矩为5.63×10¹⁹ N·m (M_W7.1).图 3对比了实际观测波形和反演合成波形，虽然其中大多数台站的观测记录得到了较好的解释，同时余震也基本分布在断层模型的深度范围内，但余震震源机制的变化可能暗示了主震震源机制在断层破裂过程中发生了变化，那么单一机制的震源模型将不足以充分解释观测记录，个别台站波形拟合不佳的情况也可能源于此.

图 2

Fig. 2

图 2 初步反演使用单段式断层模型得到的滑动分布 图中五角星标示了破裂起始点，箭头代表了滑动矢量，等值线的间隔为1.5 m. Fig. 2 Slip distribution derived from the preliminary inversion using the single-segment fault model The star indicates the rupture staring point, and the arrows indicate the slip vector. The contour interval is 1.5 m.

图 3

Fig. 3

图 3 初步反演使用单段式断层模型得到的合成波形与实际观测波形的比较 (速度时程) 在每个波形的右上方依次标出了实际观测波形和反演合成波形的最大值 (cm·s-1). Fig. 3 Comparison between the synthetic waveforms obtained from the preliminary inversion with the single-segment fault model and the observed ones (velocity) The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform respectively.

前面基于单一震源机制断层模型的初步反演验证了模型几何参数的合理性，下面使用震源机制变化的二段式断层模型再次进行了反演.图 4表示反演得到的滑动分布.断层面上出现了一大一小两个凹凸体.其中大的凹凸体位于北段断层上，在水平方向上从破裂起始点东北方向约8 km延伸到约22 km，在沿断层倾斜方向上从破裂起始点之下约4 km一直延伸到断层面最上端，这一凹凸体中滑动量最大的部分也接近于断层面最上端.较小的凹凸体出现在南段断层上，位于破裂起始点之上约8 km处，面积约4个子断层大小 (16 km²).整个断层面上最大滑动量约7.9 m，释放地震矩5.47×10¹⁹ N·m (M_w7.1).

图 4

Fig. 4

图 4 使用二段式断层模型反演得到的滑动分布 图中五角星标示了破裂起始点，箭头代表了滑动矢量，等值线的间隔为1.5 m. Fig. 4 Slip distribution derived from the inversion using the two-segment fault model The star indicates the rupture staring point, and the arrows indicate the slip vector. The contour interval is 1.5 m.

图 5表示整个断层面的破裂过程.破裂起始于南段断层，两个凹凸体均在破裂开始后约4 s时出现，其中南段断层上小凹凸体的破裂持续约4 s，北段断层上较大凹凸体的破裂持续约8 s，破裂起始后沿断层向上传播约8 s，整个断层的破裂在约18 s内结束.图 6表示断层面上各子断层的滑动速度时间函数.断层北段凹凸体区域出现了尖峰，说明在很短时间内就达到了滑动速度峰值.

图 5

Fig. 5

图 5 使用二段式断层模型反演得到的破裂过程 Fig. 5 The rupture process derived from the inversion using the two-segment fault model

图 6

Fig. 6

图 6 使用二段式断层模型反演得到的每个子断层的滑动速度时间函数，背景为断层面上的滑动量分布 Fig. 6 Slip-rate of each subfault derived from the inversion using the two-segment fault model. The background is the final slip distribution on the fault plane

图 7对比了实际观测波形和反演合成波形，相较于图 3，不仅大多数台站的波形拟合效果有所提高，而且之前个别拟合不佳的台站也得到改善，几乎所有台站都可以很好地拟合实际观测记录.在此引入波形拟合参数W_xc，作为评价波形拟合效果的标准，

(3)

图 7

Fig. 7

图 7 使用二段式断层模型反演得到的合成波形与实际观测波形的比较 (速度时程)，在每个波形的右上方依次标出了实际观测波形和反演合成波形的最大值 (cm·s-1) Fig. 7 Comparison between the synthetic waveforms obtained from the inversion with the two-segment fault model and the observed ones (velocity). The maximum values of the observed and synthetic waveforms are shown at the upper right of each waveform respectively

其中，stn为台站数目，nt为使用数据记录的长度，ncmp为记录的分量数目，D_{ist, it, icmp}为观测记录，S_{ist, it, icmp}为反演合成记录，0≤W_xc≤1，且W_xc越大说明拟合程度越好.使用单段式断层模型反演时，W_xc为0.7655，而使用二段式断层模型时，W_xc为0.8040，波形拟合效果提高了9.73%.为了解两段断层分别对各个台站的影响，我们分别合成了两段断层面反演合成波形的水平分量，如图 8所示.对于台站FKO013、KMM002、KMM005、KMMH02、KMMH03、KMMH06、MYZ001和OIT015，波形中的脉冲主要由北段断层产生.对于台站KMM012、KMMH14和NGS008，波形中的脉冲主要来自南段断层的滑动.而其余台站基本上受两段断层的同时影响，这些台站中的KMM010和KMM011，分布在断层模型南段的西南方向，此时南段断层产生的脉冲先于北段断层到达，这与断层和台站的空间位置关系相吻合；对于台站KMM003、KMM006、KMM009和KMMH16，两段断层产生的脉冲同时出现，且影响程度大约各占一半，其中KMM003和KMM009与两段断层距离几乎相同，说明这种情况也应与断层—台站的空间位置有关；而对于台站MYZ020，其处在断层北段走向的法线方向的南侧，虽然在波形上两段断层的影响同时出现，但产生的脉冲却极性相反.

图 8

Fig. 8

图 8 使用二段式断层模型反演分别得到的北段断层 (中) 和南段断层 (下) 上滑动合成波形水平分量的比较.实际观测波形也被列出 (上)，同时在每个波形的右上方依次标出了实际观测波形和反演合成波形的最大值 (cm·s-1) Fig. 8 Comparison between the horizontal component of synthetic waveforms obtained from the slips on the northern (middle) and southern (bottom) segment for the inversion using the two-segment fault model. The observed waveform is also shown (top), and the maximum values of waveforms are shown at the upper right of each waveform respectively

为了测试所使用反演方法的分辨率，我们进行了一次理论验证.首先，构建了一个与反演得到滑动分布结果相似的断层滑动模型作为输入，仍设置有一大一小两个凹凸体.在东北方向较大的凹凸体区域中，颜色由浅到深，滑动量依次为2 m、3 m、5 m和7 m；面积较小的凹凸体的滑动量为3 m.通过正演得到合成波形记录，再添加标准差为合成波形最大振幅10%的高斯噪声，得到含有噪声的合成波形记录.分别使用无噪声和有噪声的合成记录进行波形反演，反演参数与之前相同，得到的断层面滑动分布如图 9所示.结果显示，无噪声情况下几乎能完全恢复输入模型，但由于添加了光滑约束条件，反演结果相较输入模型被模糊了.在添加高斯噪声的情况下，为压制反演解的不稳定性，光滑约束条件必不可少，但反演结果依旧能在较大程度上恢复输入模型，可见，我们在本文中使用的反演方法具有较好的分辨率和稳定性.

图 9

Fig. 9

图 9 (a) 输入模型；(b) 使用无噪声合成数据反演得到的滑动分布；(c) 使用添加标准差为最大振幅10%高斯噪声的合成数据反演得到的滑动分布 Fig. 9 (a) Input model; (b) The slip distribution derived from the synthetic data without noise; (c) The slip distribution derived from the synthetic data with a Gaussian noise having a standard deviation of 10% of the maximum amplitude

在使用二段式断层模型得到的反演结果基础上，我们合成了此次地震在九州地区的峰值地面速度 (Peak Ground Velocity，简称PGV) 和峰值地面加速度 (Peak Ground Acceleration，简称PGA) 分布，并与NIED的实测烈度 (日本称震度) 和PGA分布进行了比较，如图 10所示.由于反演使用的数据频段为0.05~1.0 Hz，所以反演得到的震源破裂过程也对应这一频段，因此合成的地震动时程并不包含高频部分.而高频的地震动往往对应着更大的地面运动速度和加速度，因此造成了由反演结果合成的PGA比实测PGA要小约一个量级.但通过比较合成PGV、PGA和实测PGA，发现它们的分布特征大致相似，地震动剧烈的区域均沿着断层走向伸展，这也对应着实测地震烈度较高的区域.地震烈度描述了地震的震害影响程度，基于合成PGV、PGA分布和实测烈度、PGA较好的相似关系，可以认为使用低频数据反演得到的震源破裂过程也能够反映出整个震区的受震害程度，说明震源运动学反演得到的结果具有实际意义.

图 10

Fig. 10

图 10 (a) 使用反演结果合成的PGV分布 (0.05~1.0 Hz)；(b) 使用反演结果合成的PGA分布 (0.05~1.0 Hz)；(c) NIED的实测烈度 (日本称震度) 分布；(d) NIED的实测PGA分布 Fig. 10 (a) The distribution of PGV synthesized by the inversion result (0.05~1.0 Hz); (b) The distribution of PGA synthesized by the inversion result (0.05~1.0 Hz); (c) The distribution of Japanese seismic intensity observed by NIED; (d) The distribution of PGA observed by NIED

通过使用震源机制变化即具有两个平面断层的有限断层模型，得到了更为理想的反演结果，Kubo等 (2016)的反演过程也印证了这一现象.从另一角度看，这也说明了波形反演可以复原得到实际断层的空间位置属性 (Wu and Takeo, 2004).熊本地震序列前震和余震的震源机制都与主震有着些许差别，加之主震的规模也远大于前震和余震，都暗示着主震的发震断层更为复杂.仅使用前震、主震或余震的单一震源机制解进行反演，无法得到较为理想的结果.因此，设置断层模型时结合必要的地质调查、地球物理勘探或地表破裂的相关资料，对于反演震源机制复杂的破裂过程会有很大帮助.

本文基于余震分布、GNSS和InSAR的地表形变资料，结合地质及地球物理信息构建了熊本地震的发震构造模型，使用13个K-NET和5个KiK-net近场强震台站的波形记录，通过多时间窗函数线性波形反演法反演了2016年日本熊本县M_W7.1地震的震源破裂过程.结果显示，这是一次沿Futagawa-Hinagu断层带发生的右旋走滑破裂事件，震源机制在破裂过程中发生了变化，发震断层分为南北两段，断层面上存在一大一小两个凹凸体，均在破裂起始后约4 s出现.其中较大的凹凸体位于断层北段、破裂起始点东北方向上，破裂持续约8 s；较小的凹凸体位于断层南段，在破裂起始点上方接近断层面顶端的位置，破裂持续约4 s.断层面上滑动主要集中于断层北段，最大滑动量约7.9 m.整个断层面的破裂过程持续约18 s，释放地震矩5.47×10¹⁹N·m (M_W7.1).

衷心感谢评审专家提出的宝贵修改建议；文中所使用的强震数据源自日本强震观测台网K-NET和KiK-net (www.kyoshin.bosai.go.jp[2016-04-18])；文中所有图件均使用通用制图工具GMT (the Generic Mapping tools，Wessel and Smith, 1995) 绘制，在此一并感谢.