赵磊， 许文斌， 方楠， 刘计洪， 冯光财
中南大学地球科学与信息物理学院， 长沙 410000
2021年5月22日，在青海玛多发生MW7.3地震，据中国地震台网观测，震源深度为17 km，震中(34.59°N，98.34°E)位于巴颜喀拉块体内部的江错断裂附近(图1)，距离块体北边界的昆仑断层南侧约70 km(王未来等，2021；
潘家伟等，2021).欧亚板块和印度板块的持续碰撞，在青藏高原形成了多条大型断裂带，这些断裂带把整个高原划分为多个次级块体.巴颜喀拉地块为青藏高原中部的次级地块，是一条狭长的三角构造单元.巴颜喀拉地块南北两侧分别受鲜水河断裂和东昆仑断裂控制.块体向东挤出的过程中受到四川盆地的阻挡，形成了龙门山挤压造山带，即青藏高原东边界.在过去20多年时间里，多次中强地震围绕巴颜喀拉块体发生(Klinger et al., 2005; 王卫民等，2008；
Xu et al., 2010; Huang et al., 2013; Liu et al., 2014; 季灵运等，2017；
贺鹏超等，2018；
Zhao et al., 2018)，但鲜有对于该块体内部活动断层的研究.玛多地震为我们研究块体内部断裂的活动方式和构造应力特征提供了很好的机会，而精确的发震断层模型对于分析周缘断层的活动方式和构造应力起到关键作用.
图1 区域构造背景图(a) 紫色圆点表示余震分布(Wang et al., 2021)，白色五角星代表震中位置，红色沙滩球表示United States Geological Survey(USGS)玛多地震震源机制解，绿色实线表示断层破裂地表迹线，黑色实线表示活动断层，蓝色箭头表示震间GPS速率(Wang and Shen, 2020)；
(b) 红色虚线表示巴颜喀拉块体，黄色和绿色方框分别表示降轨和升轨数据范围，紫色沙滩球表示历史地震震源机制解(USGS).Fig.1 Regional tectonic setting(a) The purple dots denote the relocation aftershocks (from Wang et al., 2021), the white star represents the location of the epicenter, red beach balls show the focal mechanisms of the Madoi earthquake from United States Geological Survey (USGS), the green line indicates the rupture surface trace, the black lines denote the active faults, the blue arrows indicate the interseismic GPS velocity (Wang and Shen, 2020). (b) The red dotted line represents the Bayan Har block, the blue and green boxes represent the descending and ascending frame respectively, purple beach balls show the focal mechanism of historical earthquakes from the USGS.
已有研究成果表明，玛多地震的余震分布和InSAR同震形变场呈现一个复杂的发震断层几何结构，东西向延伸约150 km，发震断层在东端分支，在西端改变断层走向(Zhao et al., 2021; Wang et al., 2021; He et al., 2021).基于大地测量数据解算的同震三维形变场显示垂直形变在断层两侧正负交替，表明了断层倾向在空间分布上的复杂性(Zhao et al., 2021; 姜卫平等，2022).同震库仑应力扰动表明此次地震对周缘断层具有明显的应力加载作用(祝爱玉等，2021).静态库仑应力结果显示余震分布与应力加载区域有很好的空间相关性(王迪晋等，2022).虽然已有的研究提出了多种断层模型(He et al., 2021; Jin and Fialko, 2021; Wang et al., 2021; Zhao et al., 2021; 王阅兵等，2022；
余鹏飞等，2022He et al., 2022)，但大多基于余震和线性反演确定断层几何结构.对于余滑和余震的空间相关性及震后余滑的时空变化特征的研究较少.为了得到精确的断层几何结构，本文联合贝叶斯非线性反演方法和网格搜索反演方法估计断层几何和滑动分布模型.然后基于观测数据和精确断层模型反演同震和早期震后滑动分布并探讨两者之间的空间相关性，进而研究断层的摩擦属性.
1.1 同震和震后数据处理
本文的InSAR同震形变场和偏移量来自Liu等(2022)，而震后4.5个月内的形变场利用GAMMA软件处理得到.采用震后降轨21景和升轨17景影像，通过双轨差分干涉得到差分干涉图.在数据处理过程中，距离向和方位向分别采用10和2进行多视处理提高信噪比，采用90 m分辨率的数字高程模型(SRTM DEM)用于校正地形相位.采用自适应滤波和最小费用流法(Goldstein and Werner, 1998)得到解缠干涉图.之后掩膜相干性低于0.7的像素点和剔除垂直基线大于30 m的干涉对来提高观测区域相干性.最终分别获得77景降轨干涉图和68景升轨干涉图，使用开源软件LICSBAS(Morishita, 2021)进行时序震后形变场解算.首先，将闭环相位值设为1.5 rad并且计算每三个干涉图的闭环相位，剔除具有明显解缠误差的干涉图，然后计算形变速率的标准差及掩膜低相干性和高标准差像素点，最后采用时空滤波得到时序震后形变场(Li et al., 2022).
1.2 同震和震后滑动分布反演
根据同震干涉图和像素偏移量结果，将断层面分成4个部分：发震断层西端(断层1)、中部的主断裂段(断层2)、东段分叉两段(断层3和断层4).考虑观测与模型误差的不确定度，使用BEAT软件(Vasyura-Bathke et al. 2020)采用贝叶斯算法来确定发震断层主断裂段的震源参数.在反演过程中，升降轨数据采用相同权重.图2显示了贝叶斯搜索结果直方图，在95%的置信区间内，最优倾角为北倾83°，走向为285°的左旋走滑断层.在确定主断裂段断层参数后，我们采用SDM(Steepest Descent Method，最速下降法，Wang et al.,2013)基于最小二乘线性搜索其余三段断层最优倾向和倾角.基于模型的粗糙度和失配度搜索最优平滑因子为0.13.然后将断层面离散成2 km×2 km的断层块，将滑动角约束在[-45°，45°]范围内，假设泊松比为0.25来反演每个子断层块的滑动，并基于模拟观测值的均方根误差搜索最优倾角.在整个反演过程中，将先搜索的最优断层段作为约束来搜索之后断层段.首先，以断层2作为约束来搜索断层1的最优倾角，其次以断层1和断层2作为约束来搜索断层3的最优倾角，最后以断层1，2，3作为约束搜索断层4的最优倾角.因断层3和断层4位于发震断层东侧分叉处，反演时存在折衷关系，最后将断层3和断层4一起反演与单独反演做对比保证反演的精确度.
图2 断层几何参数的一维和二维后验概率密度分布图Fig.2 1-D and 2-D posterior probability density distribution plots of the fault geometry parameters
运动学余滑被广泛用来观测早期震后余滑特征，使用同震发震断层模型来反演震后滑动分布，将震后断层模型延伸至30 km宽度，其他断层几何参数与同震发震断层相同.
2.1 同震和震后形变场
玛多地震升降轨同震形变场整体呈东西向分布在断层两侧，断层两侧最大视线向(LOS)形变量约为1 m(图3).升降轨相反的形变信号表明玛多地震以走滑运动为主导.发震断层东端分叉处的形变信号阶跃表明发震断层几何结构的复杂性，与昆仑断裂在东端形成多条分支相似(Li et al., 2011; Ceylan et al., 2012).横贯同震形变场的剖面BB′在断层两侧信号突变，表明发震断层破裂至地表.图3g，h所示为Sentinel-1震后4.5个月升降轨震后形变场，升降轨形变信号相反表明震后形变依然以走滑运动为主.
图4和图5为震后降轨、升轨时序形变场，震后形变主要分布在断层两侧近场区域，在远场震后形变量级较小.早期震后形变主要分布在四个区域，区域2位于震中西侧同震形变量级较小的区域，在震后阶段呈现明显的震后形变，震后4.5个月断层两侧累积形变量达5 cm.在同震破裂西侧延伸区域(区域1)震后最大形变量约为3 cm.而同震时破裂至地表且同震位移最大区域(区域3)存在明显震后形变，最大LOS向形变量达5 cm.总体上，震后形变主要分布在同震形变较小的区域，而区域3震后形变与同震形变区域重叠，表明地震周期内能量释放的连续性和复杂性.
图4 玛多地震Sentinel-1降轨震后时序形变场Fig.4 Time-series of postseismic deformation after the 2021 Madoi earthquake in Sentinel-1 descending orbit
图5 玛多地震Sentinel-1升轨震后时序形变场Fig.5 Time-series of postseismic deformation after the 2021 Madoi earthquake in Sentinel-1 ascending orbit
2.2 同震滑动分布
基于非线性和线性反演，得到了玛多地震发震断层的精确几何结构.在非线性反演得到的最优主破裂段几何约束下，断层1最优倾角为北倾75°(图6a)，分开搜索结果显示东端分叉处断层3最优倾角为81°(图6b)，断层4最优倾角为84°(图6c)，两段都为南倾.断层3和断层4联合反演结果(图6d)与单独搜索结果相同，确保了断层倾角搜索结果的准确性.断层1和断层2北倾，断层3和断层4南倾，而震后短期余震分布也呈现在断层中西端位于发震断层北侧，而在东端分段处位于断层南侧.余震的空间分布特征验证了发震断层的几何特征.
图6 (a) 断层1最优倾角，蓝绿色五角星表示最优值；
(b) 断层3最优倾角；
(c) 断层4最优倾角；
(d) 断层3和断层4联合搜索最优倾角；
(e) 玛多地震同震滑动分布；
(f) 玛多地震震后滑动分布，蓝色圆点表示精定位余震，彩色直线表示断层地表迹线，红色五角星代表震中Fig.6 (a) The optimal dip angle of Seg.1, the cyan star represents the optimal value; (b) The optimal dip angle of Seg.3; (c) The optimal dip angle of Seg.4; (d) The optimal dip angle of jointly inversion of Seg.3 and Seg.4; (e) Coseismic slip distribution of the Madoi earthquake; (f) Postseismic slip distribution of the Madoi earthquake, the blue dots represent the aftershocks, the colored lines represent the rupture surface trace, the red star represents the epicenter
图6e为玛多地震同震滑动分布，整体呈现长条状分布在上地壳浅部.同震滑动分布破裂至地表且主要分为四个区域：断层2中东部和中西部的浅部滑动；
断层1和断层3的浅部滑动.最大同震滑动量约为4 m位于断层2中东部，平均滑动角为5°，同震滑动分布整体呈现左旋走滑运动.模拟残差图如图7所示，较大的残差值主要分布在断层近场相干性较差的区域.假设剪切模量为30 GPa，玛多地震同震释放能量为～1.4446×1020Nm，等同于矩震级MW7.37，与Global Centroid-Moment-Tensor(GCMT)和USGS一致.
图7 玛多地震同震形变场(a)(d)(g)(j)分别为ALOS2降轨、ALOS2升轨、Sentinel-1降轨、Sentinel-1升轨观测值，(b)(e)(h)(k)相应的模拟值，(c)(f)(i)(l)相应的残差值，白色五角星代表震中.Fig.7 Coseismic deformation field of the Madoi earthquake(a), (d), (g), (j) ALOS2 descending, ALOS2 ascending, Sentinel-1 descending and Sentinel-1 ascending data, respectively; (b), (e), (h), (k) The corresponding modeled data; (c), (f), (i), (l) The corresponding residuals, white star represents the epicenter.
同震滑动分布在上地壳浅部不连续并存在多个滑动亏损区域，一处位于断层2的两个主要同震滑动区域之间，以及东端同震滑动未破裂至地表区域.同震滑动分布周围未破裂以及滑动亏损区域通常为速度强化体或者同震时应力未达到破裂阈值的速度弱化体(Yamanakaand Kikuchi, 2004；
Miyazaki et al., 2004).因此，滑动亏损区域以及未破裂区域在震后阶段的活动特征对于区分它们的物理属性具有重要意义.
采用Jackknife重采样技术(Melgar et al., 2015)来评估断层滑动模型的稳定性.随机从数据集中采取50%的数据进行反演，一次反演完成后，将数据放回数据集再重新随机抽取50%数据进行反演，如此重复进行50次并计算每个断层块的标准差.图8a所示为同震滑动分布的子断层标准差，最大标准差大约为0.2 m，位于同震最大滑动处，约占同震滑动量级的4%，表明了数据与模型的强依赖关系.
图8 (a) 同震滑动分布模型不确定度，彩色直线表示断层地表迹线，红色五角星代表震中；
(b) 震后早期余滑分布模型不确定度；
(c) 同震静态库仑应力变化；
(d) 发震断层速率状态摩擦参数a-bFig.8 (a) Model coseismic slip parameter uncertainties, the colored lines represent the rupture surface trace, the red star represents the epicenter; (b) Model early afterslip parameter uncertainties; (c) Coseismic static Coulomb stress change; (d) The calculated rate and state friction parameter a-b of the seismogenic fault
2.3 震后滑动分布
震后滑动分布采用的升降轨影像对为2021-05-26—2021-10-05，升降轨震后模拟数据都能较好的拟合观测数据(图9).震后余滑分布在多个区域(图6f)，在上地壳浅部，余滑整体分布在同震滑动亏损区域内，最大滑动量约为20 cm.在中下地壳存在30 km×10 km的深部余滑，震后4.5个月累积余滑量约为14 cm.最大余滑位于断层2两个主要同震滑动区域之间的滑动亏损区域内.而多个余滑区域显示较大的同震滑动区域在震后阶段持续滑动.在发震断层东端和西端延伸区域，余滑分布在断层14～15 km深度范围内，在断层3表现为深部余滑.早期震后余滑释放能量约为5.0417 × 1018Nm，等同于矩震级MW6.4，震后短期内的快速余滑是同震应力扰动释放的主要方式.Heki和Tamura(1997)发现震后24 h内的余滑约占同震滑动的30%，多项研究也表明短期余滑在整个震后周期内占据较大比重(Kim and Dreger, 2008; Yano et al., 2014; Liu and Xu, 2019；
Ragon et al., 2019 ).震后第一景影像距离同震发生时间为震后4天，缺失了短期内的震后形变，因此我们得到的震后矩震级低于实际量级.
图9 玛多地震震后形变场(a)(d) 分别为Sentinel-1降轨、升轨观测值; (b)(e) 相应的模拟值；
(c)(f)相应的残差值，白色五角星代表震中.Fig.9 Postseismic deformation field of the Madoi earthquake(a) (d) represent Sentinel-1 descending and ascending data, respectively; (b) (e) The corresponding modeled data; (c) (f) The corresponding residuals, white star represents the epicenter.
采用相同的方法来评估震后数据与模型的匹配度，图8b所示为震后滑动分布的子断层的均方根误差，较大的均方根误差主要分布在断层浅部震后余滑较小的区域，整体均方根误差～1 cm，表明了震后数据与模型的强依赖关系.
2.4 余滑空间分布揭示的摩擦性质
震后余滑一般分布在同震滑动区域周围，在同震阶段未发生破裂的障碍体上，以及沿发震断层走向和倾向延伸的区域(孙权等，2022).余滑一般分布在同震库仑应力加载区域内(Dogan et al., 2014; Wallace et al., 2018).运动学余滑通常被用来观测短期震后余滑的时空分布特征.本文使用Coulomb 3.3软件计算了玛多地震断层面上的同震库仑应力变化(图8c)，余滑整体分布在同震库仑应力加载区域内.
同震滑动区域是速度弱化体成核区域，余滑一般围绕同震滑动区域分布，主要位于周边的速度强化体区域.前面所提及的，上地壳浅部震后余滑部分位于同震滑动量值最大区域.Yabe和Ida(2018)通过摩擦属性非均匀分布的断层面模拟实验表明，如果同震阶段断层未能充分释放积累的能量，残余应力会在震后阶段以余滑和余震的方式进行释放.地下介质在空间分布上是非均匀的，同震破裂时快速的滑动可以有效的穿过局部的速度强化体.此外，断裂带的摩擦属性在地震周期中也会发生变化，当摩擦属性沿断层发生变化并且同时存在速度弱化体和速度强化体，这些断层区域会同时参与同震和震后滑动.综上，我们认为上地壳浅部同震和震后滑动重叠区域主要是由于残余应力的释放驱动了震后余滑在这些区域的滑动.相似的震后余滑位于同震滑动区域在很多地震事件中可以观测到(Diao et al., 2014; Iinuma et al., 2016; Liu and Xu, 2019).震后余滑的传播也会对断层面的应力加载产生影响，进而引发余震的发生.余震与余滑在空间上形成很好的互补关系，基本填充了同震库仑应力加载区域.位于深部中下地壳的滑动分布，震后呈现线性滑动状态(图10b)，其中包含了下地壳黏塑性形变机制黏弹性松弛的影响，可能意味着中下地壳较为软弱的物质层的存在，因此之后的研究中应重点关注震后远场区域应力驱动余滑和黏弹性松弛的折衷关系，分析青藏高原东北部下地壳和上地幔的流变性.
图10 (a) 同震、震后早期余滑和余震归一化累积地震矩随深度变化；
(b) 震后早期余滑地震矩随深度的时序变化Fig.10 (a) The normalized accumulated seismic moment of coseismic, early afterslip and aftershocks with depth; (b) Temporal evolution of early afterslip seismic moment with depth
2.5 震后余滑的时空分布
时序震后余滑模型显示余滑量级随时间逐渐增加.基于震后滑动分布计算的震后时序地震矩(图10b)显示，上地壳浅部余滑在2～8 km深度范围内快速衰减，震后3个月后趋于稳定.归一化的同震和震后滑动分布显示(图10a)，震后滑动分布在同震滑动量较小的区域，两者在空间上形成良好的互补关系.Gualandi等(2014)表明长期的余滑位于距离同震破裂较远的区域.实验室模拟实验也表明浅部和深部余滑速率不同，有效正应力是与深度相关的，深部余滑速率小于浅部余滑(Ariyoshi et al., 2007).这表明同震库仑应力扰动触发了上地壳浅部的速度强化体产生快速余滑，深部余滑随着距离同震破裂区域的距离增加趋于稳定滑动.对于大地震，余滑可以持续几年到几十年(Barbot et al., 2009)，浅部余滑在震后4.5个月基本释放，深部将是未来主要余滑区域.
本文基于ALOS-2和Sentinel-1同震形变场反演获得了玛多地震主破裂段断层几何参数，并确定了发震断层在西端和东端分叉段的倾角.同震滑动模型显示断层位错主要发生在上地壳浅部并且破裂至地表，最大滑动量为4 m，位于2 km深处.利用N-SBAS技术获得震后4.5个月升降轨时序震后形变场，最大震后形变量近场约为5 cm，远场约为2 cm，上地壳浅部最大余滑量约为20 cm.震后余滑在上地壳浅部分布在多个区域，且在深部存在长条状余滑.时空分布的震后余滑表明上倾角和下倾角区域不同的摩擦属性，上地壳浅部在震后阶段快速余滑，而深部余滑可能持续更长时间.
致谢感谢European Space Agency (ESA) 提供的免费Sentinel-1雷达影像数据，日本Japan Aerospace Exploration Agency(JAXA)提供的ALOS-2雷达影像数据(ER3A2N521).本文所有图片采用GMT和MATLAB绘制.
猜你喜欢同震玛多发震另一种时间观读者(2023年3期)2023-02-28——以盈江地区为例">基于构造应力场识别震源机制解节面中发震断层面——以盈江地区为例地球物理学报(2022年11期)2022-10-312021年玛多MS7.4地震的深部构造背景地震地质(2021年4期)2021-09-27基于钻孔应变观测约束的2016年**呼图壁M6.2地震的发震断层研究地震研究(2021年1期)2021-04-132021年5月22日青海玛多MS 7.4地震总结地震地磁观测与研究(2021年6期)2021-03-19云南思茅大寨井水位地震同震响应特征分析*地震研究(2018年4期)2018-11-23芦山地震发震构造及其与汶川地震关系讨论中国学术期刊文摘(2016年1期)2016-02-13芦山地震前后介质波速变化与GPS应变场相关性研究∗地震科学进展(2015年9期)2015-05-13芦山Ms7.0地震引起的水位同震响应特征分析四川地震(2014年2期)2014-12-02川滇地区钻孔四分量应变仪记录的同震应变阶分析1震灾防御技术(2013年4期)2013-11-26