2021年青海玛多7.4级地震的同震变形分析

杨君妍， 孙文科， 洪顺英， 苑争一， 李瑜， 陈伟， 孟国杰

1 中国地震局地震预测研究所， 北京 100039 2 中国地震台网中心， 北京 100045 3 中国科学院大学地球与行星科学学院， 北京 100049

0 引言

2021年5月22日凌晨2时4分，青海省果洛州玛多县发生7.4级地震(以下简称玛多地震)，是我国继2017年九寨沟7.0级地震后的又一次震级超过7级的地震.地震位于34.59°N， 98.34°E，震源深度17 km，震中区人烟稀少，造成13人受轻伤.青海玛多地震发生于青藏高原东北缘的玛多—甘德断裂带上(詹艳等, 2021; 张喆和许力生, 2021)，该断裂带位于青藏高原东北缘巴颜喀拉块体北边界东昆仑断裂以南70 km，整体呈N105°E走向，全长约151 km(潘家伟等, 2021).新生代以来，由于太平洋板块向欧亚板块的俯冲汇聚以及印度板块向北的持续碰撞挤压，造就了中国大陆西侧挤压造山以及东侧拉张成盆的基础格局.青藏高原及其周缘地区作为活跃造山带因而形成了强震多发的态势.伴随着高原的持续扩展，青藏高原周缘地区存在明显的侧向挤出，其运动学特征存在明显的受复杂断层运动组合影响，包括沿昆仑断层的左旋走滑运动、青藏高原内部的正断层以及青藏高原撞击中国中部稳定的大陆地壳时产生的反向运动等(Yin and Harrison, 2000; 熊仁伟等, 2010).青藏高原北缘巴颜喀拉块体及其周缘地区构造变形强烈，1997年以来，先后发生2001年昆仑山口西8.1级地震、2008年汶川8.0级地震、2010年玉树7.1级地震、2013年芦山7.0级地震、2017年九寨沟7.0级地震等9次7级以上地震(图1).尤其集中在东昆仑断裂带和龙门山断裂带，该区域是中国大陆中强地震的发震地，也是科学家持续关注的地震多发地带.

图1 玛多地震地理位置及其周缘7级以上大地震分布Fig.1 The geographical location of the Madoi earthquake and distribution of large earthquakes above M 7 in the surrounding area

全球导航卫星系统GNSS和干涉合成孔径雷达InSAR是观测地震变形最重要的现代大地测量观测技术.GNSS可以精确检测出地表形变，尤其是水平方向精度非常高，广泛应用于地壳变形、地球动力学和地质构造的相关研究中.如，GNSS观测数据的累计滑移量绘制出青藏高原和中国大陆水平运动速度场，从而推动了中国大陆地壳运动与活动块体特征的研究(顾国华, 2005; Gan et al., 2007)；提供同震、震后以及震间地壳形变信息，作为经典地震学研究地震破裂的有力补充(Sabadini et al., 2009; Via et al., 2012).InSAR数据也可以提取地震形变量，在地壳变形和地震监测方面发挥了重要作用.2021年玛多M7.4地震发生以后，大量研究机构开展了灾区应急考察和震后GNSS野外测量工作.此次地震科考主要任务是调查地表破裂、监测余震序列和地壳形变、探测地下介质，研究地震发震构造和地震孕育环境.上述观测资料为了解地震机理，解释地表位移场以及反演断层滑动模型等提供了重要基础数据.为了完成这些数据解释或断层反演工作，使用合理的地震位错理论非常重要.因此，通过实际震例研究地震位错理论的适用性是一个重要的科学问题.

地震位错理论是地震断层面滑动与地球物理场变化之间的理论关系，是研究震源机制、反演断层滑动、解释测量数据等必要的理论基础，在地震学、地球动力学以及大地测量学中具有核心地位.地震位错理论相当于震源与地表变形之间的系统函数：一方面可以根据震源参数计算地表的同震变形；另一方面可以根据地表变形反演震源参数.自Steketee (1958)将位错理论引入到地震学后，地球物理学家对半无限空间均匀地球模型下的同震变形做了大量的研究(Berry and Sales, 1962; Press, 1965; Savage and Hastie, 1966; Davis, 1983).Okada(1985, 1992)总结上述工作，给出了半无限空间介质内部的同震变形计算方法.由于缺少地球曲率和层状构造的影响，半无限空间地球模型下的同震变形与真实大地测量观测数据差别较大.Sun(Sun, 1992; Sun and Okubo, 1993; Sun et al., 1996, 2009)和Pollitz (1997)基于更实际的地球模型，如SNREI(球对称、不旋转、弹性和各向同性)，给出了准静态球形位错理论.该理论可以计算任意断层在地表面任意位置产生的同震变形，包括位移、重力、大地水准面、倾斜和应变等.球形地震位错理论可以精确地解释大地震产生的地表重力和位移变化(Sun and Okubo, 1998; Wang et al., 2014; Yang et al., 2015; Zhou et al., 2018; Yang and Sun, 2020).理论研究结果表明，相对于均匀半无限空间地球模型，地球分层效应引起的偏差会达到25%，在理论模拟时不容忽视(Sun and Okubo, 2002; Dong et al., 2014).目前，球形地球模型位错理论已经发展到了考虑地球内部三维不均匀结构的研究，如利用微扰原理研究球形地球横向不均匀性对同震变形的影响(Fu and Sun, 2007, 2008)；利用大地测量观测数据反演地震位错Love数和格林函数等(Yang et al., 2015; Yang and Sun, 2020; 杨君妍和孙文科, 2020).精确的地球模型位错理论对于解释大地测量观测数据、反演断层滑动模型的科学问题至关重要.

另一方面，地震发生后，科学家们通常根据不同观测数据给出不同的断层滑动模型，然而，这些断层滑动模型往往差异很大.玛多地震发生后，很多科研团队提出了不同的断层滑动模型，例如，美国地质调查局USGS、中国地震局地质研究所InSAR工作组(简称地质所模型)和北京大学张勇教授工作组等(详见下节).上述断层滑动模型能否合理解释大地测量观测数据或者哪个模型的理论计算值更吻合观测数据，是一个值得探讨的问题.本文利用2021年玛多M7.4地震的三个断层滑动模型分别计算了半无限空间地球模型(Okada, 1985)和球形分层地球模型(Sun et al., 1996)下，巴颜喀拉地块及其周缘的同震位移场，并与GNSS观测同震位移和InSAR视线向形变量对比分析，确定出较优的断层模型和地震位错理论.同时，分析讨论2021青海玛多7.4级地震的同震变形特征和动力学机制.

1 玛多地震断层滑动模型

为了解释大地测量技术观测到的同震变形，需要选择合理的断层滑动分布模型.青海玛多7.4级地震发生以后，许多机构和研究者发布了基于不同观测数据获取的断层滑动模型.本文选取了三家不同机构提供的断层模型，其中，两个模型由地震波形数据反演得到，另一个则由InSAR数据反演得到，断层最大滑动量、走向和倾角都不同.三个断层滑动模型分别来自于USGS (https:∥earthquake.usgs.gov/realtime/product/finite-fault/us7000e54r_1/us

/1623171757136/basic_inversion.param)、地质所模型和北京大学张勇工作组(Zhang et al., 2021)(http:∥geophy.pku.edu.cn/tpxw/353689.htm)(图2).

图2 玛多地震断层滑动分布(a) USGS模型； (b) 地质所模型； (c) Zhang等(2021)模型.Fig.2 The fault slip distribution of the Madoi earthquake(a) USGS model; (b) Institute of Geology model; (c) Zhang et al. (2021) model.

USGS模型采用宽频带地震波数据，分析了24个远震P波、28个SH波和79个长周期面波资料.震中定为(34.613°N，98.246°E)，震源深度10 km，震级为MW7.3，走向为106°，倾角为76°.其断层分别沿走向方向和倾角方向均匀地分成52×9块子断层，每个子断层的尺寸为3.5 km×3.5 km(图2a).中国地震局地质研究所InSAR工作组使用欧空局(ESA) Sentinel-1 SAR数据，获取了此次地震的InSAR视线向形变量.InSAR形变约束的发震断层以左旋走滑运动为主，倾向北东，倾角约80°，主破裂深度位于10 km以上，震源位于(34.59°N，98.34°E)，震级为MW7.45.其断层分别沿走向和倾角方向均匀地分成43×7块子断层，每个子断层的尺寸为5.03 km×4.99 km(图2b).张勇团队震后5小时发布了由35个宽频带地震波数据反演得到的断层滑动模型，震中定为(34.586°N，98.255°E)，震源深度10 km，震级为MW7.4，走向为92°，倾角为67°.其断层分别沿走向方向和倾角方向均匀地分成41×5块子断层，每个子断层的尺寸为5 km×5 km(图2c).

三个模型中，由InSAR数据获取的断层滑动模型的最大滑动量最大，约5.0 m；地震波反演得到的断层模型的最大滑动量分别为3.30 m和3.86 m.地质所模型相应的震级也大于其他两个模型.由于三个模型采用的数据和反演方法不同，地震波获取的断层数据具有固定的走向和倾角，而InSAR模型则根据发震断层，走向和倾角随之变化.为了找到最佳模型，有必要比较三个断层滑动模型理论同震位移与观测值之间的差异.接下来将用玛多地震GNSS同震水平位移、InSAR视线向形变量分别与理论计算值进行比较.

2 大地测量观测的同震变形

2.1 GNSS同震水平位移

本文使用了李志才等(2021)提供的同震变形，包括近场基于PPP-AR技术处理的2 h时段9个高频(1 Hz采样)GNSS数据的同震位移和远场12个CORS站震前、震后各3天低频(30 s采样)数据的非差PPP(Precise Point Positioning)解算同震形变场.QHMD、MADU和QHAE三个站点变形在水平量级上保持着良好的一致性，均显示往西约0.24 m，往北约0.09 m的位移；QHAJ站点变形水平分量接近0.6 m，在所有CORS站中最大，移动方向朝西偏南，垂向变形约有5 mm的隆升(图3).

图3 玛多地震的GNSS同震水平位移Fig.3 Coseismic horizontal displacement of the Madoi earthquake by GNSS

图3同震形变结果显示玛多地震是一个典型的左旋走滑断层，断层上盘往西北方向运动，下盘向东南方向运动，发震断层破裂方式与青藏高原向东运动一致.据中国大陆构造环境监测网络提供资料，距离震源30多公里的玛多台站记录到东西向永久位移约25 cm；中国地质科学院地质研究所应急考察团队发现，在接近G214国道处地表破裂切过公路形成最大左旋位错量约20 cm，最小约5 cm.这些GNSS观测同震变形和地质考察资料为认识玛多地震的发生机理提供了基础资料.

2.2 InSAR同震形变量

利用Sentinel-1A/1B卫星的IW模式SLC(Single Look Complex)数据，对青海地震的同震变形场进行信息提取.数据处理采用ISCE软件(Rosen et al., 2012)进行，通过“二轨法”对干涉像对进行处理，外部地形采用SRTM 30 m数据进行拟合.经过相位解缠、地理编码以及四叉树降采样后，升、降轨的同震变形场如图4所示.在降轨形变场中，破裂带北盘呈LOS向拉伸运动，南盘呈LOS向缩短运动，而在升轨形变场中，形变特征恰相反.升降轨的同震形变场呈现截然相反的特征，反映了该地震变形场以东西方向的水平运动为主，表现出明显的左旋走滑特征.在升轨变形场上，北盘最大的LOS向缩短变形约1.10 m，南盘最大的LOS向拉伸变形约0.77 m；降轨变形场上，北盘最大的LOS向拉伸变形约1.21 m，南盘最大的LOS向缩短变形约1.11 m，断层两侧的相对LOS向运动可达2.30 m，反映了该地震的地表破裂具有明显的错动位移.

图4 玛多地震的InSAR同震形变量(a) 升轨数据； (b) 降轨数据.Fig.4 Coseismic deformation of the Madoi earthquake observed by InSAR(a) Ascending data; (b) Descending data.

3 理论同震位移的模拟计算与比较

为了利用GNSS和InSAR数据比较不同断层滑动模型的同震位移，我们选用均匀半无限空间位错理论(Okada, 1985)和球形分层地球模型位错理论(Sun et al., 1996)两种模型.Okada(1985)总结前人的工作并给出了一套同震变形的简洁计算公式.由于该公式定义在平面笛卡尔坐标系下，有必要将球坐标系转换至统一坐标系下.通过公式计算均匀半无限空间下的同震位移，泊松比选用v=0.25.最后再将结果通过走向角转换回地方坐标系下得到同震变形量.

在球形地球模型中，采用一个球形的、无旋转的、弹性和各向同性(SNREI)的地球模型PREM (Dziewonski and Anderson, 1981).PREM模型是目前较为完善的地球模型，考虑了地球径向非均匀性，但该模型仍为一维平均地球模型，无法反映研究区域的局部构造影响.而玛多地震位于巴颜喀拉块体东北缘，该块体是青藏高原到扬子克拉通的过渡带，也是青藏高原隆升、物质东流速度自西向东由快到慢的汇聚地带(徐锡伟等, 2008)，在地质构造上GNSS水平位移速度场与周边块体存在差异(陈长云等, 2013; 杨光亮等, 2020).本文通过Crust1.0模型(Laske et al., 2013)提供的地球横向三维数据，修改100 km深度以内的PREM模型参数.Sun等(2009)定义了4种独立的震源：水平走滑震源、垂直倾滑震源、水平引张震源和垂直引张震源，并给出了相应的格林函数.通过对格林函数的样条插值和组合就能求得任意地震在地表面任意位置产生的同震变形.利用这两种位错理论计算得到的三个断层滑动模型在GNSS观测站引起的同震位移显示在图5中.

图5 玛多地震的理论同震水平位移(a) 球形位错理论； (b) 半无限空间位错理论.黑色箭头表示GNSS观测值；红色箭头表示USGS模型理论同震位移；蓝色箭头表示地质所模型理论同震变形；绿色箭头表示Zhang等(2021)模型同震变形.Fig.5 The theoretical coseismic horizontal displacement of the Madoi earthquake(a) Spherical dislocation theory; (b) Half-space dislocation theory. Black arrows represent GPS observations; Red arrows represent the theoretical coseismic displacement of the USGS model; Blue arrows represent the theoretical coseismic deformation of the Institute of Geology model; Green arrows indicate the theoretical coseismic deformation of the Zhang et al. (2021) model.

利用三个不同断层滑动模型计算的理论同震水平位移显示，近场GNSS观测站有明显的水平形变量，并且大部分站点的理论值大于观测值.远场观测站同震水平位移相对较小，只有1～4 cm，说明玛多地震的同震形变场随震中距离衰减较快.地质所模型下的同震位移与GNSS观测值吻合较好；其他两个断层滑动模型在QHAJ站的理论计算值与观测值的大小和方向都有所差异，并且在不同地球模型位错理论下表现出类似的特征.从三个地球模型计算结果的比较分析可以看出，通过InSAR数据反演的断层滑动模型的理论值与GNSS观测值的符合程度总体上优于其他两个模型的计算结果.

图6和图7是三个断层滑动模型分别在球形地球模型位错理论和半无限空间地球模型位错理论下的同震变形，经处理后变为视线向形变量，左侧为升轨数据，右侧为降轨数据.从图中可以看出，球形分层地球模型下的理论同震变形升轨数据的最大相对形变量要大于观测值，降轨数据的最大相对形变量略小于观测值.尤其是利用地震波数据反演的断层滑动模型计算的LOS方向升轨数据最大形变量较观测值大0.3 m(图6).三个断层滑动模型中，地质所模型理论值与InSAR视线向符合得更好.均匀半无限空间模型和球形层状模型得到的结果基本类似(图7).

图6 利用球形位错理论计算的玛多地震视线向形变量(a)和(b)、(c)和(d)、(e)和(f)分别为USGS、地质所和Zhang等(2021)的断层滑动模型，左侧为理论模拟升轨同震形变量，右侧为理论模拟降轨同震形变量.Fig.6 The LOS deformation of the Madoi earthquake calculated by the spherical dislocation theory(a) and (b), (c) and (d), (e) and (f) are the fault slip models of USGS, Institute of Geology and Zhang et al. (2021), respectively. The left side is the theoretical simulation ascending coseismic deformation, and the right side is the theoretical simulation descending coseismic deformation.

图7 利用半无限空间位错理论计算的玛多地震视线向形变量(a)和(b)、(c)和(d)、(e)和(f)分别为USGS、地质所和Zhang等(2021)的断层滑动模型，左侧为理论模拟升轨同震形变量，右侧为理论模拟降轨同震形变量.Fig.7 The LOS deformation of the Madoi earthquake calculated by the half-space dislocation theory(a) and (b), (c) and (d), (e) and (f) are the fault slip models of USGS, Institute of Geology and Zhang et al. (2021), respectively. The left side is the theoretical simulation ascending coseismic deformation, and the right side is the theoretical simulation descending coseismic deformation.

为了定量说明三个断层滑动模型之间的差异和两种位错理论的差异，我们采用了统计均方根(RMS)的方法衡量.RMS定义如下：

(1)

GNSS数据和InSAR升轨/降轨形变量的RMS数值结果列于表1中.表1显示由球形分层位错理论得到的形变量比均匀半无限空间的形变量更符合大地测量观测数据，表明球体位错理论的层状效应不能忽略.通过比较2008年汶川地震(MW8.0)、2011年日本东北大地震(MW9.0)在不同地球模型下理论与观测数据的同震位移，也得出类似的结论(Wang et al., 2010; Pollitz et al., 2011; Zhou et al., 2012).无论是球形层状位错理论还是半无限空间位错理论，地质所的断层滑动模型为这三者中最佳，GNSS水平位移的RMS是0.04 m，InSAR升轨数据形变量的RMS是0.06 m，降轨数据形变量的RMS是0.11 m，优于其他两个断层滑动模型.

表1 球形与均匀半无限空间理论同震位移和大地测量观测之间的RMSTable 1 RMS comparison of spherical/half-space modeled and geodetic observed coseismic displacement

综合上述对GNSS同震位移和InSAR视线向形变量的理论模拟分析得出，中国地震局地质研究所通过InSAR数据反演的2021年青海玛多7.4级地震断层模型的理论结果与观测值最为吻合.玛多地震属于浅源地震，由于地质所使用的InSAR升降轨数据覆盖了震源区域，较其他两个断层滑动模型所使用的全球地震波数据，更加清晰的约束地震的断层几何形状和浅源滑动分布.

4 青海玛多7.4级地震的同震变形

经过上述研究与讨论，本文选取球形分层地球模型位错理论和中国地震局地质研究所InSAR工作组反演的断层滑动模型及参数来模拟计算2021年青海玛多7.4级地震的同震变形，主要包括位移、大地水准面、重力和应变变化，为后期解释大地测量观测数据和重力观测数据提供理论依据.结果如图8—10所示.

图8 球形位错理论计算的玛多地震同震位移场(a) 东西方向位移，向东为正； (b) 南北方向位移，向北为正； (c) 垂直位移，地面向上为正.Fig.8 Coseismic displacement field of the Madoi earthquake calculated by the spherical dislocation theory(a) In the east-west direction, east is positive; (b) In the north-south direction, north is positive; (c) In vertical displacement, ground-up is positive.

图9 球形位错理论计算的玛多地震同震大地水准面和重力变化(a) 大地水准面变化； (b) 重力变化.Fig.9 Coseismic geoid and gravity changes of the Madoi earthquake calculated by the spherical dislocation theory(a) Geoid changes; (b) Gravity changes.

图10 球形位错理论计算的玛多地震同震应变变化(a) 东西方向应变； (b) 南北方向应变； (c) 体应变； (d) 剪切应变.Fig.10 Coseismic strain changes of the Madoi earthquake calculated by the spherical dislocation theory(a) East-West strain; (b) North-South strain; (c) Dilatation strain; (d) Shear strain.

图8显示，东西向位移主要以玛多—甘德断裂为中心，两侧发生相对错动.上盘最大位移达到1.6 m，下盘最大位移为0.9 m(图8a).南北向位移在破裂带西段以北向运动为主，东段则以向南运动为主，相对位错约为1.2 m(图8b).垂直位移在断层上盘西部表现出明显的上升运动，下盘则下沉，沿断裂带东部，则表现出相反的运动，最大相对垂直位移达到0.42 m(图8c).水平位移总体以左旋走滑运动为主兼正断层运动方式，断层两侧沿断裂带向外扩张，呈现明显的南北向破裂带；垂直位移相对较小，在破裂带东部拐角处呈现最大垂直位移.

在断层附近的区域用球形层状模型位错理论和半无限空间地球模型位错理论计算的理论值与GNSS和InSAR观测值都吻合较好，而在远离断层破裂带的区域球形分层地球模型的理论值与观测值大小较为接近.其中用InSAR数据反演得到的断层参数及滑动模型的理论值和观测值更为接近，这可能与该模型使用的观测数据覆盖了整个断层面和地质调查的地震破裂资料等有关.

图9显示利用球形层状地球模型位错理论计算的2021年玛多M7.4地震同震大地水准面变化和重力变化.同震大地水准面在破裂带附近为下降，最大达0.6 mm，在破裂带北侧最大上升0.2 mm.总体上破裂带呈负变化.同震重力变化也有类似形态，断层面整体负变化，断裂带拐角南侧显示30 μGal的正变化，在震源东侧出现最大负变化255 μGal.该理论同震重力变化是基于地表面计算的，由于青海地震震级小于M8.0，很难被重力卫星GRACE检测到(Sun and Okubo, 2004)，所以本文尚未计算空间固定点的重力变化.

最后，我们计算了球形地球模型下的同震应变变化(图10).同震应变变化在断层线附近变化较大，变化范围达-5×10-5～5×10-5.由于应变是位移的微分，随着震中距的增加应变迅速衰减.上述应变理论值可以用来解释应变观测结果和计算应力分布.

5 讨论与结论

2021年青海玛多7.4级地震是印度板块和欧亚板块碰撞的结果，该地震位于巴颜喀拉地块北部边界东昆仑断裂带以南70 km，属于块体内部地震.地震发生后，中国大陆构造环境监测网络(陆态网)项目组GNSS观测数据显示，距离震源约30 km的玛多台站(QHMD)记录到该地震产生约25 cm的东西向永久位移.为了更好地拟合和解释大地测量观测数据，本研究选取三个地震滑动模型(USGS模型、中国地震局地质研究所InSAR工作组模型和Zhang等(2021)模型)和两种地球模型位错理论(球形地球位错理论和半无限空间模型位错理论)进行了比较.结果表明，地质所通过InSAR数据反演的断层滑动模型跟观测值更加吻合.另外，球形分层地球模型的理论值更接近于观测位移场，明显优于半无限空间地球模型.基于地质所InSAR数据反演的断层滑动模型，本文进一步计算了在球形地球模型下的同震变形，包括位移、大地水准面、重力和应变变化等.这些结果揭示了地壳变形的空间分布，为反演地震断层、研究地球内部结构、解释大地测量观测数据提供参考依据.值得注意的是，目前所使用位错理论是基于球对称地球模型讨论的，该模型为一维平均模型，无法考虑震源区域非均匀性结构的影响.通过大地测量观测数据反演地震位错Love数和格林函数，调节局部地球模型参数，对解释大地测量观测现象、精确反演地震断层滑动分布、研究地球内部结构等问题具有重要的理论价值和应用意义，今后需进一步研究讨论.

致谢感谢李志才教授工作组提供的2021年玛多地震的GNSS同震位移数据；中国地震局地质研究所InSAR工作组提供的2021年玛多地震的断层滑动模型；中国地质科学院地质研究所应急考察队提供的野外地质考察资料.