地球物理学报  2021, Vol. 64 Issue (12): 4378-4393   PDF    
基于区域宽频带和近震波形数据的多点源反演——以2016 MW7.1日本熊本地震为例
施其斌1,2, 韦生吉1,2     
1. 南洋理工大学亚洲环境学院, 新加坡 639798;
2. 南洋理工大学新加坡地球观测中心, 新加坡 639798
摘要:中强度(M6-7)地震可以对人类社会基础设施和生命财产造成威胁,并且其发震频率要远高于M8+地震.从震源物理的角度来说,此类地震的发震断层的几何形态和破裂过程往往较为复杂.因此,准确地确定此类地震的震源参数对地震防治和了解震源物理都具有重大意义.作为承接地震单点源解和有限破裂模型的桥梁,地震多点源解可以为有限破裂模型反演提供可靠的断层几何,并有效地帮助约束破裂的时空演化过程.本文基于马尔可夫链蒙特卡洛(MCMC)非线性反演算法,将一种多点源波形反演方法的适用范围从近震距离推广到了区域距离,并验证了利用余/前震(或历史地震)对地震波传播路径进行校正的有效性.我们将改进后的方法应用于具有代表性的2016 MW7.1熊本地震的研究中.结果发现,利用一个MW6.0和一个MW5.4的前震,可以对半径100 km范围内强震记录和1000 km内的宽频带记录进行有效的路径校正,筛选出现有速度模型所适用的地震台和震相,并且可以确定波形拟合中高质量Pnl波和面波所能利用的频率上限,即高达0.3 Hz(Pnl)和0.2 Hz(面波).研究结果还显示区域宽频带数据反演的结果和强震+高频GPS数据结果有高度的一致性,验证了将该方法推广到区域数据反演的可行性.对熊本MW7.1主震的反演结果显示,在反演所用的频段,我们需要至少四个子事件来表示该地震的运动学过程.第一个MW6.7子事件发生在地下12 km深的Hinagu断层上,并具有高倾角的纯走滑位错解;第二(MW6.7)和第三个(MW6.7)子事件的深度均为7 km,其走向和Futagawa断层一致并有一定的正断滑移成分;最后一个子事件(MW6.6)最浅(1~2 km),也具有最强(滑动角=-137°)的正断滑移成分.这些子事件揭示了发震层的厚度随着主震向东北方向的破裂传播距离增加而变薄,并且其深度随着破裂位置靠近Aso火山而变得越浅.当使用离断层最近的0701高频GPS台站时,五点源反演可以分辨出在破裂最后阶段沿一个次级断层的传播,我们认为该断层分叉现象有可能对地震的停止起到了关键性的作用.
关键词: 地震破裂过程      非线性反演      多点源      非对称震源时间函数      熊本地震     
Multiple point source inversion with regional broadband and local waveform data -a case study of the 2016 MW7.1 Kumamoto, Japan earthquake
SHI QiBin1,2, WEI ShengJi1,2     
1. Asian School of the Environment, Nanyang Technological University, Singapore 639798;
2. Earth Observatory of Singapore, Nanyang Technological University, Singapore 639798
Abstract: Moderate-size (M6-7) earthquakes possess great threats to human society,as they occur much more frequent than M8+ earthquakes and many take place inland. The rupture processes and fault geometry of this kind of earthquake are usually very complex. Precise determination of their kinematic rupture parameters is a fundamental task in seismology,which is critical to understanding the earthquake physics as well as to better prepare for seismic hazard. As an earthquake source model with a precision between point source and finite fault model,Multiple Point Source (MPS) solution of large earthquakes provide reliable fault geometry and first order rupture evolution of earthquakes,which is extremely helpful for more detailed finite fault inversion and modeling. Here we generalize a Markov-Chain-Monte-Carlo (MCMC) based MPS inversion algorithm to the regional broadband waveform inversion and apply to the 2016 MW7.1 Kumamoto mainshock. We first calibrate the path and select the frequency ranges for the strong motion (< 100 km) and regional (< 1000 km) broadband stations through the point source modeling of the MW5.4 and MW6.0 foreshocks,which results in up to 0.3 Hz for Pnl waves and 0.2 Hz for surface waves. The inversion results from regional,strong motion and high-rate GPS datasets show highly consistent solutions,validating the generalization of our MPS inversion method. Our preferred model show that the earthquake is composed of at least four sub-events. The first (MW6.7) sub-event is a high-dip-angle purely strike-slip event,taking place near the hypocenter at the depth of 12 km. The strike of the first sub-event is well consistent with the strike of the Hinagu fault. The second (MW6.7) and the third (MW6.7) sub-events all took place at the depth of 7 km,with strike more consistent with the Futagawa fault and with slight normal slip component. The last sub-event (MW6.6) is the shallowest (1~2 km) and its focal mechanism is composed of almost equal percentage of strike-slip and normal-slip components. Our five point source inversion with the closest high-rate GPS station (0701) shows that the last stage of the rupture bifurcated onto a secondly fault,probably playing an important role in stopping the rupture.
Keywords: Earthquake rupture process    Non-linear inversion    Multiple point source    Asymmetric source time function    Kumamoto earthquake    
0 引言

地震震源的运动学过程是地震学研究的一个基本内容.随着地震波观测数据的质量的提高——包括台站密度的增加和观测频段范围的拓展等——震源运动学的发展经历了从简单点源(例如地震的位置、大小、震级和矩张量解(Aki and Richards, 2002))到更接近真实地震的复杂的有限破裂模型(姚振兴和纪晨, 1997; Ji et al., 2002; Xu and Chen, 2002; 许力生和陈运泰, 2004; 张勇等, 2008; Zhang et al., 2009; Yue et al., 2012; 岳汉等, 2019)的进步.点源对于中强度(MW6-7)地震的描述往往过于简单,然而有限破裂模型又常常用过多的参数来非唯一地表示地震的破裂过程.对后者而言,研究者一般很难估计断层几何等参数的误差.地震的多点源解刚好填补了点源和有限破裂模型之间的空缺,弥补了上述两者的不足之处.通常情况下,多点源解所适用的反演频段也介于单点源和有限破裂模型的频段之间.在此频段,多点源解通常只需要用几个(例如3~8个)双力偶位错点源来描述地震的破裂过程,要远比单矩张量解更能代表地震的复杂性.同时,多点源解的参数空间远小于有限破裂模型,因而不会过度拟合观测数据,更易于对震源参数进行误差估计.近年来,多点源反演已经被多次应用于大地震研究中(Chen et al., 2014; Duputel and Rivera, 2017; Yue and Lay, 2020; Jia et al., 2020),此类研究通常使用远震体波和长周期的面波.对于那些密集强震记录到的中强度地震,多点源反演也获得了很好的效果(Shi et al., 2018; Shi and Wei, 2020).然而,对中强度地震的多点源研究还没有充分利用区域距离(< 1000 km)的台站记录.实际上,除了少数地震频发的发达地区(例如日本和美国西海岸),绝大部分中强度地震频发的地区并没有密集的强震台网.相比之下,区域距离范围内台站数量对全球绝大部分地区都更加可观.因此,是否能利用此类数据对中强度地震的破裂过程进行多点源反演非常值得探究.为此,本文通过对2016 MW7.1熊本地震的研究来测试利用区域宽频带地震数据进行多点源地震反演.我们将利用前/余震来对区域台站的不同分量进行路径校正,以得到适用于主震反演的频段.然后对比加入了更为密集的强震和高频GPS数据得到的多点源解和现有的有限破裂模型,来进一步验证多点源解.

我们改进了Shi等(2018)所提出的基于马尔可夫链蒙特卡洛(Markov Chain Monte Carlo,MCMC)的反演方法.该方法借鉴了前人所提出的“剪切黏贴”(Cut-And-Paste,CAP)理念(Zhao and Helmberger, 1994; Zhu and Helmberger, 1996),将每个台站的三分量地震波数据截成Pnl波和面波部分,在使用理论地震波拟合数据时允许对不同分量使用不同的时移进行拟合,然后计算时移后的波形拟合误差.Shi等(2018)采用该方法,利用强震数据研究了2016熊本地震序列的MW6.2前震,结果揭示了发震区断层几何形态的复杂性.在本文中,我们将该多点源反演方法应用到包含高频GPS数据的密集近震观测,并拓展到对区域地震数据的分析.在此基础上,提出了使用正则化的Yoffe震源时间函数来代替之前的三角形函数,并且还提出了基于台站和波形分量来选择滤波参数,从而最大化地利用地震波形所包含的信息.为了对每个分量选择合适的滤波频段,我们利用比研究地震震级更小的地震来进行路径校正.这些地震的震级应当小到在主震反演的滤波频段可以被合理地简化成一个点源,但同时也要大到可以被足够多的台站记录到,不管是强震仪还是宽频带地震仪.在后文将给出更多相关的技术细节.

1 2016年熊本地震序列

2016年的MW7.1熊本地震序列发生在日本九州岛的Beppu-Shimabara地堑地区,起始于Futagawa断层和Hinagu断层相交处.该断层系统的形成主要是由菲律宾板块向西北方向略微斜向俯冲所导致,板块俯冲的速度5~6 cm·a-1,斜冲投影在Futagawa断层和Hinagu断层上的滑动速率为7~8 mm·a-1 (Nishimura and Hashimoto, 2006)(图 1a).

图 1 板块俯冲区构造和地震观测网 (a) 菲律宾板块斜插俯冲在欧亚板块之下.日本的F-net宽频地震观测台网被标记为深绿色三角形. (b) 九州地区的活动断层(黑色线段)和密集分布的K-net以及Kik-net强震地震仪(黑色三角形).该序列的余震由圆点表示,深度用不同颜色显示. Fig. 1 Tectonic setting of the plate subduction zone and seismic observation networks (a) Oblique subduction of the Philippine Sea Plate beneath the Eurasian plate. The dark green triangles denote the location of the F-net broadband stations. (b) Active faults (black lines) and the K-net and Kik-net strong-motion stations (black triangles) at Kyushu. Aftershocks are denoted by dots color-coded by depth.

序列始于4月14日发生的MW6.2前震,约2.5 h后,发生了一个MW6.0的余震,MW7.1主震发生在当地时间4月16日的凌晨,震中位于前震的东北方向(Kato et al., 2016).该地震序列的余震极其活跃,产生了至少30000个Mj (JMA(Japan Meteorological Agency)震级,即日本气象局测定的地震震级)大于2.5的余震(Kato et al., 2016).地震序列还造成了复杂的地表破裂,以右旋走滑为主但有明显的正断层滑动(Shirahama et al., 2016),尤其在两个断层相交处和破裂的东北段(图 1b).

该地震序列有着高质量的密集强震观测、卫星影像观测以及实地考察记录,为震源运动学研究提供了宝贵的数据.前人利用地震学数据或者联合测地学数据对MW6.2前震(Asano and Iwata, 2016)和MW7.1主震(Asano and Iwata, 2016; Yagi et al., 2016; Hao et al., 2017; Kobayashi, 2017; Yue et al., 2017)进行了有限破裂模型的反演.但是不同研究者所得到的主震破裂模型之间仍存在不小的差异.例如,Asano和Iwata(2016)的结果显示主震的破裂主要集中在8~10 km深度,位于震中东北约20 km处;而Yue等(2017)Kobayashi(2017)的结果则显示主震主要的位错比较浅,已经破裂到了地表.这些有限破裂模型使用2~4个断层面来近似主震较为复杂的断层几何.相对的,前人对MW6.2前震进行的有限破裂模型反演采用了相对简单的断层几何形态,仅为一个矩形的断层面(Asano and Iwata, 2016).不同于有限破裂模型,Shi等(2018)提出了一种基于MCMC的地震多点源反演方法,用于在相对高频反演MW6.2前震所产生强震波形.他们使用了一个MW5.4地震的波形记录来对所使用的波形分量进行筛选并确定最佳滤波频段,在此基础上的多点源反演结果揭示了MW6.2前震较为复杂的断层几何,即发现该地震发生在沿相反方向倾斜的多个断层上.

2 方法和数据

如前所述,Shi等(2018)所提出的多点源反演方法已经在MW6.2前震中得到了成功的应用.该方法尽管对强震数据反演非常有效,但还没有在高频GPS数据和区域宽频带地震波的反演中尝试过.为此,本文将该方法拓展到高频GPS和区域宽频带地震记录来对MW7.1熊本主震的破裂过程进行研究,并将结果和强震+高频GPS记录反演的多点源解以及已经发表的有限破裂模型进行对比,旨在了解该方法在中强度地震研究中的优势和局限性.我们还对该方法做了进一步的改进.Shi等(2018)的方法中对所有观测台站相同的波形分量采用相同的滤波频段(例如用0.02~0.2 Hz对Pnl统一进行滤波),但是由于每个台站所对应的速度结构不相同,其波形分量对应的最佳拟合频段也有所不同.这里我们采用了台站相关的滤波参数,以最大限度地从波形数据中提取地震破裂过程的信息.另外,还改进了反演所用的震源时间函数,不同于之前的等腰三角形函数,这里使用正则化的Yoffe函数(Tinti et al., 2005)来作为震源时间函数(图 2),以测试所用的数据是否能进一步分辨地震位错发生的时间过程.实际上Yoffe函数更符合滑动弱化震源动力学破裂模型中断层上应力的演化过程(Yoffe, 1951; Tinti et al., 2005),这对于理解震源破裂物理过程具有重要的意义.新的震源时间函数用αβ两个参数来定义函数的形状(图 2),数学表达如下:

(1)

图 2 正则化Yoffe震源时间函数的数学定义 该震源时间函数表示为黑色,由三个对称轴相同但周期不同的cosine/sine函数组成,分别显示为红色、蓝色和绿色. Fig. 2 The mathematical composition of the reguralized Yoffe-type source time function The Yoffe function is shown in black, composed of the three cosine/sine functions with the same symmetric axis but different periods and shown in red, blue and green, respectively.

其中f1(t), f2(t)和f3(t)分别用αβ来表示:

(2)

(3)

(4)

表达式(1)具有更多的灵活性,既可以表示一个对称的函数(α=0,β=1),也可以表示非对称的“长尾”的函数(例如α=0.5,β=0.2).

所使用的波形数据包括强震、高频GPS和宽频带地震数据(图 3).三分量强震数据是从KiK-Net和K-net网站上(https://www.kyoshin.bosai.go.jp)直接下载,其原始波形为加速度记录,将其积分至速度和位移波形用于反演.宽频带波形数据是从NIED(https://www.fnet.bosai.go.jp/top.php)上下载,去仪器响应后转化成位移和速度波形两种记录.我们还从日本的Nippon GPS Data Service Corporation购买了1 Hz采样的主震高频GPS数据.除了主震(MW7.1),还从这些地震台网上下载了两个校正地震事件的波形数据.这两个校正事件分别发生在Futagawa断层(MW5.4,32.776°N,130.815°E,2016-04-14-13∶07∶35.290)和Hinagu断层上(MW6.0,32.7007°N,130.778°E,2016-04-14-15∶03∶46.450).这两个事件均为MW6.2前震的余震,其破裂过程相对简单,在相对较高(例如0.3 Hz)的频段可以看作点源.利用这两个校正事件,我们确定了所使用的台站和波形分量,以及台站相关的滤波参数(例如以80%互相关系数来选择滤波频段和分量).最终选取了110个强震台,包含276个Pnl波形分量和251个面波波形分量,36个宽频带台站,包含90个Pnl波形分量和83个面波波形分量,以及26个高频GPS台站.图 3给出了这三个地震以及MW6.2前震在几个代表性台站上的0.2 Hz低通滤波滤波波形.可以看到两个校正事件(MW5.4,MW6.0)的波形相对简单,而MW6.2前震则要复杂得多,说明其破裂过程更加复杂.值得指出的是,在KMM005台站上,MW6.2前震的波形甚至比MW7.1主震还要复杂.这一方面是因为主震是向东北方向的单侧破裂,而KMM005台站(方位角38°)正好位于破裂的传播方向上,因此多普勒效应将波形压缩显得更为简单;另一方面是因为主震的震源时间函数更长、更平滑,导致了主震的波形的频率成分更低一些,看起来没那么复杂.而背离主震破裂方向的KMM008台站其波形振幅要低很多,但是持续时间更长,波形也更为复杂.这些都符合主震单侧破裂的特征.同时宽频带台站的波形数据则显示MW7.1主震的波形要比其他地震复杂得多(图 3cd),因此将这些数据用于多点源反演可以对其破裂过程进行很好的约束.相对于强震台站,我们需要选择更低频段的滤波才能对区域台站的校正事件波形实现较高互相关系数的拟合.区域台站的路径校正主要是通过MW6.0前震实现的.

图 3 近Futagawa断层的熊本主震及前震的强震和宽频数据对比 (a)和(b)为强震仪KMM006和KMM005的观测数据低通滤波后的波形;(c)和(d)为宽频地震仪TKO和SIB记录低通滤波后的波形.7.1级主震的波形由蓝色特别显示.由上至下分别为竖直、南北和东西分量的原始地震波形.所有波形均通过了角频率为0.2 Hz的低通滤波. Fig. 3 Strong-motion and broadband observations of the Kumamoto mainshock and foreshocks near the Futagawa fault (a) and (b) show the low-pass filtered waveforms of the strong-motion stations KMM006 and KMM005, respectively. (c) and (d) show the low-pass filtered waveforms of the broadband stations TKO and SIB, respectively. The mainshock waveforms are highlighted in blue. The vertical, north-south and east-west waveforms are arranged from the top to the bottom. All waveforms are low-pass filtered with the corner frequency of 0.2 Hz.
3 反演结果

我们使用经过路径校正得到的波形分量和滤波频段,来对主震的多点源解进行反演.从单个点源反演开始逐渐增加点源的数量,直到反演得到的波形互相关系数的统计分布非常接近校正事件的波形拟合效果,这样可以避免对数据的过度拟合.图 4给出了联合使用强震数据和区域宽频带数据,和只使用区域宽频带数据得到的多点源反演结果,包括用速度波形和位移波形得到的结果.包含强震数据的反演结果显示我们需要至少4个子事件来足够准确的表示主震的运动学破裂过程.我们也将三点源反演的结果和四点源解进行比对(图 5),发现四点源解包含更多细节.如图 6所示,当只用3个子事件时,得到的互相关系数的统计直方图和用校正事件得到的直方图还有明显的差距,而当用到4个子事件时统计的结果和校正事件更为接近.我们注意到3个和4个子事件的波形拟合结果最明显的差别集中在离地震最近的几个强震台站上,把几个代表性的台站在图 7中标注了出来.对这几个台站,3个子事件明显对波形和振幅的拟合都差于4个子事件的解.总体而言,这两种数据组合的反演结果非常稳定和可靠,这点从反演的误差分布可以看出来(图 8).包含强震数据的反演所对应的误差分析和统计显示,各个子事件的走向、倾角和滑动角的误差都在5°以内,震级的误差小于0.1,深度和水平位置的误差则控制在2 km以内(图 8).而区域宽频带数据所对应的误差水平要相对更大一些.同时我们还看到用区域宽频带速度波形得到的反演结果和强震数据得到的结果更一致,也更稳定.这样的差别是合理的,因为首先区域数据所对应的滤波频段要更低,而且台站距离地震更远,对地震的时空分辨率不可避免的要更低;其次在相同的滤波频段速度记录所对应的波形要比位移波形更复杂,因此对各子事件约束更为可靠.总之,区域宽频带数据的结果和强震数据结果非常一致,所对应的误差范围也足以分辨断层几何的变化和破裂过程的时空演化,验证了可以将该方法扩展到区域台站的反演中的有效性.

图 4 矩震级7.1级的主震破裂过程以及早期余震分布 (a)强震和宽频数据联合反演所得的四点源解将主震分为四个由西南向东北传播的子事件,其水平位置用较大圆圈表示,相应的震源机制解绘制成红色沙滩球形状.余震则由较小圆圈表示,并且震源深度用不同颜色显示.黑色和红色线段分别表示地震前所测绘的活动断层和此次地震所产生的地表破裂(Shirahama et al., 2016).近地震台站用黑色三角形表示.(b)该四点源解(同(a))与测地学资料导出的有限断层模型的对比.此次地震序列的主要事件由红色五角星标记.右下角地震矩张量投影依次为四点源解的合成解、Global Centroid-Moment-Tensor (CMT)反演解和美国地质调查局的公布结果.(c)联合高采样GPS数据,强震地震数据和宽频地震数据的反演结果.左上角为该四点源解的子事件震源时间函数(黑色)以及累加震源时间函数(灰色阴影).地震仪和GPS台站分别由黑色三角形和白色方块表示.(d)只利用宽频地震仪反演的四点源解.其中黑色和灰色分别表示位移波形和速度波形的反演结果.左上角为位移波形反演得到的震源时间函数,而右下角为速度波形得到的震源时间函数. Fig. 4 Rupture process of the MW7.1 mainshock and early aftershocks (a) The four-point-source solution of the joint inversion of the strong-motion and broadband data divides the mainshock into four subevents propagating from south-west to north-east. Four large circles denote their locations and four red beachballs indicate their focal mechanisms. Aftershocks are indicated by small circles that are color-coded by event depth. Black and red lines represent the previously mapped active faults and the surface ruptures (Shirahama et al., 2016) of the MW7.1 earthquake, respectively. Near-fault strong-motion stations are denoted by black triangles. (b) The comparison between the four-point-source solution (same as (a)) and the finite-fault model derived from the geodetic data. The red stars show the epicenters of the major earthquakes that occurred in this sequence. The lower-right inset shows the total moment tensor of the four-point-source solution in comparison with the Global Centroid-Moment-Tensor (CMT) and United States Geological Survey (USGS) solutions. (c) The four-point-source solution obtained from the inversion incorporating the high-rate GPS, strong-motion and broadband seismic data. Top-left inset shows the source time functions of subevents (black) and the total sum (gray-shaded area). The seismic and GPS stations are denoted by black triangles and white squares, respectively. (d) Four-point-source solutions obtained from the inversions with only F-net broadband data. Black and gray color indicates inversion results obtained from the displacement and velocity waveforms, respectively. The top-left shows the source time functions corresponding to the solution obtained from the displacement waveform inversion. The lower-right is corresponding to the solution obtained from the velocity waveform inversion.
图 5 强震和宽频数据联合反演所得的四点源解和三点源解的对比 (a)四点源解的震源参数和地图展示.空心圆圈表示子事件的水平位置,与其相对应的震源机制解相连.空心和实心的三角形分别表示朝向和远离破裂方向的地震台站.底部小图显示四个子事件震源时间函数(黑色)以及累加震源时间函数(灰色区域).(b)与(a)相似,但为三点源解的震源参数和地图表示.底部小图中绿色虚线为(a)中的累加震源时间函数. Fig. 5 Comparison between the four-point-source and three-point-source solutions obtained from the joint inversions of the strong-motion and broadband data (a) Source parameters and the map view of the four-point-source solution. The subevent locations are denoted by open circles, connected with the beachballs that show the corresponding focal mechanisms. The open and solid triangles indicate the strong-motion stations that are toward and away from the rupture directivity. The bottom inset shows the source time functions of four subevents (black) and the sum of them (gray-shaded area). (b) Similar to (a), but for the source parameters and the map view of the three-point-source solution. The green dashed line in the bottom inset shows the cummulative source time function in (a).
图 6 路径校正事件和主震三点源以及四点源反演的波形拟合效果比较 此处只对比了利用强震和宽频数据联合反演的情况.四幅直方图分别统计不同波形分量的互相关系数. Fig. 6 Comparison of the quality of the waveform fits between the path calibration event, three-point-source inversion and four-point-source inversion Here the strong-motion and broadband waveforms are jointly used. The four panels of histograms show the statistics of the cross-correlation coefficients of different waveform components.
图 7 部分波形拟合效果 (a)和(b)分别为图 5中四点源解和三点源解的波形拟合.带通滤波后的观测数据和合成数据分别显示为黑色和红色,其拟合所需时移以及表示拟合效果的互相关系数分别标记在波形下方.台站名、台站震中距、台站方位角标记在波形左侧,波形分量名标记在顶部. 滤波频段因台站不同而有所区别,Pnl波的频段总体在0.01~0.30 Hz,面波的频段总体在0.01~0.22 Hz. Fig. 7 Selected waveform fits (a) and (b) are waveform fits of the four-point-source and three-point-source shown in Fig. 5. The bandpass filtered data and synthetic waveforms are shown in black and red, respectively. Numbers below each set of waveforms are the time shift used for alignment and cross-correlation coefficient of the two waveforms. Station names are shown on the left with the distance and azimuth indicated above and below respectively. The names of component are shown on the top. The filtering frequency is station-dependent, within 0.01~0.30 Hz for Pnl wave and 0.01~0.22 Hz for surface waves.
图 8 多组数据反演所得四点源解的不确定性分析 (a)表示强震仪波形和宽频地震仪波形数据的联合反演不确定性.(b)表示GPS,强震仪波形和宽频地震仪波形数据的联合反演不确定性.(c)和(d)各代表只使用宽频地震仪位移波形和只用宽频地震以速度波形的反演的不确定性分析.以上四幅图通过相似方式呈现MCMC采样过程来反映其不确定性.四个子事件按照发震顺序分别由红、蓝、紫、绿四色显示.地图上主要绘制了四个子事件水平位置采样和反演得到的最佳震源机制.左上角绘制四个子事件断层解抽样.左侧中间小图为子事件震源时间函数的抽样,其中粗黑线为反演最终震源时间函数,并且四点源累加震源时间函数由灰色填充区域显示.右下角小图为PP′竖直剖面上投影的子事件的深度抽样,最佳深度显示为黑色圆圈. Fig. 8 Uncertainty analysis of the four-point-source inversions using multiple datasets (a) shows the uncertainty of the joint inversion of the strong-motion and broadband waveforms. (b) shows the uncertainty of the joint inversions incorporating the GPS, strong-motion and broadband data. (c) and (d) are for the inversions using only the broadband displacement and only the velocity waveforms, respectively. The uncertainty is represented by MCMC samples of the four subevents. The four subevents are shown in red, blue, purple and green in the order of occurrence. The sampled horizontal locations and the optimal focal mechanisms is shown in map view. The top-left inset shows the sampled fault-plane solutions of the subevents. The middle-left inset shows the sampled source time functions of the subevents with the final solution denoted by the thick black lines. The gray shaded area is the total source time function of the four-point-source solution. The bottom-right inset shows sampled depths projected to the PP′ cross-section with the optimal depths denoted by black circles.

此外,我们注意到当把高频GPS数据包括在反演中时,得到的结果和强震数据得到的结果几乎完全相同(图 9),这说明数据之间的一致性非常好,也进一步验证多点源反演的稳定性和可靠性.

图 9 GPS、强震仪和宽频地震数据的联合反演所得到的五点源解 (a)五点源解的地图表示.最靠近主断层的GPS台站由白色方块表示.右下角为GPS台0701附近区域的放大视图.左上角为五个子事件对应的震源时间函数以及累积震源时间函数.(b)部分GPS台站位移波形拟合,台站位置如(a)所示.带通滤波后的观测数据和合成数据分别显示为黑色和红色,其拟合所需时移以及表示拟合效果的互相关系数分别标记在波形下方.台站名、台站震中距、台站方位角标记在波形左侧,波形分量名标记在顶部.滤波频段因台站不同而有所区别,总体在0.001~0.1 Hz.(c)0701观测波形(黑色) 与五点源解的合成波形(红色)以及图 4a中四点源解的合成波形(浅蓝)的对比. Fig. 9 Five-point-source solution obtained from the joint inversion of GPS, strong-motion and broadband data (a) The map view of the five-point-source solution. The nearest GPS stations are denoted by the white squares. The lower-right inset is a enlarged view near station 0701. The top-left inset shows the source time functions of the five subevents and the cumulative source time function. (b) Selected waveform fits of GPS stations shown in (a). The bandpass filtered data and synthetic waveforms are shown in black and red, respectively. The numbers below each set of waveforms are the time shift used for alignment and cross-correlation coefficient of the two waveforms. Station names are shown on the left with the distance and azimuth indicated above and below respectively. The names of component are shown on the top. The filtering frequency is station dependent and is within 0.001~0.1 Hz. (c) The comparison between the observation (black) of 0701, five-point-source synthetics (red), and four-point-source synthetics (cyan) that was calculated with the solution in Fig. 4a.

近震(强震+高频GPS)和区域数据联合反演的结果显示需要至少4个子事件来表示主震.第一个子事件(MW6.75)位于12 km深度,非常靠近其震中(大约在Hinagu断层和Futagawa断层的交界处),具有一个高倾角(78°)的纯走滑震源机制解,其走向(208°)和Hinagu断层的走向非常一致.第二个(MW6.67)和第三个(MW6.62)子事件的深度几乎相同(7 km),分别位于震中东北侧12 km和20 km处,它们的震源机制解(图 5a)也非常相似,值得指出的是该机制解的断层走向(230°±3°)和Futagawa断层的走向(232°)更为接近.从位置和断层几何来看,这两个子事件应当发生在Futagawa断层上,显示少量的正断层分量.最后一个子事件(MW6.54)最浅(2 km),其震源机制解显示斜向的滑动角,其走滑和正断层分量几乎相当.整个地震的震源持续时间大约20 s,最大的地震矩释放发生在~8 s.我们得到的反演结果和Yue等(2017)以及Hao等(2017)的主震有限破裂模型具有很好的一致性,尤其是和Hao等(2017).我们的多点源解所揭示的断层几何形态,包括最北段的低倾角性质也与Yue等(2017)通过余震分布约束得到的几何状态吻合.我们将多点源解相加得到的等效矩张量解和其他机构给出的单个点源矩张量解也具有高度一致性(图 4b),该一致性是对多点源解的有效验证(例如Wei et al., 2011; Avouac et al., 2014Shi et al., 2018).

4 讨论

反演的结果显示,随着地震破裂的发展,子事件(1~4)质心深度不断变浅并且越往后正断滑移的分量越大,其水平位置处的余震的深度也随之变浅.我们的子事件深度和Hao等(2017)Yue等(2017)的有限破裂模型结果非常一致.Hao等(2017)将其解释为断层发震层深度(或者温度)受附近的Aso火山影响所致,越靠近火山3~400 ℃等温线的深度越浅,Lythgoe等(2021)对2018年Lombok地震序列的解释也类似.第二和第三个子事件所在的位置(图 4a)正好对应了余震的空区,余震则分布在更深和更浅处,这和其他地震(例如Wei et al., 2011)的观测一致,可以很好地用应力变化来解释;Yue等(2017)也指出了这一点.然而,第四个子事件附近几乎没有余震,而且在沿着断层走向的方向约15 km处余震才重新出现.Yue等(2017)提出了地形应力加载对断层上正应力和剪应力的影响,其模型计算发现,尽管设置了倾角为65°的正断层,地形产生的正应力却显著升高反而不利于正断层滑动.除温度和地形应力的影响之外,不均匀分布的孔隙流体压力的潜在影响也不能忽视.Aizawa等(2017)对于熊本地震区域电阻率研究结果显示,Futagawa断层浅部近Aso火山部分存在明显高电导率,与高流体含量的特征相符.Matsumoto等(2018)利用Asano和Iwata(2016)的有限断层模型,结合历史地震活动,计算了Futagawa断层上的理论孔隙流体压力,其结果表明Aso火山附近的浅部孔隙流体压力超过了该断层的理论强度.这对断层黏滑特性的改变是比较可观的,同时也促进了主震破裂过程中该区域弹性势能的充分释放.值得注意的是,尽管本次地震序列中Aso火山口西北缘缺乏余震,但从2002年至2016年历史地震活动记录中,Uchide等(2016)发现这一区域实际上发生了大量地震.因此,这一区域的地震活动性有待未来进一步研究.

我们通过多点源反演所得到的四个子事件显示出Futagawa断层相比于Hinagu断层具有更多的正断层成分,尤其是其北段靠近Aso火山口附近的最后一个子事件拥有几乎同等大小的正断层和走滑分量.Yue等(2017)则指出该部分主震滑移包含较大比例的正断层成分.实际上,在相对远离Aso火山口的Futagawa断层部分,第二个和第三个子事件由走滑主导.如此差异可能是由于火山附近地形增加了正断层性质的相对剪切应力,而未同时增加走滑方向的剪切应力.精细的震后实地勘测(Shirahama et al., 2016)表明,Futagawa断层地表破裂主要为走滑性质,而其东南侧大约2 km处的Idenokuchi断层的地表破裂为正断层性质.结合野外测量和卫星遥感数据,Toda等(2016)提出了倾斜滑移分解模型,即Futagawa断层主要释放了长期弹性积累的走滑分量,而Idenokuchi断层主要释放了正断分量,从而形成小规模地堑.而靠近火山口的Futagawa断层可能缺乏这种走滑和正断层的分解机制,从而倾滑过程(正断+走滑)主导了同震位错.长期的地壳尺度的地表形变观测也可以为我们分析断层性质成因提供参考.通过对2000—2010年分布在九州岛的GNSS数据的反演,Mochizuki和Mitsui(2016)发现Aso火山附近尤其是其西缘的GPS台站(32.87°N,131.00°E)的下移分量会被均一结构模型远远低估;而在Aso下方添加一块与火山岩浆活动相关的席状形变源则可以充分模拟出该台站的下移.该结果表明,由于特殊的深部地质构造,Aso附近的Futagawa断层可能在长期地壳形变中积累比远离Aso区域更多的竖直方向的弹性应变,从而提高其地震活动释放的正断层分量.尽管从卫星干涉影像中并未获得最后一个子事件附近的正断层所产生的明显地表形变,但我们认为最后一个子事件所揭示的浅部倾滑的特性是可靠的,并且可能并未在地表产生容易被勘测到的正断层成分.

我们注意到位于断层上或者离断层很近的台站(GPS台0701和强震仪KMMH16)的波形拟合不是很理想(图 9c,浅蓝色波形).我们认为这可能是多点源解对有限破裂过程近似所造成的影响.因为断层上的台站更容易受到最近破裂的影响,而这些局部破裂振幅虽然不比其他破裂更大,但由于非常靠近台站,其近场项(随距离的立方衰减)波场振幅在该台站的记录更大,因此没有被反演.作为对比,我们参考了其他学者对于熊本地震震源过程的研究,发现大部分有限断层模型都会明显低估0701台的总水平位移;而Zhang等(2018)的模型则是包含了靠近0701的近地表高滑移区,从而可以充分拟合0701的静态位移,尽管其他GPS台站位移量拟合不是很好.为了更好的拟合在断层上的波形数据,必须要考虑到有限破裂过程的影响,使用有限破裂模型反演对这些近台数据的拟合将会有进一步的提高.为了测试GPS台站对于此次地震破裂究竟对破裂的矩心位置还是对最近的破裂更加敏感,我们增加了三组模拟测试(图 10).对于每组测试,先用六个点源构建一个输入震源模型.该输入震源模型包含两个滑移集中区,其中一个位于西北部,靠近熊本主震的震中,由一个MW6.7的点源表示;另一个位于东北部,靠近GPS观测站0701,由线性排列的五个MW6.3的点源组成(如图 10中红色圆圈所示).我们用构建的输入震源模型合成理论波形,并且加入相对强度为20%的高斯噪声,用于反演测试.三组模拟测试分别针对了第二个滑移集中区的矩心位置相对临近(图 10a)和远离GPS台0701(图 10b10c)的情况,其中图 10b10c分别对应该矩心位置位于0701的西侧和东侧的输入震源模型.在接下来的反演测试中,我们只用两个点源去恢复输入震源模型,并且使用了和实际数据一样的台站选择(包括近发震断层的GPS台站)和滤波频段.三组反演结果(图 10)很好地恢复了输入震源模型的矩心位置,(近似)震源时间函数以及震源机制解.因此,尽管两个点源简化了震源模型,加入接近断层的波形数据依然能有效恢复输入模型.仍应注意,在一些特殊情况下,即使使用了有限断层模型,拟合非常靠近断层的地震观测也存在挑战.比如在对断层结构复杂性缺乏充分认知的情况,或是台站所在位置的地质条件对波场传播产生巨大影响的情况.例如,Lee等(2019)对2018年MW6.4中国台湾花莲地震的有限断层反演中发现,尽管设置了三个断层来模拟近震数据,但是该模型远远低估紧靠断层的尤其是其中两个断层之间的台站的速度振幅.他们将这解释为该台站所具有的较强的场地效应以及厚沉积层的影响.对于MW7.1熊本地震,我们试图添加一个额外的子事件来模拟0701的波形.如图 9所示,通过联合反演包含0701在内的GPS数据、强震仪波形和宽频地震仪波形,获得五点源解(图 9a).和四点源解类似,五点源解依旧显示了自西北向东南方向的连续破裂过程,并且主要破裂区域的断面几何、震级大小和震源深度等都非常相似.但主要区别在于,五点源解将四点源解的最后一个子事件分为两个位于Aso火山口旁的更小的正断层为主的子事件(MW6.2-6.3),并且这两个子事件代表不同的断层走向(图 9a及其右下角小图).考虑这两个子事件的位置误差在内,它们反映出Futagawa断层在Aso火山口附近形成了两条分支断层,其中一条保持主断层的走向(大致232°,由第五个子事件表示),另一条更加近东西走向(大致250°,由第四个子事件表示).此处的断层分支很可能对地震破裂的停止起到了关键性作用.该断层分支与Aso火山活动以及火山特殊地质结构的关系值得进一步探究,包括对长期历史地震数据、火山监测数据、Aso区域地壳结构以及岩石学特征等可以提供较为全面的研究角度.上述五点源解可以非常好地拟合台站0701的波形,并且合成图和观测波形的互相关系数达到100%(图 9c).值得注意的是,尽管先前所得到的四点源解无法理想地拟合0701的观测波形,但是可以与五点源解几乎同等程度地拟合其余GPS台站(图 9b).总之,这印证了多点源反演方法的优势:在波形数据指导下,通过提高点源数量来测试观测数据对地震破裂过程的敏感性.比如在不使用紧靠Futagawa断层的0701GPS台站,四点源解可以充分拟合其余近震GPS、强震仪和区域宽频地震仪的波形数据;而更加反映局部断层破裂过程的五点源解则需要利用非常靠近该局部断层的0701台站的波形.这一经验突出了紧靠潜在的发震断层布置观测台站的重要性,这对普遍的地震破裂过程的研究非常关键(比如2018年MW6.4花莲地震).

图 10 有限破裂效应的模拟测试 (a) 输入震源模型由六个点源组成(红色),用来合成模拟波形,用于后续两点源反演.反演结果(蓝色)有效地恢复了输入震源模型中的两个滑移集中区的矩心位置、震源机制解以及震源时间函数(见右下小图).输入震源模型和反演结果所对应的GPS台站0701的波形三分量分别由左上小图中的红线和蓝线表示.滤波频段为0.001~0.1 Hz. 图(b)和(c)与(a)的绘制方法类似,对应另外两组模拟测试. (b)中第二个滑移集中区相对于(a)朝西南方向平移了大约2 km.相反,(c)中第二个滑移集中区相对于(a)朝东北方向平移了大约2 km. Fig. 10 Synthetic tests for the finite rupture effect (a) The input six-point-source model (red) are used to compute the synthetic waveforms for a two-point-source inversion. The inversion result (blue) well recovered the centroid locations, focal mechanisms and source-time functions (lower-right inset) of the two asperities of the input source model. The three-component waveforms computed for the input and inverted source models are shown in red and blue, respectively, for GPS station 0701 and filtered to 0.001~0.1 Hz. (b) and (c)Similar to (a), but for another two synthetic tests. The centroid locations of the second asperity in (b) is shifted to southwest by about 2 km relative to (a). The centroid locations of the second asperity in (c) is shifted to northeast by about 2 km relative to (a).

此外,0701台的强振幅波形有可能包含了其他同震地质事件的贡献,比如该地震序列引发的最大的山体滑坡(Hung et al., 2018; Song et al., 2019)刚好发生在0701西北方约1 km处.Hung等(2018)认为该滑坡附近的强震仪明显记录到了滑坡造成的位移波形,并由此确定滑坡对应的波形起始时间为主震发震时刻之后18.45~21 s.这一结果支持该最大滑坡是由主震动态触发.但是是否可以从0701台站的波形数据中分离出明显的对应滑坡的成分,这需要对滑坡的运动学过程进行进一步模型设置,这超出了本文的探究范围.我们用剪切位错(双力偶)源对0701波形的拟合已接近完美,并且得到的断层走向和地表破裂极其吻合,我们认为滑坡对0701波形的贡献应当比较有限.

熊本地震的宽频带台站数据质量相当高,并且主要的传播路径是在大陆上,如果地质条件更为复杂,例如地震波传播要经过大陆和海洋,尤其是经过海沟,或者是发生在海沟附近的地震,其对应的地震波记录会更加复杂.在此情况下,一维格林函数就不再适用,如果仍旧使用一维格林函数则必须要使用更长周期来对数据进行滤波,那么对震源破裂过程的分辨率就会下降.而Qian等(2019)Wu等(2018)所提出的方法和应用表明,三维格林函数对源区结构复杂的波场传播在10~20 s的周期是不可忽略的,这些研究也表明在没有很厚沉积层的情况下考虑地形和水层就可以对波形拟合大为改进.因此使用三维格林函数进行多点源反演是将来一项重要的工作.我们注意到Duputel和Rivera(2017)在地震多点源反演中使用了全球SEM(Spectral Element Method,谱元法)程序(Komatitsch和Tromp, 2002)计算的三维格林函数,但是全球版SEM程序对海水使用了应力近似边界条件,并不能模拟水波及其多次波.可能也是由于这个原因,Duputel和Rivera(2017)尽管使用了三维格林函数,依然用了50 s以及更长周期的地震波来对2016年发生的MW7.8 Kaikouro地震进行反演.

我们所使用的正则化的Yoffe函数在多点源地震反演更符合震源破裂动力学(Tinti et al., 2005).在地震破裂过程的实际研究中,由于震源时间函数和其他震源参数之间存在取舍关系,并且高频地震波的快速衰减不容忽视,因此破裂过程的震源时间函数往往极难分辨.显然,在波形拟合可靠(经过校正)的情况下,反演所用的频段越高就越能分辨破裂过程的震源时间函数.我们的结果显示前三个子事件的震源时间函数更接近于高斯函数,在速度和位移波形反演中均是如此.最后一个子事件在速度波形反演结果中更接近于Yoffe函数.对比于离该子事件最近的几个台站波形(图 7),我们认为这是一个可靠的结果.这说明浅部的地震破裂很可能受到前面子事件造成的动态应力传播的影响.在反演所用的频段,我们不需要具有强不对称性的Yoffe函数来拟合数据,这有两种解释:(1)熊本地震的主震的震源时间函数确实更接近于高斯函数;(2)需要用更高频段的地震波来进一步分析.我们提出的正则化的Yoffe函数震源时间函数可以在其他地震的破裂过程研究中使用,其灵活性将有助于我们进一步了解破裂形成和发展的物理机理.

5 结论

我们将一个基于MCMC的地震多点源波形反演方法扩展到了区域距离的地震波形数据,并对2016年MW7.1熊本地震进行了应用.结果显示区域地震数据可以较好地分辨此类中强度(M6-7)地震的时空破裂过程.利用离断层很近的0701高频GPS数据可以分辨出地震破裂最后阶段更为精细的破裂过程.我们可以将该方法进一步推广到全球中强度地震的运动学过程研究中,来更准确地确定基本的震源物理参数,以理解震源物理过程.

致谢  谨此祝贺陈颙先生从事地球物理教学科研工作60周年.感谢两位匿名审稿人宝贵的建设性意见,感谢日本Kik-net,K-net,F-net为本文研究所提供的波形数据.
References
Aizawa K, Asaue H, Koike K, et al. 2017. Seismicity controlled by resistivity structure: the 2016 Kumamoto earthquakes, Kyushu Island, Japan. Earth, Planets and Space, 69(1): 4. DOI:10.1186/s40623-016-0590-2
Aki K, Richards P G. 2002. QuantitativeSeismology. Sausalito: University Science Books.
Asano K, Iwata T. 2016. Source rupture processes of the foreshock and mainshock in the 2016 Kumamoto earthquake sequence estimated from the kinematic waveforminversion of strong motion data. Earth, Planets and Space, 68(1): 147. DOI:10.1186/s40623-016-0519-9
Avouac J P, Ayoub F, Wei S J, et al. 2014. The 2013, MW7.7 Balochistan earthquake, energetic strike-slip reactivation of a thrust fault. Earth and Planetary Science Letters, 391: 128-134.
Chen Y, Wen L X, Ji C. 2014. A cascading failure during the 24 May 2013 great Okhotsk deep earthquake. Journal of Geophysical Research: Solid Earth, 119(4): 3035-3049. DOI:10.1002/2013JB010926
Duputel Z, Rivera L. 2017. Long-period analysisof the 2016 Kaikoura earthquake. Physics of the Earth and Planetary Interiors, 265: 62-66. DOI:10.1016/j.pepi.2017.02.004
Hao J L, Ji C, Yao Z X. 2017. Slip history of the 2016 MW7.0 Kumamoto earthquake: Intraplate rupture in complex tectonic environment. Geophysical Research Letters, 44(2): 743-750.
Hung C, Lin G W, Syu H S, et al. 2018. Analysis of the Aso-bridge landslide during the 2016 Kumamoto earthquakes in Japan. Bulletin of Engineering Geology and the Environment, 77(4): 1439-1449. DOI:10.1007/s10064-017-1103-7
Ji C, Wald D J, Helmberger D V. 2002. Source description of the 1999 Hector Mine, California, earthquake, part I: Wavelet domain inversion theory and resolution analysis. Bulletin of the Seismological Society of America, 92(4): 1192-1207. DOI:10.1785/0120000916
Jia Z, Shen Z C, Zhan Z W, et al. 2020. The 2018 Fiji MW 8.2 and 7.9 deep earthquakes: One doublet in two slabs. Earth and Planetary Science Letters, 531: 115997. DOI:10.1016/j.epsl.2019.115997
Kato A, Nakamura K, Hiyama Y. 2016. The 2016 Kumamoto earthquake sequence. Proceedings of the Japan Academy, Series B, 92(8): 358-371. DOI:10.2183/pjab.92.359
Kobayashi T. 2017. Earthquake rupture properties of the 2016 Kumamoto earthquake foreshocks (Mj 6.5 and Mj 6.4) revealed by conventional and multiple-aperture InSAR. Earth, Planets and Space, 69(1): 7.
Komatitsch D, Tromp J. 2002. Spectral-element simulations of global seismic wave propagation-I.Validation. Geophysical Journal International, 149(2): 390-412. DOI:10.1046/j.1365-246X.2002.01653.x
Lee S J, Lin T C, Liu T Y, et al. 2019. Fault-to-fault jumping rupture of the 2018 MW6.4 Hualien Earthquake in eastern Taiwan. Seismological Research Letters, 90(1): 30-39.
Lythgoe K, Muzli M, Bradley K, et al. 2021. Thermal squeezing of the seismogenic zone controlled rupture of the volcano-rooted flores thrust. Science Advances, 7(5): eabe2348. DOI:10.1126/sciadv.abe2348
Matsumoto S, Yamashita Y, Nakamoto M, et al. 2018. Prestate of stress and fault behavior during the 2016 Kumamoto earthquake (M7.3). Geophysical Research Letters, 45(2): 637-645. DOI:10.1002/2017GL075725
Mochizuki K, Mitsui Y. 2016. Crustal deformation model of the Beppu-Shimabara graben area, central Kyushu, Japan, based on inversion of three-component GNSS data in 2000-2010. Earth, Planets and Space, 68(1): 177. DOI:10.1186/s40623-016-0550-x
Nishimura S, Hashimoto M. 2006. A model with rigid rotations and slip deficits for the GPS-derived velocity field in southwest Japan. Tectonophysics, 421(3-4): 187-207. DOI:10.1016/j.tecto.2006.04.017
Qian Y Y, Wei S J, Wu W, et al. 2019. Teleseismic waveform complexities caused by near trench structures and their impacts on earthquake source study: Application to the 2015 Illapel aftershocks (Central Chile). Journal of Geophysical Research: Solid Earth, 124(1): 870-889. DOI:10.1029/2018JB016143
Shi Q B, Wei S J, Chen M. 2018. An MCMC multiple point sources inversion scheme and its application to the 2016 Kumamoto MW6.2 earthquake. Geophysical Journal International, 215(2): 737-752. DOI:10.1093/gji/ggy302
Shi Q B, Wei S J. 2020. Highly heterogeneous pore fluid pressure enabled rupture of orthogonal faults during the 2019 Ridgecrest MW7.0 earthquake. Geophysical Research Letters, 47(20): e2020GL089827. DOI:10.1029/2020GL089827
Shirahama Y, Yoshimi M, Awata Y, et al. 2016. Characteristics of the surface ruptures associated with the 2016 Kumamoto earthquake sequence, central Kyushu, Japan. Earth, Planets and Space, 68(1): 191. DOI:10.1186/s40623-016-0559-1
Song K, Wang F W, Dai Z L, et al. 2019. Geological characteristics of landslides triggered by the 2016 Kumamoto earthquake in Mt.Aso volcano, Japan. Bulletin of Engineering Geology and the Environment, 78(1): 167-176. DOI:10.1007/s10064-017-1097-1
Tinti E, Fukuyama E, Piatanesi A, et al. 2005. A kinematic source-time function compatible with earthquake dynamics. Bulletin of the Seismological Society of America, 95(4): 1211-1223. DOI:10.1785/0120040177
Toda S, Kaneda H, Okada S, et al. 2016. Slip-partitioned surface ruptures for the MW7.0 16 April 2016 Kumamoto, Japan, earthquake. Earth, Planets and Space, 68(1): 188. DOI:10.1186/s40623-016-0560-8
Uchide T, Horikawa H, Nakai M, et al. 2016. The 2016 Kumamoto-Oita earthquake sequence: aftershock seismicity gap and dynamic triggering in volcanic areas. Earth, Planets and Space, 68(1): 180. DOI:10.1186/s40623-016-0556-4
Wei S J, Fielding E, Leprince S, et al. 2011. Superficial simplicity of the 2010 El Mayor-Cucapah earthquake of Baja California in Mexico. Nature Geoscience, 4(9): 615-618. DOI:10.1038/ngeo1213
Wu W B, Ni S D, Zhan Z W, et al. 2018. An SEM-DSM three-dimensional hybrid method for modelling teleseismic waves with complicated source-side structures. Geophysical Journal International, 215(1): 133-154. DOI:10.1093/gji/ggy273
Xu L S, Chen Y T. 2002. Source time function and rupture process of earthquake. Seismological and Geomagnetic Observation and Research (in Chinese), 23(6): 1-8.
Xu L S, Chen Y T. 2005. Temporal and spatial rupture process of the great Kunlun Mountain Pass earthquake of November 14, 2001 from the GDSN long period waveform data. Science in China Series D: Earth Sciences (in Chinese), 34(3): 256-264.
Yagi Y, Okuwaki R, Enescu B, et al. 2016. Rupture process of the 2016 Kumamoto earthquake in relation to the thermal structure around Aso volcano. Earth, Planets and Space, 68(1): 118. DOI:10.1186/s40623-016-0492-3
Yao Z, Ji C. 1997. The inverse problem of finite fault study in time domain. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 40(5): 691-701.
Yoffe E H. 1951. LXXV.The moving griffith crack. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 42(330): 739-750. DOI:10.1080/14786445108561302
Yue H, Lay T, Koper K D. 2012. En échelon and orthogonal fault ruptures of the 11 April 2012 great intraplate earthquakes. Nature, 490(7419): 245-249. DOI:10.1038/nature11492
Yue H, Ross Z E, Liang C R, et al. 2017. The 2016 Kumamoto MW=7.0 earthquake: A significant event in a fault-volcano system. Journal of Geophysical Research: Solid Earth, 122(11): 9166-9183.
Yue H, Lay T. 2020. Resolving complicated faulting process using multi-point-source representation: Iterative inversion algorithm improvement and application to recent complex earthquakes. Journal of Geophysical Research: Solid Earth, 125(2): e2019JB018601. DOI:10.1029/2019JB018601
Yue H, Zhang Y, Gai Z X, et al. 2020. Resolving rupture processes of great earthquakes: Reviews and perspective from fast response to joint inversion. Science China Earth Sciences, 63(4): 492-511. DOI:10.1007/s11430-019-9549-1
Zhang Y, Feng W P, Xu L S, et al. 2009. Spatio-temporal rupture process of the 2008 great Wenchuan earthquake. Science in China Series D: Earth Sciences, 52(2): 145-154. DOI:10.1007/s11430-008-0148-7
Zhang Y, Feng W P, Xu L S, et al. 2008. Spatio-temporal rupture process of the 2008 great Wenchuan earthquake. Science in China Series D: Earth Sciences, 38(10): 1186-1194.
Zhang Y F, Shan X J, Zhang G H, et al. 2018. Source model of the 2016 Kumamoto, Japan, earthquake constrained by InSAR, GPS, and strong-motion data: Fault slip under extensional stress. Bulletin of the Seismological Society of America, 108(5A): 2675-2686. DOI:10.1785/0120180023
Zhao L S, Helmberger D V. 1994. Source estimation from broadband regional seismograms. Bulletin of the Seismological Society of America, 84(1): 91-104.
Zhu L P, Helmberger D V. 1996. Advancement in source estimation techniques using broadband regional seismograms. Bulletin of the Seismological Society of America, 86(5): 1634-1641. DOI:10.1785/BSSA0860051634
许力生, 陈运泰. 2004. 震源时间函数与震源破裂过程. 地震地磁观测与研究, 23(6): 1-8. DOI:10.3969/j.issn.1003-3246.2004.06.001
许力生, 陈运泰. 2004. 从全球长周期波形资料反演2001年11月14日昆仑山口地震时空破裂过程. 中国科学D辑: 地球科学, 34(3): 256-264.
姚振兴, 纪晨. 1997. 时间域内有限地震断层的反演问题. 地球物理学报, 40(5): 691-701. DOI:10.3321/j.issn:0001-5733.1997.05.010
岳汉, 张勇, 盖增喜, 等. 2019. 大地震震源破裂模型: 从快速响应到联合反演的技术进展及展望. 中国科学: 地球科学, 50(4): 515-537.
张勇, 冯万鹏, 许力生, 等. 2008. 2008年汶川大地震的时空破裂过程. 中国科学D辑: 地球科学, 38(10): 1186-1194. DOI:10.3321/j.issn:1006-9267.2008.10.002