2022年1月8日MS6.9青海门源地震震前三维地壳变形与应变分配

李长军，郝明，2*，李煜航，宋尚武

1 中国地震局第二监测中心, 西安 710054 2 南方海洋科学与工程广东省实验室(珠海), 广东珠海 519080

0 引言

2022年01月08日01时45分在青海海北州门源县(37.77°N，101.26°E)发生MS6.9地震，震源深度10 km(图1).初步研究认为该地震发生于祁连山断裂带中东段，同时破裂了托莱山断裂东段和冷龙岭断裂西段，最大滑动量达到了3.5 m(李振洪等，2022)，USGS给出的震源机制解表现为走滑型地震(https:∥earthquake.usgs.gov/earthquakes/)，矩震级为MW6.6.在本次地震以东～50 km的冷龙岭断裂北侧的分支断裂上，曾于1986和2016年分别发生了两次MS6.4逆冲型强震(图1；胡朝忠等，2016；姜文亮等，2017；李彦宝等，2017)，呈现出该地区较强的地震活动性.

震源区发育有托莱山断裂、冷龙岭断裂、民乐—大马营断裂和皇城—双塔断裂等，共同构成了祁连山断裂带中东段的复杂活动断裂体系(图1).其中，托莱山断裂和冷龙岭断裂是祁连—海原大型活动断裂带的一部分，其左旋滑动速率达到了2～19 mm·a-1，兼具一定的逆冲分量(葛伟鹏等，2013；Guo et al., 2019；Li et al., 2009, 2018；何文贵等，2010；Duvall and Clark, 2010; Loveless and Meade, 2011; Yuan et al., 2013; Gaudemer et al., 1995).民乐—大马营断裂为晚更新世活动的逆断裂，与皇城—双塔断裂西段以左旋阶区方式相接，其垂直滑动速率为0.91±0.09 mm·a-1(雷惊昊等，2017)、0.25～0.28 mm·a-1和0.14～0.24 mm·a-1(Liu et al., 2021)，该断裂目前尚无历史地震证据(图1；李彦宝等，2017).皇城—双塔断裂以逆冲运动为主，兼有一定的左旋走滑运动，是全新世活动断裂，其垂直滑动速率为0.54～0.8 mm·a-1(陈文彬，2003)，被认为是1927古浪M8地震的发震断裂(图1；侯康明，1998).

图1 青藏高原东北缘主要活动断裂和历史地震分布

祁连山断裂带中东段位于河西走廊，是“一带一路”的重要纽带，也是诸多国家重大工程的必经之路.因此，开展该区域的现今三维地壳变形特征及其强震孕育背景的研究具有非常重要的意义.虽然，基于大地测量资料开展青藏高原东北缘祁连山断裂带及其周边区域的地壳变形特征及主要活动断裂的滑动速率，前人已开展了诸多研究.然而，上述研究主要集中在包括阿尔金断裂和海原断裂带等在内的大尺度主边界断裂，研究内容主要聚焦于大尺度地壳变形模式(Gan et al., 2007; Li et al., 2018; Loveless and Meade, 2011；葛伟鹏等，2013；Duan et al.2022)、“天祝地震空区”的蠕滑变形特征(Jolivet et al., 2012; Li et al., 2021a)、发生过1920海原地震的狭义海原断裂和存在较高地震危险性的六盘山断裂等(M7专项工作组, 2012；郝明等，2017；杜方等，2018；Su et al., 2019)，而对祁连山断裂带中东段冷龙岭和托莱山段的三维地壳变形特征和强震孕育背景则鲜有研究.

本研究首先搜集整理了祁连山断裂带中东段的GPS和水准等大地测量观测资料，通过精密的数据处理，构建了该区域的高精度三维地壳运动速度场；其次，通过计算应变率场和厘定主要活动断裂的滑动速率，分析了该区域在2022年MS6.9门源地震之前的精细三维地壳变形特征及其变形模式，探讨了该区域的强震孕育背景.

1 数据和方法

1.1 GPS数据及其处理方法

本研究使用的GPS数据主要来自中国地壳运动观测网络项目于1999—2019年观测，以及中国大陆构造环境监测网络项目于2009—2019年观测.首先利用 GAMIT 软件获得GPS测站坐标、卫星轨道及台站对流层天顶延迟等参数及其方差-协方差矩阵的单日松弛解.为保证定位结果的自洽性，采用相同的模型和方法对同时段的全球70个IGS核心站数据进行了处理.然后利用GLOBK软件将区域站和全球IGS站的单日松弛解捆绑在一起进行网平差，通过7参数转换方法将区域站坐标固定在ITRF2014框架下，得到区域站的GPS位移时间序列.由于研究区受到2001年昆仑山口西8.1级大地震的同震及震后变形影响，根据Wang和Shen(2020)的处理策略，只采用2009年之后的GPS时间序列.基于最小二乘方法拟合区域站的GPS线性速率，然后将其转换至相对华南地块(图2).

1.2 水准数据及其处理方法

本研究收集整理了1970—2014年穿过震中的一条长约330 km的水准测线(图2)，均按照国家一等水准测量规范开展了3—8期重复观测.利用水准数据获取地壳垂直运动速率的方法都是基于线性速率模型，即地壳的垂直运动速率是恒定的，得到的结果也只能是几十年时间尺度上的一种总的趋势或平均运动的定量描述.因此，我们以相邻水准点之间的测段高差作为观测值，基于动态线性平差模型，以河西走廊内的水准点(测线最北端)为参考基准，估计出各水准点相对河西走廊的垂直运动速率.平差计算得到的验后单位权中误差为0.96 mm·km-1，也就是说采用线性速率模型，平差后得到的一等水准测量每公里高差中误差为0.96 mm，这与我国精密水准测量规范规定的一等水准测量每公里高差测量全中误差1.0 mm一致，说明线性速率模型的模型误差较小，估计出的垂直运动速度是可靠的.

1.3 水平应变率场

为了获取可靠的水平应变率场，本文使用了基于中位值滤波的MELD(Median Estimation of Local Deformation)应变率场稳健计算方法(Kreemer et al., 2020, 2018).在球面坐标系下，速度场和应变率场之间的关系可以表示为

(1)

(2)

(3)

(4)

2 结果与讨论

2.1 水平形变场

相对于稳定华南地块的GPS速度场显示，以祁连山断裂带(包括托莱山断裂、冷龙岭断裂等)为界，其南侧的GPS速度矢量为NE方向，而其北侧的GPS速度矢量则偏转至N或NW方向且速率大小由～10 mm·a-1减小到～4 mm·a-1(图2)，表明祁连山断裂带地壳存在明显的左旋走滑和缩短变形，与该区域的活动构造背景一致(李彦宝等，2017；Duvall and Clark, 2010; Yuan et al., 2013; Zuza and Yin, 2016)，反映了在印度板块的N向推挤和阿拉善块体的阻挡作用下，青藏高原东北缘的地壳变形被高原地壳的整体E向挤出和祁连山构造带的挤压逆冲所转换吸收(Li et al., 2018, 2021b; Zheng et al., 2013).水平应变率场表明，祁连山断裂带地壳变形以NE-SW向的挤压和NW-SE向的拉张为主，呈现出较强的剪切和挤压特征(图3—5)，主压应变方向与该区域的主压构造应力方向近乎平行(图3b；Heidbach et al., 2010；周琳等，2022).其中，沿祁连山断裂带的最大剪应变率为(15～20)×10-9/a(图3a)；面应变率为(-10～-20)×10-9/a(图3b)；旋转应变率为(5～10)×10-9/a(图3c)，表明该区域地壳整体呈现逆时针旋转运动，与前人结果基本一致(Li et al., 2018, 2021b; Wang and Shen, 2020).位于祁连山断裂带南、北两侧的柴达木盆地和阿拉善块体的变形则明显弱于祁连山断裂带.此外，本文的应变率场揭示出，此次门源地震发生在祁连山断裂带中东段的最大剪应变高值区(图3a)、两个压性面膨胀高值区的交汇边缘位置(图3b)、以及旋转率高值区边缘(图3c)；而2016MS6.4逆冲型地震震源区的剪应变率相较于此次地震略低(图3a和4a)，面应变率和旋转率则略高(图3b和3c).为了更进一步探讨发震前震源区的地壳变形特征，本文分别计算了1999—2016和2016—2019年的应变率场(图4和图5).其中，1999—2016年的应变率场结果表明，2016MS6.4地震发生于剪应变率高值区边缘(图4a)、面应变率高值区(图4b)；2016—2019年的应变率场揭示出2016MS6.4地震之后，面应变率高值区明显向西迁移(图5b)，表明2016MS6.4地震释放了部分累积的压性面应变，此次地震则发生在剪应变率高值区和面应变率高值区边缘(图3和图5).在此需要指出的是，虽然计算2016—2019年的应变率场所用的GPS观测数据较少而导致应变率场的整体有效空间分辨率较低，但是，在此次地震震源区较高的空间分辨率有助于更好地判定震前地壳变形特征.

图3 1999—2019年的应变率场

图4 同图3，除了应变率场的时间范围为1999—2016年

图5 同图3，除了应变率场的时间范围为2016—2019年

2.2 断裂滑动速率及三维地壳变形模式

为了准确厘定祁连山断裂带中东段主要活动断裂(包括托莱山断裂、冷龙岭断裂、皇城—双塔断裂和民乐—大马营断裂等)的滑动速率，并分析该区域的应变分配特征，横跨祁连山断裂带提取2条GPS速率剖面(图6)和1条水准速率剖面(图7)，剖面位置如图2所示.在提取GPS速率剖面时，将GPS速度矢量分解到平行于主要活动断裂(在此为托莱山断裂和冷龙岭断裂)方向(图6a和6c)和垂直于主要活动断裂方向(图6b和6d).图6能够清楚地反映出断裂两侧的GPS站点分布及计算断裂滑动速率所用的GPS点数.由于断裂之间的观测点稀少，在此首先通过计算整个断裂带两侧的GPS速度分量的均值的差值，以获取整个祁连山断裂带的总体走滑和缩短速率.其中，托莱山段的总走滑速率为3.6±0.5 mm·a-1，缩短速率为2.1±0.5 mm·a-1，基本与以GPS为约束反演得到的2.5～5.8 mm·a-1的左旋和3.7 mm·a-1的缩短速率，以及活动构造研究给出的3.5～5.5 mm·a-1的缩短速率一致(葛伟鹏等, 2013; Li et al., 2018; Hu et al., 2021).此外，基于改进的断裂位错模型(邹镇宇等，2015)和图6a所示的GPS剖面数据，本文利用Levenberg-Marquardt方法反演了断层的走滑速率和闭锁深度.反演结果中，拖莱山断裂的左旋走滑速率为4.5±0.5 mm·a-1，闭锁深度为15.0±7.8 km(图6a).虽然，位于托莱山断裂和民乐—大马营断裂之间仅有一个GPS点(图6a和6b)，但该点能够有效地用于研究该区域的应变分配.通过计算断裂两侧GPS速度分量的差值，得到托莱山断裂的走滑和缩短速率分别为2.5±0.3 mm·a-1和1.3±0.4 mm·a-1，与Hu等(2021)的研究指出托莱山断裂存在1～2 mm·a-1的缩短率一致；民乐—大马营断裂的走滑和缩短速率分别为1.1±0.3 mm·a-1和0.8±0.3 mm·a-1，与地质上给出的0.24～0.36 mm·a-1的缩短基本一致(Liu et al., 2021)，远小于李彦宝等(2017)的5.9±0.8 mm·a-1.活动构造研究给出的冷龙岭断裂晚第四纪以来的滑动速率为2～19 mm·a-1(Lasserre et al., 1999；何文贵等，2010；Li et al., 2009; Guo et al., 2019).限于冷龙岭断裂北侧缺少GPS观测，本文并未利用该剖面的观测数据进行反演，只是简单地利用断裂带两侧的GPS速度分量均值的差值计算得到冷龙岭和皇城—双塔断裂的总体走滑速率为3.1±0.7 mm·a-1，缩短速率为3.0±0.6 mm·a-1，与Zheng等(2013)的3.9 mm·a-1一致.尽管位于冷龙岭断裂北侧的皇城—双塔断裂也属于全新世活动断裂，且存在0.5～0.8 mm·a-1的垂直滑动速率，但该断裂以逆冲运动为主(陈文彬，2003)，因此，GPS剖面揭示出的～3 mm·a-1的左旋和缩短速率主要反映了冷龙岭断裂的现今滑动速率.

图6 跨祁连山断裂带中东段主要活断层的GPS剖面，取3σ置信水平.剖面位置见图2

横跨祁连山断裂带中东段的水准垂直速率表明，在印度板块的N向推挤和阿拉善块体的阻挡作用下，位于青藏高原东北缘扩展变形前缘的祁连山断裂带呈现出整体隆升.其中，以托莱山断裂为界，其南侧的隆升速率为1.0～1.5 mm·a-1，略高于北侧的～0.5 mm·a-1(图7).基于断裂两侧的差异运动，计算得到托莱山断裂的垂直滑动速率为1.2±0.6 mm·a-1，与Hu等(2021)利用阶地陡坎得到的1.5 ± 0.1 mm·a-1一致；民乐—大马营断裂的垂直滑动速率为0.5±0.5 mm·a-1，与雷惊昊等(2017)利用阶地陡坎计算得到的0.91±0.09 mm·a-1一致，略大于Liu等(2021)的0.25～0.28 mm·a-1和0.14～0.24 mm·a-1.结合改进的倾滑断裂位错模型(邹镇宇等，2015)，反演得到拖莱山断裂的倾滑速率为1.8±0.5 mm·a-1.此外，在托莱山断裂两侧和民乐—大马营断裂的北侧也存在一系列的水准速率拐点，速率差异<1 mm·a-1，可能反映了其他逆冲断裂(包括门源断裂、肃南—祁连断裂等)的滑动状态，受资料所限，本文在此不做详细讨论.

图7 跨祁连山中东段的水准观测剖面(圆点，取1σ置信水平)及地形(灰色实线)

跨断裂GPS速率剖面表明托莱山—冷龙岭断裂的水平缩短速率为2～3 mm·a-1，水准速率剖面揭示该断裂具有1.2±0.6 mm·a-1的垂直滑动速率.根据单轴挤压泊松横向变形和褶皱变形两种情况来考虑，其伴生垂直位移量与水平位移缩短量之间的比例最大为1∶2(极端紧密褶皱).因此祁连山断裂带中东段的抬升速率主要以区域北东向地壳挤压缩短的形式实现.同时，GPS应变率场、三维速率剖面以及断裂缩短和垂直速率均揭示出～300 km宽的祁连山断裂带中东段的变形分配比较均匀，主断裂和次级断裂的滑动速率和变形方式基本一致，均表现为左旋走滑兼具逆冲，且滑动速率未出现明显差异，表明该区域的地壳变形可能满足连续—弥散变形模式，热动力学成因可能来自于较高温度(塑性)的中—下地壳，使脆性上地壳断裂与其周围完整岩石的强度差变小，从而更容易产生破裂，孕育地震(孙玉军等, 2013; Hu et al., 2021; Li et al., 2020; Yang et al., 2020；Wang et al., 2022).

2.3 对地震危险性的指示意义

“弹性回跳理论”和负位错理论阐明了断层震间应变积累及地震孕育和发生过程(Reid, 1910; Savage et al., 1983).Scholz(1998)结合断层摩擦定律，将其进行了延伸拓展，给出了断层面不同深度的运动随时间和深度的变化模式：即断层经过震后愈合之后进入下一次地震的震间期，其中，在震间早期的变形集中在断层附近很狭窄的范围内，而到了震间晚期断层两侧的变形则趋于平缓且变形范围逐渐扩大，而变形离开断层距离的衰减速度与断层闭锁深度有关(Savage, 1983).具体表现为，距离下一次地震时间越近，则断层附近所能观测到的相对运动与变形就越小(江在森和武艳强，2012).

本研究的GPS剪应变率场揭示出2022年MS6.9门源地震发生在与构造变形背景相一致的剪应变率高值区(图3a和5a)，但未发生在面应变率的压性高值区而是在边缘区(图3b和5b)，而2016MS6.4门源地震则发生在面应变率高值区，进一步验证了江在森和武艳强(2012)基于地壳变形的异常信息识别强震地点的判据.跨断裂GPS速率剖面表明，GPS速率分量随着到断层的距离越近而明显减小(图6)，表明2022年门源地震前，该地区处于剪切和挤压应变积累状态.水准垂直速率剖面表明，托莱山和冷龙岭断裂以南地壳垂直变形并不明显(图7)，同样揭示出门源地震震前应变在持续增加.因此，震前三维地壳变形反映出2022年MS6.9门源地震之前，托莱山—冷龙岭断裂带应力-应变积累处于较高水平，具有发生强震的背景.

1986年和2016年发生在冷龙岭断裂北侧分支断裂的两次MS6.4地震(图1；姜文亮等，2017)都表现为逆断特征，与2022年MS6.9地震的左旋错动截然不同.尽管区域三维地壳运动速度场揭示出托莱山和冷龙岭断裂兼具左旋走滑和逆冲倾滑的构造背景(图6和图7)，但此次MS6.9地震发生在这两条左旋走滑断裂的拉张阶区内，所以发震断层总体表现为走滑错动，而逆冲倾滑分量并不明显.

3 结论

本文收集整理了2022年1月8日MS6.9青海门源地震之前祁连山断裂带中东段的GPS和水准等大地测量资料，计算了应变率场和主要活动断裂的滑动速率，探讨分析了此次门源地震震间晚期的三维地壳变形特征及其强震孕育背景.主要结论包括：

(1)祁连山断裂带地壳变形以NE-SW向的挤压和NW-SE向的拉张为主，呈现出较强的剪切和挤压特征，其中，最大剪应变率为(15～20)×10-9/a、面应变率为(-10～-20)×10-9/a、旋转应变率为(5～10)×10-9/a，而位于祁连山断裂带南、北两侧的柴达木盆地和阿拉善块体的变形则明显弱于祁连山断裂带.

(2)祁连山断裂带托莱山段的总走滑速率为3.6±0.5 mm·a-1，缩短速率为2.1±0.5 mm·a-1，其中，托莱山断裂的走滑和缩短速率分别为2.5±0.3 mm·a-1和1.3±0.4 mm·a-1、垂直滑动速率为1.2±0.6 mm·a-1，民乐—大马营断裂的走滑和缩短速率分别为1.1±0.3 mm·a-1和0.8±0.3 mm·a-1、垂直滑动速率为0.5±0.5 mm·a-1；冷龙岭段的走滑速率为3.1±0.7 mm·a-1，缩短速率为3.0±0.6 mm·a-1.

(3)2022年门源MS6.9地震发生在剪应变率高值区、两个压性面膨胀高值区的交汇边缘位置、以及旋转率高值区边缘.在此次地震前，该地区处于剪切和挤压应变积累状态，而托莱山和冷龙岭断裂以南地壳垂直变形并不明显，揭示出门源地震震前应变在持续增加，同时，基于震间GPS剖面反演得到拖莱山断裂的闭锁深度为15.0±7.8 km，表明该断裂具备发生强震的背景.

(4)2016MS6.4门源地震发生于面应变率高值区，该地震之后面应变率高值区明显向西迁移，表明2016MS6.4地震释放了部分累积的压性面应变.

致谢中国地震台网中心为本研究提供了GPS观测数据，本文使用的画图软件为GMT，作者在此一并表示感谢.同时感谢所有参与GPS和水准观测的工作人员.