岳 汉
北京大学地球与空间科学学院,北京 100871
大地震依然是造成最多人类伤亡的自然灾害.一次大地震发生后,解析其破裂过程是指导震后救灾的重要内容,也是实现长期地震危险性研究的重要基础研究支撑(Yue et al., 2020). Yue 等(2020)讨论了震后快速响应以及破裂过程联合反演在震后应急中的重要性,集中分析了多数据联合反演对于破裂过程解析度的互补以及各个方法在大地震快速响应中的应用,而对于“复杂破裂过程”的破裂过程反演以及应急响应,现有研究框架尚未完全包括. 本文针对复杂破裂过程研究中的灵活性需求进行简要概述,并提出结合反投影技术和多点源反演算法研究复杂地震破裂过程的技术路线.
大地震破裂过程研究的主要思路是:
(1)假设破裂断层面为矩形,通过震源机制解获得断层面几何参数,并通过震级估计其尺度.
(2)参考震源机制以及区域构造特征选择合适的断层面.
(3)搜集相关数据反演其破裂过程.
这种平面断层的假设在很多情况下可以满足破裂过程快速响应的需求,但是在破裂过程复杂的很多大陆内部强震中,平面断层假设无法满足解析破裂过程的需求. 现代地震学和大地测量测数据显示,很多大地震地表破裂往往不局限于单一断层面而包括若干条断层的先后或者同时破裂(如 Sun et al.2018; Wang et al., 2018). 这些现象并不是因为地震变得更加复杂,而是因为观测让我们能够解析到更多的地震破裂细节. 例如2012 年印度洋底地震(Meng et al., 2012; Yue et al., 2012)、2016 年新西兰Kaikoura 地震(Duputel and Rivera, 2017; Wang et al., 2018)、2013 年巴基斯坦地震(Avouac et al., 2014)、2016 年日本熊本地震(Yue et al.,2017)、2019 年美国加州Ridgecrest 地震(Wang et al., 2020; Yue et al., 2021)等,这些地震在断层几何形态以及分段特征上都发生了明显改变,破裂面都不限于单一平面断层. 对于若干中等规模地震,如6~7 级地震,科研人员也发现了地震的复杂性,如2014 年云南鲁甸地震(刘成利,2013,2014)、2017 年四川九寨沟地震(Sun et al., 2018)、2018年台湾花莲地震等(Lo et al., 2019). 应用了点源解的单一平面来参数化断层面,往往让这些地震的有限断层反演无法准确表达地震的破裂机制,为地震危险性的快速估计带来困难.
对于大陆内部复杂地震而言,研究人员通常可以参考地质考察、卫星影响以及余震分布刻画非平面断层. 地表勘察以及卫星影像数据往往可以得到比较精细的地表断层形态. 而精定位后的余震目录往往可以在深度分布上刻画断层形态. 利用地表和余震数据,科研人员可以刻画出精细的断层形态,作为破裂过程反演的参考模型. 例如,2016 年日本熊本地震(Yue et al., 2017)、2019 年美国加州Ridgecrest 地震(Ross et al., 2019)等,都有研究人员通过上述方法获得了多段非平面断层面. 然而,卫星数据的获得以及余震目录的获取往往需要若干天的时间,无法满足快速应急的需求. 此外,发生在海中的地震,往往没有足够精细的表层观测以及余震目录实现上述断层刻画过程. 因此对于复杂地震的破裂过程,我们尚需要一套可以在全球范围内普遍适用的、可以反演破裂在多断层破裂上的发展过程的、灵活度较高的反演方法. 下文结合大地震研究中四种常用的方法进行简要介绍.
对于复杂破裂过程,基于远震地震波的多种反演技术为非单一平面的破裂过程反演提供了依据. 本文列举了四种基于远震地震波的重要方法,分别为单点源震源机制解、多点源震源机制解、反投影技术以及有限断层模型. 对于一次大地震的复杂破裂过程,上述四种方法都可以发挥作用. 如图1 所示是针对2017 年的Kaikoura 地震,应用几个不同研究方法的成果. 这些方法都反映了地震东北向的传播特征以及主要的两次能量辐射. 这些研究结果都可以通过远震地震波得到. 远震地震波是震源快速响应的重要数据. 因为其具有全球可得的特点,可以形成一套全球地震普遍适用的方法体系,为全球范围内复杂地震的快速响应提供了可能. 这些基于远震地震波的方法中,有限断层模型反映了最为准确的破裂细节,但是由于有限断层模型需要先将断层面参数化,对于复杂破裂过程而言,无法通过单点源解的最佳破裂面得到有效的断层结构. 多点源方法和反投影技术(上游技术)为有限断层模型(下游技术)的参数化提供了依据,但是这些方法之间的衔接目前依然依赖于研究人员的经验和手动调试,这一过程通常比较耗时,且无法在快速响应中应用. 因此实现复杂地震快速响应的瓶颈不在数据层面,而在技术层面.复杂地震破裂过程的研究迫切需要一种可以将上游结果作为下游输入参数的自动参数化手段,实现不同反演技术之间的衔接,实现整个流程自动化.
图 1 2016 年新西兰Kaikoura 地震的多种震源研究方法比较. (a)GCMT 和Wphase 的点源反演结果;(b)多点源反演结果(修改自Duputel and Rivera, 2017);(c)反投影成像结果(修改自Hollingsworth et al., 2017);(d)有限断层反演结果(修改自Wang et al., 2018)Fig. 1 Comparison between different inversion results for the 2016 Kaikoura earthquake. (a) GCMT and Wphase point source solution; (b) Multi-point-source solution (modified from Duputel and Rivera, 2017); (c) Back-projection (modified from Hollingsworth et al., 2017); (d) Finite fault model inversion (modified from Wang et al., 2018)
上述四种方法在学界内还处于大量单独应用的“单兵作战”模式. 而在单独应用各个方法时通常因为方法本身的缺陷而无法了解地震的全面信息.下面对于各个方法的优缺点进行简要的介绍. 由于单点源、反投影、有限断层等方法的相关文章较多,本文只做简要介绍. 而对于多点源方法,由于其满足了震源研究灵活性的需求,本文将做重点介绍.
数学上可以证明,发生在断层面上的点位错源产生的位移场可以等效表示为双力偶源产生的位移场(震源表示定理). 因此用空间中集中于一点的双力偶源描述震源模型(震源机制解)成为描述地震破裂模型的基本形式. 震源机制解将震源描述为一个无限小的点源,描述地震断层面的几何特征,不包括尺度范围,也可以称为单点源解. 目前,基于远震波形的单点源反演,如GCMT 点源解、(Ekström et al., 2012)以 及Wphase(Kanamori and Rivera, 2008; Duputel et al., 2012b),已经自动成为监测全球大地震机制的重要手段. 不同于GCMT 和Wphase 的反演原理,Zhu 和Helmberger(1996)提出了用最大化区域震相时移后的相关函数的方法,建立了搜索震源机制的Cut-and-Paste(CAP)方法,对于中小型地震的区域数据反演震源机制解, 并在大量区域地震的研究中得到应用(如Zheng et al., 2009; Zheng et al., 2010; Zhao et al., 2013). 除了以波形为基础的反演方法之外,利用P 波初动、P/S 振幅比等观测方法也一直在大地震点源解的反演中发挥着重要的作用(如Dreger and Helmberger, 1993; 胡新亮等, 2004; 胡幸平等,2008). 在时间域单点源反演中通常不将破裂过程近似为脉冲源,而一般假设震源时间函数近似为一个三角形,通过搜索三角形时长的方法估计震源持续时间,如GCMT. 在有经验格林函数的情况下,可以假设大地震和参考地震发生的位置相同,然后通过反卷积的方法提取比较精细的震源时间函数(Mueller, 1985).
N
方差(Xu Y et al., 2009)、Multi-taper(Meng et al., 2011)等技术的应用使得反投影技术可以快速自动化地成像大地震破裂过程. 不同于震源反演技术,反投影技术应用远震地震波的高频信息,不需要对破裂过程做出先验假设,可以快速得到破裂过程的时空分布. 结合理论和经验的到时校正,其对于破裂过程中能量位置的解析度可以达到10~20 km 的量级(Yue et al., 2020).反投影叠加得到的地震能量对于震源空间相关性比较敏感,并不与地震矩直接相关. 因此反投影技术具有快速成像、空间解析度高的优势,但是对于地震矩、震源机制没有解析度.N
个点源的目录就比单点源多N
倍的信息.Kikuchi 和Kanamori(1982, 1986, 1991)自1980年代开始研发基于远震体波的多点源反演算法,利用时间域迭代反演算法(剥皮法)实现了大地震的多点源反演. 本世纪的大地震记录给多点源更多的发展空间,基于不同反演算法的多点源反演工作也在多个地震中得到了应用. 2004 年苏门答腊地震后,Tsai 等(2005)借鉴GCMT 反演的方法研发出multi-CMT 反演,在限定地震位置的前提下反演其机制和地震矩. 对于这次9.0+级的地震,GCMT 使用的地震波频宽已经达到饱和,无法正确获得地震矩. GCMT 地震矩为9.0,明显小于地球简正模(Park et al., 2005)以及大地测量数据(Banerjee et al., 2005)给出 的 地震 矩(9.1~9.3). Tsai 等(2005)将苏门答腊地震划分成5 个点源,通过多次面波进行反演,反演得到的地震矩达到9.3,震源破裂时间为500 s,更符合区域数据得到的破裂过程. 2012 年印度洋底地震后, Duputel 等(2012)应用马尔科夫链-蒙特卡罗(Markov Chain Monte Carlo, MCMC)方法,将这次地震分解为两个点源,并且反演了点源的震级、位置和深度. 其得到的破裂模型与有限断层和反投影技术得到的主要断层破裂一致(如Meng et al., 2012; Yue et al., 2012).2016 年Kaikoura 地震后,Duputel 和Rivera (2017)利用相似技术,应用全球三维速度模型计算的波形格林函数,用四个点源模拟了破裂过程. 该模型表示,破裂后期的主要能量释放阶段,包含了两个时间、空间和震级相近、而机制不同的逆冲和走滑分量的地震(M
7.6),分别对应了Kekerengu 断层以及俯冲板块上的破裂. 2018 年斐济深震后,Jia 等(2019)通过MCMC 方法,解析了发生于毗连的两个板块上的两个点源解. 对于中等规模的地震(M
6~7),多点源反演同样能发挥作用. Shi 等(2018)研发了基于马尔科夫链-蒙特卡罗(MCMC)的算法,对于2016 年6.2 级熊本地震前震进行了多点源反演. 对于中等规模的地震(6~7 级),作者应用时间域迭代反演算法,并且结合远震、近震数据对于2017 年九寨沟地震(Sun et al., 2018)和2018 年花莲地震(Lo et al., 2019)进行了多点源反演. 在九寨沟地震中,Sun 等(2018)通过波形初动分析以及多点源反演的方法,发现了破裂初期包括了与主断层共轭的M
5.6 的逆冲型地震(图2).在2018 年花莲地震的研究中,Lo 等(2018)结合远震体波以及近震波形,用6 个点源解释这次地震的破裂过程(图3),得出了从低角度走滑+逆冲转化为高角度走滑、最终到垂直走滑的破裂机制,与通过大地测量学以及地震学手段得到的多断层破裂模型符合较好. 该研究显示花莲地震涉及了外海板间断层,以及板内走滑断层(美伦断层与岭顶断层)的先后破裂.
上述震例中,在断层几何形态不明朗的情况下,多点源方法作为具有较高自由度的方法,通过多个震源机制解反映破裂过程中断层面几何形态的变化. 在这些震例中,多点源算法的共同优势是可以提高反演的自由度,解析复杂的破裂过程. 自由度高的优势也带来稳定性降低的不足. Yue 和 Lay(2020)针对反演稳定性提出了基于时间域反演算法的先验信息的应用,通过限制点源时间函数、位置以及机制组合等策略增加对于单点源时空以及机制的约束,如图4 所示.
Yue 和Lay(2020) 将改进的多点源反演方法应用于过去10 年内发生的7 个大地震中, 通过与有限断层模型的对比,系统地讨论了基于远震数据的多点源反演的时空解析度以及反演结果稳定性问题. 并指出,基于远震体波的多点源反演的空间解析度上限在20 km 左右,该误差范围对于参数化断层模型依然比较粗糙. 因此,单纯基于多点源反演的结果很难给出准确的破裂面几何形态. Sun 等(2018)与Lo 等(2019)结合了大地测量数据、区域地震目录以及多点源反演结果,构造了包含多个断层面的断层破裂模型,但是这些分析强烈依赖于大地测量数据,时效性不高,无法在地震快速响应方面取得突破.
图 2 (a)2017 年九寨沟地震的发震位置,其南侧分别发生了1973 年黄龙地震、1976 年的松潘震群. (b)2017 年九寨沟地震多点源反演. GCMT 点源解表示为大的震源球,左侧4 个远震体波波形以及单点源和多点源拟合的波形表示为黑色、蓝色和红色的波形,多点源反演的结果呈现于右下方,多点源叠加的震源机制解和有限断层叠加的震源机制解画在右上方(修改自Sun et al., 2018)Fig. 2 (a) 2017 Jiuzhaigou earthquake location. 1973 Huanglong earthquake and 1976 Songpan earthquake occurred on the south Huya fault. (b) Multi-point source inversion results of the 2017 Jiuzhaigou earthquake. Example waveforms of four southern teleseismic stations are plotted. Waveforms fitted by single point and multi-point solutions are plotted in blue and red, respectively. Solutions of 4 point sources of the MPS inversion result are plotted in the bottom (modified from Sun et al., 2018)
图 3 (a)2018 年花莲地震的发生位置以及区域构造. 花莲地震发生于欧亚和菲律宾板块的交接处,是造山和俯冲构造的交界点,破裂过程复杂. (b)花莲地震多点源反演结果,不同颜色的震源球表示了远震体波和近场三分量地震波形反演的多点源解,震源球的颜色代表发震断层:紫色为板间断层,蓝色为美伦断层,绿色为岭顶断层. 每个点源发生的时间、震级以及深度在右下角呈现(修改自Lo et al., 2019)Fig. 3 (a) 2018 Hualien earthquake rupture. (b) MPS inversion results of the 2018 Hualien earthquake. Point sources on the off-shore fault, Meilun fault and Linding faults are plotted in magenta, blue and green focal mechanisms, respectively (modified from Lo et al., 2019)
不同于震源机制解,有限断层模型注重于描述断层面的空间尺度,断层面的几何形态也不限制于单一平面. 如果有限断层的反演中包含动态信息,如远震体波、面波、区域震相等,有限断层模型通常也可以描述地震破裂过程. 相对于震源机制解而言,有限断层模型更加接近于真实的地震破裂模型,与野外地质以及大地测量观测更加相关,同时也更加有助于强地面震动的模拟以及地震灾害的估计.自从有了有限断层模型的先驱性工作(Trifunac,1974),经过地震学家、大地测量学家数十年的努力,有限断层的反演方法已经日趋成熟,并逐渐成为大家接受的并且广泛使用的震源表示方法. 目前快速有限断层反演在我国以及美国都已投入使用.有限断层的快速反演通常可以在若干小时之内对于大地震的破裂过程给出第一手的模型,为救灾工作提供参考(张培震等, 2009; Hayes, 2017).
快速有限断层反演一般依赖于地震的点源解预设断层面,然后通过远震体波和面波进行有限断层的反演. 对于大多数地震,单一断层面的有限断层反演可以囊括地震的主要破裂特征. 然而,现在地震学以及大地测量学发现在一些大地震破裂过程中,可以发生多个断层的先后破裂,而破裂在多个断层之间跳跃的过程中可能发生震源机制的变化. 过去几十年的观测已经积累了一些地震在多断层之间跳跃的破裂过程的案例. 比如,2008 年汶川地震破裂了龙门山破裂带的多条断层,震源机制经历了逆冲到走滑断层的变化(Wang et al., 2008; Shen et al.,2009; Xu X et al., 2009; Zhang et al., 2009; Zheng et al., 2009);2012 年印度洋底8.6 级地震在空间上包含了5~6 个走滑型断层的先后破裂,破裂范围超过400 km(Meng et al., 2012; Yue et al., 2012);2016 年新西兰Kaikoura 地震产生了至少10 个断层的地表破裂,同时也造成了俯冲板块表面的破裂(Wang et al., 2018; Xu et al., 2018). 这些地震的复杂破裂过程无法用单一断层面的破裂模型进行解释,用单断层的有限断层反演算法将会对破裂过程造成系统偏差.
图 4 (a)第一行的左右分别为生成模拟波形的点源位置以及时间函数. 第二和第三行分别为有先验约束的点源反演以及无先验约束的点源反演的多点源机制以及时间函数. (b)模拟波形以及每个点源产生的模拟波形. 图(b)左右分别为有无先验约束的波形拟合(修改自Yue and Lay, 2020)Fig. 4 (a) Mechanisms and source time functions of synthetic earthquakes are plotted in the left and right panels of the top row,respectively. MPS inversion results with and without a priori constraints are plotted in the second and third rows. (b) Waveforms fit by either MPS inversions are plotted in the bottom panels. Synthetic and modeled waveforms are plotted in black and red curves in the top panel. Synthetic waveforms of each point source subevent are plotted as color coded waveforms at the bottom (modified from Yue and Lay, 2020)
相对于日本、智利等国家的俯冲带地震而言,我国的地震灾害主要来自于大陆内部的板内地震,更适合现代地球物理学观测的解析能力. 这一特征适应了现代地球物理学以及大地测量学的观测趋势,也给我国的地震学研究提出了挑战. 对于大陆内部的复杂板内地震,依然需要研发快速的、可以广泛推广的断层面划分算法,在地震破裂模型的快速反演、提高震灾估计以及救援的时效方面提供有效的工具.
在复杂地震破裂过程的反演工作中,大地测量数据、地震地质资料、地震波波形分析等方法通常是划分复杂断层几何的依据,而地震的快速响应几乎必须依赖波形数据. 上述三种反演方法(多点源、反投影、有限断层)都可以基于远震体波进行处理, 在解析破裂过程中各有优缺点,在功能上可以互补. 三种方法的优缺点如表1所示.
总体来说, 先验约束是在反演稳定性不高的时候通过经验或者其他信息缩小参数空间的策略,因此反演稳定性越高的算法往往需要的先验条件越少,但是这些方法对于破裂的多方面特征(如地震矩、震源机机制等)往往解析度不完备. 反投影方法需要的预设信息最少,通常只需要给出震源区的范围以及网格信息,就可以对任意复杂的破裂过程给出高时间和空间精度的成像. 但是反投影技术应用高频信息叠加而成,叠加方法也通常不是线性方法,无法恢复地震矩以及震源机制信息. 多点源反演在强的先验模型下才能得出比较稳定的破裂模型,但是对于点源的时间和空间解析度较差. 有限断层模型在所有方法中对于先验约束的依赖性最强,通常需要对断层面几何形态以及破裂运动学过程给出比较完备的描述,其反演结果的时空解析度以及地震矩的解析度也是最高的.
从表1 的方法中可以看出,从上游到下游(反投影>多点源>有限断层)遵循先验约束逐渐增加,进而模型解析度逐渐增加的特征. 在逻辑上比较合理的策略是将上游方法的结论作为先验约束输入下游方法,从而约束下游方法反演结果的稳定性. 在一些综合分析中(如Yue et al., 2012; Sun et al., 2018; Lo et al., 2019),也依赖研究人员的主观判断,将上游的信息转换成为下游方法输入参数.因此,将上述方法串行需要更加客观的标准以及算法,才能够实现快速、普适的震源模型参数化策略.
表 1 反投影、多点源以及有限断层模型反演方法的先验约束以及解析度分析Table 1 A priori constraints and resolution analysis of back projection, multi-point source and finite fault inversion methods
上文分析了算法从上游到下游在先验约束、结果稳定性上的差别,因此通过先验约束少的上游算法逐步增加对下游算法的约束,在实现算法灵活性的同时,增加其稳定性的可行性. 本文提出一种通过反投影结果约束多点源先验参数,再通过多点源反演结果构造多断层几何模型,从而增加破裂过程反演自由度的串行方法. 下文就串行方法的若干思路进行简要描述. 反投影算法可以有两种呈现方式:(a)每一时刻的集中辐射点源的离散表示(如图2 反投影结果下图);(b)聚束能量场的时空表示(如图2 反投影结果上图). 而每一种表示方法与现有多点源算法有比较容易的串行方式. 本文提出两种从反投影到多点源串行算法的思路,并针对其各自的优势进行分析.
(1)离散表示结果与时间域多点源算法的串联
反投影结果的离散表示与时间域迭代算法的输入比较契合. 在离散表示中,反投影的能量呈现为空间离散的若干个点源,这与时间域迭代算法需要提前设定的空间离散的震源位置情形类似. 而每个离散点源的能量释放时间函数(包括起始时间以及震源时长)可以用反投影聚束结果给出比较好的约束. 因此,通过反投影聚束结果,我们可以比较容易给出各个点源的位置、时间信息. 这一方式极大地减少了时间域多点源反演中对于预设断层面位置、网格划分、点源起始时间以及时间函数的需求. 而多点源反演可以在目标点源位置和时间获得其震源机制信息. 在点源个数方面,也可以根据反投影点源的能量大小,从大到小筛选若干个点源进行反演.
(2)聚束能量场时空分布与MCMC 多点源算法串联
反投影的能量聚束直接体现为聚束能量的时空分布. 如果直接将聚束能量场的时空函数转化为具有时空分布的先验概率分布,并通过MCMC 方法进行多点源反演. 可以通过拟合实际波形获得电源位置以及发震时间的后验概率. 该思路的难点体现在两个方面:(a)如何构建聚束能量到先验概率的数学表达. (b)如何通过有先验概率的采样以及震源机制的采样提高反演效率. 如果不对MCMC 的参数采样做技术改进,MCMC 应该会直接搜索全部参数空间,包括每个点源的走向、倾角、滑动角以及地震矩. 而如果震源时间以及位置的采样已经通过先验概率给出,则基于该时空采样的震源机制可以通过线性算法获得. 这种结合采样方法和线性方法的联合策略有望提升MCMC 的反演效率.
类似将反投影结果应用于多点源反演的思路,张喆和许力生(2020) 将其应用于2013 年南斯科舍海岭M
7.8 的地震中,并且取得了较好的反演结果. 而形成自动化的算法以及具有严格概率含义的先验概率转化方法,尚需要进行更多的探讨. 同时,反投影成像技术应用不同的台网可能获得不同的时空分布. 这与三维介质下远震体波的射线矫正相关. 通过三维射线矫正可以较好地解决该问题(Liu et al., 2018). 如果反投影体现的高频辐射与多点源的低频辐射存在差别,也可以通过增加先验概率函数时空范围的方法获得二者兼顾的反演效果. 同时多点源反演的不确定性往往高于单点源反演,如GCMT 的角度误差约15°. 如果用MCMC 方法实现多点源反演,往往可以对其不确定性给出较好的估计. 在后续工作中,结合地质构造信息,有可能进一步缩小震源机制的不确定性.实现反投影结果与多点源数据的串行方法,不仅为我们稳定地解析大地震破裂过程提供了策略,也为震后快速响应提供了可行性方案. 基于远震体波的反投影,多点源和破裂过程反演对于数据的依赖较弱. 通常可以在震后几十分钟内获得. 三种方法的计算效率都较高,在参数稳定的情况下,完成一次串行反演通常可以在1 小时之内完成. 这为地震破裂过程速报以及救援快速部署提供了重要依据.在震源破裂过程快速反演之后,通过数值模拟的方法进一步计算震源区的强地面运动,也是一套可能为震后救灾提供重要依据的模型驱动思路. 数值模拟方法一般基于有限断层模型,对于破裂过程以及震源机制准确度要求较高. 因此通过多点源方法获得更加准确的断层面,以及在多点源机制约束下的有限断层模型也可以进一步增加用数值模拟计算震源区强地面运动的准确性.