青藏高原东北缘构造应力场及动力学特征

孟文, 郭祥云, 李永华, 韩立波, 张重远

1 中国地震局地球物理研究所, 北京 100081 2 中国地震局震源物理重点实验室, 北京 100081 3 中国地质科学院地质力学研究所, 北京 100081 4 自然资源部活动构造与地质安全重点实验室, 北京 100081

0 引言

新生代以来，印度板块与欧亚板块的碰撞和持续汇聚作用导致喜马拉雅造山带和青藏高原的快速崛起与不断向外扩展(Molnar and Tapponnier, 1975).前人围绕青藏高原深浅地质过程和变形模式提出了多种动力学模型，主要包括大陆侧向逃逸模型(Tapponnier et al., 1982, 2001)、地壳连续变形模型(England and Houseman, 1986; Dewey et al., 1988)和下地壳流模型(Royden et al., 1997, 2008; Clark and Royden, 2000).侧向逃逸模型认为青藏高原内部块体沿深大断裂发生侧向挤出，其变形主要局限在块体的边界断裂带上；地壳连续变形模型以持续的挤压动力、均匀的黏性物质和垂直连贯变形为特征，认为地表的隆升主要依靠地壳均匀的黏性增厚；下地壳流模型以上地壳高强度和无明显缩短，以及上下地壳发生解耦为特征.各种变形机制和动力学假说都获得了来自地震地质、地球物理与形变测量等学科的证据(Shen et al., 2001a; Chen et al., 2004; 王新胜等, 2013; Wang et al., 2016; Liu et al., 2014; Sun et al., 2012)，有关青藏高原隆升和动力学机制的争论仍在继续.实际上，整个青藏高原大陆变形很难用某一种变形模型去解释，不同地区或块体可能以某一种变形方式为主，同时兼具其他变形方式(Replumaz and Tapponnier, 2003; 许志琴等, 2010; Liu et al., 2014; 李海兵等, 2021)，高原形成机制与构造模式以及变形方式主次仍存在争议的一个重要原因就是应力应变状态没有得到很好的限定(潘正洋等, 2020).

青藏高原东北缘是青藏高原扩张的前沿地带.青藏高原快速隆升并向北东方向推挤，在其东北缘受到阿拉善块体和鄂尔多斯块体的阻挡，导致区内地壳发生缩短、抬升和左旋剪切，构造变形强烈且复杂(黄汲清等, 1980; 臧绍先等, 1992; 邓起东等, 2002; 张培震等, 2003)，发育北西西向和北北西向两组主导性活动构造，包括以左旋走滑为主的阿尔金断裂带和海原断裂带，以逆冲-左旋走滑为主的祁连山北缘断裂带、西秦岭北缘断裂带和东昆仑断裂带，以左旋逆冲为主的六盘山断裂，以及以右旋走滑为主的日月山断裂带和鄂拉山断裂带等边界断裂带(图1，Yin and Harrison, 2000; Tapponnier et al., 2001; 虢顺民等, 2000; 张培震等, 2013; 邓起东等, 2002).这些断裂直接制约和影响了东北缘晚第四纪活动构造变形，形成了以挤压推覆构造、次级剪切构造、剪切压扁构造和弧形挤出构造为主的4种构造变形方式，且由活动断裂控制的菱形块体旋转运动明显(Duvall and Clark, 2010; 虢顺民等, 2000; 袁道阳等, 2004)，形成了青藏高原东北缘强烈而独特的构造环境，是我国西部挤压构造向东部拉张构造转换的过渡区域(Lasserre et al., 1999)，也是新生代构造变形最为强烈的地区之一.历史上发生过1920年海原8.5级地震、1927年古浪8.0级地震、1932年祁连山北西端昌马7.6级地震和2001年昆仑山口西8.1级地震.

地壳构造变形和地震的孕育均与地壳应力状态密切相关，地震是在区域构造应力场作用下，岩体沿某一构造面发生破裂的结果(Scholz, 2002)，因而地壳构造应力场研究能够帮助深入理解青藏高原东北缘的构造变形及孕震环境.地壳构造应力场的研究方法诸多，除水压致裂法、应力解除法等可直接获得应力状态的原地应力测量方法外，通过地质学、地球物理学等相关方法反演也可提取应力信息(Zoback, 2007).地震震源机制不仅可以确定地震断层类型，也反映了区域构造应力状态，是研究区域构造应力场，特别是地壳深部应力状态的有效方法之一(Gephart and Forsyth, 1984; Hardebeck and Michael, 2006; Wan et al., 2016； Li et al., 2018; Tian et al., 2019).地震各向异性与应力应变场密切相关，也是揭示岩石圈变形最直接的有效手段(Boness and Zoback, 2006; Gao et al., 2011; Crampin and Peacock, 2005; 高原等, 2018).利用近震剪切波分裂开展的上地壳各向异研究(张辉等, 2012; 郭桂红等, 2015; Hu et al., 2020)认为青藏高原东北缘块体内部地壳各向异性与研究区主压应力有关，而深大断裂带附近的各向异性则与断裂活动变形有关.部分学者基于XKS波分裂研究认为青藏高原东北缘壳幔变形具有强烈的耦合性，呈现垂直连贯变形的特征(常利军等, 2016; Chang et al., 2017; Wang et al., 2016).也有学者利用接收函数确定地壳各向异性特征的分析方法，认为地壳各向异性主要来源于中下地壳，且受控于构造应力场及下地壳介质逃逸作用，同时揭示了存在局部壳幔解耦的可能(邵若潼等, 2019).利用震源机制解反演应力场的研究同样认为东北缘部分区域壳幔以及上下地壳变形可能是解耦的(Han et al., 2019).前人通过地震各向异性或震源机制研究对青藏高原东北缘变形和构造应力场分布特征取得了一定认识，但关于研究区壳幔变形方式仍然存在争议，并且缺乏联合不同研究方法和研究成果的综合性分析.

本文利用CAP全波形反演方法及HASH初动极性和振幅比方法，系统反演了32°N—41°N，96°E—107°E区域内ML≥3.0地震震源机制解781组，同时收集了Global CMT给出的MW≥4.6中强地震震源机制解96组，并利用STASI(阻尼区域应力场反演)算法对该区的构造应力场分布进行分析，同时联合基于地应力实测获得的地壳浅表层地应力状态、基于地震各向异性获得的地壳和上地幔变形特征以及基于GPS观测获得的地壳浅表层形变特征，讨论了地壳浅表层、上地壳、中下地壳和上地幔的垂向变形特征，以期能够对青藏高原东北缘的变形特征和动力学机制获得较为深入的认识.

1 数据资料与计算方法

1.1 数据资料

本研究收集了研究区及周边103个宽频带固定台站2008年10月—2019年6月记录的ML≥3.0以上地震波形数据(台站分布见图1b).为了保证研究结果的可靠性，仅从中选择信噪比较高的地震波形记录进行分析.同时，本研究补充采用了全球矩心矩张量目录(GCMT)中1976年1月—2020年12月96组MW≥4.6地震事件的震源机制解数据，共同参与应力场反演.

图1 青藏高原东北缘及区域构造背景图(a) 青藏高原东北缘构造背景. 浅灰色圆形符号代表本文收集并进行反演的1976—2020年地震事件(ML≥3.0); 黑色曲线代表全新世活动断裂, 灰色曲线代表更新世活动断裂(邓起东等, 2002); 紫色虚线表示图12中地形剖面位置. (b) 青藏高原东北缘区域构造背景. 紫红色三角符号代表本研究使用的宽频带固定台站，浅紫色粗实线代表块体边界.Fig.1 Tectonic settings in and around the NE margin of the Tibetan Plateau(a) Geological structure in the NE margin of the Tibetan Plateau. Light gray circles denote local earthquakes with magnitudes equal to or greater than 3.0 occurred between 1976 and 2020 collected and inverted in this study; The black curves denote Holocene active faults, and the gray curves denote Pleistocene active faults (Deng et al., 2002); The purple dashed line indicates the position of the terrain profiles in Fig.12. (b) Regional tectonic settings in the NE margin of the Tibetan Plateau. The fuchsia triangular symbols represent 103 permanent seismic stations deployed in the NE margin of the Tibetan Plateau and the light purple thick solid lines represent block boundaries.

1.2 震源机制解计算

1.2.1 CAP全波形反演法及数据处理

本研究中ML≥3.5地震事件的震源机制解采用“裁剪-粘贴”(Cut-And-Paste，CAP)全波形反演方法(Zhao and Helmberger, 1994; Zhu and Helmberger, 1996).其主要思想是将宽频带地震记录分成Pnl波(包含P波及其部分后续波形)和面波波段两个部分进行反演.通过分别对Pnl和面波两部分波形赋予不同的权重，计算其理论合成波形和真实观测波形的互相关系数及误差函数，并运用网格搜索方法搜索出理论合成波形与真实观测波形差异最小时的地震震源机制解及最佳震源矩心深度.CAP方法具有对速度模型和台站分布依赖性较小、计算结果相对稳定等特点.

采用CAP方法反演之前，需要对实测波形数据和理论波形数据做如下准备工作：首先对原始台站波形数据进行预处理，主要进行去仪器响应并校正参考时间为发震时刻；而后将N-E-U坐标系下的三分量波形旋转至R-T-Z坐标系下，并将单位转换为cm·s-1；利用F-K频率-波速法计算格林函数时(Zhu and Rivera, 2002)，采样间隔为0.05 s，采样点数为1024个.在震源机制反演过程中，对Pnl和面波设置不同的带通滤波范围，Pnl波为0.02～0.15 Hz，面波为0.02～0.12 Hz，这样可以尽量减少背景噪声对波形的干扰，亦能压制地下小尺度结构散射造成的影响(陈伟文等，2012).为确保反演结果的可靠性，选取震中距500 km以内、位置分布均匀、三分量波形完整并且信噪比较高的台站参与反演，去除理论波形与实际波形的拟合相关系数较低的波形.最终获得研究区域287组3.5级以上地震事件的震源机制解.

1.2.2 HASH反演法及数据处理

对于研究区内一些地震震级较小(3.0≤ML<3.5)的地震事件，本文利用基于P波初动极性和S/P振幅比的HASH方法(Hardebeck and Shearer, 2002, 2003)开展了震源机制解反演.HASH方法的基本原理是基于不同方位、不同震中距的台站的地震波形记录，利用从震源向上射出的直达P波(Pg)和直达S波(Sg)引起的地动位移振幅比求解震源机制.其主要优点是可在一定程度上克服震源位置、速度模型和极性观测误差等，对于中小震地震有良好的适用性.在求解过程中为保证结果的可靠性，数据处理时仅选择震中距200 km以内的、信噪比大于3的数字波形资料.利用SAC软件人工读取了每个台站垂直向的初至波极性及S波到时，对于振幅的度量，将原始波形旋转至R-T-Z坐标系下.P波振幅的度量：第一个半周期的峰值，量取Z和R方向的振幅，S波振幅的度量：三分量Sg波到时2 s内的最大值.在震源机制反演过程中，台站分布范围也是十分关键.在数据处理过程中，每个地震事件的相邻台站最大空隙角小于150°；离源角小于90°.HASH方法还从断层不确定度、台站分布比、振幅比的数量等参数对断层面解的结果进行质量评价.在后文的讨论中只采用了A、B、C三类的结果，共计494组.

1.2.3 速度模型的选取

现有结果表明，青藏高原存在着明显的横向不均匀性(李敏娟等, 2018; 许英才等, 2018).如果在研究区域采用单一的速度模型计算离源角，可能会造成较大的计算误差.HASH方法能够最多采用10个分区的速度模型，因此本文基于可靠的地壳速度结构测量，采用了多速度模型(青海省地震局日常定位用速度模型; 王周元, 1984; 张少泉等, 1985; 赵珠和张润生, 1987; 周民都等, 2012; 金春华等, 2015)参与离源角计算并反演震源机制解(图2).CAP方法在震源机制反演过程中，允许体波和面波进行时间平移，因此反演结果对速度模型的依赖性较小(韦生吉等，2009)，本文直接采用每个地震震中位置处Crust2.0一维速度模型计算格林函数(Bassin et al., 2000).

图2 HASH反演法采用的地壳速度结构模型Fig.2 Velocity model used by HASH inversion method in this study

1.3 应力场反演方法

本研究采用Hardebeck和Michael(2006)提出的、基于震源机制解的MSATSI程序(Martínez-Garzón et al., 2014)进行构造应力场反演.该方法主要是基于线性反演方法，与网格搜索方法相比反演速度更快，能够最大程度上拟合数据，消除人为影响.主要工作原理是首先将震源机制解按照分布区域划分到若干个网格，然后采用阻尼最小二乘法进行反演，得到每个网格内的应力张量，并在最小二乘反演的过程中通过加入最优阻尼系数进行约束，以消除传统应力反演方法中人为因素带来的应力旋转，最大限度减少相邻区域间的应力差异.反演过程中，阻尼系数的确定是通过设定一系列阻尼值，得到应力场反演模型长度与数据拟合误差关系的折中曲线.曲线拐点即为最佳的阻尼系数值，低于该值，模型复杂程度得到提高，但对反演残差的改善作用不大，且观测值与理论值匹配度急剧降低；高于该值，随着模型的简化，反演误差会急剧增加(Hardebeck and Michael, 2006; 王晓山等, 2015; 郭祥云等, 2017).本研究将研究区划分为1°×1°的网格进行构造应力场反演.根据双力偶震源模型，存在两个相同可能性的断层面，即发震节面和辅助节面，因此MSATSI程序通过多次重采样的方法，随机选择其中的一个节面进行反演，统计在一定置信度下的反演结果，从而确定最优解和不确定范围.本研究反演过程中，在95%的置信区间内对原始数据进行2000次bootstrap重采样.

1.4 R值反演方法

应力张量是对称张量，有6个独立参数，完整的应力张量需要三个主应力的量值及P、T和B轴方向，但是确定绝对应力大小是非常困难的(万永革等, 2006).通常采用三个主应力的比例关系归一化应力大小，引入应力形因子R表征相对应力大小(Gephart and Forsyth, 1984; 万永革等, 2008)：

(1)

其中σ1、σ2、σ3分别表示最大、中间、最小主应力.当R值为0.5时，三个主应力成等差排列，σ1和σ3明确；当R值接近1时，σ2与σ3表现为一致的张应力状态，并在与稳定的压应力轴σ1垂直的平面内自由旋转呈双轴拉张状态，两轴不易区分；当R值接近0时，σ1与σ2表现为一致的压应力状态，并在与稳定的张应力轴σ3垂直的平面内自由旋转呈双轴压缩状态，两轴不易区分(Guiraud et al., 1989; Wan et al., 2016; 刘兆才等, 2019).

2 结果与分析

2.1 反演计算结果

图3给出了2012年6月3日18时49分青海共和MS4.2地震CAP反演结果，图3a为在最佳拟合深度处各台站的理论波形与实际波形拟合图，其中约78%的分量拟合相关系数均高于80%，表明理论波形与实际波形具有较好的拟合关系.图3b为该地震CAP波形反演误差随深度的变化图，具有最小误差的最佳拟合深度为6 km.

图3 利用CAP方法反演2012年6月3日18时49分青海共和MS4.2级地震示例(a) 地震震源机制和理论波形与观测波形拟合图. 红线代表理论地震图, 黑线代表观测地震图, 波形下方数字分别代表理论地震图相对观测地震图的偏移时间(单位: s)及二者的相关系数(%), 台站下方数字分别代表震中距(单位: km)和方位角(单位: °); (b) 地震震源机制解反演残差随深度分布图.Fig.3 Example of the focal mechanism determination by the cut-and-paste (CAP) method for an earthquake that occurred at 18∶49 on 3 June 2012(a) The corresponding waveform comparisons of synthetics (red) and records (black) at 8 stations. The two numbers under each segment are the time shift in second (upper) between the synthetic and record (positive means a delayed record) and the waveform correlation coefficient (lower, 100 indicates the best fit). Epicentral distances (in km) and azimuths (in degree) are given below the station codes; (b) Variation of fitting error with depth during the focal mechanism inversion of the earthquake occurred at 18∶49 on 3 June 2012.

为了验证震源机制解反演结果的可靠性，选取研究区部分震级大于3.8级的中小地震，同时采用CAP方法和HASH方法反演其震源机制解，并对结果进行对比；选取5.0级以上中强地震，将CAP反演结果与GCMT结果进行对比，结果如表1和图4a所示.不同方法的反演结果具有高度的一致性，说明采用的方法和反演结果可靠.

表1 不同方法求取震源机制结果比较Table 1 Comparison of different methods to obtain focal mechanism solutions

图4 (a) 本研究获得的781个地震事件的震源机制解分布, 其中绿色填充的膨胀区表示走滑型, 蓝色填充的膨胀区表示逆断型和逆走滑型(统称为逆断型), 红色填充的膨胀区表示正断型和正走滑型(统称为正断型), 灰色矩形内展示了分别利用不同方法得到的同一地震事件的震源机制解; (b) 从GCMT收集的96个地震事件的震源机制解分布, 沙滩球颜色标注同(a); (c) 本研究获得的及从GCMT收集的共计877个地震事件的震源机制解P轴水平分量投影分布; (d) 877个地震事件的震源机制解T轴水平分量投影分布. 其他标注符号同图1Fig.4 (a) 781 FMSs determined in this study. the focal mechanisms with green, blue, and red filled in expansion quadrants represent strike-slip, thrust, and normal fault earthquakes, respectively; the grey boxes show the comparison of the same events with GCMT, CAP and HASH solutions, respectively; (b) Same as (a) but for 96 FMSs collected from GCMT catalogue; (c) Horizontal projections of the P axes (blue line segments) of the 877 FMSs, including 781 events in this study and 96 events from the GCMT catalogue; (d) Same as (c) but for T axes (red line segments). The other labelling are the same as those in Fig.1

2.2 震源机制解总体特征

基于本文获得的781个震源机制解，参照世界应力图WSM(World Stress Map)的划分原则(Zoback and Healy, 1992; Heidbach et al., 2018)，根据震源机制解P、T、B轴倾伏角的大小，将震源机制解类型划分为6种，分别为正断型(NF)、正走滑型(NS)、走滑型(SS)、逆走滑型(TS)、逆断型(TF)和不确定型(U)，如表2所示.结果表明：研究区内走滑型地震为433个，占55.44%，正走滑型和正断型地震为177个，占22.66%，逆走滑型和逆断型地震为171个，占21.90%.研究区整体以走滑型地震为主，其次为正断型和具有一定走滑分量的正断型地震.

表2 震源机制解分类表Table 2 Categories of tectonic stress regime for focal mechanism

图4a给出了本研究获得的781个地震的震源机制解，图4b为收集的GCMT中96个地震的震源机制解，图4c和4d分别为本研究反演获得的以及收集的共计877个震源机制解的主压应力P轴和主张应力T轴的水平投影分布.整体上P轴呈NNE向～ENE向分布，受控于印度板块北东向欧亚板块持续俯冲和青藏高原物质东向扩散作用.P轴、T轴大多倾伏角较小(其中，P轴62.26%<30°，79.13%<45°；T轴52.57%<30°，74.91%<45°)，体现了近水平的挤压应力作用，且分布特征与区内走滑型断裂活动性质相对应，揭示了研究区所处的青藏高原东北缘块体间强烈的剪切作用方式.

2.3 构造应力场特征

利用MSATSI程序对本研究获得的781组震源机制解和收集的96组震源机制解数据进行应力场反演，最佳阻尼值取为1.1(图5)，网格内震源机制解的最小个数为5个，得到各应力单元应力张量反演结果及最佳拟合构造应力张量(图6、图7、表3)，在此过程中，剔除了不确定性较大的数据.反演结果表明，研究区最大主压应力轴σ1方向以NE向为主，且由巴颜喀拉块体西南部向外呈扇形辐射，至祁连造山带逆时针偏转为NNE向，至秦岭造山带顺时针偏转为ENE向—ESE向，与鄂拉山断裂、日月山断裂、祁连山北缘断裂、海源断裂、东昆仑断裂、西秦岭北缘断裂等边界断裂活动性质一致，且与青藏高原东北缘地形陡变带迹线近垂直，表明构造变形和地壳隆升与构造应力场关系密切.最小主压应力轴σ3方向以NW向—NNW向为主.σ1轴和σ3轴倾伏角较小，大部分近水平，表明研究区整体处于近水平的挤压或剪切应力环境，局部最大主压应力轴近垂直，处于拉张应力环境.中间主应力轴σ2倾伏角较大，多数近垂直.应力形因子R表征相对应力大小，为了更直观的获得研究区相对应力分布特征，采用插值法给出了连续的R值分布(图7、图9).揭示了研究区内R值具有明显的空间差异分布特征，表明块体间相互作用强烈.R值较小区域表明中间主应力值与主压应力值较为接近，中间主应力轴表现为压应力特征，R值较大区域表明中间主应力值与主张应力值较为接近，中间主应力轴表现为张应力特征.

图5 模型长度与数据拟合误差之间的折中曲线图空心圆旁所标数字是阻尼参数, “+”符号表示所选择的最佳阻尼系数.Fig.5 Trade-off between model length and misfit calculated from the corresponding damping parametersDamping values are shown beside each point and the best selected damping parameter is shown with a cross.

图6 主应力轴投影分布图黑色圆点表示主应力轴σ1、σ2、σ3; 深灰色扇形区表示主应力轴方向玫瑰图; 等值线云图表示主应力轴分布密度, 由冷色调蓝到暖色调红表示主应力轴投影分布密度由低到高.Fig.6 Stereonet map of principal stress axesThe black dots indicate the principal stress axes σ1, σ2, and σ3; The dark gray sector represents rose diagram of the principal stress direction; The contour map shows the distribution density of the principal stress axes, from cool blue to warm red, indicating the distribution density of the principal stress axes from low to high.

图7 应力场反演结果红色、绿色和蓝色圆点分别表示95%置信度下最大主压应力轴σ1、σ2和σ3的不确定范围; 黑色实线加号表示最优解; 底图颜色表示反演的应力形因子R值分布; 灰色曲线表示活动断裂; 蓝色曲线表示河流; 网格内数字代表该网格内参与反演的震源机制解个数.Fig.7 Stereomap of stress field inversion resultsThe coloured dots represent the stress axes in the 95 percent confidence interval by bootstrap, while the red, green and blue dots denote the orientations of the σ1, σ2 and σ3, respectively; the black crosses denote the best solutions for the principal stress axes in each grid; the color of base map represents distribution of the inverted stress shape factor R value; the gray curves represent active faults; the blue curves represent rivers; the number in the box at each grid shows the number of FMSs assigned to that grid in the stress inversion.

MSATSI程序采用Lund和Townend(2007)的方法计算了各应力单元的水平最大主应力SHmax方向，并基于WSM分类标准对其应力状态进行了划分(Zoback and Healy, 1992; Heidbach et al., 2018).笔者基于由MSATSI程序反演得到的最佳拟合构造应力张量，同样采用Lund和Townend(2007)的计算方法，得到各应力单元水平最大主应力最优方向，结果如图8所示.SHmax方向与最大主压应力轴基本一致，同样表现为由研究区西南部向外呈扇形扩展.研究区内应力结构类型以走滑型为主，祁连造山带和柴达木盆地西部区域、陇西块体东部区域和秦岭造山带则以逆断型或无法明确应力结构类型为主，这些区域地震活动较为频发，且地震震源机制类型多样(图4a)，地壳同时存在挤压、剪切和拉张变形.此外，从《中国及邻区现代构造应力场图》中搜集研究区内基于水压致裂和应力解除等原位实测方法获得的地壳浅表层(<2 km)SHmax方向(谢富仁和崔效锋, 2015)，并按照本研究应力单元划分标准，将地壳浅表层SHmax方向进行分区平均，同样绘制到图8中.在研究区中部区域，地壳浅表层SHmax方向总体以NE向为优势方位，同时存在近NS向和近EW向分布.西南部以NE向为主，东南部以NW向为主，自西向东呈现顺时针偏转，与震源深度SHmax方向分布特征一致.地壳浅表层实测SHmax方向与所在应力单元内震源层SHmax方向近平行或呈一定角度斜交，两者夹角大多数在30°以内，表明地壳浅表层应力场与震源深度应力场具有较好的一致性，与深部运动即青藏高原物质东北向扩散相耦合.但同时地壳浅表层应力场影响因素较多，受地形地貌、构造、岩性等影响易形成局部应力场，因而与深部应力场存在一定的差异.在研究区东南部，即柴达木—陇西块体、巴颜喀拉块体与秦岭造山带、四川盆地北缘过渡区域，地壳浅表层实测SHmax方向由NE向偏转为NW向，且与震源深度SHmax方向夹角较大，该区地壳深浅构造应力场以及变形特征可能存在较大的差异性.

3 讨论

3.1 区域构造应力场与动力学过程

据中国构造应力分区图(谢富仁和崔效锋, 2015)，研究区及周边区域被划分为祁连山—河西走廊应力区、海原—六盘山应力区、柴达木盆地应力区、西秦岭应力区和巴颜喀拉应力区，其中，除祁连山—河西走廊应力区为挤压区外，其余应力区均为剪切区(图9).该应力分区主要考虑块体和断裂相互作用的影响，将应力状态具有较好一致性的区域划分为同一分区.本研究表明，构造应力在空间上的非均匀性分布特征较为显著，在现有应力分区基础上，可进行次一级或局部应力分区.根据应力形因子R值的定义(式(1))，图9中暖色调区域表示双轴拉张区，冷色调区域表示双轴压缩区.双轴压缩区以祁连造山带，陇西块体东部，以及由鄂拉山断裂带、西秦岭北缘断裂带、贵德断裂和东昆仑断裂带围限区域为主；双轴拉张区以柴达木盆地、秦岭造山带、陇西块体西部和巴颜喀拉块体为主.双轴拉张区内主应力轴方向以σ1明确为特征，表现为NE-SW向压缩和NW-SE向拉张，且挤压应力大于拉张应力.双轴压缩区内主应力轴方向以σ3明确为特征，表现为NE-SW向挤压和NW-SE向拉张，拉张应力可能远小于挤压应力，也可能与挤压应力相差不大.Wang和Shen(2020)基于1991至2016年GPS观测结果，得到了中国大陆地壳形变场.其研究结果表明，在祁连山地区以NE-SW向缩短和NW-SE扩展为特征，并且缩短速率明显高于扩展速率，而在东昆仑断裂带显示出不低于20×10-9/yr的左旋应变率.本文研究揭示的应力机制与基于GPS观测揭示的地壳形变特征基本吻合，表明研究区应力场与应变率场基本耦合，地壳形变主要受控于背景构造应力场(图9).综合区内R值、主应力分布特征，表明研究区总体主要受控于印欧板块碰撞造成的NE向挤压，并在这种挤压背景下，存在高原物质向北东推挤受周边稳定刚性块体阻挡后由重力位能垂向差异产生的拉张力，并向外扩展，最终形成研究区内现今走滑为主的应力结构类型(图8).同时区内构造受控于区域应力场，吸收拉张和剪切构造变形，在地质演化时期内存在由挤压向走滑构造转变的过程，如鄂拉山断裂在第四纪初期由挤压逆冲转换为右行走滑(袁道阳，2004).

表3 1°×1°网格划分反演得到的应力场参数Table 3 The inverted stress field parameters on grid of 1°×1°

续表3

图8 (a) 水平最大主应力SHmax分布. 中心为圆形符号短线表示震源机制解反演结果，其中绿色、蓝色、红色和黑色短线分别代表走滑型、逆断型、正断型和过渡型应力状态; 中心为菱形符号短线表示由水压致裂和应力解除等原位实测方法得到的最大水平主应力方向(谢富仁和崔效锋, 2015), 应力单元左下角或左上角数字角度代表该单元内震源深度SHmax方向与地壳浅表层SHmax方向的夹角; 32°N—33°N纬度行及96°E—97°E经度列左上角圆圈内数字代表应力单元坐标, 纬度行对应y坐标, 经度列对应x坐标. (b) 震源深度SHmax方向与地壳浅表层SHmax方向夹角统计图Fig.8 (a) Distribution of the maximum horizontal stress orientations. Short lines with a circular symbol in the center represent SHmax orientations obtained by inversion of focal mechanism solutions, and the green, blue, red and black symbols represent strike-slip, reverse, normal and transitional stress states, respectively; Short lines with a diamond symbol in the center represent SHmax orientations obtained by in-situ measurement methods such as hydraulic fracturing and stress relief (Xie and Cui, 2015); The numbers in the lower left or upper left corner of the stress unit represent the angles between SHmax orientation of the focal depths and those of the shallow crust; The number in the circle in the upper left corner of 32°N—33°N latitude row and 96°E—97°E longitude column stand for the coordinate of the stress unit, the latitude row corresponds to the y coordinate, and the longitude column corresponds to the x coordinate. (b) Statistical diagram of the angle between SHmax orientation of the focal depths and those of the shallow crust

图9 青藏高原东北缘动力模型图彩色底图为由震源机制解反演得到的R值分布, 其中暖色调区域表示双轴拉张区, 冷色调区域表示双轴压缩区; 应力分区引自中国构造应力分区图(谢富仁和崔效锋, 2015); 紫红色互相垂直箭头代表基于GPS得到的主应变率及其方向(Wang and Shen, 2020); 浅灰色宽箭头代表块体运动方向.Fig.9 The dynamic model of the NE margin of the Tibetan PlateauColor base map is the distribution of R value obtained by inversion of focal mechanism solutions, the warm color area represents the biaxial stretching area, and the cool color area represents the biaxial compression area. Stress zones is quoted from China's tectonic stress field map (Xie and Cui, 2015). The fuchsia mutually perpendicular arrows represent the principal strain rate and its direction based on GPS observations (Wang and Shen, 2020); The light gray wide arrows represent the direction of movement of the blocks.

图10 地壳不同深度应力场分布特征(a,c,e) 应力场反演结果(95%置信区间), 其他符号代表意义同图7; (b,d,f) 水平最大主应力分布(95%置信区间, 应力状态分类采用WSM标准: 绿色代表走滑型, 蓝色代表逆断型, 红色代表正断型, 黑色代表过渡型), FD: 震源深度.Fig.10 Stress field of different depths of the crust(a,c,e) Stereomap for the principal stress field (within 95 per cent confidence region) of the shallow and deep crust (see Fig.7 for the other labelling); (b,d,f) Maximum horizontal stress orientations (SHmax, within 95 per cent confidence region with stress regime categories according to the WSM standard: Blue, green, red and black lines denote thrust, strike-slip, normal and mixed environment, respectively) in different depths of the crust; FD: Focal depth.

图11 青藏高原东北缘构造应力场及各向异性分布特征(a) 不同深度地壳构造应力场分布, 虚线矩形区指示(b)范围; (b) 青藏高原东北缘东南部区域深浅地壳构造应力场分布(<2 km, ≤10 km, 和>20 km), 其中2km以浅SHmax方向由原位地应力实测获得, 震源深度SHmax方向由震源机制解反演获得, 背景图为该区20 km深度地壳横波速度(Shen et al., 2016); (c) 近震直达S波分裂(Hu et al., 2020; 钱旗伟等, 2017; 张辉等, 2012); (d) Ps波分裂(Xu et al., 2018; Wang et al., 2016; 邵若潼等, 2019); (e) XKS波分裂(Chang et al., 2017); (f) 青藏高原东北缘深部动力学机制示意图, 修自Clark and Royden (2000).Fig.11 Stress field and seismic anisotropy of the NE margin of the Tibetan Plateau(a) Characteristics of stress fields in different crustal depths, the dotted rectangular area indicates the range of (b); (b) Stress fields of the shallow and deep crust (<2 km, ≤10 km, and >20 km) in southeastern of the NE margin of the Tibetan Plateau, SHmax within 2 km are obtained by in situ stress measurement methods, SHmax in seismogenic depths are obtained by focal mechanism inversion, the background colours denote the average S-wave velocity at depth of 20km (Shen et al., 2016). (c) Local S wave splitting (Hu et al., 2020; Qian et al., 2017; Zhang et al., 2012); (d) Ps wave splitting (Xu et al., 2018; Wang et al., 2016; Shao et al., 2019); (e) XKS wave splitting (Chang et al., 2017); (f) Deep dynamic mechanism model of the NE margin of the Tibetan Plateau, modified from Clark and Royden (2000).

利用本研究CAP反演法得到的最佳震源矩心深度，给出了不同地壳深度(≤10 km, 10～20 km和>20 km)的应力场分布特征(图10, 图11a).由于按地壳深度分层后地震活动空间分布差异较大，在该部分研究中每个网格采用的震源机制至少为1个，尽管部分网格数据量相对较少，但程序采用阻尼最小二乘法反演同时得到每个网格点的应力张量并进行平滑，依然可以得到较为可靠的结果(Menke, 1989; 刘兆才等, 2019).结果表明不同深度的SHmax方位均以NE向为主，总体上地壳构造应力场分布具有垂向一致性，且与GPS速度场总体特征一致，受控于印欧板块持续碰撞导致的青藏高原快速隆升并向北东方向推挤作用.Han等(2019)通过对比深(>20 km)浅(<20 km)地壳应力场，得到与本文较为一致的应力场分布特征.图11a明显揭示出SHmax方向由祁连造山带和柴达木—陇西块体向秦岭造山带呈现由NE→ENE→ESE向偏转的特征，大致与地形差异边界或地块碰撞边界垂直.已有研究表明，地震震源深度的垂向分布受到岩石圈热流变性质的控制，地震主要发生于岩石圈的脆性域内(Ranalli and Murphy, 1987; 张国民等, 2002)，因此震源深度可以反映地壳岩石圈脆性变形层的厚度.从图11a可以看出，地震多集中于20 km以浅，中下地壳多以塑性为主，青藏高原东北缘构造应力场分布特征似乎是深部软弱物质北东向挤出并向周缘扩散作用在浅部的响应.在柴达木—陇西块体和巴颜喀拉块体东南部与秦岭造山带和四川盆地的过渡区域，中下地壳应力场与上地壳应力场存在差异，深部地壳变形可能与上地壳变形解耦.考虑到震源深度的不确定性，为了更为清楚的表示深浅地壳应力场分布特征，图11b给出了10 km以浅和20 km以深由震源机制解反演得到的SHmax方向，同时叠加了2 km以浅由原位地应力测量方法实测得到的SHmax方向.可以看出，陇西块体10 km以浅地壳深度SHmax方向耦合性较好，而秦岭造山带和巴颜喀拉块体东南部不同地壳深度SHmax方向差异较大(图11b中虚线椭圆范围内).2 km以浅SHmax方向以WNW-ESE为主，10 km以浅SHmax方向以WSW-ENE为主，而20 km以深SHmax方向同时具有WNW-ESE和WSW-ENE两种优势方位.Han等(2019)认为深部地壳SHmax方向一般垂直于低到高的横波速度边界.图11b显示，在横波速度高低过渡区域地壳SHmax方向出现明显的解耦现象，深部SHmax方向耦合于下地壳流方向(Clark and Royden, 2000; Zhao et al., 2021)，与青藏高原东南缘共同构成下地壳弱物质东向流出的两个分支通道(Clark and Royden, 2000; Sun et al., 2012).

3.2 壳幔变形特征与变形机制

青藏高原东北缘近震各向异性研究表明上地壳各向异性同时受到区域应力与断裂构造的双重约束(图11c).快波偏振以NE向为主，与SHmax方向一致，同时耦合于地表形变场(张辉等, 2012; 郭桂红等, 2015; 钱旗伟等, 2017; Wang and Shen, 2020).该区与阿拉善块体交接区受到NE向强烈挤压作用，特别是祁连造山带受到南北两侧阿拉善块体和柴达木块体的俯冲以双向逆冲形式挤出，汇聚速率达16 mm·a-1(许志琴等, 1999)，表明该区经历着NE向地壳缩短.阿拉善块体、祁连块体、柴达木—陇西块体交汇区内，同时表现出与断裂构造一致的NW向优势偏振方向，且与上地幔XKS偏振方向和板块绝对运动方向(APM)一致(图11c, e)，壳幔变形方式同时存在差异性和关联性，与块体交汇区复杂的深部动力学机制相关.柴达木—陇西块体东南部和秦岭造山带快波偏振方向则以WNW向—NW向为主，与上地幔XKS偏振方向以及中下地壳SHmax方向近似平行(张辉等, 2012; 郭桂红等, 2015; Chang et al., 2017)，与上地壳SHmax方向呈大角度斜交，因此主要受断裂构造影响.

基于远震接收函数数据的研究揭示了青藏高原东北缘中下地壳显著的各向异性特征，其快剪切波方向整体以WNW-ESE或NW-SE向为主(图11d)，大致与断裂构造及造山带走向一致且垂直于水平最大主压应力方向，从而将青藏高原东北缘中下地壳各向异性归因于挤压碰撞相关的纯剪切引起的角闪石晶格优势排列(Wang et al., 2016; Xu et al., 2018; 邵若潼等, 2019).但在共和盆地和青海湖盆地交汇区域，以及鄂尔多斯块体西缘则呈现以NE向—ENE向为主的快波偏振方向(图11d椭圆形区域)，与SHmax方向、GPS速度场、上地壳各向异性等均具有较好的一致性，而与近NW向分布的XKS偏振方向呈大角度斜交或近于垂直(图11f)，表明壳幔具有不同的变形机制.在秦岭造山带及其周边区域，Ps快波方向与中下地壳SHmax方向一致，而与上地壳SHmax方向解耦，地壳通道流的存在可能是造成该区深浅地壳应力场解耦以及中下地壳各向异性的主要原因(Han et al., 2019; Zhao et al., 2021).

相对于青藏高原东北缘大多小于0.4 s的地壳各向异性时间延迟(张辉等, 2012; 郭桂红等, 2015; 钱旗伟等, 2017; Wang and Shen, 2020; Wang et al., 2016; Xu et al., 2018; 邵若潼等, 2019)，XKS波分裂结果估算的时间延迟平均为1.3 s，因此主要反映了上地幔各向异性的特征(常利军等, 2021).Silver(1996)认为大陆上地幔各向异性受构造应力作用引起的岩石圈变形影响显著，构造活动区快波方向通常与大型走滑断裂和造山带构造方向一致.XKS波分裂研究表明青藏高原东北缘上地幔快波方向以NW-SE为主(图11e)，与区域走滑断裂和造山带构造走向一致，揭示了强烈的NE-SW挤压应力环境，与由震源机制解揭示的NE-SW向主压应力一致.而在秦岭造山带及周边区域，XKS偏振方向偏转为E-W或WNW-ESE向，且与中下地壳SHmax方向一致，前人根据该区较大的XKS慢波时间延迟和较薄的岩石圈，认为秦岭造山带是青藏高原岩石圈物质的挤出通道，同时也是软流圈物质东流的地幔流通道(王琼等, 2013; Chang et al., 2017; 常利军等, 2021).根据地壳应力场和上地幔各向异性特征，不能排除该区存在壳幔解耦变形的可能.

3.3 青藏高原东北缘动力学特征

综合以上研究，本文初步提出了青藏高原东北缘岩石圈多种变形机制共存的模式(图11f, 图12).高原物质向北东推移过程中受到阿拉善和鄂尔多斯两个稳定刚性块体的阻挡作用，总体上沿NW-SE向发生了伸展变形，形成了一系列NW-SE方向的造山带和断裂带，在上地幔产生了NW-SE方向的快剪切波(王琼等, 2013; Chang et al., 2017; 常利军等, 2021).在阿拉善块体与柴达木盆地之间的祁连造山带北部，地壳各向异性也以NW-SE向分布的快波偏振方向为主要特征(张辉等, 2012; 钱旗伟等, 2017; Wang et al., 2016; Xu et al., 2018)，SHmax方向则与快波偏振方向近垂直，体现了以岩石圈连续增厚变形模式为主.在柴达木盆地东南部地区以及鄂尔多斯块体西缘，地壳剪切波快波方向以NE-SW向为主，各向异性来源于由SHmax作用导致的沿应力方向定向排列的微裂隙(张辉等, 2012; 钱旗伟等, 2017; Wang et al., 2016; Xu et al., 2018)，地壳和地幔的变形机制可能不同.柴达木盆地强度低于阿拉善块体和鄂尔多斯块体，对下地壳弱物质流没有形成有效的阻挡，物质流可能继续向北东方向推移直至受到阿拉善块体阻挡，继而产生伸展变形以及地壳增厚.鄂尔多斯块体西缘是构造活动强烈的青藏高原向稳定的鄂尔多斯块体的过渡地带，深部运动和构造复杂，壳幔可能存在解耦现象(许英才等, 2019).柴达木—陇西块体和巴颜喀拉块体东南部与秦岭造山带和四川盆地过渡区域壳幔剪切波快波偏振方向均以近EW向或WNW向为主(张辉等, 2012; 钱旗伟等, 2017; Hu et al., 2020; Wang et al., 2016; Xu et al., 2018; 邵若潼等, 2019; Chang et al., 2017; Han et al., 2019)，中下地壳SHmax方向与壳幔剪切波快波偏振方向一致(图8, 图11a,b)，耦合于下地壳流方向，与青藏高原东南缘共同构成下地壳弱物质东向流出的两个分支通道(Clark and Royden, 2000; Sun et al., 2012).上地壳SHmax方向与快波偏振方向斜交，深浅地壳SHmax方向存在一定差异，深部地壳变形与浅部发生解耦.

总体而言，青藏高原东北缘壳幔变形机制和深部动力学特征极为复杂，很难用单一模式解释，可能同时存在多种变形方式.根据地壳构造应力场分布特征并结合壳幔各向异性，青藏高原东北缘可能以岩石圈连续增厚变形与下地壳通道流的共同作用模式为主，共同影响并导致地壳增厚和高原隆升(图12).青藏高原北东向持续推挤，向北受到具有古老克拉通性质的阿拉善块体阻挡，发生地壳均匀增厚导致的地壳缩短，从而形成陡峭的地形边缘，如图12所示A-A′剖面在约50 km范围内地形陡变近3000 m.同时，受青藏高原北东向挤出作用以及上地幔物质加热作用，巴颜喀拉块体深部地壳物质发生部分熔融并产生北东向扩展的下地壳通道流(Liu et al., 2014; Zhao et al., 2021)，导致下地壳增厚并伴随着地壳流的向外扩展形成较为平缓的地形，如图12所示B-B′和C-C′剖面.而D-D′剖面位于两种变形方式的过渡区域，其地形起伏特征介于两种典型地形之间.

图12 (a) 青藏高原东北缘深部动力学机制三维示意图, 蓝色箭头代表高原物质逃逸方向, 紫色箭头代表下地壳变形方向, 白色实线代表(b)中地形剖面; (b) 青藏高原东北缘典型地形剖面, 剖面位置如(a)及图1所示Fig.12 (a) 3D perspectives of deep dynamic mechanism model of the NE margin of the Tibetan Plateau, the blue arrows indicate orientation of material escape, the purple arrows indicate the lower crustal deformation, the white solid lines represent the terrain profile in (b); (b) Typical topographic profiles for the NE margin of the Tibetan Plateau, the locations of profiles are shown in (a) and Fig.1

4 结论

本研究利用CAP全波形反演方法及HASH初动极性和振幅比方法，系统反演了研究区内ML≥3.0地震震源机制解781个，同时收集了MW≥4.6中强地震震源机制解96个，并使用MSTASI程序进行了应力场反演，获得了研究区应力张量、R值及水平最大主应力优势方向，详细分析了该地区震源机制解和构造应力场特征，最后根据计算和分析结果，联合地表GPS观测和壳幔各向异性研究，综合讨论了青藏高原东北缘变形机制和动力学特征.分析结果显示：

(1)研究区震源机制解类型以走滑型为主，P轴、T轴大多倾伏角较小，最大主压应力轴σ1方向以NE向为主，且由巴颜喀拉块体西南部向外呈扇形辐射，至祁连造山带逆时针偏转为NNE向，至秦岭造山带顺时针偏转为ENE向—ESE向，分布特征与区内走滑型边界断裂活动性质相对应.σ1轴和σ3轴倾伏角较小，大多近水平，表明研究区整体处于近水平的挤压或剪切应力环境.SHmax方向与最大主压应力轴分布特征基本一致，应力结构类型以走滑型为主.

(2)青藏高原东北缘最大主压应力方向与最大主应变方向具有很好的一致性，应力机制与基于GPS观测揭示的地壳形变特征基本吻合，表明研究区应力场与应变场基本耦合，地表形变主要受控于区域构造应力场.大部分地区不同地壳深度的SHmax方向一致性较好，而柴达木—陇西块体和巴颜喀拉块体东南部与秦岭造山带和四川盆地的过渡区域，深浅地壳SHmax方向具有差异性，深部地壳变形与浅部可能发生解耦.地壳构造应力场可能是深部软弱物质北东向挤出并向周缘扩散作用在浅部的响应.

(3)青藏高原东北缘壳幔变形机制和深部动力学特征极为复杂，很难用单一模式解释.阿拉善块体与柴达木盆地之间的祁连造山带北部地区，壳幔变形垂向连贯，体现了以岩石圈连续增厚变形模式为主.柴达木盆地东南部及鄂尔多斯块体西缘，地壳各向异性主要来源于由SHmax作用导致的沿应力方向定向排列的微裂隙，壳幔变形机制可能不同.青藏高原东北缘东南部区域，深部地壳变形可能与浅部发生解耦，至鄂尔多斯块体与四川盆地之间的秦岭造山带区域，中下地壳SHmax方向与壳幔剪切波快波偏振方向一致，并耦合于下地壳流方向.青藏高原东北缘可能以岩石圈连续增厚变形与下地壳通道流的共同作用模式为主.