利用空间观测技术研究青海玛多7.4级地震孕育发生变形时空特征

姜卫平，许才军，李志伟，武艳强，谭凯，耿江辉，屈春燕，郑刚，温扬茂，贺克锋，周晓慧，刘计洪，陈长云，冯光财，省天琛

1 武汉大学卫星导航定位技术研究中心，武汉 430079 2 武汉大学测绘学院，武汉 430079 3 中南大学地球科学与信息物理学院，长沙 410012 4 中国地震局第一监测中心，天津 300180 5 中国地震局地震研究所，武汉 430071 6 中国地震局地质研究所，北京 100029 7 青海省基础测绘院，西宁 810101

0 引言

2021年5月22日，青海玛多发生了7.4级地震，震中位于北纬34.59°，东经98.34°，此次地震延续了近20年中国大陆7级以上地震活动格局，即围绕巴颜喀拉地块的分布特征.作为青藏高原内部活动地块之一的巴颜喀拉地块(图1)，是一个东西方向长近2000 km，而南北方向最窄处仅200 km的长条形地块，东西横跨青藏高原腹地(李平恩等，2019).在欧亚参考框架下，南部的羌塘地块和川滇地块的东-南东向运动速率明显快于巴颜喀拉地块，而巴颜喀拉地块东向运动又明显快于其北侧的柴达木地块，因此巴颜喀拉地块南北边界以及内部近东西向断裂均表现出显著的左旋走滑特征(张培震等，2003).巴颜喀拉地块北边界的东昆仑断裂带的滑动速率呈现明显的自西向东衰减特征，其中西大滩—东大滩段全新世晚期滑动速率10～12 mm·a-1(任金卫等，1993；青海省地震局和中国地震局地壳应力研究所，1999；Van Der Woerd et al., 2000)；托索湖段滑动速率6～7 mm·a-1(Harkins et al., 2010；Kirby et al., 2007)；下大武段滑动速率为6～7 mm·a-1(青海省地震局和中国地震局地壳应力研究所，1999)；玛沁—玛曲段全新世以来的滑动速率3～5 mm·a-1(青海省地震局和中国地震局地壳应力研究所，1999；何文贵等，2006；Kirby et al., 2007);若尔盖段的滑动速率为1～3 mm·a-1(付俊东，2012)；塔藏段全新世右旋滑动速率0.8 mm·a-1，逆冲速率约为0.3 mm·a-1(Ren et al., 2013).

自20世纪90年代中叶起，青藏高原内部的一系列7级及以上强震都主要发生在巴颜喀拉地块及周边，包括2001年昆仑山口西8.1级、2008年汶川8.0级共2次8级以上地震，以及1997年玛尼7.5级、2008年新疆于田7.3级、2010年玉树7.1级、2013年芦山7.0级、2014年新疆于田7.3级、2017年九寨沟7.0级、2021年玛多7.4级等共7次7级以上地震(图1).2021年青海玛多7.4级地震的发震断裂为昆仑山口—江错断裂，是一条走向北西西的高倾角断裂，全长约370 km；该断裂发育于三叠纪砂岩、板岩中，形成于印支期，断裂活动强烈，航片卫片显示断裂在空间展布呈雁列状分布，左旋位错现象较明显，震间GNSS(Global Navigation Satellite System)观测结果给出江错断裂的滑动速率约为1.2±0.8 mm·a-1(朱亚戈等，2021)，古地震研究结果显示该断裂的地震复发周期约3500年(张裕明等，1996).

2021年玛多7.4级地震是中国大陆继1947年达日7.7级地震后(刘雷，2021)，又一次发生在活动地块内部的地震.地震发生后，一般测震学最先可以得到其震源机制、余震分布(王未来等，2021)、以及简单的地震破裂滑动模型.GNSS和InSAR(Interferometric Synthetic Aperture Radar)资料自20世纪90年代以来被广泛应用于中国大陆地壳形变研究，可以有力推进强震变形与机理研究(单新建等，2004，2017；万永革等，2008；张培震等，2009；江在森等，2009；Shen et al., 2009；Wang et al., 2011).因此有必要针对玛多地震震前形变、同震位移分布和发震断层位错开展深入的研究工作.为此，本文系统收集了青海地区的GNSS观测数据，以及Sentinel-1和ALOS-2(Advanced Land Observation Satellite-2)SAR数据，研究获取了玛多7.4级地震前GNSS形变场、同震InSAR和GNSS形变场，在此基础上反演得到了发震断层的位错滑动分布.研究结果为认识板内7级以上强震(特别是活动地块内部7级以上强震)的孕育发生机理提供了可靠的大地测量结果与模型基础.

图1 区域构造背景与公元前780以来7级以上强震分布一、二级活动地块边界、绿色箭头(数字代表地块运动速率)依据张培震等(2003)，活动断裂数据依据邓起东等(2002).Fig.1 Regional tectonic background and distribution of strong earthquakes at/above magnitude 7 since 780 B.C.Boundaries of first- and second-order active blocks and green arrows (numbers represent movement rates of blocks) are from Zhang et al. (2003), and data of active faults are from Deng et al. (2002).

1 数据与方法

1.1 GNSS数据与处理方法

为了获取高空间分辨率的GNSS形变结果，在前期获取青海地区139个GNSS测站的基础上(Zheng et al., 2017)，本文还收集了青海省卫星导航定位基准站网(以下简称青海CORS)的GNSS观测数据，共包括67个测站2017—2021年的观测资料.本文采用高精度GNSS数据处理开源软件PRIDE PPP-AR处理了这部分GNSS观测数据，表1给出了具体的数据处理策略(Geng et al., 2019；Geng and Mao, 2021)，解算得到了67个测站在ITRF2014参考框架下的单天解坐标时间序列.

表1 青海CORS GNSS数据处理策略(Geng et al., 2019; Geng and Mao, 2021)Table 1 The processing strategies of GNSS data from Qinghai CORS (Geng et al., 2019; Geng and Mao, 2021)

为了提高青海CORS站坐标时间序列的可信度，本文对环境荷载非构造影响进行扣除.具体处理包括如下步骤：利用陆地水储量模型(包括积雪深度和土壤湿度；https:∥hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH10_3H.2.1/；Rodell et al., 2004；Beaudoing and Rodell, 2020)计算了水文负载形变，时间分辨率为3 h，空间分辨率为1°×1°；采用美国国家环境预测中心(National Centers for Environmental Prediction, NCEP)发布的全球地表气压产品(https:∥psl.noaa.gov/data/gridded/data.ncep.reanalysis.surface.html；Kalnay et al., 1996)计算了非潮汐大气压负载形变，时间分辨率为6 h，空间分辨率为2.5°×2.5°；采用美国国家海洋合作计划(National Oceanographic Partnership Program, NOPP)制定的海洋环流及气候估值(Estimating the Circulation & Climate of the Ocean, ECCO)的全球海底压力格网数据(https:∥ecco-group.org/products.htm；Gross et al., 2005；Gross, 2009)计算了非潮汐海洋负载形变，时间分辨率为12 h，空间分辨率为1°×(0.3°～1°).模型分析表明，负载形变周年信号振幅平均约为GNSS的E、N、U三分量时间序列周年振幅的36.5%、47.6%、63.7%.利用经过模型处理后的测站坐标时间序列，本文计算了ITRF2014参考框架下的测站速度，然后采用Altamimi等(2017)给出的ITRF2014欧亚板块欧拉矢量(Wx=-0.085,Wy=-0.531,Wz=0.770 mas·a-1)将测站速度转换到稳定欧亚参考框架下，图2为基于青海CORS站的青海省现今地壳运动GNSS速度场.

同时，本文将青海CORS站在玛多地震震前、震后各3天的数据均值做差，解算得到了地震在测站处产生的同震位移，结果如图8蓝色和红色箭头所示.结果表明：1)最大同震位移出现在距震中约47 km的JDUO站，约为28.5 cm；东方向最大同震位移也位于JDUO站，约为28.0 cm；北方向最大同震位移则位于距离震中约86 km的KANQ站，约为21.4 cm.2)地震同震影响范围较大，在距离震中超过350 km的测站处，仍然可以探测到约3 mm的同震位移量.3)同震位移场显示玛多地震是一个典型的左旋走滑型地震，断层上盘从东向西由西南向运动逐渐转为西北向运动，下盘从西向东由东北向运动逐渐转为东南向运动.

1.2 InSAR数据与处理方法

为了获取玛多7.4级地震的高分辨率InSAR同震形变结果，本文使用了升降轨Sentinel-1卫星和升降轨ALOS-2卫星的4对SAR影像(见表2).通常情况下，Sentinel-1卫星以TOPS模式获取C波段的SAR数据，单景SAR数据幅宽可达250 km，ALOS-2卫星通常可获取条带模式(70 km幅宽)和宽幅扫描模式(350 km幅宽)的L波段SAR数据.相较于C波段SAR数据，ALOS-2卫星SAR数据的波长更长，因此在抵御失相关噪声和测量大梯度形变方面具有较为明显的优势.

在SAR数据处理方面，本文综合利用差分干涉合成孔径雷达测量方法(differential interferometric synthetic aperture radar，DInSAR；Gabriel et al., 1989)、像素偏移量追踪方法(pixel offset-tracking，POT；Michel et al., 1999)、多孔径干涉测量方法(multiple aperture interferometry，MAI；Bechor and Zebker, 2006)和子带重叠区域干涉测量方法(burst overlap interferometry，BOI；Grandin et al., 2016)处理上述4对SAR影像.具体而言，针对Sentinel-1数据，升轨ALOS-2数据和降轨ALOS-2数据分别采用8×30、8×20和32×5(方位向×距离向)的多视因子，得到了形变观测值的空间分辨率约为100 m×100 m的形变结果.在DInSAR数据处理过程中，采用了一种改进的Goldstein滤波方法(Li et al., 2008)来降低失相干噪声对结果的影响，然后将相干性小于0.8的像素进行掩膜，利用最小费用流方法(Chen and Zebker, 2002)进行相位解缠.在POT数据处理过程中，采用了128×128的匹配窗口进行像素偏移量估计，其中过采样因子为2；同时，在进行POT处理之前对Sentinel-1数据进行了去斜处理(Wegnüller et al., 2016).在ALOS-2数据的MAI处理过程中，对全孔径的SAR数据进行前后视分频后，对前视和后视孔径的SAR数据分别进行标准的DInSAR数据处理，然后再对前后视的DInSAR干涉图进行差分即可得到方位向形变(Bechor and Zebker, 2006).在Sentinel-1数据的BOI处理过程中，针对已配准的两景影像，根据方位向每一行的时间即可提取出子带重叠区域的相邻子带的两个SAR数据，这两个SAR数据可认为是MAI中的前后视SAR数据，进而类似MAI的数据处理过程即可得到BOI结果(Grandin et al., 2016).

图2 基于青海CORS GNSS观测数据的青海省现今地壳水平运动速度场(a) 负载改正前的速度场; (b) 负载改正后的速度场.Fig.2 Horizontal velocity field of present-day crustal movement of Qinghai Province derived from Qinghai CORS GNSS data(a) shows the velocity field with loading uncorrected and (b) shows the velocity field with loading corrected.

表2 本文所用的Sentinel-1卫星和ALOS-2卫星SAR数据的基本信息Table 2 The outline of the SAR data from Sentinel-1 and ALOS-2 satellites used in this study

图3为本文获取的青海玛多地震16个SAR形变场结果，可以发现成像几何相似的SAR形变场的形变空间特征较为相似.例如，降轨Sentinel-1数据的DInSAR形变场(图3e)、降轨Sentinel-1数据的POT方法LOS向形变场(图3f)和降轨ALOS-2数据的DInSAR形变场(图3m)较为相似.由于ALOS-2数据比Sentinel-1数据的波长更长，因此，相对于Sentinel-1数据(图3a, e)，基于ALOS-2数据的DInSAR方法(图3i,m)可测量的形变梯度较高，即ALOS-2数据可获得更为完整的近场同震形变结果.然而，从ALOS-2数据的方位向形变结果(图3k,l,o,p)可以看出，ALOS-2数据受到了严重的电离层扰动影响，以至于相应的方位向观测数据中的形变信号被电离层信号完全掩盖，无法用于后续的形变分析和地震解译.另外，电离层扰动也会影响LOS向的观测结果.本文采用了距离向谱分解(range spectrum split，RSS)方法(Gomba et al., 2016)用来估计并去除ALOS-2数据中LOS向观测值的电离层相位.对于POT形变观测值而言，其测量精度明显低于DInSAR和BOI方法.尤其对于Sentinel-1数据POT方法的方位向形变(图3c,g)和降轨ALOS-2数据POT方法的LOS向形变(图3n)，其相应的空间分辨率均较低(见表2)，因此相应的形变场包含了较为严重的噪声信号.尽管BOI方法仅能获取子带重叠区域的方位向形变结果(图3d,h)，但对于青海玛多地震而言，BOI方法获取的形变结果很好地揭示了断层两侧的相对运动，为后续分析解译提供了重要的数据基础.

2 玛多地震前地壳变形特征分析

在本文处理的67个青海CORS站速度场结果(图2)的基础上，进一步融合了前期得到的青海省139个GNSS测站速度(Zheng et al., 2017).青海CORS站东向和北向速度平均精度分别约为0.62 mm·a-1和0.75 mm·a-1，而前期GNSS测站东向和北向速度平均精度分别约为1.02 mm·a-1和0.86 mm·a-1，均略低于青海CORS站的精度.虽然前期青海地区GNSS数据的时间跨度为1998—2015年(Zheng et al., 2017)，跨度较长，而本文青海CORS数据的时间跨度为2017—2021年，跨度较短，但是前期GNSS数据大部分为流动观测，平均间隔1～2年才有数天观测，而本文青海CORS数据均为连续观测，所以得到的测站速度精度会相对高一些.本文青海CORS站速度场和前期GNSS速度场均转换到了稳定欧亚参考框架下，但是不同数据源之间可能会存在参考框架差异，如果两个数据源有适量的公共数据，则可以利用公共数据对差异进行拟合，但是本文两个数据源之间没有公共数据，而考虑到这种差异对速度场带来的影响通常很小，一般小于速度场本身精度(Zheng et al., 2017)，所以本文忽略了参考框架差异对速度场融合带来的影响.图4给出了最终的稳定欧亚参考框架下总共包含206个测站的玛多地震震前GNSS速度场.

图3 基于升降轨Sentinel-1卫星/ALOS-2卫星SAR数据以及DInSAR/POT/MAI/BOI获取的青海玛多地震沿卫星视线向(line-of-sight，LOS)和方位向(azimuth，AZI)的地表形变观测值紫色折线代表断层线，黑色长箭头代表对应卫星飞行方向，红色短箭头代表LOS方向.Fig.3 The surface deformation fields of the Qinghai Madoi earthquake in the line-of-sight and azimuth directions derived from the ascending- and descending-track Sentinel-1 and ALOS-2 SAR data by DInSAR/POT/MAI/BOI methods, respectivelyThe purple line represents the fault; the black long arrow shows the azimuth direction; the red short arrow shows the line-of-sight direction.

图4 稳定欧亚参考框架下的青海省现今地壳运动GNSS速度场红色箭头表示本文处理的青海CORS 67个GNSS测站的速度，蓝色箭头表示融合的Zheng 等(2017)中139个位于青海省的GNSS测站速度，误差椭圆表示95%的置信区间，绿色粗线表示昆仑断裂，红色粗线表示昆仑山口—江错断裂.Fig.4 GNSS velocity field of present-day crustal movement of Qinghai Province with respect to stable EurasiaThe red arrows show the velocities at 67 Qinghai CORS GNSS stations processed in this study, and the blue ones show the velocities at 139 GNSS stations taken from Zheng et al. (2017). Error ellipses are 95% confidence interval. The green thick line shows the Kunlun fault and the red thick line shows the KunlunshanPass-Jiangcuo fault.

图5 青海省现今地壳应变率场(a) 应变率张量第二不变量分布; (b) 最大剪应变率场.红白沙滩球表示近20年发生在青海省的3次大地震震源机制解，昆仑断裂和昆仑山口—江错断裂加粗显示.Fig.5 Present-day strain rate fields in Qinghai Province(a) Distribution of second invariant of strain rate tensor; (b) Maximum shear strain rate field. The beach-balls show the focal mechanisms of three large earthquakes occurred in Qinghai Province in the recent 20 years. The Kunlun and KunlunshanPass-Jiangcuo faults are shown in bold.

为了直观给出玛多震源区附近的变形模式和量值，基于图4所示GNSS速度场，应用球面近似方程(Savage et al., 2001；Wang and Wright, 2012)得到了应变率张量第二不变量分布(图5a)和最大剪应变率场(图5b).结果表明：1)应变率张量第二不变量分布和最大剪应变率场的高值区基本都位于2001年可可西里地震的发震断层——昆仑断裂西段及附近；2)近20年青海省发生的3次大地震(2001年可可西里7.8级地震、2010年玉树6.9级地震、2021年玛多7.4级地震)均位于应变率高值区与低值区的转换地带.这里需要特别指出，由于本文采用的GNSS资料综合了近20年观测数据的平均结果，尽管在数据解算中已经对2001年可可西里7.8级地震和2010年玉树6.9级地震等强震影响做了扣除，但难免有残留影响，因此关于这两个地震震源区的震前变形特征仅供参考.另一方面，由于图4和图5涉及的全部GNSS观测数据均来源于玛多7.4级地震前，因此该地震的震前变形特征对强震孕育研究具有重要参考意义.

3 玛多地震三维同震形变场与同震破裂模型

3.1 玛多地震同震三维形变场分析

本文用InSAR资料获取了丰富的视线向(LOS)同震形变场，但只是一维形变.GNSS同震形变较少，且垂直形变精度较低.为了直观，也为了获得较丰富的三维同震形变，我们将玛多地震InSAR和GNSS同震形变融合成高精度同震三维形变场.本文使用了一种基于地表应变模型(strain model，SM)和方差分量估计(variance component estimation，VCE)的InSAR三维地表形变分解方法(SM-VCE；Liu et al., 2018).鉴于ALOS-2数据的方位向形变观测结果(图3k,l,o,p)和降轨Sentinel-1数据POT方法的方位向形变观测结果(图3g)受电离层扰动和失相干噪声的影响较为严重，本文仅利用剩余的11个SAR形变观测数据进行分析，图6给出了基于SM-VCE方法获取的青海玛多地震的同震三维形变场.结果显示，此次地震的整体形变符合左旋走滑断层的形变特征，最大水平形变发生在破裂断层的东部，量值高达2.4 m.在震中附近的东西向形变量级较小，沿断层轨迹从震中向东和向西区域的东西向形变量级较大，并且在断层东端北侧附近，水平形变方向从西南变为西北.与水平形变相比，垂直向形变主要发生在近断层区域，且形变量级较小(-0.52～0.64 m)，并且垂直向形变沿断层走向在断层两侧出现了正负交替的现象，显示此次地震破裂的断层倾角存在非均一性特征.总体上，玛多7.4级地震的三维形变场结果显示，大部分形变主要发生在巴颜喀拉地块内部，地块北边界的东昆仑断裂带、南边界的甘孜—玉树以及鲜水河断裂带的形变不明显.

对比青海玛多地震近场GNSS的水平同震位移结果(图6b—c的圆圈处)，InSAR与GNSS东西向和南北向形变差值的均方根误差(root mean square error，RMSE)分别为5.2 cm和12.2 cm，其中InSAR和GNSS在震中最近处测点南北向形变之差高达35 cm，直接导致了南北向形变的RMSE值较大.原因可能是在断层附近破裂严重，导致GNSS站点发生了局部剧烈形变，也可能是因为InSAR结果的南北向形变几乎全部由精度较低的方位向形变观测数据得到所致，SAR南北向形变本身的精度较低.在剔除此GNSS站点之后，重新计算的东西向和南北向RMSE数值分别为5.2 cm和6.4 cm，说明本文获取的青海玛多地震三维形变场是可靠的.

3.2 三维同震形变约束下的同震破裂模型

为获得玛多地震同震破裂的精细几何结构和同震滑动分布，本文首先基于InSAR测量、现场勘察(潘家伟等，2021)和重定位余震(王未来等，2021)构建了一个倾向沿走向变化的三维同震破裂几何模型，断层整体为高倾角特征，从西北端的西南倾沿走向向东逐渐转为高倾角的东北倾.随后，将断层划分为3 km×3 km的均匀三角单元，以InSAR和GNSS观测合成的三维同震形变数据为约束，进行了玛多地震同震破裂的滑动分布反演.三维同震形变具有直观、水平形变和垂直形变精度都很高的优点，为减少冗余和提高效率，采用基于断层模型的自适应降采样方法(Lohman and Simons, 2005)降采样InSAR三维形变.同时在反演过程中为避免得到一个力学不可信且过于粗糙的滑动分布结果，对滑动分布反演施加了一个改进的平滑约束(Maerten et al., 2005)，既考虑了主-分支断层的几何关系，又有效避免了断层浅部由于传统平滑约束所导致的滑动亏损.反演在剪切模量为30 GPa、泊松比为0.25的均匀弹性半空间内进行(Meade, 2007).此外，考虑到玛多地震走滑为主的破裂特征(潘家伟等，2021)以及余震分布特点(王未来等，2021)，反演中仅允许滑动向量的滑动角在-50°～50°范围内变化，整个反演过程依据模拟残余和模型粗糙度的L曲线来确定一个合理的平滑因子，以至于使滑动分布模型在不至于过分粗糙的情况下，能够很好地拟合观测数据.

图6 青海玛多地震的同震三维形变场(a) 代表此次地震的同震三维形变场，其中底图代表垂直向形变，箭头代表降采样后的水平形变； (b)和(c)分别代表此次地震的东西向形变和南北向形变，其中圆圈代表GNSS观测站点.Fig.6 The 3-D co-seismic deformation field of the Qinghai Madoi earthquake(a)shows the 3-D co-seismic deformation field, in which the base map shows the vertical deformation, and the arrows show the down-sampling horizontal deformation. (b) and (c) show the deformation fields in the E-W and N-S directions, respectively. The circles in (b) and (c) represent GNSS stations.

模拟结果显示玛多地震破裂至地表，滑动主要分布在地表至15 km深度，最大滑动位于近地表且在地表没有明显的弱化(图7)，最大滑动量值约4 m，释放累积矩1.78×1020Nm，对应矩震级MW7.4，滑动分布沿断层呈现峰-谷交替的特征，显示断层破裂具有明显的横向不均匀性.倾滑分量主要发生在断层走向转换区，且这些转换区主要对应于拉分阶区，表明断层几何和滑动分布有很好的空间对应关系.余震主要分布于同震滑动分布高值区周缘，和高滑动分布位置呈现空间互补的特征.正演GNSS和InSAR三维形变(图8和图9)，可以发现模拟值和观测值有很好的对应性，对应的加权残余误差分别为4.9 mm和15.9 cm.表明滑动分布结果可以有效再现GNSS和InSAR两个数据集的主形变特征.两个台站GPS形变信号拟合较差的主要原因可能是由于这两个站点建立在非基岩上，因此信号受到场地效应的影响，致使其偏离了同震形变主趋势.InSAR三维形变拟合相对更差的原因可能是三维形变导出误差和InSAR观测原始数据对应解缠和大气扰动等因素共同作用所致.InSAR的残差较大测点主要集中在近断层位置，推测可能是由于弹性模型无法解释近断层的非弹性形变或主断层形变模式无法模拟近断层形变分布所致(潘家伟等，2021).此外，对比原始InSAR LOS观测和GPS联合反演结果，发现两者具有相似的滑动分布模式，意味着这两个模型具有很好的一致性，均能较好的解释相应的观测资料.相较于原始InSAR LOS观测直接用于滑动分布反演，InSAR三维形变数据由于引入了物理约束(应变)，可以有效避免原始观测数据集中局部噪声对滑动分布结果的影响，虽然这可能会导致部分区域信号被平滑.此外，由于InSAR三维形变数据有效避免了多数据集之间的定权问题，因此也就降低了不合理定权对反演结果的影响.

图7 青海玛多地震同震滑动分布Fig.7 Co-seismic slip distribution of the Qinghai Madoi earthquake

图8 玛多地震GNSS同震形变观测值和模拟值Fig.8 Observed and modeled values of GNSS co-seismic deformation of the Madoi earthquake

图9 InSAR导出和同震模拟的地表三维形变以及三方向残差第一列到第三列分别为东西、南北和垂直方向形变分量；其中(a—c)为观测值，(d—f)为模拟值； (g—i)为观测值和模拟值之差.Fig.9 3-D InSAR-derived and modeled surface deformation fields and residuals in three directionsThe subfigures in the first, second and third columns show the deformation in EW, NS and UP directions, respectively; (a—c) show the observed values; (d—f) show the modeled values; (g—i) show the differences between the observed and modeled values.

4 讨论

自1997年玛尼7.5级地震以来，巴颜喀拉地块及周边发生了一系列7级以上地震，其中2001年昆仑山口西8.1级和2008年汶川8级地震的相继发生，巴颜喀拉地块东向运动增强(江在森等，2009)，导致地块北边界以及与之邻近的近东西向走滑断裂(比如发生玛多地震的昆仑山口—江错断裂)表现出显著的左旋运动加载效应，该动力效应有利于玛多地震的孕育发生.

从长期地震危险性研究角度，玛多7.4级地震发生在中国地震局M7专项工作组给出的昆仑山—玛多危险区内(M7专项工作组，2012)，潜在的发震断裂主要为东昆仑断裂带西大滩段和昆仑山口—达日断裂等，判定依据主要为大型活动断裂带上的地震空区、2001年昆仑山口西地震破裂演化以及历史大地震的时间-序次关系等.从形变场特征与强震发生地点角度来看，此次玛多7.4级地震与2001年以来其他4次7级以上强震(2001年昆仑山口西8.1级、2010年玉树7.1级、2014年新疆于田7.3级、2017年九寨沟7.0级)具有一致性，均发生在最大剪应变率高值区边缘(武艳强等，2020).具体到发震断层尺度，GNSS速度剖面结果显示玛多地震发震断层的滑动速率明显小于其北侧东昆仑断裂带托索湖段(朱亚戈等，2021)，同时玛多地震发震断层两侧连续变形更缓慢，反映了一定程度的孕震晚期特征(Savage and Burford,1973；Meade and Hager, 2005).

地震发生后，测震学最先用地震波数据约束反演获得震源机制以及破裂滑动分布模型，为地震危险性评估提供快速参考.USGS基于震源机制的走向106°、倾角76°构建了一个平板断层面，用地震波约束反演了有限断层滑动分布，矩震级MW7.4、最大滑动量4 m的数值与本文结果一致，但其最大滑动量位于本文最大滑动量西侧约50 km处，从最大滑动处向两侧分散.Wang等(2021)虽然在地震波的基础加了Sentinel InSAR数据约束，但还是没有考虑断裂东侧的破裂情况，其最大滑动达7 m，位于本文最大滑动东侧几十公里处，埋深较大，可能更多的受到远场地震波的影响.Chen等(2021)用升降轨Sentinel数据约束的模型最大滑动5 m，最大滑动位置和深度与本文比较接近.本文用近场更多的升降轨Sentinel和ALOS InSAR数据和GNSS数据分辨出5个以上滑动峰值区，揭示了更多的破裂细节和断裂东侧的分叉破裂，与余震活动分布一致，可能更多的受到近场大地测量形变的影响.

本文的滑动分布反演结果显示2021年青海玛多MW7.4地震发生在昆仑山口—江错断裂带，破裂长度～160 km，最大破裂位于近地表，且滑动分布在地表未见明显的滑动衰减.这一滑动特征在合成3D形变中也有所反映，形变最大值位于毗连断层的两侧，显示断层滑动分布最大值至少位于近地表，否则由于浅地表的空间滤波效应，地表形变将更加连续平滑且形变峰值将偏离断层一定距离.这一特征区别于常见的伴随非成熟断层破裂产生的浅部滑动亏损(shallow slip deficiency，SSD；例如：2019年Ridgecrest地震)，表明玛多地震发震断层(昆仑山口—江错断裂)是一条发育成熟的断裂带.由于合成3D形变可能导致形变在空间上相较于多数据集原始形变更加连续平滑，但基本的破裂滑动分布特征仍然可以假定其稳定性，如最大凹凸体位置和大小，滑动分布沿走向变化的整体空间特征.结合地质现场勘察(潘家伟等，2021)和动态破裂反演(USGS, 2021)，我们可以发现玛多地震成核于近阶区位置，而进一步沿直且连续的简单断层几何传播，且在该简单几何区间形成最大破裂，最后破裂终止于两个释放阶区，符合典型的走滑地震成核-破裂-传播-终止特征(如：1999年Hector Mine地震).此外，玛多地震在断层两侧释放阶区产生大的滑动集中区可能是由于释放阶区断层一般相对更弱，更容易产生分支断裂和发生中小能量的集中释放(Scholz, 2019).

玛多地震同震变形结果显示，此次地震的应变释放主要集中在巴颜喀拉地块内部的狭长区域.地震断层破裂滑动最大达4 m，若以断层滑动1.2 mm·a-1计算，则地震复发周期是3333 a，与古地震研究结果接近.而据GNSS速度场剖面，巴颜喀拉地块的边界断层滑动速率一般在3～12 mm·a-1，地震复发周期更短.发震断裂北侧的东昆仑断裂带托索湖段和玛沁—玛曲段对地震的形变响应不明显，但玛多地震对该断裂、龙日坝断裂带以及龙门山断裂带南段呈现明显的应力加载效应.另一方面，库仑应力计算结果同样表明，东昆仑断裂带的托索湖段与玛沁—玛曲段、龙日坝断裂带为库仑应力增加段(Li et al., 2021).因此，综合断裂带对玛多地震的同震响应特征以及库仑应力影响，需重点关注上述断裂带的强震紧迫性.

5 结论

GNSS和多源InSAR资料显示2021年玛多7.4级地震导致了昆仑山口—江错断裂显著的左旋破裂，最大水平形变发生在断层东部、量值达到2.4 m，同时InSAR资料还表明垂直形变主要发生在断层近场，在断层两侧呈现正负交替现象，表明地震破裂的断层倾角在空间上呈现非均一特征.GNSS和InSAR联合反演结果显示，此次玛多地震的断层滑动主要分布在15 km以上，破裂延续到地表，最大滑动量约4 m，矩震级为MW7.4，破裂沿断层走向呈现明显的分段特征，形成了多个显著破裂区.结果表明，发震的江错断裂地震复发周期大约3000 a，青藏高原地块内部断裂与边界断裂一起，分摊消化印度—欧亚板块的汇聚变形.

致谢感谢“中国地壳运动观测网络”和青海省基础测绘院提供GNSS数据.感谢日本宇航局提供的ALOS-2数据(PER2A2N038)，感谢欧空局提供的Sentinel-1数据.