0 区域构造与历史地震

据中国地震台网中心正式测定，2022年1月8日1时45分27秒青海海北藏族自治州门源县发生M_S6.9地震，震中位于(37.77°N，101.26°E)，震源深度10km ^①。此次地震发生在青藏高原东北缘的昌马堡—古浪—海原构造带(李振洪等，2022)，该区域经受着强烈的地壳缩短进而左旋剪切，构造复杂且活动强烈(张培震等，2003)。宏观震中位于冷龙岭断裂西段，中国地震台网中心、中国地震局地球物理研究所和美国地质调查局(USGS)等多家单位给出的震源机制解均显示，此次地震为一次走滑型为主的破裂，与冷龙岭断裂的构造运动背景一致(图 1)。发震断裂长约120km，西接托莱山断裂，东连金强河和天桥沟—黄羊川断裂断裂，与托莱山、金强河、毛毛山、老虎山和海原断裂等一起组成祁连—海原断裂(陈文彬，2003；Zheng et al，2013)，晚第四纪时期主要表现为左旋走滑运动，局部兼具倾滑分量(Gaudemer et al，1995；何文贵等，2000、2010；Lasserre et al，2002)。

① https://www.cenc.ac.cn/cenc/dzxx/396391/index.html

② https://www.globalcmt.org/CMTsearch.html

冷龙岭断裂带所在区域是中国大陆地壳运动最强烈、地震活动频度最高、强度最大的地区之一(李智敏等，2022)。据统计，1900年以来震中100km范围内共发生过6.0级以上地震6次，6.0~6.9级地震5次，8.0~8.9级地震1次(1927年5月23日甘肃古浪8.0级地震)。时间最近的一次强震为2016年青海门源M_S6.4地震，空间距离约为33km，震源机制显示为一次逆冲型破裂。研究表明，广义海原断裂带存在两个明显的挤压弯曲带(或挤压阶区)，分别为“海原挤压弯曲带”和“冷龙岭挤压弯曲带”，类似区域易于聚集应力形成孕震构造(Lozos et al，2011；Nabavi et al，2018；Wang et al，2020b)，2016年门源逆冲型地震加速了冷龙岭断裂的应力积累，对2022年门源地震的发生有一定的促进作用(潘家伟等，2022)。库仑应力触发的相关研究结果也显示，2022年门源地震的发震断层处在2016年门源地震同震引起的库仑应力加载区，加之该区域高闭锁和应变积累较强的背景(赵静等，2016；石富强等，2018)，二者具有显著的触发作用(李振洪等，2022；朱琳等，2022)。

震后野外考察显示出此次地震具有南北两条地表破裂带，主破裂位于冷龙岭断裂西段，南支的次级破裂向西延伸到了托莱山断裂东段局部段上(李智敏等，2022；潘家伟等，2022；Yang et al，2022)。分析研究两大断裂带震间期的深、浅部构造变形特征，对于深入理解该区域孕震过程，以及结合同震库仑应力变化评估断裂带未来地震危险性，均具有重要的理论和实际意义。

1 研究区域的运动学特征 1.1 速度场空间展布

图 2为1991—2016期相对于稳定的欧亚板块的中国大陆GNSS速度场空间分布(Wang et al，2020a)，以此作为背景场，分析研究区水平向地壳变形特征。从图 2可以看出，在受印度板块北东向推挤和阿拉善、鄂尔多斯块体阻挡的共同作用下，青藏高原东北缘的上地壳水平运动在跨过一系列北西向展布的断裂带后锐减，并发生顺时针旋转；到达稳定的阿拉善和鄂尔多斯块体后，北东向推挤速率减小到原来的1/3左右，沿主要断裂带形成速率变化的高梯度带，从而导致海原断裂带呈现出左旋走滑为主、兼具挤压逆冲的构造运动特征，并在断裂附近积累弹性应变能(李彦川，2016)。1976年以来M_W4.0以上地震的震源机制解显示，研究区域以逆冲和走滑型地震为主，与GNSS速度场显示的左旋挤压为主的变形模式基本一致。

板块内部各次级块体在周围块体的作用下发生刚性旋转和弹塑性变形，因此，扣除欧拉运动模型和速度场拟合的区域整体旋转分量后，残余速度分量的空间展布能够凸显局部精细化形变特征，从而反映不同构造段落的差异化运动特征，计算方法如下

$ \left[\begin{array}{l} V_e \\ V_n \end{array}\right]=r\left[\begin{array}{ccc} -\sin \phi \cos \lambda & -\sin \phi \sin \lambda & \cos \phi \\ \sin \lambda & -\cos \lambda & 0 \end{array}\right]\left[\begin{array}{c} \omega_x \\ \omega_y \\ \omega_z \end{array}\right] $

(1)

其中，(λ，ω)表示块体内部各点的经纬度坐标，r表示研究块体的平均曲率半径，(V_e，V_n)表示东、北向速度矢量，(ω_x，ω_y，ω_z)表示块体运动的欧拉矢量。从扣除区域整体旋转运动后的速度场空间展布来看(图 2)，青藏高原东北缘不同断裂带的运动学状态存在显著差异。

祁连西段南侧的青藏高原内部地区呈现出逆时针旋转为主的运动特征，伴随一定的地壳缩短变形；托莱山至毛毛山、老虎山一带的青藏高原内部地区以垂直于断层的地壳缩短运动为主，挤压变形突出；东部六盘山地区存在显著的顺时针旋转运动。阿拉善块体整体以顺时针旋转运动为主，断层附近各点运动方向与断层走向近乎平行，挤压运动不明显。青藏高原北东向推挤受到稳定的阿拉善和鄂尔多斯块体阻挡后，表现出地壳缩短和侧向扩展为主的变形特征，与前人基于其他数据得到的研究结果类似(葛伟鹏等，2013)。

1.2 跨断层速度剖面

选取门源M_S6.9地震震中附近的断裂两侧一定区间范围内的GNSS台站，将各台站的速度矢量在断裂上投影，得到平行和垂直于断裂走向的速度分量随距离断裂的分布，剖面位置见图 2中黑色矩形框所示，平行和垂直速度的剖面投影结果如图 2(b)、2(c)所示。平行断层的速度投影显示，震中附近的冷龙岭和托莱山断裂左旋运动显著，变形宽度约为58km(变形总量的80%)，具备反正切函数的分布特征，表明上述断裂存在明显的应变积累现象；垂直断层的速度投影显示，断裂带两侧伴随近似均匀的地壳缩短变形。

2 GNSS应变场及其动态演化特征

强震孕育发生过程始终伴随着地壳应变的积累释放现象，研究GNSS应变场的空间分布、时间演化特征，对于强震地点预测具有重要意义；另一方面，GNSS速度场能够直观表达不同地区地壳运动的构造差异特征，但易受参考框架的影响，而应变分布不受参考框架影响，能够全面表达地壳变形的性质和强度，更有利于震前变形特征识别和异常演化过程研究(中国地震局监测预报司，2020；江在森等，2020)。常用的计算应变场的方法包括局部法和整体法，局部法是利用应变公式获取块体或网格单元的局部应变场，进而将局部应变场合并得到整体应变场，常用方法有REHSM和RELSM等(李延兴等，2001)；整体法则通过建立位移与位置的函数关系，进而对位置求导获得应变场分布，常用方法有多面函数法(杨国华等，2002)、球谐函数法(石耀霖等，2006)、最小二乘配置(江在森等，2010)、球面小波方法(Tape et al，2009)等。

同时，为探索GNSS观测资料在中强地震预测中的应用，很多研究者应用静态或动态GNSS应变参数，论证了其时空分布和演化特征与中强地震发生时间、地点对应关系。例如王启欣等(2020)总结了汶川M_S8.0、玉树M_S7.1、芦山M_S7.0和岷县—漳县M_S6.6地震周边剪应变的空间分布特征，发现上述地震均发生在剪应变高低值转换区域；江在森等(2013)研究表明，云南地区的M≥5.0地震多为走滑型地震，地震的破裂类型多含有剪切破裂，震前应变积累包含一定的剪切应变积累；Wu等(2015)总结了中国大陆地区7级以上强震与GNSS应变率分布的关系，指出倾滑型地震多发生在应变积累高值区边缘或过渡区与倾滑断层交会区，走滑型地震多发生在最大剪应变高值区边缘与大型走滑断裂带的交会地带；Wu等(2022)利用1999—2007期速度场和最小二乘配置算法计算的“变形总强度”，与1990—2017年期间三个时段的地震活动累积率作统计分析发现，二者对应率良好，显著优于自然概率，即“变形总强度”与地震发生地点具有较好的相关性。

本文研究区域GNSS站点分布不均匀，最小二乘配置方法通过建立区域水平速度经验协方差函数，采用位移-应变微分关系获得断裂带附近的应变率场，能够较大限度地削弱应变空间量值的大小受观测点分布密度的影响。因此采取上述方法，基于1991—2016期(Wang et al，2020b)和2017—2020期中国大陆及周边GNSS速度场，在扣除区域欧拉旋转分量后，计算了门源M_S6.9地震前GNSS面膨胀率、最大剪应变率的空间分布，并对比分析了其时间演化特征。

整个祁连—海原断裂带呈现出压剪性的应变特征，主应变以北东向的挤压和北西向的拉张为主，图 3面膨胀率计算结果表明：①以托莱山、肃南—祁连断裂一带为中心，形成挤压高值区，其西南侧的巴颜喀拉东北边界呈弱张性变化特征；②以1991—2016期的面膨胀率作为背景场，2017—2020期震源区附近压性应变略有增强，由1.51×10^-8/a增加到1.55×10^-8/a，挤压应变高值中心沿祁连带向震中方向迁移。

图 4最大剪应变率计算结果表明：①高值区与块体边界具有较好的吻合性，并沿区域内两大型走滑断裂祁连—海原和东昆仑断裂带分布；②对比两期计算结果，震源区附近最大剪应变率变化不大，其东西两侧断裂带，特别是海原断裂带东段有所增强。

综上，以1991—2016期GNSS应变场作为背景场，2017—2020期应变场的空间分布演化特征显示，震中附近应变特征变化不大，表明断层处于孕震晚期，变形较为微弱；同时，此次地震处在最大剪应变高值区边缘与大型走滑断裂带的交会地带，与前人关于强震发生地点的认识一致。

3 断裂带滑动速率与闭锁深度

断层的运动学参数对深入理解构造运动和孕震过程具有重要意义，也能为地震危险性评估提供定量依据(Burchfiel et al，2008；Zhang et al，2010；Diao et al，2019；朱亚戈等，2021)。作为运动学参数求解的常用方法之一，GNSS速度场的跨断层段剖面分析广泛应用于断裂带形变特征的定量研究(Le Beon et al，2008；张培震等，2008；江在森等，2009；武艳强等，2012；魏文薪等，2012)。利用基于位错理论的经典走滑位移公式拟合震间形变走滑分量，可以得到断层的滑动速率和闭锁深度信息(Savage et al，1973)。由于经典走滑公式通常不考虑断层倾角对地表的影响，邹镇宇等(2015、2018)在经典走滑位移公式的基础上进行了改进，引入断层倾角信息，使其能够反映断层倾角对地表形变的影响，并将其应用到汶川地震前后川滇块体东边界主要走滑断裂形变特征的剖面分析研究中。基于GNSS观测对海原断裂带滑动速率的相关研究表明，从西向东呈现出逐步减弱的趋势，平均滑动速率在5~9mm/a，毛毛山—老虎山段约为8mm/a，中、西段速率为5~7.6mm/a，东段为2.9~5mm/a(胡亚轩等，2009；葛伟鹏等，2013)

为研究门源M_S6.9地震周边断层的震前运动状态，本研究选取震源区周边跨托莱山和冷龙岭断裂的两个剖面，剖面位置详见图 1。基于剖面投影法和1991—2016期GNSS速度场，得到平行于断层的水平速度投影分量，进而利用改进走滑位移公式计算两个断层的长期滑动速率和闭锁深度，具体计算方法如下

$ V_{\text {Paralell }}=V_0+\frac{V}{\pi} \tan ^{-1}\left(\frac{x-d_{\text {offset }}}{D}\right) $

(2)

式中，V₀表示断层两侧块体的共同运动速度，V_Parallel表示平行于断层的GNSS水平速度投影分量，V表示断层的长期滑动速率，x表示垂直于断层走向的距离，d_offset表示实际断层的闭锁与滑动分界点相对于断层出露地表点的偏移量，D表示断层的闭锁深度。V₀、V、d_offset、D为拟合参数，通常采用网格搜索或非线性拟合算法求解，本文采用后者进行各参数的反演得到参数最佳估值。

如图 5所示，红色曲线表示模型给出的最佳拟合速度，灰色条带表示不确定度区间(加减一倍均方差)，每个点表示带有误差棒的平行于断层走向的水平速度分量。结果显示两个剖面均满足反正切的震间变形特征，存在不同程度的闭锁现象。长期滑动速率方面，托莱山和冷龙岭断裂带分别为(6.0±0.7)mm/a、(3.8±0.5)mm/a，与前人(何文贵等，2000)研究基本一致；闭锁深度方面，托莱山段远大于冷龙岭段，二者分别为46km和9km；另外托莱山段的变形宽度远大于冷龙岭段，更加符合震间变形期的孕震晚期阶段特征。

4 讨论与分析

本文以1991—2016期中国大陆GNSS速度场作为背景场，扣除区域刚性旋转后分析了门源M_S6.9地震及其周边区域精细化的形变特征，结合2017—2020期速度场分析了应变率场的时空演化特征，同时定量分析了震中附近托莱山和冷龙岭断裂的震间变形特征。取得的认识和存在问题如下：

(1) 研究区内不同段落的运动学特征差别较大，青藏高原内部以旋转和地壳缩短变形为主，托莱山至毛毛山、老虎山一带地壳缩短变形显著，以北区域具有逆时针旋转特征，以南区域具有顺时针旋转特征，阿拉善块体整体以顺时针旋转运动为主，区域变形受控于周边大型块体间的推挤与阻挡，呈现出局部旋转和整体弥散的特征，与前人利用其他数据得到的认识基本一致。

(2) 此次门源M_S6.9地震处于面膨胀率和最大剪应变率高值区的边缘，与前人关于强震地点与应变率场空间分布对应关系的认识一致；相比于1991—2016期背景场，近年来震源区附近压应变略有增强，剪应变持续弱化，存在闭锁可能；震中附近应变特征整体变化不大，断层变形较为微弱，符合孕震晚期的变化特征。

(3) 跨冷龙岭断裂的速度剖面中近断层点位稀疏，缺少数据约束会对反演结果产生一定影响。目前两个剖面投影和走滑速率最佳拟合结果显示，托莱山和冷龙岭断裂均存在不同程度的闭锁现象，托莱山的滑动速率、变形宽度及闭锁深度均大于冷龙岭断裂，更加符合孕震晚期的变形特征；结合朱琳等(2022)关于青藏高原北部1900年以来强震对托莱山断裂带库仑应力加载的研究结果，托莱山断裂未来的地震危险性仍值得关注。

致谢: 感谢中国地震局第一监测中心提供的GNSS速度场，感谢武艳强研究员提供的最小二乘应变率计算软件包，感谢匿名审稿专家对本文提出的宝贵修改意见。