2022, Vol. 42
2. 成都理工大学地质灾害防治与地质环境保护国家重点实验室，成都市二仙桥东三路1号，610059;
3. 西南石油大学土木工程与测绘学院，成都市新都大道8号，610500;
4. 西南交通大学地球科学与环境工程学院，成都市犀安路999号，611756
2019-10-11菲律宾棉兰老岛南部区域发生MW5.1地震。此后2个月，在50 km范围内共发生4次MW＞6.0强震，分别为2019-10-16 MW6.4地震、2019-10-29 MW6.6地震、2019-10-31 MW6.5地震以及2019-12-15 MW6.8地震，下文将4次地震按发震时间分别简写为地震1、地震2、地震3、地震4。此次地震序列均发生于菲律宾棉兰老岛，该区域地质构造较为复杂，位于巽他板块、菲律宾海板块、太平洋板块和澳大利亚板块的大规模俯冲交会边界带上[1]。复杂的板块构造关系导致整个菲律宾群岛地震频发，在最近的50 a里，距离此次地震序列250 km范围内共发生34次6.5级以上强震。
科研人员对此次地震序列的同震断层破裂模型展开研究，并讨论相关地震间的关联性[2-3]。但当前研究成果中，针对地震1的模型正演InSAR形变场与原始观测值之间存在不可忽视的形态分布差异[2-3]，且USGS、GCMT等研究机构给出的震源机制解相差较大。此外，在同震断层平面维度内，研究先发地震对后发地震的定量影响，特别是先发序列地震对后续地震震源位置与显著滑动区库仑应力传输与累积的影响，有助于揭示地震序列之间的触发关系，而此次地震序列缺少这方面的系统研究。
为此，本文搜集ALOS-2和Sentinel-1卫星升降轨雷达影像数据，利用InSAR技术提取此次地震序列不同事件的精细同震地表形变场，并在弹性半空间位错理论框架下，利用地表形变数据反演获得各地震事件的同震断层破裂模型，确定4次地震事件的发震断层，揭示其运动模式、空间分布和滑动量级。并进一步利用反演同震断层运动模型计算不同地壳深度和同震断层面上的库仑应力变化，分析讨论此次地震序列之间以及主震对余震的触发关系。
1 InSAR数据收集与处理收集ALOS-2卫星和Sentinel-1卫星的震前、震后L波段和C波段SAR影像，并对其进行干涉处理，最终得到4次地震的InSAR同震形变场(图 1)。由于地震2和地震3发震时间过于接近，两次地震造成的地表形变无法直接分离提取(图 1(b)、(c))，后续将利用大地测量数据对这两次地震的断层运动模型进行统一解译。
|
红色五角星为地震震中位置, 黑色线段为推测的孕震断层位置 图 1 4次地震的InSAR同震形变场 Fig. 1 The coseismic InSAR deformation filed of the four earthquakes |
已有研究成果显示[2-3]，目前地震1发震断层仍存在较大争议。由图 1的同震InSAR形变场可知，AB和CD所代表的断层均有可能是此次地震的发震断层，而USGS给出的双节面解显示2个可能的发震断层走向角分别为42°和132°。若断层走向为132°，则位于断层CD西南一侧的断层上盘形变量应显著大于东北一侧的断层下盘形变量，而这与InSAR观测结果明显不符。因此，初步推断走向角为42°的这组解与实际情况更为接近。同时，其他已有研究以约300°的走向角构建了发震断层模型[2-3]。本文将在断层运动模型反演部分对比分析上述2个模型的合理性。
图 1(b)、(c)中同时包含地震2和地震3的同震形变场，根据其空间分布可知，形变场整体沿东北-西南方向分布，因此，初步判定2次地震孕震断层可能均为东北-西南走向(图 1(b)、(c)中黑色线段)。结合USGS给出的震源机制解可知，这2次地震发震断层不同的走向角意味着这2次地震事件的发震断层具有不同的倾向，即两者为相互独立的断层。由图 1(d)、(e)可知，地震4断层走向为西北-东南(图 1(d)、(e)黑色线段)，且向东北方向倾斜。
3 同震断层运动模型反演与分析基于4次地震的InSAR同震形变场开展断层运动模型的反演与分析。首先利用四叉树算法对原始InSAR形变场进行降采样处理，获得数量较少且质量可靠的InSAR形变数据，用于同震断层运动模型的约束反演。然后根据前文InSAR形变场与USGS震源机制解耦合分析确定的最佳断层面解参数信息(表 1)，构建初始的断层几何模型，其中断层走向和倾角分别设定为USGS震源机制解给定的初值。
|
表 1 4次地震断层参数 Tab. 1 Fault parameters of the four earthquakes |
在上述断层初始几何模型基础上，利用Yang等[4-5]提出的基于大地测量的断层运动模型反演算法，获取最优的断层几何参数(表 1)。在上述参数基础上，进一步反演获得精细的同震断层滑动分布，结果如图 2~4所示(五角星为震源位置，黑色箭头为断层滑动的方向和量级，黑色实线框为显著的同震滑动区)。为验证滑动反演结果的准确性，根据上述反演得到的断层运动模型，正演获得模拟的同震InSAR形变场，并根据模拟数据与观测数据的残差分布(图 5)分析反演断层模型的可靠性。
|
图 2 地震1同震断层滑动模型 Fig. 2 Coseismic fault slip model of earthquake 1 |
|
图 3 地震2、地震3同震断层滑动模型 Fig. 3 Coseismic fault slip model of earthquake 2 and earthquake 3 |
|
图 4 地震4同震断层滑动模型 Fig. 4 Coseismic fault slip model of earthquake 4 |
|
图 5 模拟形变场及其对应残差 Fig. 5 Simulation deformation field and its corresponding residual |
对于地震1，正演结果(图 5(a))与原始InSAR观测数据(图 1(a))不管是在形变分布还是量级上均具有较高的一致性，同时，断层近场也未发现明显的残差分布(图 5(b))，仅在远离断层的西北区域存在一个相对显著的残差区，但并未发生明显的地表变形。造成上述残差的主要原因是该区域被密集的植被覆盖，导致InSAR干涉失相干，因此，原始InSAR观测数据在此区域并不可靠。此外，由表 1可知本文与USGS的估计震级规模基本一致，但本文反演获得的断层几何参数与USGS的反演结果存在一定差异。这主要是由于USGS是利用远场地震波数据开展震源机制参数估计，但受限于精细地壳模型和近场观测台站缺失，反演结果往往具有较大的不确定性[6]；而本文基于InSAR技术，在断层近场和远场均提取到了可靠的同震形变数据，可为同震断层运动模型的反演提供较好的约束。因此相比于USGS，本文的反演成果具有更高的可靠性。
此外，为讨论约300°走向角发震断层模型的可靠性，参考Li等[2]的研究成果，增加304.5°断层走向的反演(图 2(c)、(d))，并利用反演断层运动模型正演计算该模型对应的形变场(图 5(c))及其与InSAR观测值之间的残差(图 5(d))。通过与本文计算的正演形变场(图 5(a)、(b))对比可知，虽然二者在形变量级上均与InSAR观测结果有较高的一致性，但本文结果在空间分布上与InSAR观测结果的一致性更高且残差更小；图 5(d)中还存在3处显著的残差区(图中椭圆区域)。因此本文结果具有更高的可靠性。
进一步分析图 2(a)所示的地震1发震断层模型可知，在此次地震震源区域(图 2(a)中黑色矩形所示A区域)存在一个运动量级不大的滑动密集区，其平均滑动量级约为0.15 m，呈近纯右旋走滑破裂模式。显著滑动区域(B区域)断层运动以右旋走滑为主，兼具正断层破裂，最大滑动量约为0.5 m，发生在震源东北向约20 km、地下约15 km处。B区域的断层运动释放了此次地震超过95%的能量，并在同震区域造成cm级的地表形变(图 1(a))。
3.2 地震2、地震3的同震断层运动模型分析地震2、地震3的断层运动模型正演地表形变(图 5(e)、(g))与原始InSAR观测值(图 1(b)、(c))间具有极好的一致性，升降轨均仅在断层近场区域发现少量形变残余信号。这主要是因为断层近场存在剧烈地表形变，易造成InSAR干涉失相干，进而导致InSAR观测值与正演结果的不一致，而远场存在的量级不显著的残余信号则主要由大气延迟分量贡献。相较于Li等[2]对这2次地震的联合反演，Sentinel-1升轨SAR数据干涉结果质量较差，并且受大气及干涉失相干噪声影响较为显著，导致其结果残差较大。综上所述，本文联合ALOS-2及Sentinel-1数据的正演形变结果与观测数据间的高度一致性说明反演断层运动模型具有较高的可信度。
由图 3(a)、(c)可知，地震2事件断层破裂主要包括A、B、C三个显著滑动区域。其中区域A滑动量最大，最大滑动量约为2.0 m, 并以右旋走滑运动为主；区域B最大滑动量约为1.1 m，滑动形式同样以右旋走滑为主；区域C滑动量级相对较小，平均滑动量约为0.3 m，但滑动模式以正断层破裂为主。
地震3反演破裂模型如图 3(b)、(d)所示，相比于地震2，地震3的断层滑动分布更为集中，仅在震源正上方区域存在单一的滑动中心D。D区域滑动以右旋走滑为主，最大滑动量约为1.1 m，在0~8.0 km深度范围内断层运动除右旋走滑外还伴随有轻微的逆冲运动分量，而在8.0 km以下深度则转变为轻微的正断层破裂。此外，观察图 5所示的各孕震断层空间分布可知，F1(地震1孕震断层)和F3(地震3孕震断层)断层的显著滑动区在空间上相邻，表 1中的断层几何参数也显示F1和F3具有相似的断层走向和倾角，同时2次地震的破裂机制也基本一致(图 2(a)、图 3(c))。因此，综合推断地震1与地震3应该是发生在同一孕震断层相邻的2个区段上。
3.3 地震4的同震断层运动模型分析图 4所示的断层滑动模型显示，地震4为典型的走滑地震，与USGS给定的震源机制解结论一致。显著的断层滑动主要分布在地下2.8~11.3 km处，断层破裂以左旋走滑为主并兼具微量的逆冲和正断层运动，其中逆冲分量主要集中在震源北侧以及西北侧区域。通过InSAR升降轨干涉结果(图 1(d)、(e))可以看出，在震源位置北侧以及西北侧区域，升降轨结果均呈现LOS向缩短信号，因此可以推断出该区域地表形变以抬升为主，即断层在此区域除了整体的走滑运动外还兼具逆冲运动。此外，断层滑动模型中显示出的微量正断层运动在滑动模型中仅占有2~3个格网区域，且量级较小，整体以走滑为主。造成这种现象的原因可能是此次地震震级较大，且伴有多次余震事件，极易受到地震震后效应的影响，同时由于此次地震发震断层的高倾角特性，受自身重力影响，易出现正断层运动。上述地震断层破裂以走滑运动为主，同时兼具一定分量的逆冲和正断层破裂，此现象与该区域2017年MW6.5地震[4]类似。
地震4的断层最大滑动量约为3.0 m，显著大于之前的3次地震事件，这也造成断层近场存在剧烈的地表形变，并导致InSAR观测的干涉失相干现象(图 1(d)、(e))。断层的滑动主要分布在地表 3.0 km以下区域，而在接近地表 0~3.0 km深度的滑动量逐渐衰减为0，所以此次地震并未造成明显的断层地表出露。
对比图 5(i)、(k)与图 1(d)、(e)可知，升降轨模拟InSAR形变与原始观测在形变分布、量级和趋势上均保持较高的一致性，这说明反演同震断层破裂模型具有较高的可靠性。但需要注意的是，在图 5(j)、(l)所示的残差图中，升降轨断层近场均存在一定的残余形变分量，这主要是由于断层近场区域被密集植被覆盖，InSAR相干性不易维持所致；同时，约3.0 m的同震断层破裂在断层近场也造成显著的地表形变，区域形变梯度可能会超过C波段的监测能力，导致该区域InSAR干涉成果可靠性较差。此外，强震后的余震序列也是造成上述形变不一致的因素之一。
4 棉兰老岛地震序列库仑应力传输与分析图 6为地震序列后余震密集的20 km深度处的库仑应力变化[7]，选择与主震事件一致的接收断层参数(表 1)，摩擦系数为0.4[8]。观察图 6(a)、(b)可知，地震1、地震2的空间分布与地震造成的正应力变化在空间上有较好的一致性，这意味着此次地震后的余震可能与主震有较为一致的破裂机制。然而，图 6(c)、(d)中余震与正、负应力变化的对应关系差异不大，表明这些余震除了受到同震破裂事件影响外，区域构造背景应力、区域断裂和褶皱发育分布以及动态库仑应力传输等其他因素也对余震事件的发生造成了不可忽视的影响[8]。
|
图 6 地震事件在20 km深度处造成的库仑应力变化 Fig. 6 The Coulomb failure stress change at the depth of 20 km caused by the earthquakes |
为进一步定量分析此次地震序列4次强震间的触发关系，利用上述反演同震断层破裂模型计算每次强震对其后序列地震孕震断层的静态库仑应力传输[9-10]，其中接收断层参数根据前文反演结果确定，具体参见图 7。图 7(g)、(h)、(i)为先前地震事件造成的孕震断层面库仑应力变化的总和。
|
图 7 地震序列间库仑应力传输 Fig. 7 The Coulomb failure stress transfer between the earthquakes sequence |
由图 7可见，地震1对其后3次强震事件的显著同震破裂区域(A、D、E)均造成明显的正向库仑应力传输，即促进了这些区域的破裂。地震2的显著滑动区A(该区域释放了此次地震超过90%的能量)库仑应力平均增加约0.3 bar，显著超过0.1 bar的最小触发量级。但需要说明的是，地震1在地震2的震源区域造成的库仑应力变化较为微弱，且对B、C显著滑动区的应力传输为负值。综上所述，地震2的初始破裂为区域构造应力持续累积的结果，此次地震事件最显著滑动区的破裂受地震1的正向触发影响显著。
由图 7(b)可见，地震1对地震3的震源区域存在负的库仑应力传输，量级约为-0.34 bar；但地震1对地震3的显著滑动区B超过70%的区域存在正向的库仑应力传输，平均量级约为0.3 bar。同时，由图 7(d)可见，地震2显著增加了地震3震源区域的库仑应力，量级约为1.0 bar；但地震2对显著滑动区D的浅部约50%的区域造成不可忽视的负值库仑应力传输，平均量级约为-4.7 bar，抵消了先前地震1对该滑动区域的正向库仑应力传输(图 7(h))。综上可知，地震2显著触发了2 d后地震3的初始破裂，之后在背景构造应力和动态库仑应力传输等因素的综合作用下，造成地震3孕震断层面的破裂。
地震4受到前3次地震事件的综合影响，图 7(c)、(e)和(f)显示地震1对地震4震源区域和显著滑动区均造成正向的库仑应力传输；地震2未对地震4的震源区域和显著滑动区造成影响；地震3对地震4的震源和显著滑动区造成负的库仑应力传输，平均量级约为-0.35 bar。综上可知，前3次地震对地震4的库仑应力传输总和为负值，即静态库仑应力传输可能并不是直接触发地震4的主要诱因。
从时间上看，地震4发生的时间与前3次地震间隔较远，很有可能与前3次地震不属于同一地震序列；从断层走向上看，前3次地震都是西南-东北走向的断层，而第4次地震则是西北-东南走向断层，这意味着此次地震与前3次地震并不属于同一断裂带。综上可知，前3次地震可能未直接触发第4次地震。同时，由于其他研究[2-3]计算的地震1的滑动模型采用的断层几何走向为302°，与地震4的断层走向一致，导致地震1对地震4的正向触发显著，使得前3次地震的总库仑应力表现为正向触发。而本文反演获得的地震1的断层走向与地震4的断层走向近似垂直，虽然也是正向触发关系，但由于地震1的量级较小，不能抵消地震3的负值库仑应力，因此前3次的总和为负值。综合对正演结果与原始观测形变的对比、余震与主震断层空间关系的分析、余震与主震库仑应力传输关系的分析以及反演模型残差等几个方面可知，本文反演的地震1断层模型更加可靠。因此，本文认为地震4可能不是由前3次地震触发的，而是在区域岩石特性的条件下，由累积构造应力、动态库仑应力传输、区域火山活动以及其他未知的构造活动共同导致的。
5 结语本文使用InSAR技术对2019-10~12发生在菲律宾棉兰老岛的4次震级大于MW6.0的强震序列展开研究。结果表明，此次地震序列主要由3条断裂的破裂引起。其中，地震1和地震3由同一条东南倾向断裂不同区段的先后破裂导致，同时该断裂运动以右旋走滑为主；地震2的发震断层与地震1、地震3的发震断层在空间上邻近，但反演结果显示地震2的发震断层为一条独立的西北倾向断裂，断层运动整体以右旋走滑为主并兼具一定分量的正断层破裂。地震4位于先前3次地震的东南一侧，为独立的东北倾向断裂，断层运动以左旋走滑为主。
静态库仑应力传输结果显示，在余震密集的20 km深度处，仅地震3余震与正应力变化存在较好的对应关系，其他几次地震对应关系较差，说明余震不仅受主震影响，还可能受构造背景应力、地质状况和动态库仑应力等影响。同时，地震序列间的库仑应力传输表明，地震2最显著滑动区的破裂受地震1的正向触发影响，地震2又显著增加了地震3震源区域的库仑应力；前3次地震对地震4的库仑应力传输总和为负值，即静态库仑应力传输可能并不是地震4的直接发震原因。
但上述结论受限于观测数据的丰富度和可靠性，特别是地震1仅观测获得单轨InSAR结果，导致其反演结果仍存在一定的不确定性，后续需进一步结合其他数据对模型展开精细化处理。同时，本文推测动态库仑应力的延迟触发可能影响了地震3的发生，这也有待进一步研究与论证。
| [1] |
谭皓原, 王志. 菲律宾群岛深部速度结构成像与双向俯冲板片构造特征[J]. 地球物理学报, 2018, 61(12): 4 887-4 900 (Tan Haoyuan, Wang Zhi. Seismic Velocity Imaging and Tectonic Characteristics of the Bidirectional Subducting Slabs under Philippine Archipelago[J]. Chinese Journal of Geophysics, 2018, 61(12): 4 887-4 900)
(  0) |
| [2] |
Li B Q, Li Y S, Jiang W L, et al. Conjugate Ruptures and Seismotectonic Implications of the 2019 Mindanao Earthquake Sequence Inferred from Sentinel-1 InSAR Data[J]. International Journal of Applied Earth Observation and Geoinformation, 2020, 90
( 0) |
| [3] |
Zhao L, Qu C Y, Shan X J, et al. Coseismic Deformation and Multi-Fault Slip Model of the 2019 Mindanao Earthquake Sequence Derived From Sentinel-1 and ALOS-2 Data[J]. Tectonophysics, 2021, 799
( 0) |
| [4] |
Yang Y H, Tsai M C, Hu J C, et al. Coseismic Slip Deficit of the 2017 MW6.5 Ormoc Earthquake that Occurred along a Creeping Segment and Geothermal Field of the Philippine Fault[J]. Geophysical Research Letters, 2018, 45(6): 2 659-2 668 DOI:10.1002/2017GL076417
( 0) |
| [5] |
Yang Y H, Hu J C, Yassaghi A, et al. Midcrustal Thrusting and Vertical Deformation Partitioning Constraint by 2017 MW7.3 Sarpol Zahab Earthquake in Zagros Mountain Belt, Iran[J]. Seismological Research Letters, 2018, 89(6): 2 204-2 213 DOI:10.1785/0220180022
( 0) |
| [6] |
Weston J, Ferreira A M G, Funning G J. Systematic Comparisons of Earthquake Source Models Determined Using InSAR and Seismic Data[J]. Tectonophysics, 2012, 532-535: 61-81 DOI:10.1016/j.tecto.2012.02.001
( 0) |
| [7] |
Toda S, Stein R S, Beroza G C, et al. Aftershocks Halted by Static Stress Shadows[J]. Nature Geoscience, 2012, 5(6): 410-413 DOI:10.1038/ngeo1465
( 0) |
| [8] |
Reasenberg P A, Simpson R W. Response of Regional Seismicity to the Static Stress Change Produced by the Loma Prieta Earthquake[J]. Science, 1992, 255(5 052): 1 687-1 690
( 0) |
| [9] |
Okada Y. Internal Deformation Due to Shear and Tensile Faults in a Half-Space[J]. Bulletin of the Seismological Society of America, 1992, 82(2): 1 018-1 040 DOI:10.1785/BSSA0820021018
( 0) |
| [10] |
Hsu Y J, Rivera L, Wu Y M, et al. Spatial Heterogeneity of Tectonic Stress and Friction in the Crust: New Evidence from Earthquake Focal Mechanisms in Taiwan[J]. Geophysical Journal International, 2010, 182(1): 329-342
( 0) |
2. State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, 1 Dongsan Road, Erxianqiao, Chengdu 610059, China;
3. School of Civil Engineering and Geomatics, Southwest Petroleum University, 8 Xindu Road, Chengdu 610500, China;
4. Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, 999 Xi'an Road, Chengdu 611756, China