(1)
摘要 利用远震 P 波观测数据, 研究 2012 年印度洋 Mw 7.2 地震是否近场触发海沟附近余震的问题。根据区域化水波直接触发机制和间接激发机制的假设, 提出线性断层反演模型和双点源反演模型。远震 P 波数据反演结果表明, 基于区域化水波直接触发机制的线性断层模型对波形的拟合效果更好, 并且 2012 年印度洋地震没有触发海沟附近的余震。
关键词 动态触发; 双点源反演; 线性断层反演; 反投影
在全球的活跃断层带, 广泛地观测到地震的动态触发现象。地震的触发机制大致可以分为静态应力触发和动态应力触发两种[1‒3]。静态应力触发指地震发生后, 位错区域的永久性形变对该地区发生地震事件的影响。静态应力随着距离的增加而快速衰减, 因此在大地震破裂区域或临近区域, 静态应力作用更明显。当与主震破裂区域有一定距离时, 静态应力的作用几乎可以忽略。动态应力触发指大地震激发的地震波传播到某些处于临界压力状态的区域而触发地震事件的现象。远震的地震波传播会引起动态应力变化, 此时静态应力又可以忽略不计, 因此远震的动态应力触发现象比较容易判断。对于中近地震, 区分动态应力触发与静态应力触发机制相对困难, 两种机制都有触发区域地震的能力[4‒5]。
反投影(back projection)技术通常用于确定大地震破裂过程中的高频特征[6‒8], 因其利用远震波的相干性且不需要设置断层参数, 反投影方法在解析复杂破裂过程时表现出高效。反投影方法解析得到的破裂点对应当前时窗内能量释放最大值所在位置,并不总是与真实破裂位置对应, 因此反投影得到的破裂图像中经常存在假象。如果此类辐射点源位于远离主震破裂区域, 则容易将其错误地识别为主震动态触发的余震[9‒10]。2016 年, Fan 等[9‒10]用反投影技术解析 2012 年苏门答腊岛附近 Mw7.2 级地震的破裂过程, 发现在距离主震区域约 150km 的海沟附近, 反投影图像中存在与主震破裂区域在空间上分割开的能量辐射点源, 并且紧随主震发生之后, 他们将其解释为 Mw7.2 主震的动态触发事件。2017 年, Yue 等[11]用波形分析和经验格林函数(empirical Green’s function, EGF)校准等方法, 重新解析 2012年 Mw7.2 级地震, 选取主震区域附近发生的一个Mw 6.2 地震作为校准事件, 对比该事件和主震的反投影图像, 发现在海沟附近出现与主震类似的辐射点源(校准事件与主震共享传播路径), 因此认为海沟附近的能量更有可能来源于海沟附近的特殊地形, 同时发现地震尾波中海沟附近辐射点源的共振频率接近局部混响水波的特征频率; 他们还利用该区域的真实地形进行二维数值波形模拟, 证明地震波会被海沟处特殊地形放大, 形成区域化混响水波, 并且该区域水波的位置、持续时间和间隔与反投影图像中辐射点源的特征一致。
1996 年, Ihmlé 等[12]最早提出大地震激发区域化水波的机制, 以此解释大地震主震 P 波之后的长时间尾波。近期, 此激发机制得到其他研究工作的验证, 地震尾波中的区域化水波可用于确定近海大地震中是否发生近海沟破裂[13‒17]。然而, 有关板内地震对区域化水波的激发机制还存在争议。例如, Fan 等[18]认为 2012 年 Mw 7.2 主震反投影图像中的辐射点源是由局部混响水波产生, 但他们提出一种“间接激发”机制, 即主震动态地触发了海沟附近的余震, 余震激发了海沟附近的区域化混响水波, 这种机制更容易解释主震尾波的振幅, 但其到时间隔与主震尾波不相符, 因此引入间接激发机制的必要性需要验证。Yue 等[11]的研究结果表明, 主震地震波的能量直接触发海沟附近的区域化水波, 即“直接激发机制”。
本文根据两种激发机制的条件, 设计双点源反演(间接激发机制)和线性源反演(直接激发机制)两种波形反演方法。根据两种方法波形拟合的结果, 可以判断直接触发机制和间接触发机制的合理性, 从而证明在 2012 年 Mw7.2 主震期间, 是否存在近海沟的动态触发现象。
本文选取两个事件( 2007‒10‒04 Mw 6.2 和 2012‒ 04‒20Mw5.8), 作为主震区域和余震区域的经验格林函数(分别记为 EGF1 和 EGF2)进行双点源反演。主震震中和经验格林函数震中如图 1 所示, 震源参数见表 1。从美国地震学研究联合会(Incorporated Research Institutions for Seismology, IRIS)网站下载同时记录主震、EGF1 和 EGF2 的远震波形数据, 包括 30 个台站(位置和方位角范围见图 2)。对主震、EFG1 和 EGF2 的波形进行带通滤波, 滤波频带均为0.02～0.5 Hz。
将主震和海沟附近强余震的震源时间函数分别记为 S1(t)和 S2(t), 经验格林函数分别记为 F1(t)和F2(t)。利用三角子波参数化主震和海沟强余震的震源事件函数, 构建反演矩阵进行线性反演。其中, 主震和强余震三角子波分别为 20 个和 40 个, 三角子波的半持续时间分别为 2 s 和 1 s (表 2)。
对于间接触发机制, 海沟强余震是由主震触发的, 因此其到时存在推迟时间。主震是走滑型地震, 因此我们假设海沟余震是由主震S波触发的, 则推迟时间可表示为
(1)
其中, d 为主震震中到 EGF2 的距离, θ 为 EGF2 震中相对于主震震中的方位角, azin 为主震震中相对于第 n 个台站的方位角, par 为研究区域的射线参数, Vs为剪切波波速。
设第 n 个台站接收到的波形为 Un(t):
白色圆点标记国际地震中心(International Seismological Centre, ISC)目录(http://www.isc.ac.uk)中2012‒01‒10 Mw 7.2主震发生24小时后余震的震中位置, 小圆点、中圆点、大圆点代表震级分别为3, 4, 5; 海沟处灰色圆点标记ISC目录中2005—2012年Mw>5地震的震中位置。左上角插图显示地质构造背景, 红色方框内为研究区域
图1 研究区域地质背景和地震分布
Fig. 1 Geological background and seismic distribution in study area
表1 主震和经验格林函数的震源参数
Table 1 Mainshock and EGFs parameters
事件时间USGS目录给出的震中GCMT目录给出的震中MwGCMT目录给出的震源机制解 主震2012‒01‒10 18:37:132.543/92.903/202.59/92.98/247.2 103/81/−173 EGF12007‒10‒04 12:40:302.433/93.210/192.47/92.83/126.2110/63/170 EGF22012‒04‒20 22:19:523.230/93.760/213.28/93.65/125.8167/43/113
说明: 数据来自全球地震矩张量目录(global centroid moment tensor, GCMT)(https://www.globalcmt.org/CMTsearch.html)和美国地震调查局(United States Geological Survey, USGS)的地震目录(https://earthquake.usgs.gov/earthquakes/search/); “震中”的数据格式为“纬度/经度/深度”, “震源机制解”的数据格式为“走向/倾角/滑动角”。
三角形代表地震台站, 五角星代表主震震中位置
图2 台站分布
Fig. 2 Distribution of seismic stations
表2 线性断层模型和双点源模型的震源反演参数
Table 2 Inversion parameters of the line-source and two-point source models
模型长度子断层个数子断层三角子波个数三角子波的半持续时间/s 线性断层模型1602062 模型经验格林函数位置子断层个数子断层三角子波个数三角子波的半持续时间/s 双点源模型近海1202 近海沟1401
(2)
其中, 
代表卷积运算。接下来, 我们利用不同台站的经验格林函数构建反演矩阵, 联合观测数据进行线性反演, 得到双点源模型中主震和海沟强余震的震源时间函数。
Fan 等[18]指出, 校准事件 EGF1 的波形不能很好地解释主震尾波, 需要额外添加近海沟的动态触发事件, 才能更好地拟合波形。我们知道, 时空的有限性导致主震的波形与 EGF 波形之间的主要差异。主震破裂的方向性和持续时间会影响东北方向传播的地震波, 并且这些波会激发水波震相[11]。由于单一的校准事件不能完全解释主震激发水波震相的振幅, 因此考虑震源的空间有限性非常必要。
根据“直接触发”机制, 主震的波形直接在海沟附近触发局部混响水波。完整格林函数的计算需要考虑海沟地形的散射效应, 进行 3D 数值模拟[15,19]。本文依靠一种经验格林函数校正技术, 利用 2007年的 EGF1 事件, 为每个子断层构造格林函数。主震和 EGF1 的远震观测能量来源于破裂的 P 波能量和海沟附近局部混响水波能量两个方面。该算法包括两个步骤。
1) 计算 EGF1 波形中 P 波与局部混响水波的到时差, 按照到时差将 EGF1 的波形分为 F1 和 F2 两部分。F1 代表与 EGF 事件破裂有关的波形, 包括直达 P 波、pP 波、sP 波和近源混响水波震相, 将其作为主震破裂的经验格林函数; F2 代表与近海沟局部混响水波有关的相位, 将其作为海沟局部混响水波震相的经验格林函数。同时, 沿着主震断层的走向, 设置东南‒西北走向的线性断层。基于点源假设, 可以计算出线性断层中每个子断层的 P 波与局部混响水波的到时差, 利用 F1 和 F2 构建该子断层的格林函数。
2)考虑到主震震中与每个子断层的位置差异, 对子断层的格林函数进行 P 波走时矫正。
使用与双点源反演相同的台站波形数据及滤波频段进行线性断层反演, 设置子断层的个数为 20, 每个子断层的震源函数由 6 个半持续时间为 2s 的三角子波进行参数化(见表 2)。构建好反演矩阵后, 利用非负最小二乘算法, 反演得到每个子断层上的震源时间函数, 将不同的到时累加, 得到整个地震的震源时间函数。
如图 3 所示, 用水波相对到时 tdly0 将EGF1 分为 F1 和 F2 两部分, 其中 F2 主要由横向传播的 P 和S 波激发。2012 年 Mw7.2 地震属于走滑型地震, S波对局部混响水波的激发效率更高, 因此我们使用S 波到海沟散射点的传播时间来确定局部混响水波的开始时间。按照 P 波到时, 将每个台站的 EGF 波形对齐, 同时将所有台站的 P 波初动时间设置为 0 s, 则局部混响水波的相对到时为
(3)
对于远震 P 波的有限断层反演, 子断层格林函数的构建分为两个步骤: 1)将格林函数的 P 波初动时间设置为 0 s; 2)利用震源区域的射线参数, 计算震源与子断层之间的相对到时, 并将相关的时移应用零时移的格林函数去构造子断层的格林函数矩阵, 进行有限断层反演。
本文使用类似的方法构造线性断层反演的经验格林函数矩阵。第一步, 首先将子断层的经验格林函数进行零时移校正, 然后根据子断层与主震震中的相对位置, 对子断层的格林函数进行相对走时校正。第二步与理论格林函数的计算相同, 因此这里主要介绍每个子断层零时移的经验格林函数的构造。由于 EGF1 将各个台站按照 P 波到时进行对齐, 因此每个台站的 F1 波形也是按照 P 波到时对齐的, 我们只需要根据 EGF1 与每个子断层之间的位置差异, 对 F2 进行相对到时校正。
EGF1 波形中局部混响水波的延迟时间由式(3)给出。可以使用相同的方法, 计算目标子断层处局部混响水波的到时:
(4)
其中, d1 表示散射点到子断层震中的距离; θ1 表示子断层震中相对于散射点的方位角。然后, 利用tdly0和tdly1, 可以计算目标子断层的局部混响水波相对于EGF1的相对时移:
(5)
其中, Δτ表示目标子断层相对于EGF1局部混响水波的走时校正量, d′表示EGF1震中到目标子断层震中的距离; θ′表示EGF1震中相对于目标子断层震中的方位角。
从式(5)可以看出, 每个子断层的局部混响水波走时校正包含两项: 第一项表示从震源到水波散射点激发阶段的传播时间, 第二项表示EGF1与目标子断层的相对位置差异。图4示意走时校正的原理。
F2波形的振幅可以通过局部混响水波的传播特征进行校正。局部混响水波是被上行的体波激发后在水层混响散射增强的, 其传播特征类似二维传播, 在波前圆周上的动能积分是守恒的(即A2∝d −1, 其中A表示局部混响水波的振幅, d表示波前圆周的半径)。因此, 在联合F1构造经验格林函数之前, 需要对F2的振幅进行校正。振幅校准系数Acor的表
tdly0 表示EGF1的P波与局部混响水波的到时差, d0表示从散射点到EGF1震中的距离, azi表示台站相对于主震震中的方位角, θ0 表示EGF1震中相对于海沟散射点的方位角, par表示对应台站的射线参数
图3 P波与局部混响水波到时差计算示意图
Fig. 3 Schematic diagram of calculating the arrival time difference between P wave and localized water reverberation phases
图4 子断层与EGF1局部混响水波到时校正示意图
Fig. 4 Schematic diagram of calculating arrival time difference of localized water reverberation phase between the subfault and EGF1
达式如下:
(6)
针对目标子断层格林函数GF通过下式给出:
(7)
本文通过双点源模型和线性断层模型反演得到震源时间函数(图 5)。双点源模型得到的主震矩震级为 7.0, 海沟强余震的矩震级为 6.8; 线性模型反演得到的主震矩震级为 7.14, 更接近 GCMT 给出的结果(Mw=7.2); 线性断层模型和双点源模型波形拟合的方差缩减幅度分别为 18%和 11%。双点源模型和线性源模型使用相同数量的自由参数, 可见线性源模型的波形拟合结果优于双点源模型。从波形拟合图像(图 6)中也可以得到相同的结果。
阴影区为线性断层模型反演得到的震源时间函数, 黑色折线为双点源模型反演得到的主震附近的震源时间函数, 灰色折线为海沟散射点处的震源时间函数
图5 双点源模型和线性断层模型反演得到的震源时间函数
Fig. 5 Source time function extracted by line-source model and two-point model
综上所述, 基于局部混响水波的直接激发机制的线性断层模型在地震矩恢复和模型拟合两个方面的表现均优于基于间接激发机制的双点源模型。但是, 线性断层模型的方差缩减幅度仍然小于一般的有限断层反演。影响线性断层模型波形拟合结果的因素主要有 3 个: 1)在构造子断层的经验格林函数时, 我们没有考虑子断层深度震相的影响; 2)本文设置的反演参数为 120 个, 少于一般的有限断层反演; 3)EGF1 与主震的震源机制解存在差异。由于EGF1 的信噪比较低, 线性断层模型的波形拟合结果还可以通过三维理论格林函数去进一步提高。线性断层模型也给出滑移量的分布, 从图 7 可以看出, 线性断层模型给出的主震是沿西北方向破裂的, 与反投影图像给出的方向[6,11]一致。
Fan 等[18]提出海沟附近余震的间接触发机制假说, 但未给出触发余震的震级。本文通过经验格林函数, 依据间接触发机制假说设置反演模型, 得到的海沟余震为 Mw 6.8。该震级的地震将在海沟附近产生明显的余震序列, 最大震级在 Mw 5.0 以上。我们在另一项地震识别工作中发现, 主震过后24 小时内海沟附近没有余震发生[20]。我们此前利用区域台网对海沟附近的地震进行扫描, 获得的最小震级约为 Mw 2.7, 远远低于 Mw 6.8 海沟地震应有的余震震级[20]。因此, 用 Fan 等[18]的间接触发机制获得的海沟附近地震, 无法用海沟附近没有余震分布解释, 可以认为其间接触发机制与现有观测结果相悖。
每条波形左上方标注台站名, 左下方标注相对于主震震中的方位角; 黑色曲线为台站观测波形, 红色曲线为线性断层模型的拟合波形, 蓝色曲线为双点源模型的拟合波形
图6 15个远震台站的拟合波形
Fig. 6 Waveform fits of 15 seismic stations
图7 线性断层模型的滑移量分布
Fig. 7 Slip distribution of line source model
此外, 对于间接触发机制, 本文也验证了线性断层加海沟点源的反演模型。选取 4 组不同的参数组合去拟合波形, 参数的数量和反演结果列于表3。参数组合“120+60”即是在原有线性断层模型参数的基础上, 额外增添 60 个海沟点源的反演参数, 波形拟合结果(方差缩减幅度为 19%)相对于线性断层模型(方差缩减幅度为 18%)无明显提高。将线性断层参数的数量固定为 60 个, 改变海沟点源参数的数量, 发现其波形拟合结果均差于单一线性断层的拟合结果。
表3 线性断层+海沟点源的震源反演参数
Table 3 Inversion parameters of the line-source plus trench-point-source models
次数子断层个数子断层三角子波个数海沟经验格林函数三角子波个数参数个数方差缩减幅度/% 12066018019.48 220330 9016.63 32034010016.71 42036012016.78
本文基于直接触发机制和间接触发机制的假设条件, 构建双点源反演和线性断层反演模型。反演结果表明, 线性断层模型可以更好地拟合观测波形和恢复主震震级; 线性断层模型还能给出主震的破裂方向, 与反投影图像给出的破裂方向一致。因此, 我们认为在 2012 年苏门答腊岛附近 Mw 7.2 地震发生期间, 不存在海沟附近近实时的动态触发现象。
参考文献
[1] Kilb D, Gomberg J, Bodin P. Triggering of earthquake aftershocks by dynamic stresses. Nature, 2000, 408: 570‒574
[2] Belardinelli M E, Bizzarri A, Cocco M. Earthquake triggering by static and dynamic stress changes. Journal of Geophysical Research, 2003, 108: 2135
[3] Freed A M. Earthquake triggering by static, dynamic, and postseismic stress transfer. Annual Review of Earth and Planetary Sciences, 2005, 33(1): 335‒367
[4] Gomberg J, Felzer K, Brodsky E. Earthquake dyna-mic triggering and ground motion scaling // Proc 4th Int Workshop on Statistical Seismology. Kanagawa, 2006: 45–51
[5] Richards-Dinger K, Stein R S, Toda S. Decay of after-shock density with distance does not indicate trigge-ring by dynamic stress. Nature, 2010, 467: 583‒586
[6] Ishii M, Shearer P M, Houston H, et al. Extent, du-ration and speed of the 2004 Sumatra-Andaman earth-quake imaged by the Hi-Net array. Nature, 2005, 435: 933‒936
[7] Xu Y, Koper K D, Sufri O, et al. Rupture imaging of the Mw 7.9 12 May 2008 Wenchuan earthquake from back projection of teleseismic P waves. Geochemistry, Geophysics, Geosystems, 2009, 10(4): Q04006
[8] Meng L, Inbal A, Ampuero J P. A window into the complexity of the dynamic rupture of the 2011 Mw 9 Tohoku-Oki earthquake. Geophysical Research Let-ters, 2011, 38(7): L00G07
[9] Fan W, Shearer P M. Fault interactions and triggering during the 10 January 2012 Mw 7.2 Sumatra earth-quake. Geophysical Research Letters, 2016, 43(5): 1934‒1942
[10] Fan W, Shearer P M. Local near instantaneously dyna-mically triggered aftershocks of large earthquakes. Science, 2016, 353: 1133‒1136
[11] Yue H, Castellanos J C, Yu C, et al. Localized water reverberation phases and its impact on back projec-tion images. Geophysical Research Letters, 2017, 44 (19): 9573– 9580
[12] Ihmlé P F, Madariaga R. Monochromatic body waves excited by great subduction zone earthquakes. Geo-physical Research Letters, 1996, 23(21): 2999‒3002
[13] An C, Yue H, Sun J, et al. The 2015 Illapel Chile (Mw=8.3) earthquake: direction-reversed along-dip rupture with localized water reverberation. Bull Seis-mol Soc Am, 2017, 107(5): doi: 10.1785/0120160393
[14] Lay T, Rhode A. Evaluating the updip extent of large megathrust ruptures using P coda levels. Geophysical Research Letters, 2019, 46(10): 5198‒5206
[15] Qian Y, Wei S, Wu W, et al. Teleseismic waveform complexities caused by near trench structures and their impacts on earthquake source study: application to the 2015 Illapel aftershocks (Central Chile). Jour-nal of Geophysical Research: Solid Earth, 2019, 124 (1): 870‒889
[16] Lay T, Liu C, Kanamori H. Enhancing tsunami war-ning using P wave coda. Journal of Geophysical Re-search: Solid Earth, 2019, 124: 10583‒10609
[17] Wu Z, Lay T, Ye L. Shallow megathrust slip during large earthquakes that have high P coda levels. Jour-nal of Geophysical Research: Solid Earth, 2020, 124: e2019JB018709
[18] Fan W, Shearer P M. Coherent seismic arrivals in the p wave coda of the 2012 Mw 7.2 Sumatra earthquake: water reverberations or an early aftershock?. Journal of Geophysical Research: Solid Earth, 2018, 123(4): 3147–3159
[19] Fernando B, Leng K, Nissen-Meyer T. Oceanic high-frequency global seismic wave propagation with rea-listic bathymetry. Geophysical Journal International, 2020, 222(2): 1178‒1194
[20] Jia B, Yue H, Muzli M. Absence of near-trench early triggering during the 2012 Mw 7.2 Indian Ocean strike-slip earthquake: evidence from one-day after-shocks. Seismological Research Letters, 2021, 93 (2A): 719–726
Discussion about Near Trench Triggering during 2012 Mw 7.2 Indian-Ocean Strike-Slip Earthquake
Abstract Whether the aftershocks near the trench were triggered by 2012 Indian Ocean Mw 7.2 earthquake was discussed by using teleseismic P-wave records. According to the assumption of the direct and indirect trigger mechanism of regional water waves, the line-source model and the two-point source model were proposed. The inversion results show that the line-source model based on the direct trigger mechanism was better than the two-point model on the waveform fitting, and the aftershocks near the trench were not triggered during 2012 Indian Ocean earthquake.
Key words dynamic trigger; two-point source inversion; line-source inversion; back-projection
doi: 10.13209/j.0479-8023.2022.023
国家自然科学基金(42021003)资助
收稿日期: 2021-04-07;
修回日期: 2021-05-13