巴颜喀拉块体位于青藏高原北部，是中国大陆现今地震活动最为强烈的地区之一.块体周缘发育了东昆仑断裂、鲜水河断裂、甘孜—玉树断裂和龙门山断裂等断裂带(邓起东等, 2010；Wu and Gao, 2019).自1900年以来，在巴颜喀拉块体内部及周缘共发生了20次7.0级以上地震，特别是自1997年玛尼M7.5地震以来，中国大陆所有的7级以上地震均发生在该区域.因此，研究巴颜喀拉块体边界断裂带应力演化及其与周缘强震的关系，对于理解强震间相互作用机制进而探讨强震的孕育及发生机理具有重要的意义.

强震的孕育、发生与断裂带应力状态密切相关.前人基于弹性位错理论、数值模拟等方法，开展了巴颜喀拉块体强震间相互作用及断裂带应力演化等方面的研究工作.单斌等(2012)利用黏弹性分层模型，研究松潘—甘孜块体东北端历史强震之间的相互作用以及地震活动对周边主要活动断裂的影响.李平恩等(2019)模拟给出1900年以来巴颜喀拉块体边界断裂带上7级以上强震所引起的断裂应力演化与强震之间的相互作用关系.从已有的针对块体边界断裂系统或单条断裂带上地震序列间相互作用的研究结果来看，多数强震之间存在着应力加载、甚至是触发的关系(万永革等, 2007; Xiong et al., 2010; 程佳等, 2011; 单斌等, 2012; 杨兴悦, 2012; Jia et al., 2012; Wang et al., 2016; 徐晶等, 2017; 李平恩等, 2019).同时，强震的发生也会造成周边断裂上库仑应力的明显改变(Toda et al., 2008; Parsons et al., 2008; Shan et al., 2009; Wan and Shen, 2010; Li et al., 2015; 徐杜远等, 2020).上述研究帮助我们深入认识块体边界断裂带上强震间的相互作用机制、断裂带应力演化，甚至是潜在强震危险地点的判定.然而，对于发生在块体内部的1947年达日M7.7强震在块体边界断裂带应力演化中的作用，在以往的研究中多被忽略.

达日断裂是巴颜喀拉块体内部的一条第四纪以来活动性较强的走滑断裂带(梁明剑等, 2020)，1947年M7.7达日地震就发生在该断裂带上.20世纪80年代，前人曾对该次地震地表破裂带和宏观破坏情况进行过初步调查研究(戴华光, 1983)，发现震中烈度为度，等震线长轴方向为北50°西，地震形变带在地貌上由一系列的陡坎、鼓包、错断冲沟和地裂缝等所组成，并在不同的地段表现出不同的破裂特点.

图 1

Fig. 1

图 1 巴颜喀拉块体断裂带展布与历史强震分布图(MS≥7.0, 1900—2020) 图中MS≥7.0历史强震(1900—2020年)震源机制解资料来自Molnar和Deng(1984), Global CMT Catalog (http://www.globalcmt.org/CMTsearch.html).EKLF:东昆仑断裂; DRF:达日断裂; GZ-YSF:甘孜—玉树断裂; XSHF:鲜水河断裂; LMSF:龙门山断裂; LRBF:龙日坝断裂; BF:马尔康断裂. Fig. 1 Fault configuration in the Bayan Har block and the distribution of historical strong earthquakes In which, focal mechanism solutions of the MS≥7.0 historic strong earthquakes during 1900—2020 are derived from Molnar and Deng (1984) and Global CMT Catalog (http://www.globalcmt.org/CMTsearch.html). EKLF: Eastern Kunlun fault; DRF: Dari fault; GZ-YSF: Garzê-Yushu fault; XSHF: Xianshuihe fault; LMSF: Longmen Shan fault; LRBF: Longriba fault; BF: Barkam block.

本文在前人研究的基础上，收集最新的地质、地球物理方面的基础资料，构建巴颜喀拉块体东部区域三维精细黏弹性有限元模型；以GNSS观测资料作为模型的初始边界条件，考虑同震位错与震后黏弹性松弛效应的影响，研究达日地震发生所引起的块体边界断裂带上的库仑应力变化，并探讨其与边界断裂带上历史强震间的相互作用，为区域地震危险性分析、尤其是潜在震源位置的判定提供力学参考.

结合巴颜喀拉块体东部区域活动断裂空间展布和达日地震的震中位置，确定研究区范围为91°—107°E，27°—37°N.模型中包含的断裂有：龙门山断裂、鲜水河断裂、甘孜—玉树断裂、东昆仑断裂等块体边界断裂以及达日断裂、马尔康断裂、龙日坝断裂等块体内部断裂(Burchfiel et al., 1995; Deng et al., 2003; 张培震等, 2008; 梁明剑等, 2020).其中，在模型中将倾向NW的龙门山断裂和龙日坝断裂分别设置成倾角为60°和70°的倾斜断层，作为简化，其余断裂带均设为直立断层，断层深度为20 km，各条断裂带的几何参数见表 1.同时，为了更好地模拟断裂活动，模型中所有断层均设置为接触摩擦单元，这样既可以模拟间断面，也能够模拟断层面上存在的接触摩擦机制.

表 1 (Table 1)

表 1 模型中主要活动断裂几何参数 Table 1 Geometric parameters of the major active faults in the model 断裂名称 断裂走向 断裂倾角 断裂性质 东昆仑断裂 NWW 90° 左旋走滑 达日断裂 NW 90° 左旋走滑 甘孜—玉树断裂 NW 90° 左旋走滑 鲜水河断裂 NW 90° 左旋走滑 龙门山断裂 N45°E 60° 逆冲兼右旋走滑 龙日坝断裂 N60°E 70° 右旋走滑 马尔康断裂 NW 90° 左旋走滑

表 1 模型中主要活动断裂几何参数 Table 1 Geometric parameters of the major active faults in the model

模型深度设置为80 km，包含上地壳(20 km)、中-下地壳(30 km)和上地幔(30 km)(图 2).模型单元尺寸在断裂附近控制在3 km内，其他区域控制在10 km内，共划分单元1064938个，节点190372个.

图 2

Fig. 2

图 2 巴颜喀拉块体东部三维有限元模型 Fig. 2 Three-dimensional finite element model of the eastern Bayan Har block and its adjacent areas

本文采用线性Maxwell流变学模型来模拟岩石圈的黏弹性行为，其本构关系为

其中，ε_Total为总应变，ε_D为黏性应变，ε_S为弹性应变，σ为应力，η为黏滞系数，E为杨氏模量.

基于岩石圈速度结构与密度结构的结果(https://igppweb.ucsd.edu/~gabi/litho1.0.html；Guo et al.，2019)，根据弹性模量与速度、密度及泊松比之间的关系，确定弹性模量等介质物性参数.其中，地壳内弹性模量的变化范围为(0.71~1.2) ×10¹¹ Pa，上地幔内变化范围为(1.54~1.76) ×10¹¹ Pa，华南块体弹性模量最高；泊松比为0.24~0.29.在模型中，上地壳设置为弹性介质，中-下地壳以及上地幔设置为黏弹性介质.同时，参考前人基于形变资料反演的岩石圈介质流变性质方面的研究成果(邵志刚等, 2008; He et al., 2013; Huang et al., 2014)，确定用于模拟震间构造应力变化的长期稳态黏滞系数和计算同震、震后应力变化的短期瞬态黏滞系数.并以地质体为基础单元，分别赋予对应的黏滞系数，具体见表 2.

表 2 (Table 2)

表 2 模型中各块体的黏滞系数 Table 2 Viscosity coefficient of each block in the model 柴达木块体/(Pa·s) 巴颜喀拉块体/(Pa·s) 羌塘块体/(Pa·s) 华南块体/(Pa·s) η1 η2 η1 η2 η1 η2 η1 η2 上地壳 - - - - - - - - 中-下地壳 1.0×1019 9.0×1018 1.0×1019 5.0×1017 1.0×1019 5.0×1017 1.0×1021 1.0×1020 上地幔 1.0×1020 1.0×1019 1.0×1020 1.0×1019 1.0×1020 1.0×1019 1.0×1022 1.0×1020 注：η1和η2分别表示稳态黏滞系数和瞬态黏滞系数.

表 2 模型中各块体的黏滞系数 Table 2 Viscosity coefficient of each block in the model

本文使用中国大陆1991—2015年的GNSS速度场作为有限元模型的边界约束(Zheng et al., 2017)，参照石富强等(2018)在青藏高原东北缘的研究工作，利用最小二乘配置方法(Wu et al., 2017)将离散的GNSS速度场插值到模型边界上(图 3a).由于地壳运动速度在垂向的差异尚未明确定论，作为一种近似，本文所施加的边界条件不随深度变化，且模型底部水平向自由，法向约束.模型计算使用Pylith有限元软件(https://geodynamics.org/cig/software/pylith/)进行计算.

图 3

Fig. 3

图 3 GNSS实测值与速度场模拟值结果对比(a)及残差分布(b) Fig. 3 Comparison between the observed GNSS and the modeled velocities (a) and residual distribution (b)

为了检验数值模型所用参数的合理性，将上述边界条件加载下的模型运行至平衡状态，分析平衡状态下速度场特征，并将模拟给出的速度场与GNSS地表观测结果进行对比分析.可以看出，在现有的GNSS边界约束和介质物性参数取值下，在大多数观测点，模拟值和观测值在方向和大小上拟合较好(图 3a)，残差小于2 mm·a^-1(图 3b).对于部分点存在较大差异的原因，分析认为可能与局部区域的特殊构造有关.模拟结果基本反映了研究区域的现今构造变形和应力状态，将此应力场作为研究区背景初始应力场，通过加载地震的同震位错量，模拟达日地震的发生.

本文计算库仑应力变化ΔCFS采用如下公式(King et al., 1994)：

其中，Δτ为断层面上剪应力变化，与断层滑动方向一致时为正；Δσ_n为断层面上正应力变化，张性为正.μ′为断层面上的有效摩擦系数.有效摩擦系数在不同研究中取值有所差异，通常在0.2~0.8之间，本文中有效摩擦系数取值为常用取值0.4(Stein et al., 1992; King et al., 1994; Freed, 2005).

由于各断裂带的几何参数并不相同，在计算不同断裂带上的库仑破裂应力时，需要将应力张量投影到对应的断层面上.同时，在计算库仑应力时除了需要知道断裂的走向和倾角外，还需要知道断裂的滑动角，本文中将东昆仑断裂、甘孜—玉树断裂、鲜水河断裂、达日断裂和马尔康断裂的滑动角设为0°，龙门山断裂的滑动角设为90°，龙日坝断裂滑动角设为180°.利用Pylith模拟可以得到任意点应力张量的6个独立分量，然后根据断层参数进行应力张量投影，这样就可以计算任一走向，任一倾向和任一滑动方向断层面上的正应力σ_n和剪应力τ，从而计算断层面上的库仑应力变化，具体转换方法参考石富强等(2018)的计算方法.

震间构造加载作用可能从更长时间尺度影响断层的应力积累，是地震孕育、发生的重要因素，对断层面上的库仑应力分布产生影响，使得地震发生后断层面上的应力状态更为复杂(Stein, 1999; Freed, 2005; Luo and Liu, 2010).为了研究震间构造加载方面的影响，基于图 3边界约束下的模型，我们计算了震间构造加载作用下各断裂带上的库仑应力年累积速率(图 4).可以看出，东昆仑断裂中西段、鲜水河断裂的库仑应力年累积速率分别为286 Pa·a^-1和267 Pa·a^-1，这种高水平的库仑应力年累积速率与断裂带历史强震活动频发，即特征地震复发周期较短是一致的；东昆仑断裂中段到东段的年累积速率分别为194 Pa·a^-1、162 Pa·a^-1和103 Pa·a^-1；达日断裂和甘孜—玉树断裂年累积速率分别为206 Pa·a^-1和181 Pa·a^-1；而龙日坝断裂、马尔康断裂以及龙门山断裂的库仑应力年累积速率相对较小，基本在120~160 Pa·a^-1之间，较小的累积速率说明这三条断裂所在的区域应力积累水平低，发生强震需要较长时间的应力积累.

图 4

Fig. 4

图 4 断裂带震间期构造应力年累积速率 Fig. 4 Annual accumulation rate of interseismic tectonic stress on faults

根据断层破裂长度、滑动量与地震震级之间的经验公式(邓起东等, 1992)，确定达日地震的破裂长度、水平错动量分别为252 km和2.91 m.通过对达日断裂发震段施加同震位移载荷，计算了1947年达日地震引起块体周边断裂带上的库仑应力变化(图 5a).可以看出，达日地震的发生除造成破裂段两端相邻未破裂段显著的库仑应力增加外(＞1.0×10⁶ Pa)，同时也造成了龙日坝断裂中-北段、东昆仑断裂中西段同震库仑应力明显增加，应力增加量多为(3~8)×10³ Pa，最大达到1.4×10⁴ Pa；而在龙门山断裂中-北段、马尔康断裂、鲜水河断裂中-南段和甘孜—玉树断裂西段上库仑应力增加较小.相反，东昆仑断裂中-东段、甘孜—玉树断裂中-东段、龙日坝断裂南段、鲜水河断裂北西段和龙门山断裂南段则表现为库仑应力减小，其中在东昆仑断裂中段上库仑应力减小最为明显，最大达到4.4×10⁴ Pa.

图 5

Fig. 5

图 5 达日地震引起的周围主要断裂带库仑应力变化 (a)同震应力变化；(b)震后应力变化；(c)同震与震后累积应力变化(截止2021年). Fig. 5 Coulomb stress change on the surrounding major faults caused by the Dari earthquake (a) Co-seismic stress change; (b) Post-seismic stress change; (c) The combined (co- and post-seismic) stress change by 2021.

同时，考虑震后黏弹松弛效应的影响，进一步计算了达日地震发生后到2021年由于中-下地壳和上地幔黏弹性松弛效应所引起的周边断裂的库仑应力变化(图 5b).可以看出，相比同震库仑应力变化，由震后黏弹性松弛效应引起的周边断裂带应力增加较大，大部分断裂带上库仑应力变化都大于1.0×10⁴ Pa.综合同震位错与震后黏弹性松弛效应的影响(图 5c)，龙日坝断裂中-北段、东昆仑断裂中西段、龙门山断裂中-北段、马尔康断裂、鲜水河断裂中-南段和甘孜—玉树断裂西段均表现为库仑应力的增加，其中龙日坝断裂中-北段和东昆仑断裂中西段的部分区域库仑应力变化量超过了1.0×10⁴ Pa，在龙门山断裂中-北段库仑应力变化量接近1.0×10⁴ Pa；而在东昆仑断裂中-东段、甘孜—玉树断裂中-东段、龙日坝断裂南段、鲜水河断裂北西段和龙门山断裂南段存在不同程度的库仑应力减小，在东昆仑断裂中段减小量值最大，达到1.8×10⁵ Pa.

库仑应力变化是研究强震间是否存在触发关系的重要参考，但其大小受断层有效摩擦系数取值的影响.为分析有效摩擦系数选取对库仑应力计算结果的影响，我们取有效摩擦系数为0.2、0.4、0.6和0.8，分别计算达日地震造成后续地震震源处的库仑应力变化(图 6).可以看出，对于达日地震而言，摩擦系数的取值会在一定程度上改变库仑应力的大小，但不会影响计算结果的极性.

图 6

Fig. 6

图 6 有效摩擦系数取不同值时后续地震震源处库仑应力变化 Fig. 6 Coulomb stress change at the epicenter of subsequent earthquakes with different effective friction coefficients

基于上述研究，我们选取有效摩擦系数为0.4的结果，进一步讨论达日地震与后续6次强震间的可能关系.可以看出，达日地震造成了1955年康定、1963年阿兰湖和2008年汶川3次强震震源处累积(同震+震后)库仑应力的增加(图 6)，分别增加了3.58×10³ Pa、1.16×10⁴ Pa、7.68×10³ Pa，其中阿兰湖地震前库仑应力变化达到了传统上所认识的地震触发的阈值0.1 bar(10⁵ Pa)(King et al., 1994; Harris, 1998)，表明达日地震可能触发了阿兰湖地震的提前发生，而康定地震和汶川地震库仑应力变化虽未达到触发的阈值，但达日地震对这些强震的发生仍起到了促进的作用.结合震间构造应力年累积速率结果(图 4)，达日地震造成1955年康定、1963年阿兰湖和2008年汶川地震的发震时间分别提前了约13.4年、40.6年和49.5年.另外，达日地震也造成了1976年炉霍、2001年玉树和2013年芦山3次强震震源处库仑应力的减小，分别减小了2.44×10⁴ Pa、3.06×10⁴ Pa、1.06×10⁴ Pa，说明达日地震推迟了这三次地震的发震时间.

库仑应力变化能够定量描述地震对邻区断层系统的影响，可以作为潜在危险性地点判定的重要依据，而被广泛用来研究强震时空迁移机理与断裂强震危险性分析(Tabrez et al., 2008; Shao et al., 2016; Xiong et al., 2017).例如，Stein等(1997)和Nalbant等(1998)指出土耳其北安纳托利亚断裂北西段明显的应力增加，具有较高的地震危险性.1999年Izmit7.4级地震的发生印证了这一认识(Hubert-Ferrari et al., 2000)，进一步推动了库仑应力变化在潜在强震危险地点判定方面的应用.

基于地震地质、大地测量、地球物理等多种基础资料，中国大陆地震大形势工作组(简称：大形势工作组)给出2021—2030中国大陆地震重点危险区的判定结果(引自2021—2030年中国大陆地震重点监视防御区确定工作报告)，其中有四个危险区位于研究区内(图 7)，分别是东昆仑断裂带西大滩—东大滩段、东昆仑断裂带东段—龙日坝断裂带中北段、风火山—玉树当江段和鲜水河断裂带中南段—龙门山断裂带南段.针对这些危险区内断裂的应力演化，大形势工作组利用PSGRN/PSCMP程序(Wang et al., 2006)计算了断裂的应力演化过程，可以看出危险区内多数断裂的应力变化大于1.0×10⁴ Pa，断裂应力变化的结果对于危险区有很好的支持作用.同时，从已有的巴颜喀拉块体边界断裂带上库仑应力变化的研究结果来看(Xiong et al., 2010; 杨兴悦, 2012; Shan et al., 2015; 徐晶等, 2017)，四个危险区与对应断裂的库仑应力变化有很好的对应关系，基本都表现为危险区内库仑应力加载，而在危险区外的断裂库仑应力变化则大多表现为卸载.但是，需要说明的是，前人的结果多数没有考虑块体内部达日地震对周边断裂带库仑应力变化的影响.而从本文中达日地震对相邻断裂库仑应力变化的结果来看(图 7)，东昆仑断裂中西段库仑应力增加了1.34×10⁴ Pa，与大形势工作组认为应力显著增加的东昆仑断裂西大滩—东大滩段相对应；鲜水河断裂中-南段库仑应力增加了4.37×10³ Pa，与应力显著增加的鲜水河断裂磨西段相对应；甘孜—玉树断裂西段库仑应力增加了2.13×10³ Pa，与应力显著增加的甘孜—玉树断裂当江段相对应；龙日坝断裂中-北段库仑应力增加了4.78×10⁴ Pa，与应力显著增加的龙日坝中-北段相对应，以上库仑应力存在显著增加的断裂与前人的研究结果基本一致，也与危险区有很好的对应关系，表明达日地震的发生使得这些断裂的危险性程度进一步提高.

图 7

Fig. 7

图 7 断裂带库仑应力变化与地震危险区分布 Fig. 7 Coulomb stress change on faults and distribution of seismic hazard zones

但是，也存在与大形势工作组给出的应力结果不一致的情况，比如东昆仑断裂中-东段表现为应力卸载效应.为深入认识达日地震造成东昆仑断裂中-东段库仑应力卸载的原因，我们给出了东昆仑断裂中-东段断层面上的正应力与剪应力变化.从图 8可以看出，断裂上剪应力变化以卸载为主，最大可达1.0×10⁵ Pa；正应力变化则显示断裂中东段由西边的挤压逐渐向东变为拉张作用，而断裂东段主要是以拉张为主，应力变化大部分都小于3.0×10⁴ Pa.对比断层面上剪应力与有效正应力的变化可以看出，东昆仑断裂中-东段表现为应力卸载的原因主要是由断层面上剪应力卸载造成的.对于滑动性质相同的两条近似平行断层，达日地震的发生会造成东昆仑断裂中-东段应力的释放.这与马瑾等(2002)通过物理模型试验得到的结果一致，即：滑动方向相同的平行断层之间、同一断层带上型阶相同的雁列式断层之间的关系可视作串联关系，其中一条断层的失稳会使另一条断层上的剪应力减小，从而起到减震作用.同时，也与基于数值模拟方法给出的认识一致(Wang et al., 2017; Li et al., 2020).Shan等(2015)对东昆仑断裂带玛沁—玛曲段的应力演化及地震危险性研究时，认为达日地震释放了东昆仑断裂上的应力积累，而前人的研究结果中认为该段库仑应力增加的多是没有考虑达日地震的影响，造成了东昆仑断裂库仑应力的偏高.因此，在研究断裂带库仑应力演化与开展地震危险性分析时，块体内部强震的影响不可忽略.

图 8

Fig. 8

图 8 东昆仑断裂中-东段正应力和剪应力变化 Fig. 8 Normal and shear stress changes in the middle-east section of Eastern Kunlun fault

本文综合考虑地质、地球物理等基础资料，构建巴颜喀拉块体东部区域三维精细黏弹性有限元模型.以1991—2015年GNSS观测资料作为模型边界约束条件，考虑达日地震同震位错与震后黏弹性松弛效应的影响，模拟给出达日地震发生所引起的巴颜喀拉块体东部周缘断裂带上的库仑应力变化，并探讨其与块体边界带上历史强震间的相互作用，为区域地震危险性分析、尤其是潜在震源位置的判定提供力学参考.研究结果显示：

(1) 达日地震除造成达日断裂破裂段两端相邻未破裂段显著的库仑应力增加外(＞1.0×10⁶ Pa)，同时造成龙日坝断裂中-北段、东昆仑断裂中西段等同震库仑应力的明显增加，应力增加量多为(3~8)×10³ Pa，最大达1.4×10⁴ Pa；而在龙门山断裂中-北段、马尔康断裂、鲜水河断裂中-南段和甘孜—玉树断裂西段上库仑应力增加较小.相反，在东昆仑断裂中-东段、甘孜—玉树断裂中-东段、龙日坝断裂南段、鲜水河断裂北西段和龙门山断裂南段则表现为库仑应力减小，其中东昆仑断裂中段库仑应力减小最为明显，卸载量最大达到4.4×10⁴ Pa.因此，在研究断裂带库仑应力演化与开展地震危险性分析时，块体内部强震的影响不可忽略.

(2) 考虑同震和震后黏弹性松弛效应的影响，发现达日地震导致1955年康定7.5级地震、1963年阿兰湖7.0级地震和2008年汶川8.0级地震震源处库仑应力的增加，促进了后续3次强震的提前发生.结合震间构造加载应力年累积速率的结果，认为达日地震引起的库仑应力加载效应使得这3个地震分别提前约13.4年、40.6年和49.5年.

(3) 达日地震引起的库仑应力增加的断裂带与大形势工作组认为应力显著增加的东昆仑断裂西大滩—东大滩段、鲜水河断裂磨西段、甘孜—玉树断裂当江段和龙日坝中-北段具有较好的对应关系，使得这些断裂的危险性程度进一步提高.