2021年青海玛多MW7.4地震GNSS同震形变场及其断层滑动分布

王迪晋，王东振，赵斌*，李瑜，赵利江，王阅兵，聂兆生，乔学军，王琪

1 中国地震局地震研究所，武汉 430071 2 中国地震局地震大地测量重点实验室，武汉 430071 3 中国地震台网中心，北京 100045 4 青海省地理信息与自然资源综合调查中心，西宁 810001 5 中国地质大学(武汉)地球物理与空间信息学院，武汉 430074

0 引言

据中国地震台网中心测定，北京时间2021年5月22日2时4分在我国青海省玛多县境内发生MS7.4地震，震中位于34.59°N，98.34°E，震源深度17 km.美国地质调查局(USGS)给出的震中位于玛多县黄河乡东侧(34.598°N，98.251°E)，震级为MW7.3，震源深度10 km，可能的发震断层为一条近EW走向、带有正断分量的左旋走滑断层(https:∥earthquake.usgs.gov/earthquakes/eventpage/us7000e54r/).地表破裂考察、余震分布、InSAR干涉影像等结果显示，玛多地震的发震断层为昆仑山口—江错断裂，地表破裂长约160 km，可分成四个段落，自西向东分别为鄂陵湖南段、黄河乡段、冬草阿龙湖段和昌麻河乡段(李智敏等，2021；王未来等，2021；华俊等，2021).昆仑山口—江错断裂是巴颜喀拉块体的内部断层，也是东昆仑断裂带的次级分支构造之一(李陈侠等，2011).近30年来，中国大陆7级以上的地震事件都发生于巴颜喀拉块体周缘断裂带上(邓起东，2012；邓起东等，2014)(图1).此次地震是中国大陆地区继2008年汶川地震以来震级最大的一次构造地震，也是2017年九寨沟MS7.0地震后又一次震级大于7级的板内地震.该地震的发生再次引起人们对巴颜喀拉块体内部及其边界断裂带上的地震活动性的强烈关注.

巴颜喀拉块体位于青藏高原北部，是青藏高原内部物质东向挤出的重要构造单元之一.该块体的北边界为东昆仑断裂；南边界为甘孜—玉树断裂；东边界为龙门山断裂；西边界与阿尔金断裂带交汇，这些边界的构造活动共同控制着块体的现今变形.现今地壳速度场显示青藏高原东北缘整体向北东方向运动(图1)，巴颜喀拉块体内部区域形变相对较小.伴随着巴颜喀拉块体东向运移速率自西向东逐渐衰减，北边界东昆仑断裂带的左旋走滑速率也随之衰减(Kirby et al., 2007; Lin and Guo, 2008)，在该断裂带的东段发育了多条次级分支断裂，如甘徳—玛多断裂、昆仑山口—江错断裂和达日断裂等(图1)，其中1937年在主断裂托索湖段发生M7.5地震(Molnar and Deng, 1984；沈正康等，2003)、1947年在达日断裂上发生MS7.7地震(Molnar and Deng, 1984；万永革等，2007).地理位置上，玛多地震同样发生在靠近东昆仑断裂带的区域，发震断层与相邻的其他次级分支断裂一起构成了巴颜喀拉块体北部宽阔的弥散型边界带(潘家伟等，2021).震间GNSS速度剖面的结果显示巴颜喀拉块体北部边界断裂带两侧出现明显的速度梯度变化，而本次地震发震断层的滑动速率仅为1.2 ± 0.8 mm·a-1，比主边界带东昆仑断裂托索湖段的滑动速率小(朱亚戈等，2021)，表明发震断层昆仑山口—江错断裂应力积累较慢，复发周期长，而这次地震很可能是长期震间构造应力积累下的自发破裂事件.

玛多地震在巴颜喀拉块体内部形成了怎样的形变特征？发震断层的同震滑动分布如何？对周缘断裂会造成哪些潜在的影响？这些都是亟待回答的科学问题.李智敏等(2021)通过详细的野外地震考察，获得了地震地表破裂组合特征和运动学性质、破裂分段性、同震位移特征及分布规律，并根据冲沟、道路和拉张阶区裂隙宽度确定地表同震位移量为1～2 m.潘家伟等(2021)通过经验公式估算该次地震的最大位移量约为4 m，并对发震构造进行了分析和探讨，同时论述了发震断裂与东昆仑断裂带的关系、巴颜喀拉地块构造变形的动力学机制以及区域未来强震活动性等问题.华俊等(2021)基于欧洲航空局Sentinel-1卫星C波段的升、降轨数据获取了玛多地震InSAR同震形变场，确定了地震的地表破裂迹线，并结合库仑应力扰动分析了玛多地震引起的周边区域的地震危险性，反演得到的矩震级为MW7.45，断层最大滑动量约为6 m，主体破裂集中在0～10 km深度.但由于近场失相干效应明显，利用InSAR资料得到的同震形变场缺少极震区的近场观测，这也导致了由此反演得到的震级和最大滑动量均略大于地震波的结果.

图1 (a) 青藏高原背景速度场； (b) 2021年玛多地震的构造背景图(a)中的黑色粗实线给出了巴颜喀拉块体的主边界，黑色细实线给出了东昆仑断裂附近发育的多条次级分支断裂，红色实线表示玛多地震发震断层的地表破裂轨迹，虚线矩形框指示了图(b)的位置； 图(b)中绿色圆点给出了余震序列的精定位结果(王未来等，2021).图中的速度场结果引自文献(Wang and Shen, 2020).Fig.1 (a) Background velocity field of the Tibetan plateau; (b) Tectonic setting of the 2021 Madoi earthquakeIn figure (a), the thick black solid line denotes the main boundary of the Bayan Har massif, the thin black solid lines denote secondary branch faults around the East Kunlun fault, and the red solid line denotes surface track of the rupture of the Madoi earthquake. The dashed rectangle box indicates the position of figure (b). The green circles in figure (b) represent relocated aftershocks (Wang et al., 2021). The velocity field in the figures is cited from (Wang and Shen, 2020).

随着GNSS技术的迅速发展，该项技术已经成为现代地震研究领域中的一个重要观测手段(Bock and Melgar, 2016).尤其是近场GNSS观测可以弥补InSAR影像在近场失相干区域的观测不足，直观诠释发震机制，同时为地震动态破裂过程和静态破裂分布模型提供重要的约束.李志才等(2021)展示了破裂带附近9个GNSS基准站的高频(1 Hz采样)动态形变波形，及12个GNSS连续站的同震形变场，观测到的最大水平同震位移0.6 m，距离断层20 km处.该研究虽然初步揭示了本次地震的形变特征，但由于点位稀疏，对断层滑动分布的约束不足，破裂模型稍显粗糙.

本文利用多家单位共同完成的GNSS观测数据，在给出更为全面的同震形变场和高频GNSS动态波形的基础上，分析了本次地震的同震形变特征及其对周边断裂的影响，比较了GNSS观测与InSAR卫星影像的同震形变结果，反演得到精细的地震破裂模型，最后计算了地震在发震断层面上产生的库仑应力变化，与余震精定位结果进行了比较，并分析了巴颜喀拉块体相关区域在未来的地震危险性.

1 变形观测与数据处理

1.1 观测资料

玛多地震发生以后，中国地震局地震研究所第一时间派工作组赶赴震区开展GNSS应急流动观测.利用前期国内不同机构、不同研究团队在震区布设的地壳形变观测站点，在地震发生后十天内完成了53个站点的流动观测(每站观测24～48 h).这些流动站主要包含陆态网络工程的区域站，原国家测绘局B级网点，自然资源部于2019年新建的国家大地控制点和少量的三角点.此外，还有相关科研机构在震中附近建立的流动GNSS站点，主要包括了东昆仑断裂托索湖段近南北向GNSS剖线(正好位于发震断层北侧，Diao et al., 2019)和中国地震局地震研究所在该断裂玛沁—玛曲段布设的一条GNSS剖线.

我们还收集了震区300 km范围内的53个GNSS连续观测站，主要包括陆态网络工程的连续观测站和青海省的CORS站(卫星导航与位置服务基准站).其中部分连续站点同时还具有高频采样数据，因此我们也收集到了其中17个台站1 Hz采样的高频GNSS观测数据.图2展示了所有站点的分布图.

1.2 静态形变场处理

本文所有的GNSS数据处理均采用GAMIT/GLOBK10.6软件 (Herring et al., 2015)，数据处理策略与流程与Zhao等(2015)一致，对固体潮、极潮、海潮进行改正，并应用最新的卫星、天线绝对相位中心改正模型和全球气压温度模型GPT2(Lagler et al., 2013)，估计得到包括测站坐标、卫星轨道、天顶对流层延迟的单日松弛解.接着用GLOBK将区域松弛解与SOPAC(Scripps Orbital and Permanent Array Center)解算的全球IGS(International GNSS Service)站的单日松弛解合并，得到一个包含全球IGS站和本文GNSS站的单日松弛解.最后在全球范围内选择～60个用于实现参考框架转换的参考站，以全球单日松弛解作为准观测值，利用GLOBK通过7参数的相似变换得到ITRF2014(Altamimi et al., 2016)下的单日坐标解.

图2 玛多地震震区周围GNSS点位分布图中红色实线表示发震断层的地表破裂轨迹，黑色实线表示震中附近主要活动断层.Fig.2 The distribution of GNSS sites around the Madoi earthquakeThe red solid line denotes surface track of the rupture of the Madoi earthquake and the black solid lines denote the main active faults around the epicenter in the figure.

贝叶斯后验概率密度统计方法已被广泛应用于震间断层运动参数、同震断层滑动分布、震后形变估计等研究(Sun et al., 2013；周新, 2017；Ingleby et al., 2020).本文引入一个基于MCMC(Markov Chain Monte Carlo)采样的贝叶斯方法来估计同震位移及其误差.该算法可以充分顾及震间速率估计的不确定性，尤其适用于处理流动资料.

在地壳运动研究中，GNSS坐标时间序列可表示为(Nikolaidis et al., 2001; Bock and Melgar, 2016)：

y(t)=a+bt+csin(2πt+φ1)+dsin(4πt+φ2)

(1)

式中a为初始形变量，b为长期线性速度，两个正弦函数表示周年和半周年季节性变形，g为在T时刻的同震形变量，H为阶跃函数(Heaviside function)，ε为观测误差.对于流动站点，由于采样数据不足，不考虑周期项，同时忽略早期震后效应；而对连续观测站而言，由于计算同震形变只用到了地震前后几天的数据，亦可忽略震后效应及周期形变.因此，待估参数仅剩初始形变量a、长期速度b和同震形变量g.

依照贝叶斯理论，模型参数的后验概率密度分布可表示为(Elliott et al., 2016):

(2)

式中d为GNSS观测得到震前、震后两期观测结果，m为模型参数，P(m)为模型参数的先验概率密度分布，Cd为数据的协方差矩阵，G(m)为将采样的模型参数代入公式(1)得到的预测值.

首先利用沈正康等(Shen et al., 1996)提出的应变率模拟方法(strain rate modeling method)进行空间插值，得到长期速率初始值bmod.插值得到的长期速率必然存在误差，我们利用均一分布(uniform distributions)给出其先验概率密度分布，取值范围为[bmod-2,bmod+2] mm·a-1.同震形变量g的取值范围定为 [-2,2]m，同样采用均一的先验概率密度分布.采用Metropolis-Hastings采样算法(Mosegaard and Tarantola, 1995)生成Markov链，包含6000个采样，在计算时先燃烧(burn-in)前3000个采样，以保证后一段进入收敛区域, 最终求得估计模型参数并给出其误差.由此得到玛多地震106个GNSS测站的同震位移场，图3给出水平方向的同震位移场(数值结果见电子附表1)，同震位移平均误差为4 mm.

1.3 高频GNSS动态形变数据处理

震时1 Hz高频GNSS数据处理采用武汉大学PRIDE课题组发布的开源GNSS处理软件PRIDE PPP-AR(Geng et al., 2019)进行动态解算.首先，通过对数据的预处理剔除粗差和周跳，再利用最小二乘法进行参数估计，以获取浮点解，最后再利用相位钟和相位偏差产品来完成整周模糊度的解算，获得固定解.

本文获取了陆态网络连续站和青海CORS站部分站点的高频波形数据，图4给出了其中17个台站(9个陆态网络连续站、8个青海CORS站)在玛多地震期间的三分量高频波形记录.陆态网络连续站中距离断层最近的青海玛多站(QHMD，距离断层34 km)形成了30 cm的永久形变(最大振幅南北向达到约16 cm，东西向达到约28 cm)，而青海CORS站中有多个台站均形成了永久形变，包括靠近断层的JDUO站(距离断层41 km)，该站也形成了近30 cm的永久形变.另外距离断层不到80 km的KANQ站和HSHX站均形成了明显的永久形变，其他台站也均记录到了明显的动态波形.

2 同震变形特征分析

在图3中我们利用三种比例尺的箭头，给出了不同震中距范围内的观测站点的同震形变结果.从图3中可以看出，远场连续GNSS观测站表现出明显的四象限分布特征(图3红色箭头所示)，断层西北、东南两侧测站同震位移均向远离断层的方向运动，而东北、西南两侧测站的同震位移则指向断层方向.而在近场范围内(图3蓝色箭头所示)，GNSS观测到的最大水平同震形变为1.2 m，位于野马滩附近(大野马岭站，编号：4454)，距离断层约100 m.另一个同震水平形变大于1.1 m的站点位于昌麻河乡附近(昌麻河乡站，编号：4499)，距离断层近5 km.

观测结果显示，沿地表破裂两侧区域50 km范围内有20 cm以上的地表永久变形，震中范围300 km仍可观测到同震形变，影响范围最远波及至东南方向的四川甘孜炉霍一带.在地表破裂迹线西北角距离断层20 km的鄂陵湖沿岸，有多个GNSS测站显示出有20～30 cm的NW向运动，而更靠近震中的玛多县城附近的站点则表现为NWW向运动.地表破裂迹线东南角距断层25 km处的陆态区域站J406有幅度约30 cm的SE向运动.

图3 玛多地震水平方向同震位移场不同颜色、不同比例尺箭头表示同震水平位移，红色实线表示发震断层的地表破裂轨迹，浅黑色实线表示震中附近主要活动断层.两个灰色虚线框标示了两条跨断层的GNSS剖面：P1，跨江错断裂剖面；P2，跨东昆仑断裂玛沁—玛曲段的剖面.Fig.3 Horizontal coseismic displacements of the Madoi earthquakeThe arrows with different color and different scale denote the horizontal coseismic displacements. Red solid line denotes surface track of the rupture and shallow black solid lines denote main active faults around the epicenter. Two gray dashed boxes denote two GNSS profiles, which cross the Jiangcuo fault (P1) and the Maqên-Maqu segment of East Kunlun fault (P2), respectively.

图4 连续GNSS站记录到的1 Hz高频波形数据上行给出了8个青海CORS站的记录，下行给出了9个陆态网络工程的连续站的记录.左图右侧给出的公里数表示台站的震中距.灰色虚线标示了发震时刻.Fig.4 1-Hz high-frequency waveforms recorded by the continuous GNSS sitesThe top row plots the waveforms recorded by 8 Qinghai CORS sites and the bottom row plots the waveforms recorded by 9 CMONOC continuous sites.The kilometers at the right sides of two left subfigures denotes epicentral distance of every sites.The gray dashed line in every subfigure denotes the original time of the event.

图5 两条GNSS剖面水平同震形变的投影结果(a)和(b)为跨发震断层江错断裂的GNSS剖面结果，(c)和(d)为跨东昆仑断裂玛沁—玛曲段的GNSS剖面结果.上行为平行于断层位移，下行为垂直于断层位移.F1: 东昆仑断裂，F2：久治断裂，F3: 玛多—甘德断裂，F4: 达日断裂.Fig.5 The projection of horizontal coseismic deformation of two GNSS profiles(a) and (b) are the projection of the GNSS profile which cross the seismogenic fault—the Jiangcuo fault; (c) and (d) are the projection of the GNSS profile which cross the Maqên-Maqu segment of East Kunlun fault.The top row is the projection parallel to the strike of faults and the bottom row is the projection perpendicular to the strike of faults.F1: East Kunlun fault, F2: Jiuzhi fault, F3: Madoi-Gande fault, F4: Dari fault.

近场GNSS观测结果较为清晰地勾勒出了发震断层两侧的运动模式.断层北侧台站位移方向由东到西，从西南向逐渐转为北西西方向、北西方向；南侧台站位移则呈现出一种顺时针变化的变形形态.近场形变在西北—东南方向上表现出明显拉张变形特征，在东北—西南方向表现出挤压变形特征，近场形变模式与远场的四象限变形特征相符.发震断层两侧的近场同震形变表现为明显的左旋走滑特征，与震源机制一致.

利用近场GNSS测站，本文构建了两条GNSS剖面(图3中的两个灰色虚线框)，分别为跨发震断层的剖面P1(长宽分别为400 km和180 km)和一条跨东昆仑断裂玛沁—玛曲段的剖面P2(长宽分别为300 km和80 km).图5给出了这两条GNSS剖面水平同震形变在平行于断层和垂直于断层方向上的投影结果.

跨江错断裂剖面P1平行于断层方向上的位移呈现出明显的双曲线特征(如图5a所示)，近断层同震形变较大，平行于断层的同震形变随断层距离的增加而逐渐衰减.GNSS观测到的同震形变场直观地展示了发震断层的左旋走滑性质.垂直于断层方向的同震位移相对较小(如图5b所示).

在复测的GNSS站点中，距断层最近的两个站点为位于野马滩附近的4454(距离断层～100 m)和位于黄河乡附近的2836(距离断层～1 km)，分别位于同震破裂的北侧和南侧.上述两个站点观测到的平行于断层的同震位移分别为～1.2 m和～0.6 m，因此点位间位移差为～1.8 m.这一量值与地表调查(潘家伟等，2021)获得的～2 m的地表破裂量值相当.

东昆仑断裂托索湖段距离玛多地震震中仅～70 km.跨东昆仑断裂的同震GNSS位移剖面显示(如图5a、b蓝色粗实线指示了东昆仑断裂F1的位置)，东昆仑断裂两侧站点同震响应显著.这些站均向西运动且东昆仑断裂南侧的站点形变量较大，这与东昆仑断裂左旋走滑的震间形变特征相反，因此玛多地震的发生对该断裂段的震间应力积累有缓解作用.

跨东昆仑断裂玛沁—玛曲段的剖面P2显示，断层两侧水平运动差异同样较为明显(见图5c、d)，南侧站点向远离地震震中的东南方向变形，最大值可达6 cm，断层北侧点位朝西南方向运动，最大值不到2 cm.东昆仑断裂玛沁—玛曲段两侧的地壳差异运动将造成该断层的库仑应力加载，未来地震风险性将进一步增强.

3 与InSAR图像的对比

基于合成孔径雷达差分干涉测量(D-InSAR)技术，利用欧空局(European Space Agency, ESA)C波段升、降轨哨兵1号卫星数据可获取青海玛多地震的InSAR同震形变场.InSAR资料显示玛多地震产生的同震形变影响空间范围广，形变场的长轴呈NWW向，升、降轨观测到的形变量符号相反(见图6).总体而言，玛多地震所在区域干涉条件很好，InSAR观测形成的干涉条纹连续，所获的同震形变清晰，为与GNSS观测相互验证提供了良好的资料.然而，由于受到黄河水系、湖泊、沼泽等地貌影响，在断层近处InSAR影像存在失相干区域.因此，近场的GNSS观测为失相关区域的形变特征提供了关键约束.

本文计算了4条垂直于断层InSAR剖线(图6中绿色实线)的同震观测结果，连同在图3中剖面P1的GNSS位移观测一并展示在图7中.这4条剖线分别处于四段地表破裂的中间位置.以图6中每条InSAR剖线为中轴线，东西各展开25 km的距离，形成4条GNSS剖面.结合InSAR升、降轨方向矢量，将各GNSS剖面中的观测站点的三维同震位移转化为LOS(Line-of-Sight)向形变，将其与对应的InSAR剖线升、降轨LOS向位移进行比较(见图8).利用升、降轨InSAR数据得到的最大LOS向形变量约为0.8 m，与由GNSS三维位移导出的LOS向形变量大体一致(见图7、图8).

为了拟合玛多地震同震InSAR观测的升、降轨LOS向形变和GNSS观测的位移，本文采用一种半无限空间浅表位错模型(shallow dislocation model)来表征地表位移与同震滑动的关系.在同震阶段，地表位移与发震断层的宽度(倾角方向)和平均滑动量可以简单的表示为(Segall, 2010)：

(3)

式中u3为平行断层走向的位移，s为断层平均滑动量，d1为断层宽度，x1为测站离开断层的距离，其中s和d1为待求模型参数.假设走滑断层同震破裂中平行于断层走向的位移远大于垂向位移，即忽略垂向位移，则可将平行断层走向的同震位移u3转换成对应的InSAR观测LOS向位移.

图6 欧空局哨兵1号卫星给出的玛多地震InSAR同震形变场(a)为升轨观测结果，(b)为降轨观测结果.两者形变量符号正好相反.Fig.6 InSAR coseismic deformation field caused by the Madoi earthquake from Sentinal-1A satellite of the ESA(a) is LOS displacements along the ascending orbit and (b) is LOS displacements along the descending orbit.The sign of the two displacement values are opposite.

图7 4条InSAR跨断层剖线的LOS向位移观测上图为升轨观测，下图为降轨观测.图中同时给出了前文跨江错断裂的GNSS剖面的位移结果，为方便与InSAR观测比较，GNSS位移已经通过方向矢量转化为LOS方向的投影.灰色粗实线标示了发震断层的位置.Fig.7 The LOS displacements of 4 InSAR profiles crossing the seismogenic faultThe top figure is LOS displacements obtained from the ascending orbit and the bottom figure is LOS displacements obtained from the descending orbit.The coseismic deformation of the GNSS profile crossing the Jiangcuo fault is also shown in two subfigures and the displacements had been converted to the projection along LOS with the direction vector due to the comparison with InSAR observation.The thick gray solid line in every subfigure denotes the seismogenic fault.

图8 4条InSAR剖线与对应的GNSS剖面结果的比较圆点为InSAR剖线结果，三角形为GNSS剖面的结果(已利用升、降轨方向矢量转化为对应的LOS向位移).上行为升轨结果，下行为降轨结果.红色实线为拟合结果，灰色竖直粗实线标示了发震断层的位置.Fig.8 The comparison between 4 InSAR profiles and the corresponding GNSS displacementsThe solid dots denote the LOS displacements of the InSAR profiles and the triangles denote the corresponding GNSS displacements in the LOS direction converted by the direction vectors of the ascending and descending orbit.The top row is the results along the ascending orbit and the bottom row is the results along the descending orbit.The red solid line is the fitting result and the vertical gray solid line denotes the seismogenic fault in every subfigure.

将图8中各条InSAR剖线的升、降轨LOS向位移，以及对应的GNSS剖面位移的LOS向投影作为观测数据，估计模型参数.考虑到模型参数与观测数据之间的非线性关系，采用格网搜索方法来确定s和d1的最佳取值.通过调节两个参数的取值使得观测数据与模拟结果之间的残差最小，以达到最佳的数据拟合效果.图8每列图形的行首给出了各条剖线的模型参数最佳取值，拟合结果由每个子图中的红色实线表示.参数估计结果显示，最大平均滑动量出现在冬草阿龙湖段，达到了3.765 m，其断层破裂宽度为6.8 km；而最大的断层破裂宽度则出现在昌麻河段，达到13.5 km，平均滑动量为1.565 m.

拟合曲线总体上较好地拟合了各条剖线的LOS向位移，尤其在断层近场附近，多数剖线上的拟合结果与实际观测的LOS向位移均有较好的符合.而昌麻河段InSAR剖线P4在断层北侧近场附近的升、降轨LOS向位移均与拟合曲线有相对较大的偏差.昌麻河段真实断层分布较为复杂，地表破裂野外考察(潘家伟等，2021)和InSAR影像的结果，均证实了在地表破裂昌麻河段存在分支断层.在该段的拟合过程中，未考虑分支断层，因此可能造成了拟合曲线在该段与实际观测出现相对较大的偏差.此外，走滑断裂端部的垂向位移相对较大，与假设不符，这是出现较大偏差的另一因素.

4 同震滑动分布及其应力扰动

利用本文获取的远近场GNSS同震形变场结果，我们构建了简单的断层几何模型，采用约束最小二乘算法反演了断层滑动分布.InSAR成像及野外地质考察(李智敏等，2021；潘佳伟等，2021)显示玛多地震地表破裂呈现明显的分段特征，断层走向总体为285°，但每段存在差异.如前文所述，在东侧昌麻河段还存在分支断层破裂.余震精定位结果(王未来等，2021)也显示出较为明显的分段性，但两者在地表并不重合，可能说明断层倾角的变化，也可能与余震定位的偏差有关.为此我们依据InSAR影像识别的破裂迹线构建断层面，并假设断层倾角为90°，宽度为30 km.断层倾角的微弱变化并不会改变采用GNSS资料反演的滑动分布基本特征.我们将带有分支断层的曲面离散化为1574个三角形子块，三角形的平均边长由近地表2 km随深度逐渐增加至5.5 km.

基于GNSS同震形变，本文反演得到的玛多地震破裂模型见图9.结果显示，同震破裂具有明显的分段性，并且同震滑动出露地表，与野外地表破裂考察和余震分布吻合.最大的破裂区域处于冬草阿龙湖段，黄河乡段和鄂陵湖南段均有凹凸体存在，滑动量略小.反演结果显示地震主体破裂位于断层面0～10 km的浅部区域，主破裂面最大滑动量达到了4.6 m，地震矩1.63×1020N·m，矩震级为MW7.4.昌麻河段分支断层最大滑动量可达2.3 m.

图9所示的玛多地震破裂模型可以很好地拟合GNSS观测的同震形变场(见图10)，拟合残差为3.9 cm.通过单一手段给出的断层滑动分布对原始观测量的恢复程度较好，说明断层几何参数的选取相对可靠.近场的GNSS同震形变结果可以有效填补InSAR观测在极震区由于失相干造成的观测空白，对断层滑动破裂提供有效约束.

华俊等(2021)基于InSAR观测反演得到的玛多地震断层滑动分布结果显示，断层最大滑动量约为6 m，矩震级为MW7.45.与之相比，本文基于GNSS观测反演得到的最大滑动量和矩震级都略小，但同震破裂的分段性特征一致，且均破裂到地表，主体破裂的深度范围也较一致，均集中于断层面0～10 km的浅部区域.

根据图9所示的滑动分布模型，假设摩擦系数为0.4，计算了断层面上的静态库仑应力变化(见图11).结果显示库仑应力增加区域主要集中在同震滑移下倾方向，在黄河乡段中部和冬草阿龙湖的西段存在浅部应力增强的区域.静态库仑应力增加区域与余震分布具有一致性，说明余震主要是由静态库仑应力加载而触发的.

5 讨论

观测到的同震形变显示黄河乡段和昌麻河段同震滑动显著，这一结论与InSAR结果相符(华俊等，2021).但地表调查结果显示鄂陵湖南段的地表破裂最为显著，黄河乡段和昌麻河段地表破裂相对较小(潘家伟等，2021).我们推测造成上述差异的原因可能有两方面.其一是鄂陵湖南段水系发达，地表土层含水量较高，在同震破裂中可能发生较强的非弹性形变，导致相对显著的地表破裂.其二是鄂陵湖南段的GNSS站点较少，且缺乏近断层站点约束，导致我们的同震模型对该区域的同震滑动分辨能力不足.

图9 基于GNSS同震形变反演得到的玛多地震破裂模型Fig.9 The rupture model of Madoi earthquake inverted from the GNSS observed coseismic displacements

图10 (a) GNSS观测的玛多地震水平方向同震位移与破裂模型模拟结果的比较； (b) 模拟值和观测值的拟合残差Fig.10 (a) The comparison of GNSS observed horizontal coseismic displacements of Madoi earthquake with simulated results by the rupture model; (b) The residuals between the observation and the simulations

图11 玛多地震在发震断层上产生的库仑应力变化图中灰色圆圈标示了精定位得到的余震位置(王未来等, 2021).Fig.11 The Coulomb stress changes on the modeled fault imposed by the Madoi earthquake Open gray circles represent relocated aftershocks (Wang et al., 2021).

根据震间GNSS观测约束的滑动速率，朱亚戈等(2021)估计得到相同破裂规模的地震复发周期(T)为1100～5500年.震后大地测量约束的青藏高原不同区域的下地壳稳态黏滞系数η≈1×1019Pa·s(Huang et al., 2014; Zhao et al., 2017; Diao et al., 2018; Liu et al., 2019)，由此计算的震后松弛时间(Tr=2μ/η,μ为剪切模量)远小于地震复发周期.在此条件下，震间的滑动速率强烈依赖于上次地震的离逝时间(te q), 震后早期滑动速率大于整个地震周期平均速率，而在晚期则低于平均速率.由此推测，朱亚戈等(2021)估计的结果一定程度上低估了滑动速率，高估了地震复发周期.为此，可以结合玛多地震的震后资料开展更深入的研究.

巴颜喀拉块体内部的GNSS观测站点分布较为稀疏，难以约束块体内部断裂的震间闭锁深度(朱亚戈等，2021).2021年玛多地震的同震破裂集中于断层面0～10 km浅部，本文第3节给出的4条剖线的平均断层宽度为9.6 km，亦小于10 km，表明本次地震的主要破裂并未向深部扩展，因此孕育本次地震的昆仑山口—江错断裂的震间闭锁深度很可能在10 km左右.前人利用震间大地测量观测反演的青藏高原活动断层的平均闭锁深度多为15～20 km(如Wang et al., 2011; 李煜航等，2014).昆仑山口—江错断裂“低于预期”的闭锁深度表明，震间大地测量观测获取的地表形变与常用的弹性位错模型(Savage and Burford，1973)可能容易高估活动断层(特别是滑动速率较低的块体内部断层)的闭锁深度.

长期以来，我们对类似的板块内次级分支断裂的关注度一向偏低，如龙日坝断裂、发生2017年九寨沟地震的虎牙断裂都属于这种类型.玛多地震的发生引起了我们对次级断裂活动性及其孕震能力的重新审视和思考.

6 结论

本文使用GNSS观测技术获取了玛多地震较大范围的同震形变场，揭示了玛多地震在巴颜喀拉地块内部的形变特征及对周边断裂造成的潜在影响.

(1)同震形变结果显示，玛多地震为左旋走滑型地震，同震变形场呈明显的四象限分布，观测到的最大同震位移达到1.2 m.GNSS和InSAR数据较为一致地展示了玛多地震同震形变的表现方式和变化幅度，近场GNSS为InSAR近断层失相干区域提供了必要的补充和检核，直观诠释了昆仑山口—江错断裂左旋走滑的发震机制.

(2)采用GNSS数据反演同震滑动分布的结果显示，发震断层的滑动破裂存在多个凹凸体，破裂分段特征明显且出露地表，与野外地表破裂考察和余震精定位结果一致，主体破裂位于断层面0～10 km的浅部区域，最大滑动量达到4.6 m，地震矩1.63×1020N·m，矩震级为MW7.4.同震破裂引起发震断层面上的应力扰动，余震主要集中在同震滑移下倾方向库仑应力增加的区域.

(3)从玛多地震的同震形变场空间特征来看，形变主要集中在巴颜喀拉块体内部，玛多地震的发生对东昆仑断裂托索湖段的震间应力积累有缓解作用；而东端玛沁—玛曲段的同震响应形成了明显的断层南北两侧的差异运动，因此玛多地震对东昆仑断裂玛沁—玛曲段有应力加载作用，未来地震风险性将进一步增强，需要加强监测.巴颜喀拉地块周缘强震具有明显的跳跃性特征，长期的大区域GNSS观测将为今后监测地块边界及其内部断层活动性、划定地震危险性区域提供可行的技术途径.

致谢感谢中国地质大学(武汉)熊熊教授提供相关点位的站点信息.