作为控制断层两盘相对运动的重要因素，断裂带介质力学性质与断层面上的滑动速率及应力状态(Hergert and Heidbach, 2010)、区域地壳运动速度场图像及其动力来源(Lei et al., 2010; He et al., 2013)、地震复发(Zhu and Zhang, 2010)以及地质构造活动(Gao and Wang, 2014; Dielforder et al., 2015)等密切相关.从物理上，断层力学性质可以理解为断层面两盘之间的摩擦系数(Hergert and Heidbach, 2010)或断裂带的剪切模量(Fialko et al., 2002)，但这一参数很难从野外观测的角度定量给出(Bird and Kong, 1994).研究中通常通过以下两种确定方式，一种是根据地质上的认识做定性或半定量的判断给出.如：弱化带模型中剪切模量取周边介质的1/5(祝爱玉等, 2015a)、1/10(胡勐乾等, 2014)，接触摩擦模型中取摩擦系数0.4(李玉江等, 2014)、0.65(祝爱玉等, 2015b).另一方面，随着对地壳介质认识的不断深入，特别是近年来随着GPS等观测资料的丰富，有限元模拟技术在国内外地壳动力学研究中得到广泛应用(陈连旺等, 2008; 邵志刚等, 2008; 朱守彪和张培震, 2009; Hergert and Heidbach, 2010; Luo and Liu, 2010; Moreno et al., 2010; Zhu and Zhang, 2010; Shi et al., 2014; Li et al., 2015a; 李玉江等, 2017).同时也为断层介质力学性质的定量研究提供了新的方法.朱守彪和张培震(2009)以龙门山断裂地震复发周期为约束，模拟给出汶川地震前龙门山断裂摩擦系数可达0.6.He等(2013)分析青藏高原东北缘的动力学环境时，也关注了研究区域断裂带介质力学属性对模拟结果的影响，认为青藏高原东北缘当前地壳运动速度场图像受控于海原断裂带的弱力学属性，并给出海原断裂带当前断层力学参数(摩擦系数)的最优取值可低至0.05，远低于针对该断裂有效摩擦系数的经验取值0.6~0.8(邵志刚等, 2010).Hergert和Heidbach(2010)以GPS观测数据为约束，通过有限元数值模拟手段在地球物理观测资料稀少的马尔马拉海域开展海底断裂滑动速率及地震危险性研究的模拟结果表明，区域地壳运动速度场图像可以反映出断层力学性质强弱，当前马尔马拉海海底断裂有效摩擦系数高，断层处于强闭锁状态且应力积累显著，附近的伊斯坦布尔市地震危险性值得注意.

青藏高原东北缘是印度—欧亚板块碰撞向北东方向扩展的前缘，是我国大陆地震活动频繁，地壳变形强烈的区域之一.20世纪以来该区域发生了1920年海原8.5级地震、1927年古浪8级地震以及2001年昆仑山口西8.1级地震等多次强震.区域内发育有阿尔金断裂带、广义海原断裂带以及西秦岭北缘断裂带等大型断裂带.其中广义海原断裂带是该区域一条重要的大型左旋走滑断裂带，是青藏高原东北部的地貌边界与构造边界(Tapponnier et al., 2001).东起六盘山，向西经海原、南西华山、哈思山、景泰、天祝、冷龙岭、祁连、沿大通河，到哈拉湖止，全长约980 km.构造变形上，该断裂带西段表现为左旋逆冲(Meyer et al., 1998; Hetzel et al., 2004)，中段老虎山断裂、海原断裂以左旋走滑为主(Gaudemer et al., 1995)，东段六盘山断裂表现为左旋逆冲(Burchfiel et al., 1991)，现代GPS地壳运动速度场观测也表现出相同的结果(Zhang et al., 2004; 郑文俊等, 2009).地震活动图像以及地壳形变观测揭示的断层运动状态显示：该断裂中部狭义海原断裂曾发生了1920年海原大地震，地壳形变观测结果显示该段断层活动闭锁程度较低，同时可能存在一定程度的蠕滑(Li et al., 2016)；金强河—毛毛山—老虎山断裂是天祝空段(Gaudemer et al., 1995)的主要活动断裂，断层活动闭锁程度较高(李强等, 2014; Li et al., 2016)，未来强震危险值得关注，但也有研究指出老虎山断裂存在长约35 km的无震滑移段(Jolivet et al., 2012)；东段六盘山—宝鸡断裂现代地震活动弱(M7专项工作组, 2012)，基于GPS观测的断层闭锁研究(郝明等, 2014；李强等, 2014; Li et al., 2016)指出跨六盘山断裂西侧存在约300 km的地壳缩短区域，断层闭锁程度高，未来强震危险同样值得关注.这些研究均从不同角度表明广义海原断裂带各段活动水平及应力状态差异较大.此外GPS观测显示西秦岭北缘断裂当前也处于强闭锁状态(M7专项工作组, 2012).但从构造变形角度，青藏高原东北缘是一个由北东东向的阿尔金断裂带、北西西向的广义海原断裂以及近东西向的东昆仑断裂带这三条大型边界断裂围成的相对独立的活动地块，块体内部发育有鄂拉山断裂、日月山断裂等右旋剪切构造带以及逆冲为主的祁连山汇聚带和六盘山汇聚带(袁道阳等，2004).以海原断裂、西秦岭断裂、日月山断裂以及六盘山断裂围成的陇西块体为例，构造变形研究表明六盘山断裂存在一次由东西向挤压到向南逃逸的构造转换过程(袁道阳等，2004).因此，GPS观测显示的300 km地壳缩短(郝明等, 2014)还可能是东西向挤压形成的，其应力积累可能在构造转化过程中释放.近来，又有文献(Chen and Bürgmann, 2017；Harris, 2017)报道指出闭锁率较低的蠕滑断层同样具有孕育强震的能力.

综上，从探索区域当前地壳运动的动力学机制和地震危险性分析的角度出发，有必要结合现有GPS观测数据，探索青藏高原各主要活动断裂的应力状态.为此本文在考虑区域活动断裂空间展布以及下地壳黏滞性结构的基础上构建了青藏高原东北缘有限元动力学模型，以GPS观测数据为约束，通过大量模拟试验分析青藏高原东北缘主要活动断裂的断层力学参数以及当前应力状态.

本文主要关注祁连地块，考虑青藏高原动力学环境和块体划分(邓起东等, 2003)，为避免GPS观测点空间分布稀疏对模型边界约束的可能影响(Hergert and Heidbach, 2010)，本文尽力使得模型边界通过GPS观测点较密集的区域(图 1和图 2)；为尽量减少阿拉善地块、鄂尔多斯地块与青藏高原东北缘下地壳结构差异对模拟结果的影响(He et al., 2013)，模型范围尽量保持在青藏块体内部，结合断层展布以及GPS点空间分布，在北侧边界和东侧边界少量跨入阿拉善地块和鄂尔多斯地块；本文关注的区域主要在祁连地块，因此模型的南边界尽可能避开东昆仑断裂(图 1).考虑青藏高原东北缘的动力学环境，本文模型中包含了阿尔金断裂、广义海原断裂、西秦岭北缘断裂、鄂拉山断裂以及日月山断裂等活动断裂，断裂数据采自郑文俊等(2009)和袁道阳等(2004)的研究.

图 1

Fig. 1

图 1 青藏高原东北缘主要活动断裂及历史地震破裂区分布(据M7专项工作组, 2012) Fig. 1 Distribution of major active faults and historical earthquake ruptures in the northeastern margin of the Tibetan plateau(Working Group of M7, 2012)

图 2

Fig. 2

图 2 研究区域GPS观测速度场(1999—2007) 红色误差椭圆标记的GPS观测速度用于插值给出模型的边界条件；模型内部黑色空心误差椭圆标记的GPS观测速度用于检验数值模拟结果的可靠性；洋红色矢量为本文模型边界约束. Fig. 2 GPS velocity field (1999—2007) in the study area The red error bar marked data near the model border are used to give the boundary condition for interpolation; the black error bar marked data inside the model border are used to test the reliability of the simulation results; the magenta arrows are the boundary conditions in this paper.

在上述基础上，通过对所调研的地质资料进行假设构建青藏高原东北缘动力学模型：根据青藏高原东北缘地壳速度结构(陈九辉等, 2005)，模型纵向厚度近似为50 km；考虑孕震层深度(段星北, 1997; 张国民等, 2002)，假定模型上地壳为20 km厚弹性层，弹性力学参数参考Crust1.0(Laske et al., 2013)以及嘉世旭和张先康(2008)，张洪双等(2015)关于青藏高原东北缘速度结构研究结果给出，上地壳弹性模量84.48 GPa、泊松比0.24，下地壳弹性模量102.57 GPa、泊松比0.25；考虑到本文关注重点为长期地壳形变，所以模型下地壳采用更加符合地壳长期变形模拟(邵志刚等, 2007)的Maxwell黏弹性体，黏滞系数1.0×10¹⁹Pa·s，参考He等(2013)研究结果给出.研究区域阿尔金断裂、广义海原断裂、西秦岭北缘断裂等均以走滑运动为主，将断裂带简化为10 km宽的垂直剪切变形带，其剪切模量将是本文重点讨论的断层介质力学参数.模型网格尺寸控制在10 km内，断层带网格尺寸控制在5 km内，模型共包含355482个单元，281282个节点.

由于模型涉及的断层较多，非线性摩擦接触断层模型计算耗时长，收敛困难(李玉江等, 2009)；为提高计算速度并保证得到与摩擦接触单元等效的计算结果，本文作者经过大量试验，引入复合材料力学中正交各向异性本构关系(沈观林和胡更开, 2006)，表征断层的变形行为(垂直断层方向的力学参数与周边介质一致，而走向和倾向方向剪切模量可以根据情况调整，实现走向和倾向方向的自由剪切滑动.对比研究表明(石富强等，2018)：相比直接降低弹性模量的介质弱化模型，基于正交各向异性本构关系的断层模型在地壳运动位移场以及应力场方面均与非线性摩擦接触模型的模拟结果吻合良好).假定断层带为夹在块体之间的类似于层合板的一类特殊介质：其在垂直于断层方向力学性质与周边块体相近，不易产生大挤压或拉张变形；而在平行于断层运动方向上的剪切模量较小，能够产生较大的剪切变形.则根据复合材料力学(沈观林和胡更开, 2006)基础知识，断层内应力应变关系在其三个主方向上满足关系为

(1)

其中ε=[ε₁₁ ε₂₂ ε₃₃ ε₁₂ ε₂₃ ε₁₃]^T，σ=[σ₁₁ σ₂₂ σ₃₃ σ₁₂ σ₂₃ σ₁₃]^T分别为三个主方向的应变和应力分量，K为柔度矩阵.

断层介质的变形主要是走向和倾向的剪切行为，其他方向上的介质力学性质与周边介质一致.为此正交各向异性本构关系的柔度矩阵可以表述为

(2)

其中E、v分别为周边介质弹性模量和泊松比，G_slip为断层在剪切方向(走向和倾向)上的剪切模量，独立于周边介质参数E、v.为便于表达，我们引入一个无量纲参数λ=G_slip/G₀，其中G₀=E/2(1+ν)为周边地壳介质剪切模量.不同断裂、断裂不同分段的走向和倾角不同，考虑模型坐标系和断裂走向、倾角关系将上述柔度矩阵旋转即可得到模型坐标系下断层介质的本构关系(石富强等，2018).

对于黏弹性材料的应力响应包括弹性部分和黏性部分，在小应变理论下，各向同性的黏弹性本构方程可以写为(李玉江等, 2017)：

(3)

其中σ为柯西应力，G(t)、K(t)分别为柯西应力的剪切部分和体积应力部分，利用Prony级数对其展开的：松弛时间为为介质黏弹性系数.一般认为地壳介质为体积不可压缩体(朱守彪和张培震, 2009; 高中等, 2012)，因此进一步假定K(t)=K_∞.

地壳动力学数值模拟研究中，模型初、边值条件是模型的重要组成部分，但又很难深入地壳通过测量来定量给出.GPS地壳运动速度观测直观地给出了不同区域地壳运动速度差异，是深部地壳运动在地表的直接反应.本文主要讨论的是上地壳断层剪切模量的强度问题，忽略地壳运动速度沿深度方向的差异对模拟结果的影响(刘峡等, 2010)，认为深部地壳运动速度与地表观测的一致.不考虑深部物质上涌等对地壳的作用，约束模型底部的垂向位移为0.为进一步体现研究区域断层力学参数对地壳运动速度场的影响，本文选取模型边界附近以及模型边界外的GPS观测数据充当“运动员角色”(图 2，红色误差椭圆标记的GPS观测速度)，通过线性内插给出模型边界条件(图 2，洋红色箭头矢量)；而模型内部的GPS观测数据充当“裁判员角色”来检验模拟结果(图 2，黑色空心误差椭圆标记的GPS观测速度).

关于地壳动力学数值模拟结果检验，一般追求的是模拟结果与观测结果一致，认为模拟结果与观测结果之差最小时的结果可靠(Hergert and Heidbach, 2010; 刘峡等, 2010; He et al., 2013).但实际操作过程中几乎不可能甚至做不到模拟结果与观测结果完全一致(残差为0).因此，很难评价该最小值是不是最合理.为了更好的提高模拟结果可靠性评价的可操作性，本文引入GPS观测误差椭圆，对He等人的检验方法(He et al., 2013)做了一定补充，认为只要当模拟值与观测值之差位于观测结果误差椭圆内部则模拟结果可靠，反之，可信度降低(图 3).定义δ_i为第i点观测点的评价参数，则有：

(4)

图 3

Fig. 3

图 3 模拟结果评价示意图 Fig. 3 Schematic of evaluation on simulation results

其中r为误差椭圆极半径，θ为极半径对应的幅角，则：

U_{i, x}、V_{i, x}、U_{i, y}、V_{i, y}分别为x(E)方向模拟值和观测值以及y(N)方向的模拟值和观测值，Δ_{i, x}、Δ_{i, y}分别为GPS在x(E)和y(N)方向的观测误差.在物理意义上，δ_i实为观测结果与模拟结果差矢量的模与其对应方位误差椭圆极半径之比，当δ_i≤1的时候，i点的模拟结果落于该点GPS观测误差椭圆内部.

对于整个模型的模拟效果可以采用对多个观测点的算数平均来评价，表达式为

(5)

由于本文构建的模型涉及断层较多，本文主要针对阿尔金断裂、广义海原断裂以及西秦岭北缘断裂这三条大型断裂带展开.其他断裂的剪切模量则根据前人研究(张东宁等, 2007; 郑勇等, 2007)，近似取为周边地壳介质剪切模量的1/10.为尽可能快速的取得较好的效果，本文首先从宏观角度分析区域主要活动断裂的剪切模量，重点考察阿尔金断裂(AEJ)、广义海原断裂(G_HY)和西秦岭北缘断裂(XQL)的剪切模量对区域模拟速度场的影响.对AEJ和XQL分别取λ=[0.05, 0.1, 0.5, 1]四个值，考虑He等(2013)研究指出海原断裂摩擦系数低，对G_HY取λ=[0.005, 0.01, 0.05, 0.1, 0.5, 1]六个值，将他们分别组合，共构建4×4×6=96组工况模拟计算各工况下的δ₀.结合模拟结果和GPS观测结果，利用式(5)择优选取最佳参数组合给出XQL、AEJ以及G_HY等效剪切模量.

其次考虑现有研究给出的广义海原断裂各段复杂性，继续保持其他断裂剪切模量为周边地壳介质剪切模量的1/10，AEJ和XQL有效剪切模量以前一步模拟结果为准.进一步根据历史地震、断层闭锁等研究，将广义海原断裂中的六盘山—宝鸡断裂划为六盘山段(LPS)，将天祝空段内的主要活动断裂(金强河断裂、毛毛山断裂和老虎山断裂)简化为金强河—毛毛山—老虎山段(JMLD)，将1920年海原大地震破裂的狭义海原断裂划为海原段(HYF)、冷龙岭及其以西的木里—江仓断裂简化为海原—冷龙岭—木里段(HLMD).基于GPS的断层闭锁研究显示JMLD段以及LPS段均处于强闭锁状态，因此在模拟中首先对JMLD和LPS分别取λ=[1, 0.1, 0.01, 0.001, 0.0001]，同时对HLMD取λ=[1, 0.1, 0.01, 0.001, 0.0001]，这样组合参数共构建5×5×5=125组工况，继续通过计算机模拟给出这125组模型的δ₀，结合GPS观测数据和模拟结果择优给出广义海原断裂不同分段有效剪切模量.

以阿尔金断裂、西秦岭北缘断裂以及广义海原断裂为重点考察对象设计的96组工况模拟结果如图 4a所示.阿尔金断裂(AEJ)和广义海原断裂(G_HY)力学参数变化对区域模拟速度场影响显著，而西秦岭断裂对区域模拟地壳运动速度场影响较小，但当前断层剪切模量也较高.保持图 4a中工况N1366的其他参数不变，降低阿尔金断裂剪切模量、提高广义海原断裂剪切模量，进而给出工况N0363.对比N1366和N0363模拟结果残差可知阿尔金断裂和广义海原断裂对青藏高原东北缘地壳运动速度场图像控制作用强烈(图 5a、d).这表明青藏高原东北缘地壳运动速度场图像强烈受控于阿尔金断裂以及广义海原断裂这两条大型走滑断裂，这一结果与地质研究结论一致，地质研究表明阿尔金断裂和祁连—海原断裂这两条大型左旋走滑断裂在区域构造运动转换过程中起着控制性的作用(苏琦等, 2016).

图 4

Fig. 4

图 4 断层剪切模量变化对模拟残差的影响 (a)阿尔金断裂、西秦岭北缘断裂以及广义海原断裂；(b)广义海原断裂不同段落断. Fig. 4 Simulated residuals of different cases that incorporate different fault shear modulus (a) AEJ, XQL and G_HY faults; (b) Different segments of G_HY fault.

图 5

Fig. 5

图 5 阿尔金断裂、广义海原断裂以及西秦岭北缘断裂不同剪切模量组合下的模拟残差空间分布 Fig. 5 Spatial distribution of simulation residuals from the cases that incorporate different fault shear modulus for AEJ, XQL and G_HY faults

阿尔金断裂是青藏高原的一条主控断裂，其左旋走滑及邻近区域的挤压褶皱作用吸收着印度板块向北推挤的能量(Molnar et al., 1987)，He和Chéry(2008)模拟结果显示阿尔金断裂摩擦系数低至0.06可以获得与地质结果一致的滑动速率，这一结论似乎与本文图 4a的认识相悖：本文模拟结果显示阿尔金断裂剪切模量与周边地壳介质剪切模量一致时，模拟残差最小.这是与阿尔金断裂左旋走滑速率在空间上呈现出不均匀分布有关，其0.06的有效摩擦系数主要针对阿尔金断裂的西南段.本文模型中涉及的96°E以东的阿尔金断裂滑动速率从95°E左右的2~6 mm·a^-1逐渐下降到金塔附近的1~2 mm·a^-1(郑文俊等, 2009).He和Chéry(2008)的研究还显示当摩擦系数选取0.2~0.6时在该段可以获得较好的模拟结果，表明阿尔金断裂96°E以东断层力学参数较强，与本文结果一致.保持残差最小的参数组合(图 5d)的N1366工况中其他断裂剪切模量不变，降低阿尔金断裂力学参数至λ=0.05，模拟结果与GPS观测结果残差如图 5c所示.相比图 5d最佳工况的残差分布，可以看出阿尔金断裂断层剪切模量对祁连带的地壳运动速度场图像影响显著，降低阿尔金断裂断层剪切模量使得祁连带的地壳缩短速度增大而偏离实际的GPS观测结果.因此祁连山带的GPS速度场图像表明阿尔金断裂96°E以东段落断层力学参数较高.

作为青藏高原东北缘另一条重要的大型走滑断裂带，广义海原断裂断层剪切模量的强弱对区域地壳运动速度场图像的控制作用同样非常强烈(图 4a)，模拟结果显示广义海原断裂的断层剪切模量极低，这一结果与He等(2013)的模拟结果一致，进一步证实广义海原断裂的弱介质力学性质对青藏高原东北缘地壳运动速度场图像有非常强烈的控制作用.保留图 5d的断层剪切模量不变，增大广义海原断裂断层剪切模量至λ=1.0，模拟结果与GPS观测结果残差如图 5b所示，结果显示广义海原断裂的断层力学属性不但影响着其附近的地壳运动速度场图像，甚至对远处的龙首山断裂等附近的地壳运动速度场图像也有一定程度的影响(图 5b).

考虑现有研究给出的海原断裂各段复杂性，继续保持其他断裂剪切模量为周边地壳介质剪切模量的1/10，根据上述分析AEJ和XQL断层力学参数取λ=1.0.以LPS、JMLD以及HLMD为讨论对象的125个工况模拟结果如图 4b所示.模拟结果总体上依然支持前面给出的广义海原断裂断层力学参数较弱.但六盘山断裂(LPS)却表现不同，降低LPS断层力学参数，模拟结果残差有增大的趋势，显示LPS当前断层力学性质较强.而同为高闭锁段的JMLD段的断层力学性质则比较弱，模拟结果残差随着JMLD断层剪切模量的降低而显著减小.综上取AEJ、XQL以及LPS的断层力学参数λ=1.0，取广义海原断裂金强河—毛毛山—老虎山段(JMLD)、狭义海原断裂(HYF)、冷龙岭段(LLF)以及木里—江仓段(MJD)断层力学参数λ=0.0001为最佳参数组合，给出模拟地壳运动速度与GPS观测速度对比结果如图 7所示.结果显示：整体上模拟地壳运动速度与GPS观测速度吻合良好.但狭义海原断裂(HYF)附近，GPS观测值减去模拟值的残差较大且在空间上表现出非常规律的顺时针旋转特征(王庆良等也给出类似的结果¹⁾)，集中区域与1920年海原大地震极震区(郭增建等, 1976)基本一致.会不会是狭义海原断裂断层剪切模量还不够低？在前面给出的最佳参数组合的基础上继续降低HYF段的剪切模量至λ=1.0×10^-7，测试结果如图 6a所示：继续降低HYF段断层剪切模量对模拟结果的影响并不显著，δ₀基本保持在1.486左右而不再降低.表明HYF附近规律的残差分布可能是真实的地壳运动信息，但与附近断层力学参数以及区域整体的块体运动关系不大.从剖面AB(图 2)GPS速度分解结果看(图 6b)，跨海原断裂(HYF)显著的左旋形态且靠近断层速度大、远离断层速度小，与震后短期地壳形变特征(Meade and Hager, 2005)的理论曲线非常相似，表明这可能与1920年海原大地震有关(Li et al., 2015b).具体机理还需要进一步深入研究.本文后续分析中HYF段断层剪切模量依然取前面给出的λ=0.0001.

图 6

Fig. 6

图 6 断裂剖面(AB)上最佳模型预测结果与观测结果对比(b) Fig. 6 Effect of HYF′s shear modulus on simulation residuals (a) and comparison of best fitted model prediction and the observed data (b) along the profile (AB) across the Haiyuan fault

图 7

Fig. 7

图 7 模拟速度场和GPS观测结果(a)及其整体无旋转场(b)对比 Fig. 7 Comparison of observed GPS data and simulation results (a) and their no-net-rotation fields (b)

1) 王庆良. 2015.青藏高原东北缘现今三维地壳运动监测研究.鄂尔多斯西南缘地震科学研究成果交流会.

作为印度—欧亚板块碰撞向北东方向扩展的最前缘，青藏高原东北缘受印度板块北推作用以及阿拉善地块和鄂尔多斯地块的阻挡作用, 区域内构造变形复杂.选择上述的最佳参数组合模型给出青藏高原东北缘模拟地壳运动速度场图像和实际观测GPS速度场对比(图 7a)，结果显示：除个别点之外，本文模拟结果在运动速度大小和方向上与观测结果基本一致.为更好凸显模拟结果和观测结果之间的差异，选取模型区域内GPS观测速度计算区域整体旋转欧拉极(ω=[-2.004525×10^-10, -3.898092×10^-9, 2.390402×10^-9]^T)，进而根据整体无旋转理论(杨国华等, 2005)，分别从模拟速度场和GPS观测速度场扣除模型区域的整体旋转，给出整体无旋转速度场对比结果(图 7b).结果显示除陇西地块内靠近狭义海原断裂的部分点之外，模拟结果与观测结果依然吻合良好.整体无旋转速度场图像清晰刻画了广义海原断裂从西至东“左旋逆冲—左旋走滑—左旋逆冲”的运动方式转化、祁连山带的挤压以及陇西地块的顺时针旋转.模拟结果与观测结果的高度吻合进一步证实了本文模拟结果的可靠性，保证了后续应力状态分析的可信度.

基于上述大量模拟给出的青藏高原主要活动断裂断层剪切模量，计算给出了青藏高原东北缘各主要活动断裂最大剪应力率以及区域主应力率场(图 8).模拟主应力率的空间分布显示祁连山带、六盘山应力环境表现为垂直挤压，阿尔金断裂、广义海原断裂以及西秦岭北缘断裂均表现为走滑为主、兼有少量挤压逆冲，与地震学(张辉等, 2012)方法以及He等(2013)数值模拟结果一致.由于青藏高原东北缘各主要活动断裂活动特征均以走滑为主，为此这里计算了各主要活动断裂上的最大剪应力率，结果如图 8所示.模拟结果显示：阿尔金断裂、西秦岭北缘断裂天水以西、广义海原断裂西段的木里—江仓段、中段的金强河—毛毛山—老虎山段以及东段的六盘山断裂均处于剪应力变化率高值状态.基于GPS给出的断层震间闭锁研究指出西秦岭北缘断裂西段和阿尔金断裂(M7专项工作组, 2012)、六盘山断裂和天祝空段(Li et al., 2016)当前滑动亏损较大.最大剪应力率与震间闭锁在空间上基本吻合，表明这些断裂段当前强闭锁程度高，且应力积累显著.同时这些高剪应力率区域与基于地震地质、地球物理、地震活动，大地测量等多种方法综合给出的潜在地震危险区(M7专项工作组, 2012)(图 1)在空间上有一定的对应关系，进一步表明这些区域未来地震危险性值得关注.此外德令哈附近的鄂拉山断裂、民勤附近的龙首山断裂以及岷县附近的临潭—宕昌断裂等段剪应力率也相对较高.

图 8

Fig. 8

图 8 最佳模型模拟主应力率场及断层最大剪应力率 灰色小圆是1970年以来青藏高原东北缘M≥3级以上地震，沙滩球是1970年以来青藏高原部分5级以上地震震源机制(张辉, 私人通讯). Fig. 8 Principal stress rates and maximum shear stress rates from the best fitted model Gray circles donate the M≥3 earthquakes in the northeastern Tibetan plateau since 1970; beach balls donate M≥5 earthquakes since 1970 and their focal mechanism solutions (Zhang Hui, personal communication).

本文在综合考虑青藏高原东北缘动力学环境、区域活动断裂空间展布以及下地壳黏滞性结构的基础上构建了青藏高原东北缘有限元动力学模型.通过模型边界附近的GPS观测速度内插给出模型边界条件，通过组合不同断裂的断层剪切模量共构建324组工况，模拟分析了不同断裂的断层力学参数对模型区域内地壳运动速度场图像的影响，对比不同工况下模型内模拟地壳运动速度与GPS观测地壳运动速度给出青藏高原东北缘各主要活动断裂断层剪切模量的最优取值.在此基础上分析了青藏高原东北缘各主要活动断裂应力状态.模拟结果表明：

(1) 阿尔金断裂和广义海原断裂断层控制着青藏高原东北缘地壳运动速度场图像，是青藏高原构造变形的主控断裂.但在本文所涉及的区域，二者的断层剪切力学性质相反.本文所涉及的阿尔金断裂95°E以东当前断层剪切模量高，对祁连山带地壳运动速度场图像有显著的影响；广义海原断裂带除六盘山断裂外，整体上断层剪切模量非常低，极低的断层介质力学性质强烈地控制着断裂带附近及周边的地壳运动速度场图像.六盘山断裂和西秦岭北缘断裂对模型所在区域地壳运动速度场图像影响不甚强烈，但根据本文模拟结果，其断层力学性质较强，对断层两盘块体相对运动有一定的阻碍作用.

(2) 断裂带较高的剪切模量表明其对断层两盘块体相对运动具有一定的阻碍作用，模拟结果显示：阿尔金北缘断裂东段、六盘山断裂以及西秦岭北缘断裂这些剪切模量较强的断裂当前剪应力率水平也相对较高，而剪切模量极低的狭义海原断裂剪应力率水平非常低.广义海原断裂中部的金强河—毛毛山—老虎山断裂以及西部的木里江仓段剪切模量较低，而剪应力率水平较高.剪应力率高值区与前人给出的青藏高原东北缘高闭锁区域吻合良好.动力学上的高应力率与运动学上的强闭锁良好吻合，预示着这些断裂当前滑动受阻，速率亏损较大，应力积累显著，是地震危险分析值得关注的区域.