川滇地区地震危险性数值分析

尹迪， 董培育， 曹建玲， 石耀霖

1 中国科学院大学地球与行星科学学院， 北京 100049 2 中国地震局地震研究所地震大地测量研究室， 武汉 430071 3 中国地震局地震预测研究所， 北京 100049

0 引言

川滇地区位于青藏高原东南缘，受印度板块与欧亚板块碰撞挤压作用，青藏高原内部物质被侧向挤出，而在高原东缘受古老而坚硬的华南板块的阻挡，其东向运动发生偏转，物质逐渐向南“流出”，形成围绕Assam构造结顺时针旋转的运动形态，GPS观测资料也揭示了这一运动特征(Tapponnier et al., 1982; Royden et al., 1997; 邓起东等,2002; Zheng et al., 2017; Jin et al., 2019)，如图1所示. 由此造成了川滇地区复杂的构造环境，发育了大量活动断裂带，主要包括高原东边界的龙门山逆冲推覆断裂带，构成川滇菱形块体东北边界的甘孜—玉树—鲜水河断裂带和安宁河—则木河—小江走滑断裂带，以及菱形块体西南边界的巴塘断裂带、金沙江断裂带和红河断裂带，滇西南的龙陵—瑞丽断裂带，澜沧江断裂带等(徐锡伟等,2003；吴中海等,2015).该区域构隶属于中国大陆南北地震带的南段，造活动强烈，地震活动水平高，为地震活动性科学研究提供了丰富的资料，2014年中国地震局在该区域建立了川滇地震预报实验场(吴熙彦等,2018).

据历史地震目录不完全统计，该区域1900—1999年以来共记录到M>6.0地震117次(苏有锦和秦嘉政,2001).其中主要的大震有1923年炉霍MS7.3地震，1970年通海MS7.8地震，1973年炉霍MS7.7地震，1976年5月29日在龙陵地区连发两次震级分别为MS7.3、MS7.4的强震，随后1988年11月6日在澜沧—耿马地区又连发两次强震，震级分别为MS7.2 、MS7.6.21世纪以来，该区域又陆续发生了数次强震，如2008年汶川MS8.0大地震，是罕见的内陆大震，造成了龙门山断裂带东北段长达360 km的破裂(徐锡伟等,2010)，随后又发生了2010年玉树MS7.1地震，2013年龙门山断裂带南段的MS7.0芦山地震，以及2017年MS7.0九寨沟地震.详细历史地震信息如图1和表1所示.

该地区地震危险性数值分析工作一直是地震学领域内研究热点和难点问题之一.前人的相关研究已取得了一定有意义的成果，但对地震孕育发生的物理基础和动力学过程的考虑尚不充分(陈连旺等,2008；王芃等,2019；程佳等,2020；李玉江等,2020).张培震等(2009)指出地震的孕育和发生是一个非常复杂的地球动力学过程，需要对已有的有限观测数据信息展开综合性研究才能获取更加科学可靠的结论.石耀霖等(2018)提出目前应当充分利用现有观测资料，发展地震数值预报工作，虽然地震孕育和发生的物理机制复杂，但在现有条件下，可以将地震视为地壳岩石圈孕震层岩石在构造活动作用下，岩石应力持续积累并储存能量直至破裂极限，造成断层发生错动释放应力的过程.该过程中不仅需要了解区域构造应力增长速率及岩层破裂极限，还需要对岩层的初始应力场有大致的估计.董培育等(2019)提出了一种简单的初始应力场估算方法，主要根据岩石库仑-摩尔破裂准则反推，在某一大地震发生前，震源区的应力状态必然达到了岩石承受极限值，震源区岩石强度值可通过实验室观测及其他相关资料分析得到，据此可以得到震源区相对准确的应力场状态.对于历史上没有地震的地方，推测其初始应力场尚未达到岩石应力承受极限值，据此可以对震前应力场有所估计.但是该方法仅仅考虑水平最大主应力和最小主应力，静水压力采用均匀的静岩压力值，结果显示对走滑型应力场约束较好，而对逆冲和拉张型应力场约束较差.

图1 区域构造图 其中插图中黄色箭头表示GPS观测速度场，数据来源于(Zheng et al., 2017)；白色箭头表示块体运动方向； 红色沙滩球表示1904—2017年间的MS>6.5强震，同表1.Fig.1 Tectonic setting The yellow arrows show the GPS velocity field (Zheng et al., 2017); The white arrows represent the direction of block movement; The red beach balls display the MS>6.5 earthquakes occurred between 1904 and 2017.

表1 川滇地区自1900年以来30次MS6.5级以上强震目录和震源参数Table 1 Earthquake catalogue and source parameter of earthquakes with MS≥6.5 in Sichuan-Yunnan region since 1900

基于此考虑，本文以构造应力场时空变化特征与强震的关系为切入点，为了模拟出更加准确的初始应力场，优化了应力场的反演方法.根据王辉等(2019)川滇地区应力场的应力形因子反演结果，约束区域应力场类型，建立二维弹性有限元平面模型，结合非均匀分布的静水压力，即非均匀的垂向应力，通过反演给定区域特定时刻合理的初始准三维应力场.在此基础上，利用董培育等(2020)提出的方法，综合考虑地震孕育阶段和震后调整阶段的动力学过程，以孕震层区域应力状态时空变化特征为研究对象，根据岩石破裂的物理机制——库仑-摩尔破裂准则作为判断地震发生的条件，模拟单次地震过程和多次地震循环过程.同时，对于数值模拟过程中的不确定性成分(区域初始应力场)，通过大量随机试验计算得到5000种不同的应力场结果，对不同模型的结果进行概率统计得到较为可靠的结果.本文是对董培育等(2019,2020)方法的进一步完善，以及在川滇地区开展的又一次有科学意义的探索实践.

1 有限元模型建立

考虑川滇地区地震多发生在上地壳脆性域内，同时考虑节约计算量，本文仍选择二维弹性有限元模型，暂且忽略中下地壳粘滞性影响.计算区域范围为：96°E—106°E，21.5°N—35°N.研究区域西南角位于Assam构造结附近，一方面该区域的GPS观测点较稀疏，另一方面有限的观测速度值表明该区域的运动方向存在南北向的大偏转(Zheng et al., 2017)，会对模型边界条件的选取造成很大误差，因此为了避免该问题，模型的西南角边界大致选择以国境线为界，如图2所示.模型中断裂带假定为具有一定宽度的软弱条状带，其中龙门山断裂带设为10 km宽，其他断裂带均为5 km宽.模型共剖分119038个三角形网格，235056个节点，断裂带上加密，网格边长平均约0.5 km，三角形面积约0.2 km2.Wang等(2020)研究成果认为数10年时间尺度的大地测量资料能够很好反应岩石圈区域百年尺度的地壳运动特征，因此选择2009—2015年间GPS观测平均速度场数据(Zheng et al., 2017)插值到模型边界上，作为模型的边界条件，且假设百年尺度范围内边界值不变.由于川滇地区是青藏高原物质向东-东南向流动的重要通道，研究表明其软弱的中下地壳物质流动速率比脆性的上地壳更快，对上地壳有一定的拖曳作用(Royden et al., 1997; 朱守彪和石耀霖,2004;董培育等,2016)，因此仅考虑弹性变形是不够的，还需要考虑下地壳的拖曳作用，拖曳力的选取采用尹迪等(2021)文中的差值拟合方法得到，模型材料参数(杨氏模量和泊松比)也同该文一致，在此不予赘述，具体请参考该文.在此基础上求得区域的构造应力加载速率如图2b所示.

图2 川滇地区有限元模型(修改自尹迪等，2021) (a) GPS位移边界条件; (b) 年构造应力加载速率.Fig.2 Finite element model of Sichuan-Yunnan region (Modified from Yin et al., 2021) (a) GPS boundary condition; (b) Tectonic stress loading rate.

图3 应力场演化示意图(改自：石耀霖等， 2018; 董培育等，2020) σr 历史无震区的初始应力上限值; σs 断层强度； σ0 震前 初始应力值; 应力增长速率; Δσ 应力扰动.Fig.3 Schematic diagram of stress field evolution (Modified from: Shi et al., 2018; Dong et al., 2020) σrThe upper limit of the initial stress of the historical earthquake-free zone; σs The strength of the fault; σ0 The initial stress value before the earthquake; The stress growth rate; Δσ The stress disturbance.

2 计算方法

初始应力状态的计算是复杂且重要的，本文参考董培育等(2019，2020)提出的方法，在此基础上进行方法的优化然后进行初始应力场估算.主要依据岩石库仑-摩尔破裂准则，将区域准三维应力场状态转换为最大主应力Pmax与最小主应力Pmin之间的比值关系来描述(根据不同构造区的应力状态，赋予不同方向的主应力值)，如式(1)所示，其中μ为岩石内摩擦系数，r0为理论破裂临界主应力比值，r为实际主应力比值，r≥r0，即rd≥0，表示发震临界条件；rd<0，表示安全，设定rd为危险度系数(riskdegree, 简写rd).以该准则作为发震条件的判断依据，那么反推在某一大地震发生前，震源区的应力状态必然达到了岩石承受极限，即r≥r0，r可由区域岩石内摩擦系数计算得到，并据此可约束Pmax和Pmin的大小关系，再利用垂向应力约束，可估算震前应力场的大致状态.对于历史上没有地震的地方，推测其初始应力场尚未达到岩石应力承受极限值，且断层活动区域高于内部稳定块体，因此在合理的范围内利用Monte Carlo方法随机给出：

(1)

(2)

通过该方法反演的初始应力场仅能与走滑型应力场拟合较好，但未能反映出拉张和逆冲型应力场的特征，因此本文在此基础上进行改进，优化了初始应力场的反演方法，具体如下.

2.1 根据应力形因子约束应力状态

R=(σ1-σ2)/(σ1-σ3).

(3)

不同的构造应力场对应的差应力状态不同，因此本文首先根据应力形因子来约束研究区域逆冲和拉张构造区的应力状态.根据库仑-摩尔破裂准则，控制剪切破裂的主要因素为最大主应力与最小主应力之间的差应力，即(σ1-σ3).忽略中间主应力，仅考虑最大主应力(σ1)与最小主应力(σ3)之间的关系，那么在σ1保持不变的情况下，σ3越小，差应力越大，R越小，越接近逆冲型应力状态；反之，R越大，越接近拉张型应力状态.

王辉等(2019)根据川滇地区1936年以来的数百个历史中小震震源机制解反演了区域地壳应力场状态，高原东南缘主要以剪切应力状态为主，龙门山地区以挤压型应力为主，川西高原则主要表现出拉张型应力状态(图4a)，将该结果插值即可达到整个研究区域的应力形因子分布(图4b).需要说明的是，在模型西北部(巴颜喀拉块体东部)和东南部(滇东南)区域，历史地震震源机制解相对川滇地区较少，没有该区域的应力形因子数据，因此该区域插值结果可靠性相对较低.尽管如此，结合现今研究考虑这两个地区均以剪切走滑型运动为主，因此本文将该处设定为走滑型应力状态.

2.2 有效垂向应力

董培育等(2019，2020)在反演模拟区域应力场时，将模拟的二维平面应力场叠加上岩石自重(静岩压力)Pv视为准三维应力场，设区域中的静岩压力为均匀恒定值.在理想状态下，忽略地壳介质中的孔隙流体压力，可认为任意深度处垂向正应力(压应力)等于上覆岩石的重力，也视为静岩压力.在地壳浅层地应力测量中，通常也采用该原则，将垂向正应力视为岩石自重(陈群策等,2004,2019；王成虎等,2012,2014).尽管董培育等(2019，2020)将此作为一级近似的计算，反演的应力场在走滑型构造带中与实际应力状态较为吻合，结果较为可靠，但在逆冲与拉张构造带中应力场与实际情况误差较大.

据Anderson(1951)断层理论可知，在压缩变形带中，最大主压应力为水平向，最小主压应力为垂向σv；拉张变形带中，最大主压应力为垂向σv，最小主压应力为水平向，由此可知垂向应力在逆冲和拉张变形带中起主要作用(图5).因此前期研究中给定均衡的垂向应力，势必会导致反演的逆冲和拉张变形带中的应力状态与实际应力状态不符.实际上，垂向应力不仅是上覆岩石的自重，地壳岩石圈层流体压力不可忽略，因此有效垂向应力应为岩石自重减去流体压力作用，如公式(4)所示：

σv=(ρr-ρw)gh，

(4)

其中ρw为流体密度.不同区域岩层内流体含量不同，流体压力大小也不同，尤其是逆冲型和拉张型构造带中，逆冲型构造带中的流体压力最大可等于上覆岩石重力(刘贵，2020)，与走滑型构造带中的流体压力差异较大.流体压力增大，可促进岩石发生破裂，在局部地区6 km深度处流体压力甚至能达到岩石重力(易立新等，2003)，因此流体压力不可忽略.

图4 应力形因子分布图 (a) 区域震源机制解反演得到的应力形因子分布图; (b) 应力形因子插值云图.Fig.4 Distribution of stress shape ratio (a) The distribution map of the stress shape ratio obtained by the inversion of the regional focal mechanism solution; (b) Interpolation contour of stress shape ratio.

根据库仑-摩尔破裂准则，挤压变形带产生破裂所需要的应力状态和差应力均大于拉张变形带，在有流体存在的情况下，如图5所示，摩尔圆整体向左偏移，逆冲变形带中发生破裂所需要的流体压力远大于拉张变形带(Sibson and Scott,1998).因此本文将利用应力形因子约束区域相对应力状态，考虑可能的流体压力效应，给定非均匀分布的有效垂向应力值σv，据此可给出区域准三维应力状态，并能更加准确地反映出研究区域不同构造区的主应力状态特征.

图5 存在流体状态下主应力分布示意图 (a)、(c) 拉张应力状态下的摩尔圆，σ1=σv即最大主应力为垂向应力; (b)、(d) 压缩应力状态下的摩尔圆， σ3=σv即最小主应力为垂向应力; (a)、(b) 完整岩石发生破裂摩尔圆; (c)、(d) 已有断层的重新滑动摩尔圆.Fig.5 Schematic diagram of principal stress distribution in the presence of fluid (a,c) The Mohr circle under tensile stress, σ1=σv, the maximum principal stress is the vertical stress; (b,d) The Mohr circle under compressive stress, σ3=σv, the minimum principal stress is the vertical stress; (a,b) The Mohr circle of complete rock fracture; (c,d) The Mohr circle of re-sliding of existing fault.

图6 反演得到的初始区域应力场沙滩球及关键点位 5000种主应力直方图 其中，σ1、σ2为本模型估算的最大、最小水平主应力; σv为有效垂向应力.Fig.6 The one possible initial regional stress value represented by beach balls and the histograms of 5000 kinds of stress value at three key points Among them, σ1, σ2 are the maximum and minimum horizontal principal stresses respectively; σv is the effective vertical stresses estimated by this study.

图7 历史地震序列断层单元平均危险度系数演化示意图Fig.7 Schematic diagram of the evolution of the average risk of fault units in historical earthquake sequence

3 计算结果

3.1 应力场反演结果

反演得到川滇地区15 km深度面上1904年震前初始应力场的偏应力张量，利用震源球表示应力状态，如图6所示.模拟得到的地壳应力状态(灰色震源球)与研究区历史地震(红色震源球)的震源机制解所揭示的应力状态基本一致，整体表现为川滇地区NW-SE方向的挤压以及菱形块体中部部分NE-SW拉张优势性质.为了进一步验证结果的可靠性，本文选取三个关键点位进行分析，其中1号点位位于2008汶川地震及2013芦山地震的地震空区，2号点位位于金沙江断裂带与丽江小金河断裂带交界处，3号点位位于龙陵瑞丽断裂带拐角处.图6a、6b、6c分别为1号点位、2号点位和3号点位5000次初始应力场反演结果(σ1、σ2、σv)分布直方图，可以看出，其中1号点位的最大最小水平主应力及有效垂向应力均值分别为294 MPa、115 MPa和95 MPa，即σ1>σ2>σv，表现出挤压应力场的特征；2号点位的最大最小水平主应力及有效垂向应力均值分别为120 MPa、86 MPa和265 MPa，σv>σ1>σ2表现出拉张的应力特征；3号点位最大最小水平主应力及有效垂向应力均值分别为179 MPa、70 MPa和106 MPa，σ1>σv>σ2表现为走滑应力场的特征，由此得知选取的三个关键点位反演得到的应力性质与区域实际地震表现出的应力状态一致.据此认为反演的应力场结果在局部地区与实际状态吻合较好，能够反映出不同构造区的应力特征，尤其是不同构造区垂向应力的非均匀分布特征更为明显.

尽管该方法的计算结果较之前的方法有所改进，但仍然存在误差.如龙门山断裂带具有明显分段构造差异，西南段以逆冲挤压为主，东北段以右旋走滑为主，应力形因子如图4所示，表现出相同的应力场特征.然而大范围区域模拟结果则显示龙门山断裂带均为走滑应力状态，未能体现出西南段的逆冲应力状态.主要原因在于龙门山断裂带位于巴颜喀拉块体及华南块体分界，两大块体内部均体现为走滑应力状态，因此在龙门山地区可能受插值影响(图6)，未能显示局部的逆冲特征.同时阿萨姆构造结处应力场模拟与实际也存在一些差异，主要是受计算边界效应的影响.基于二维模型的局限性，认为这些差异的存在是基本合理的，本文结果与已有的区域构造应力场的研究结果(许忠淮等,1987;Zhao et al.,2013;王晓山等,2015;Tian et al.,2019)基本一致.

3.2 地震序列发展过程

从反推得到的研究区域1904年初始应力场开始正演计算，该过程中同时考虑震间构造应力加载效应及前序地震的同震应力扰动.历史地震的发震条件，采用库仑应力破裂准则为约束，即当断层单元超过断层破裂极限时(转换为rd≥0)视为 “地震”发生的条件.同震应力扰动采用横向各项同性“杀伤”单元法计算(董培育和石耀霖，2013)，即通过弱化破裂单元的垂直断层方向的剪切模量，实现同震变化的计算(暂不考虑震后长期应力松弛效应).依次完成30次历史强震序列发展过程的模拟，通过展现历次大地震震源区的平均危险度值随时间变化的曲线，来表示历史地震的发展过程.在前述5000种不同的初始应力场条件下均一一展开正演计算，图7为其中一个模型的30次历史地震破裂区平均危险度系数值随时间变化曲线，正演模拟基本复现了川滇地区历史强震的有序发生.对于所有历史地震，在其各自地震发生前，破裂区平均危险度系数基本处于零线以下，但维持较高水平.随时间推移危险度系数逐渐增加，反映出随着构造应力的持续加载效应，应力逐渐积累接近破裂极值，直至地震发生临界时间应力累积达到破裂极值(即平均危险度系数达到零值线)发生地震.随后应力释放危险度系数大幅下降.其中在震间期，危险度系数演化曲线还存在明显的上下波动阶跃变化(危险度瞬间提高或者下降)，主要反应的是周边地震的同震应力扰动效应.

如下图所示(图7c)，1981道孚地震破裂区的发展曲线(黑色)在1923年时存在明显阶跃抬升现象，主要是受1923炉霍地震的影响.由图1可知，1923年炉霍地震与1981年道孚地震均发生在鲜水河断裂带西北段，且距离相近，前者发生在西段，后者在东段.1923年发生在炉霍的左旋走滑型破裂事件在其震源区东部尾端(即1981年道孚地区)产生典型的应力加载效应，即导致了1981年道孚震源区曲线的阶跃抬升.另外，受2008年汶川地震的影响，2013年芦山震源区的危险度曲线(深蓝色)在2008年也有抬升现象，而2017年九寨沟震源区曲线(红色)在2008年有下降现象，主要是受到汶川地震的应力卸载效应.

3.3 现今地震危险性概率

完成历史地震序列的演化计算之后则得到区域最后一次地震的震后应力场状态，据此可以计算危险度系数分布(公式(4)).显然单次的演化结果有很多不确定性，因此将所有反演给出的5000种初始应力场模型，均按照上述方法完成独立正演计算，对其中不能完成真实历史强震有序发生的算例进行修正或剔除，保证每一个算例演化的准确性.该过程中每一个算例都是确定性的，但是由于初始条件的不确定性，在合理范围内调整模型参数，确保足够多的算例可以实现与真实历史地震序列吻合.对所有算例完成正演计算之后，即得到最后一个地震——2017年8月8日九寨沟MS7.0地震震后的危险性状态，对所有的模型结果进行概率统计得到川滇地区现今地震危险性概率如图8所示 .

图8 川滇地区现今地震危险性概率图Fig.8 Probability map of current earthquake hazard in Sichuan-Yunnan area

计算结果表明，现今危险性概率较高的地区主要集中在以下区域：龙门山断裂带东北段，发震概率高达30%，主要是受到汶川地震震后应力调整的影响.苏生瑞等(2012)对汶川地震前后龙门山地区区域构造应力场演化研究发现，地震后剪应力在龙门山中央断裂上集中分布的区域逐渐向东北方向移动，造成龙门山东北段库仑应力增加，余震频发且沿着断层向东北方向单向延伸(单斌等，2009；Liu et al.,2019).本文在该区域的应力场演化结果与地震破裂过程的反演结果基本相符(陈九辉等，2009)，地震高危险度概率的变化规律与震后余震的分布规律一致，解释为汶川地震的余震作用.其原因是计算模型中仅考虑了主震破裂造成的应力释放，忽略了余震的应力释放，造成了该段显示较高的危险概率.

龙门山断裂带西南段两段地震破裂空区(汶川地震破裂区与芦山地震破裂区中间区段及芦山地震破裂区最南段与鲜水河断裂带交界处地震空区)，发震概率约为20%～25%, 主要是受汶川地震与芦山地震应力扰动共同作用的影响.单斌等(2013)及黄禄渊等(2019)通过计算指出汶川地震对芦山地震破裂区及龙门山西南段造成的库仑应力加载均为正值，且随时间推移震后库仑应力增加，造成区域地震危险性增加.而Liu等(2020)通过计算区域过去300年完整地震(M≥6.0)目录在该区域引起的库仑应力变化指出：两段地震破裂空区均受到汶川地震及芦山地震的应力扰动的加载，造成了区域危险度系数的提升.但龙门山断裂西南端因受到鲜水河断裂带前序多次历史地震的应力卸载，整体区域处于低应力状态，汶川地震及芦山地震的影响虽造成了该段的库仑应力增加，但该段整体仍处于相对低应力的安全状态；而汶川地震与芦山地震间的破裂空区缺少前序地震的应力卸载，整体处于高应力状态，受汶川地震及芦山地震的影响不容忽视，该地震空区的危险性更值得注意.

鲜水河断裂带西北段(包括与龙门山断裂带交叉地区)，金沙江断裂带南段区域，未来强震概率约为15%～20%；滇西南地区的龙陵—瑞丽、南汀河和澜沧江断裂带附近地区，红河断裂带西北段和中段与小江断裂带交汇区域，马边—延津断裂中南段，楚雄—建水断裂带东南段等也表现出了较高的强震危险性，约为10%～15%.徐锡伟等(2014)及赵静等(2015)根据地震地质资料、块体构造运动资料、历史地震活动性资料及GPS观测数据等定性或半定量研究认为，鲜水河西北段和东南段、安宁河断裂中段、大凉山断裂、则木河断裂中北段、小江断裂北段(东川附近)和南段(开远附近)以及红河断裂带中南段等地区均处于最新的地震活跃期，有发生强震的可能性；廖思佩等(2016)关于川滇地区应力场地震数值模拟研究中指出，川滇地区应力的不均匀分布主要集中在龙门山断裂带、鲜水河断裂带和安宁河断裂带交界区域南北两侧、小金河—丽江断裂带西段与澜沧江断裂带交界处南北侧、南汀河断裂带西段南北侧及红河断裂带中段东西侧；Wang 等(2015)指出红河断裂带西北段、红河断裂带与小江断裂带交界段以及滇西南龙陵瑞丽断裂带及澜沧江断裂带周缘均表现出较高的地震危险性.综上所述本文研究结果与前人结果具有较好的一致性.

4 讨论及结论

本文选取川滇地区自1900年以来震级大于6.5的30次历史地震进行初始应力场的约束，基于库仑-摩尔岩石破裂准则，推测发生过地震的破裂区应力水平至少达到了断层临近破裂的状态或者半临界状态；未发生过地震的区域应力水平低于临界破裂值.通过上千次反演重建川滇地区可能的初始应力场，在此基础上依次重现30次历史地震的发展过程，得到研究区域现今应力场并据此求得区域现今地震危险度概率分布.本文结果显示龙门山断裂带东北段和西南段局部地区存在较高的地震危险性概率，此外小江断裂带、金沙江断裂带、红河断裂带西南段、丽江—小金河断裂带也有部分高发震概率区段，与前人研究结果具有较好的一致性.另外需要指出，本文研究范围和计算量均较大，为了提高计算效率，仅考虑主要因素为一级近似，忽略了一些次要因素，导致局部计算结果存在误差如下.

(1)局部地区的构造异常

首先，滇西南云南腾冲地区(24.5°N—25.5°N，98.2°E—98.7°E)位于青藏高原东南边缘，也是本模型边界地区，一方面受到模型边界效应影响存在误差；另一方面其特殊的构造环境(印度与欧亚大陆碰撞、走滑逃逸作用和汇聚作用集中的地区)，使其成为我国最年轻的火山活动区之一，形成了区域复杂的应力场特征和频繁的地震活动(皇甫岗，2009).其次，本文未能考虑介质孔隙压力变化对地震活动性的影响，如四川盆地南部2018年12月16日MS5.7以及2019年1月3日MS5.3地震，Lei等(2019)认为与区域注水造成地下岩石介质孔隙压力变化有关，且2019年6月份以来该地区陆续发生多次MS>5.0的破坏性地震，包括MS6.0的四川长宁地震，均被认为是流体注入到地下开采引起的断层中，导致断层活化而发生的(Zhou et al.,2021).而水力压裂导致发震的断层大多是未知的，而且在目前的构造应力场的作用下，其产状不利于发生滑动破裂甚至与发震应力场状态相异(Yang et al.,2021;胡幸平等，2021).

(2)模型的设置

本文假设地震均发生在已有的断裂带上，但实际上，由于区域构造活动和地震错动的复杂性，有不少强震发生在尚未探明断层的块体内部，且难以判定其发震断层(邵志刚等，2010).如2014年8月3日，在川滇块体的东侧，莲峰—昭通断裂带以北的鲁甸地区发生MS6.5地震，野外地质调查认为其发震断层为大凉山断裂带的南端分支构造，且这一区域在近十年内发生数十次MS>5.0地震；同年10月7日在川滇块体西南边界(红河断裂带)西侧的景谷地区，发生了MS6.6地震，此次地震发震断裂带为隐伏的景谷—云仙断裂(也称之为无量山断裂)(徐锡伟等，2014).2017年8月8日九寨沟地震发生之后，地震地质学勘查结果判定此次地震的发震断层为隐伏的虎牙断层北段，而在这之前，均认为虎牙断层没有贯穿到北部(唐文清等，2004).在本文写作之际，2021年5月22日，青海玛多地区发生MS7.4地震，该地震同样发生于块体内部次级分支断裂——江错断裂，而在本文模型中并未包含该断裂.此外，川滇地区尤其滇西地区断层更为发育，很多断层存在较多的分支结构.但在本文计算模型中，仅考虑其主干断裂带，忽略分支构造，均简化成一条断层参与计算.忽略了块体内部具有强震潜能的断裂以及部分断裂带的分支断层等，这些简化均会影响到局部计算结果与实际结果不符.

此外，本文考虑计算的复杂性，在模拟过程中对块体划分及断层的参数仅做简化处理，地震的发生仅考虑应力释放，其断层摩擦强度没有随之变化，因此应力水平降低而强度却未有效变化，该计算缺陷需要在后续深入分析中有所体现.另外本文研究问题为三维，但实际为了提高计算效率，采用二维简化模型计算.利用二维简化模型的计算结果来描述分析三维问题，虽然可取得不错的结果，但必然存在一定的偏离.例如，应力增长速率计算等与实际的偏差，应力形因子约束中垂向应力的偏差等.在后续模型优化中，需要克服更多困难建立三维模型，更加接近实际情况.

(3)历史地震的参数及震后应力扰动

断层错动的同时会对周边区域应力场产生扰动，且不同的错动方式会造成不同的应力变化特征，因此对于历史地震破裂参数的选取，也会影响到计算结果.本文仅选取研究区30个MS>6.5历史地震参与计算.一方面有些历史地震的研究不充分，断层的错动方式及滑动量存在争议，如1933年叠溪地震，发生在东昆仑断裂带东端分支的岷江断裂南部，其震源机制存在多种不同的观点，如左旋走滑型(Wang and Shen, 2011)，右旋逆冲型(Chen et al.,1994)，拉张型破裂(Ren et al., 2018)等，在本文工作中暂未考虑此次地震.另外目前应力演化过程中我们只计算了前序地震的同震影响，忽略了下地壳粘滞性震后长期应力扰动效应，同样会对结果造成一定的影响.因此在后续细化工作中，应尽可能考虑更多的地震信息，对断层的应力状态进行更加精确的约束.

综上，本文是在一定的假设条件和简化模型的基础上，将模型的设置和参数选择采取一级近似对区域地震活动性趋势展开的概率性预报的探索性工作，所得的结果为一级近似的结果，与实际地质情况及前人的研究成果整体上存在一致性.尽管本文模型存在缺陷，局部地区结果存在误差，但该结果仍可以为区域中长期地震趋势研判做出一定的定量化科学参考.在未来的进一步深入研究工作中，将对本文计算方法逐步进行改善.