鄂尔多斯活动地块西缘强震间库仑应力作用

张瑞, 张竹琪*, 郑德文, 刘兴旺, 雷启云, 邵延秀

1 中国地震局地质研究所 地震动力学国家重点实验室，北京 100029 2 中国科学院广州地球化学研究所 同位素地球化学国家重点实验室，广州 510640 3 中国地震局兰州地震研究所，兰州 730000 4 宁夏回族自治区地震局，银川 750001 5 天津大学 表层地球系统科学研究院，天津 300072

0 引言

鄂尔多斯活动地块西缘位于南北地震带北段，是鄂尔多斯活动地块与阿拉善地块、柴达木—陇西地块的交汇区，区内包含祁连—海原、三关口—牛首山、罗山东麓、云雾山、六盘山等多条地块边界断裂带，具备发生8级以上强震的能力(图1).自公元876年到1920年近千年的时间里，该地区记录到的M≥61/2历史强震事件多达11个.韩竹军等(2008)和刘方斌等(2014)以同震弹性模型为基础，选择包含鄂尔多斯活动地块西缘部分事件的不同的强震序列，分别计算了1561后强震在前序地震作用下的累积同震库仑应力变化.结果显示1622—1920年这一活跃期的强震震中均位于前序强震引起的应力增加区，研究认为这些强震间普遍存在触发关系(韩竹军等, 2008; 刘方斌等, 2014).

库仑应力变化一般应用于主震对余震活动的触发作用，余震在空间上的集中发生或者主震前后地震活动性的变化可以认为是库仑应力变化发挥了作用(Das and Scholz, 1981; Stein and Lisowski, 1983; King et al., 1994; Toda et al., 2011; Toda and Stein, 2020; 缪淼和朱守彪, 2012; 李振月等，2020).部分研究利用库仑应力变化研究强震序列中各震源或断层之间的触发作用(Stein et al., 1997; Freed and Lin, 2001; Papadimitriou et al., 2004; Toda and Stein, 2008; 韩竹军等，2008；刘方斌等，2014; Shao et al., 2016; 孟令媛等, 2016).强震记录往往较短或不完整，难以准确衡量强震活动是否因为库仑应力作用发生了典型的变化，因此这些研究一般简单地将震源处或断层中心处库仑应力的增加与触发作用相对应.正是由于强震活动性的这种不确定性，当探讨强震序列中的触发关系时，应在历史记录比较完整的基础上，尽可能选择更长时间序列的强震活动来作对比和分析.很明显，由于地震数量、震级、震源相对位置和震源类型等不同，选择不同的地震序列，相应的累积应力变化计算结果会有所不同(Luo and Liu, 2010, 2018; Wang et al., 2017).

由于非仪器观测的历史地震参数多由定性描述获得，这些参数一般存在较大不确定性，在构建库仑应力模型以及计算分析时，应考虑这些参数的不确定性.库仑应力变化的分布对断层的几何与分段特征比较敏感(King et al., 1994; Lin and Stein, 2004; Mildon et al., 2016, 2019)，其特征或绝对大小还受介质模型、等效摩擦因数等参数选择的影响(Shan et al., 2013; Xiong et al., 2017).对于间隔几十年以上的强震事件，下地壳和上地幔震后黏性松弛效应导致的震后应力调整也可能显著影响模型结果(Freed and Lin, 2001; 万永革等, 2007; 雷兴林等, 2013).模型参数的不确定性也可能导致结论截然相反.比如，关于2008年汶川8级地震对2017年MW6.5九寨沟地震的库仑应力作用，汶川地震触发(汪建军和许才军, 2017; 单斌等, 2017; Liu et al., 2019; 黄禄渊等, 2019; 徐杜远等, 2020)和延迟(徐晶等, 2017; Jia et al., 2018)九寨沟地震的研究结果同时存在.除以上因素外，讨论强震间的应力作用时，由于断层几何尺度较大，如果用点取样的方法对应力变化进行分析，相应的结论还可能受到应力变化空间分布不均匀性的影响，进行三维断层面上库仑应力分布研究显得尤为重要(Mildon et al., 2016, 2019; Pino et al., 2019).

为探讨鄂尔多斯活动地块西缘强震间触发关系，我们选择了该区域公元1219—1920年间10次M≥61/2的强震，基于黏弹性分层模型计算分析了前序强震引起的后续强震发震断层面上同震和震后库仑应力变化.针对历史地震参数的不确定性，我们综合更新的历史地震和活动构造研究资料(宁夏地震工程研究院和宁夏回族自治区地震局, 2016)，重新梳理了相关历史地震的发震构造和震源参数，从而建立了考虑更多构造细节的断层模型.在此基础上，获得了沿断层面的应力变化分布，以及累积应力变化的时间演化特征，使模型分析更为全面.

1 研究区构造背景与历史强震活动

鄂尔多斯活动地块西缘处于柴达木—陇西、鄂尔多斯和阿拉善等多个活动地块相互作用的交接部位，也是南北地震带与祁连—秦岭断裂带等多组深大断裂的构造交汇复合部位，发育了不同规模、不同性质的断裂，构造形式复杂多样(图1).广泛分布的活动构造不仅控制着区域构造格局和形态，同时也控制着该地区强震的发生，晚第四纪以来该区域构造变形强烈、地震活动频繁(郑文俊等, 2016).

本研究中1219—1920年的10次历史强震主要分布在六盘山东麓断裂、会宁—义岗断裂、清水河断裂、烟筒山断裂、罗山东麓断裂、香山—天景山断裂、贺兰山东麓断裂、老虎山断裂和海原断裂上(图1).除会宁—义岗断裂位于地块内部，这些活动断裂带基本涵盖了鄂尔多斯活动地块西缘边界带的多个主要构造变形区，包括由海原断裂、香山—天景山断裂和六盘山断裂等组成的弧形构造带，由三关口—牛首山断裂、罗山断裂、云雾山断裂等组成的转换过渡带和由贺兰山东麓断裂、银川—平罗断裂、芦花台断裂、黄河断裂组成的银川地堑(国家地震局, 1988).海原断裂和香山—天景山断裂在西部以左旋走滑为主，受鄂尔多斯地块阻挡，在东部表现出明显的“端部效应”，逆冲分量增加，到六盘山南段以逆断性质为主(郑文俊等, 2016).北部银川盆地内部，处于拉张的构造环境，盆地内发育有多条正断为主的大型断裂带(国家地震局, 1988).以下简要介绍上述断裂的几何与变形性质，及相关的强震活动.

图1 鄂尔多斯活动地块西缘主要活动构造及历史地震分布

六盘山东麓断裂：六盘山东麓断裂北段在硝口附近与海原断裂相接，至马家辛庄附近与陇县—宝鸡断裂相接.整个六盘山东麓断裂带全长约90 km，走向315°～350°，倾向SW，倾角50°～70°.该断裂在开城以北运动性质以左旋为主，古地震平均复发周期约为4062±1885 a，开城以南断层运动性质以逆冲为主(向宏发等, 1998, 1999; 史志刚等, 2011).对历史地震资料的考证和实地调查表明，该断裂为1219年M7强震的主要发震断裂，重度破坏区位于固原、平凉和德隆一带(袁道阳等, 2008)，对应的发震构造为S1-1和S1-2(图3)，破裂长度约60 km.

小关山东麓断裂：小关山东麓断裂北起步湾以北，向南到泾源县的沙南以南，长约60 km，断面西倾，倾角50°～70°，整个断层由西向东逆冲现象明显，是一条典型的逆断层(宁夏地震工程研究院和宁夏回族自治区地震局, 2016).袁道阳等(2008)通过等震线分布认为六盘山东麓断裂北段与小关山断裂共同作用导致了1219年M7强震，但是也有研究表明小关山断裂自身活动性较弱，在地质地貌和断层剖面上其活动证据不明显(Zhang et al., 1991).本研究考虑了位于1219年地震重破坏区的小关山断裂北段(S1-3)(袁道阳等, 2008)，长度约40 km.

清水河断裂：宁夏地震工程研究院和宁夏回族自治区地震局(2016)地质调查研究表明该断层沿清水河盆地东岸依山而走，主体是一条东倾逆冲为主的隐伏断层，其晚第四纪以来活动强烈，均错断全新世以来的地质、地貌.宁夏地震工程研究院和宁夏回族自治区地震局(2016)综合断层活动特征、发震能力、震害范围和地震重现周期分析认为，1306年固原M7地震的发震构造为从固原南到开城乡附近的南小段(S2)，全长约20 km，该次地震具有破坏影响范围小，但极震区烈度大的特点.

会宁—义岗断裂：会宁—义岗断裂总长约50 km，走向310°～315°，倾向NE，倾角50°～70°，是一条逆冲兼左旋性质的断层.中国地震局兰州地震研究所(1993)野外考察时在会宁县城西北向东南经炭山村、新添、止于义岗附近发现一条长50 km的形变带，与会宁—义岗断裂吻合.何新社等(2007)对地震历史资料和地震等震线图分析表明破坏区的分布也与会宁—义岗断裂展布吻合，这些证据表明该断裂(S3)为1352年会宁M7强震的发震构造，破裂长度约50 km(何新社等, 2007).

烟筒山断裂：烟筒山断裂带是发育在青藏高原东北缘外围第三条弧形构造，处在香山—天景山断裂与牛首山—罗山断裂之间，整条断层总长超过150 km，北西起余丁金沙牙石沟，向南穿窑山东麓及炭山西麓，在云雾山一带可能与牛首山—罗山断裂带合二为一，沿途可见山脊、河流左旋位错.北段榆树沟以北，走向近310°，倾向北东，倾角70°左右；中段窑山东麓及炭山西麓断裂走向转为340°左右，倾向南西，倾角40°～70°(张维歧等, 2015; 宁夏地震工程研究院和宁夏回族自治区地震局, 2016).对于1495年强震的研究较少，本研究结合历史地震震中位置，利用经验关系得到地震破裂长度(S4)约25 km.1622年强震的重度破坏区位于固原以北，且探槽测年确定了该断裂南段(S6-1)为1622年M7地震的发震构造，破裂长度约60 km(宁夏地震工程研究院和宁夏回族自治区地震局, 2016).此前，对于1622年强震缺少发震构造研究结果，韩竹军等(2008)和刘方斌等(2014)将云雾山断裂(S6-2)作为了发震断层.

罗山东麓断裂：罗山东麓断裂北起西泉，与牛首山断裂相连，南到庙山，与云雾山断裂相连，全长60 km，总体走向345°～350°，倾向西南，倾角约62°～85°.断裂带在晚第四纪有显著活动，在全新世以来以右旋走滑为主兼有逆冲分量，古地震平均复发周期为2046±94 a(国家地震局地质研究所和宁夏回族自治区地震局, 1990; 闵伟等, 1993; 柴炽章等, 1999; Li et al., 2013).沿整个罗山东麓断裂大量分布最新破裂遗迹，断层陡坎、水系、山脊等地貌标志的右旋位移以及地震崩积楔等，显示出该断裂(S5)是1561年M71/4强震的发震构造，平均右旋位错为3.6 m(闵伟等, 1992).

香山—天景山断裂：该断裂带整体呈现向北东方向凸起的弧形，西梁头为东西两段分界，双井子为东段、东南段分界，全长近200 km.西段整体为走滑兼正断性质，倾向和倾角存在明显差异，东段和东南段以左旋走滑为主兼具逆冲分量，倾向西南，倾角60°～75°(国家地震局, 1988; Li et al., 2013; 李新男, 2014).整条断裂的古地震复发周期东西段有所差异，复发周期较长，约5000～6000 a(汪一鹏等, 1990; 李新男, 2014).张维歧等(2015)综合地震沟槽、陡坎、滑坡和土堆、水系断错、裂缝与鼓包等地表破裂标志的分布范围，认为1709年M71/2强震主要发生在断裂的东段，破裂带总的展布范围西起大堆堆附近延伸至小红湾附近(S7-1，S7-2)，全长约60 km，冲沟错动最大6 m左右，垂直位错多在0.5～1 m(国家地震局地质研究所和宁夏回族自治区地震局, 1990; 张维歧等, 2015).

贺兰山东麓断裂：贺兰山东麓断裂位于银川盆地，大致在套门沟位置处分南北两段，断裂全长120 km，整体走向约40°，倾向南东，倾角约60°～80°(国家地震局, 1988; 杜鹏等, 2009).晚更新世该断裂中北段较南段活跃，古地震平均复发周期中北段为2713±141 a，南段为5400±1327 a(Deng et al., 1996; 杜鹏等, 2009).杜鹏等(2009)通过地形地貌及古地震研究表明1739年M8强震没有延伸到套门沟以南，地震破裂(S8)集中在套门沟以北88 km范围内(国家地震局, 1988; 杜鹏等, 2009).其中，插旗口—苏峪口附近破裂带最为清晰，垂直断距在1.3～3.2 m之间，水平断距约1.45 m(国家地震局, 1988; 雷启云等, 2015)，也有研究认为该强震的发震断层可能为银川隐伏断裂(李梦銮和万自成, 1984; 郭增建等, 1988)或黄河断裂(包国栋等, 2019).

老虎山断裂：老虎山断裂东起景泰县东南的老庄沟，西至黑马圈河沟脑的独山子，约78 km，整体走向N70°～80°W，中更新世以来以左旋走滑为主，古地震平均复发周期为1080±92 a.从东到西依次为喜集水盆地段、老虎山段、草峡段、黑马圈河段(刘百篪等, 2013; 袁道阳等, 1997).何文贵等(1994)在老虎山段发现至今保存完好的形变带，大量的断层陡坎、左旋冲沟及洪积扇位错、地震鼓包、裂缝、陷坑以及喷砂孔等，表明该段(S9)为1888年景泰63/4～7级强震发震构造，长约30 km，倾角70°～80°，倾向SW.周俊喜等(1992)对老虎山段地震破裂带上的断距统计结果表明水平平均位错为3.6 m，垂直位错1～1.5 m.

海原断裂带：总体走向北西西，倾向南南西，倾角较大，自早更新世中晚期至中更新世初以来，断裂以左旋走滑为主(国家地震局地质研究所和宁夏回族自治区地震局, 1990).以大营水附近为西(S10-1)、中(S10-2)段分界，以南华山与月亮山交界处为东段(S10-3)和中段的分界(国家地震局地质研究所和宁夏回族自治区地震局, 1990; Ren et al., 2016).西段、中段古地震平均复发周期约1000～2000 a(张培震等, 2003; Ren et al., 2016)，东段古地震复发周期4196±1286 a(张培震等, 2003).1920年海原M81/2大地震形成了长达237 km的地表破裂带，宏观震源深度17.8 km，最大左旋位错达10 m，垂直位错可达数米(国家地震局地质研究所和宁夏回族自治区地震局, 1990).

2 方法与模型

2.1 库仑应力作用模型

库仑应力作用模型描述源断层地震活动对接收断层的应力作用.当源断层发生地震位错时，将引起周围介质发生弹性变形，相应地，介质中应力状态发生变化；另外，震后下地壳黏性物质发生黏性松弛时，由于偏应力的迁移，也会导致上地壳孕震层应力发生变化.源断层地震活动导致的弹性或者黏性应力变化，通过应力变化张量的投影，受影响的接收断层面上库仑应力变化都可表示为

ΔCFS=Δτ+μ′Δσn,

(1)

其中，Δτ表示接收断层上剪应力的变化，Δσn表示接收断层上正应力的变化.μ′为等效摩擦因数，它既包括了孔隙流体的影响又包括了断层区介质性质的影响，一般模型中取经验值.在没有特别说明的情况下，本文取较等效摩擦因数为中等的常用值0.4(King et al., 1994).一般认为，当地震引起接收断层上库仑应力增加，会促进下一次地震的发生，反之，将抑制或延迟地震活动.

2.2 黏弹性介质模型

在地震应力变化与迁移研究中，当时间尺度超过几十年、几百年，震后的黏弹性松弛引发的地壳应力场的持续调整不容忽视.这种黏弹性松弛是由弹性上地壳和黏弹性中下地壳的耦合形变所导致(Freed and Lin, 2001; 沈正康等, 2003).邵志刚等(2007)指出以Burgers体为黏弹介质模型可以解释地震引起的瞬时弹性响应、呈指数衰减的短期响应及线性增加的稳态长期响应，也可以解决Maxwell体和Kelvin体在模拟震后短期和长期形变时的不协调问题.本文基于Burgers体的分层黏弹性介质模型，使用基于PSGRN/PSCMP(Wang et al., 2006)的Geotaos进行同震和震后库仑应力变化的计算.

本文分层黏弹性介质模型中的黏度、P波波速和密度直接引用了孟令媛等(2016)对邻近的祁连—海原断裂带上1920年M81/2强震后的6次M≥7强震的同震和震后库仑应力演化研究中的参数，见表1.模型中S波速度参考并简化了陈九辉等(2005)利用远震体波波形资料和接收函数方法获得的青藏高原东北缘—鄂尔多斯地块地壳上地幔S波速度结构.宁夏地震工程研究院和宁夏回族自治区地震局(2016)结合地震深部探测与重力反演结果，获得固原地区莫霍面埋藏深度大约为47～52 km.同时，对1970年至2014年固原附近2107次ML5.9以下地震的精定位结果显示，除部分地震震源较浅，70.4%以上震源深度集中在10～25 km(图2).本文综合以上莫霍面深度和中小地震深度分布，设置黏弹性分层模型介质参数如表1所示.

表1 本研究采用的地壳与上地幔介质参数(陈九辉等, 2005; 宁夏地震工程研究院和宁夏回族自治区地震局, 2016; 孟令媛等, 2016)

图2 中小地震深度分布(改自宁夏地震工程研究院和宁夏回族自治区地震局，2016)

2.3 断层位错模型

本研究中断层位错模型采用简化的平直断层面，根据断层的活动构造分段或走向、倾角等几何特征的变化，模型中各地震的发震断层由单个或多个分段组成(图3).断层位错模型参数包括各断层面的长度、宽度、走向、倾角、同震位错的滑动角以及源断层的同震位错量(表2).这些参数的确定主要依据活动构造和历史地震等研究资料.

表2 断层位错模型参数

图3 发震断层分段图

如本文第1部分介绍，模型中同震断层的走向和倾角直接依据已公开的活动断层资料(闵伟等, 1992, 1993; Deng and Liao, 1996; 史志刚, 2011; 李新男, 2014; 张维歧等, 2015)确定.由于不同地震相关的资料详细程度不一，同震破裂长度和位错量的判定方法分为两种情况.一种情况是资料中的数据比较详细，我们直接使用资料中的研究数据，如1709—1920年及1561年强震的同震破裂长度和同震位错量均为野外调查结果(国家地震局地质研究所和宁夏回族自治区地震局, 1990; 闵伟等, 1992, 1993, 1994; 周俊喜等, 1992; 何文贵等, 1994; Deng and Liao, 1996; 张维歧等, 2015; 刘百篪等, 2013).另外一种情况，由于受到野外工作条件限制，资料无法提供足够详细的结果.这种情况下，我们在参考地震地表破裂调查、震害和烈度分布研究(何新社等, 2007; 袁道阳等, 2008; 宁夏地震工程研究院和宁夏回族自治区地震局, 2016)的基础上，还综合了破裂长度和同震位错量与震级的经验关系(邓起东等, 1992; Wells et al., 1994).这种情况涉及1219—1495年及1622年强震.模型中，大部分断层同震位错的滑动角由位错的水平和垂直分量的比值及相关的三角函数计算得到.由于缺少直接观测结果，历史地震的断层宽度缺少准确约束，根据该地区小震震源深度的分布范围(宁夏地震工程研究院和宁夏回族自治区地震局, 2016)，模型中统一设置断层深度为25 km.

3 计算结果

为了探讨强震间库仑应力的综合作用，我们计算了1306年至1920年9次M≥61/2地震临震时各自断层面上由前序强震活动引起的同震和震后累积库仑应力变化.计算结果如图4所示，其中1219—1495年和1495—1920年两个阶段强震间的库仑应力作用有明显区别：除了1306年7级地震发震断层自7～8 km深度至地表部分有明显应力增加(图4a)，1306—1495年强震发震断层面上库仑应力增加均不明显，甚至库仑应力变化为负值(图4a—4c)；与此相反，1561—1920年强震的发震断层面上库仑应力均显著增加(图4d—4i)，至临震前，这些库仑应力变化均接近或超过0.01 MPa.

地震引起的库仑应力变化一般比构造应力和同震的应力降小得多，通常0.01 MPa被认为是能触发地震的库仑应力变化的下限(Harris, 1998).为了方便定量地分析强震间的应力作用，我们将强震断层面上典型部位的累积库仑应力变化值与触发“阈值”0.01 MPa进行对比.如表3所示，“中心值”为发震断层面在17.5 km深度中点位置的库仑应力变化，该深度为固原地区地震精定位所得震源深度分布的优势深度(宁夏地震工程研究院和宁夏回族自治区地震局, 2016)(图2)，其中1709和1920年强震对应多段发震断层，我们选择S7-2和S10-2进行计算.“最大值”、“平均值”表示震源优势深度处，沿断层走向分布的库仑应力变化的最大值、平均值.“±”表示以上各库仑应力变化值是否高于触发阈值，其中“+”表示更高，“-”表示更低.“高值比例”表示库仑应力变化高于阈值部分所占断层长度的比例.从表3可以看出，除了1561年M71/4地震断层面上中心值略微低于阈值，1561—1920年其它强震库仑应力变化可达阈值1.2～13倍，优势深度处的最大值均达到阈值，可达阈值的2～300倍.另外，除了1561年强震优势深度处的平均值未达到阈值，其它事件的平均值可达阈值的1.2～6.4倍，高于阈值应力变化的断层段长度比例高达40%～100%.

表3 断层面上典型部位累积库仑应力变化

以上结果显示出临震时断层面上库仑应力变化的分布存在空间不均匀性.为了更加全面地反映强震间的库仑应力作用，我们对震前多个时间点的震源优势深度处沿各断层面走向的库仑应力变化进行采样分析.如图5所示，随着前序强震的发生，1306—1495年强震断层整段的库仑应力变化始终为负值(图5a—5c)，而1561—1920年强震断层上主要区段的库仑应力变化为正值(图5d—5j).另外，除了1561年和1920年海原地震断层(S5,S10-2)靠近端部出现应力下降(图5d,5i)，其它地震断层上库仑应力变化随时间发生了整体升高(图5e—5h,5j).

图4 各强震临震时断层面上同震和震后累积库仑应力变化

图5 优势深度断层面上同震和震后库仑应力变化随时间演化过程

为了更好地反映断层面上库仑应力变化的时间演化过程，我们将1306—1920年强震的发震断层段均匀分为9(3×3)个子断层面，并选取各子断层面中心点作为采样点分析了库仑应力变化随时间的变化(图6).结果表明，1306年强震断层面(图6a,S2)上自1219年以后库仑应力变化的平均值保持在大约-0.1～-0.01 MPa的水平，但是平均值加标准差为正值，约为0.07 MPa，说明该断层面存在局部的库仑应力变化为正值.1352年强震断层面(图6b, S3)上库仑应力变化的平均值为-0.010～-0.012 MPa(图4 a)，1495年强震断层面(图6c, S4)上库仑应力变化的平均值为-0.0016～-0.0005 MPa，即使考虑标准差，断层面上库仑应力变化也保持为负值.整体上看，1306—1495年3次强震发生前，各自断层面上库仑应力变化一直处于负值水平.

图6 各强震断层面上震前同震和震后累积库仑应力变化的时间演化过程

与以上3个地震断层面不同，1561年、1709年、1739年，以及1888年强震的断层面上，库仑应力变化的平均值在地震发生以前发生快速升高(图6d,f,g,h).如图6d所示，1561年强震断层面因受1495年强震的应力作用，库仑应力变化由零附近升高到0.01 MPa左右，并在1561地震前的66年时间里维持在这个水平附近，而在1495年强震发生之前276年时间里，1561年强震断层面上库仑应力变化的幅度约在±0.002 MPa以内(图6d).如果简单地用变化幅度除以时间长短来计算库仑应力变化的速率，1495年前后的速率分别为1.45×10-5MPa·a-1和1.52×10-4MPa·a-1，也就是说1561年强震发生前相对较短时间内，其断层面上库仑应力变化的速率增加了大约10倍.与1561年强震断层面相似，1709年强震断层面上库仑应力在1561年以后约148年时间里平均增加了0.04 MPa，而在此前约342年时间里，平均库仑应力变化幅度约为-0.005～0.001 MPa(图6f)，说明1709年强震前的490年里，在距其发生约30%的时间段内，该断层面上强震引起的库仑应力变化速度增加了近15倍.在1739年强震断层面上，也存在库仑应力快速变化的现象，在1561年以前三百多年里，该断层面上强震引起的平均库仑应力变化幅度小于0.001 MPa，但随之在约178年里增加了近0.03 MPa(图6g)，该断层面上强震引起的库仑应力变化速度增加超过50倍.1888年强震断层面上平均库仑应力变化也是在1561年强震以后发生快速增加，这个倍数相对略低，但也达到约5倍(图6h).

在1622年和1920年强震断层面上，也可以看到地震前强震引起的库仑应力发生快速增加的现象，但只是局部时间特征，而且幅度并不明显(图6e,i).如图6e所示，1219年强震对1622年发震断层南端的端部效应导致1622年强震断层面上库仑应力平均值增加明显，超过0.02 MPa，如果考虑标准差，可以看到±0.025 MPa的变化范围，1219年强震后S6-1断层面上应力分布极不均匀，局部库仑应力增加值将达到0.05 MPa左右，此后因1306年和1352年强震作用，库仑应力略微降低，此后200多年时间里维持在稳定的水平.直至该事件发生61年前，1561年强震的发生导致平均库仑应力增加约0.005～0.01 MPa，且标准差有所减小，表明断层上库仑应力的不均匀性有所减小(图6e).在1920年强震断层面上，平均库仑应力变化整体趋势平稳，在1219—1561年时间段内，库仑应力增加水平大约在0.01 MPa附近，1561年后降低为负值(图6i).之后，1709—1888年强震导致平均库仑应力增加了大约0.02 MPa，总体库仑应力达0.01 MPa，但整体上库仑应力变化的相对幅度较小，趋势平缓(图6i).

以上结果表明，断层面上库仑应力变化分布的不均匀性有可能影响对1306年、1561年、1622年和1920年强震震源受到的应力作用过程的解释(图6a,6d,6e,6i).因此，我们进一步分析了这4次强震断层面上的库仑应力分布特点.如图7所示，1306年强震临震前其断层面S2上的库仑应力变化在近地表约7～8 km深度的范围内为正，达0.05～0.16 MPa，在8 km以下的深度库仑应力变化为负值(图7a)，从库仑应力变化时间曲线可以看出，除浅部1、4、7号区域库仑应力变化为正，其他大部分区域的库仑应力变化保持为负值(图7e).

1561、1622和1920年强震断层S5、S6-1、S10-2上1—6号区域均表现出临震前库仑应力变化整体快速增加(图7b—d, f—h)，其中S5在1495年之前的276年里库仑应力变化幅度为±0.0025 MPa，1495年之后到震前库仑应力变化可达0.018 MPa，库仑应力变化速率最大增加了15倍(图7f).类似地，S6-1在1561年之后的库仑应力变化速率约为之前的5倍(图7g)，S10-2在1709年前后库仑应力变化速率增加了近7倍(图7h).可以看出，1561、1622和1920年3次强震断层面上大部分区域经历了明显的库仑应力变化快速增加，只是在断层的端部出现局部的库仑应力集中或者降低.

总体上，断层面上库仑应力变化出现大的标准差，以及出现库仑应力变化平均值与其标准偏差的时间演化趋势特征明显不同的情况(图6a,d,e,i)，主要是由于断层面上局部区域发生与断层面多数区域相反的库仑应力变化(图7a,b,d)，或者虽然局部具有与整体相同极性的库仑应力变化，但由于端部效应等因素，这种局部变化特别大(图7c)，因而导致断层面上平均库仑应力变化时间曲线无法真实反映断层面的整体情况.

图7 1306、1561、1622和1920年断层面同震和震后库仑应力变化随时间演化过程

4 讨论

4.1 海原强震库仑应力演化

1920年海原强震对应的发震构造由三段组成(S10-1、S10-2、S10-3)，前面的分析中我们选取了活动构造资料给出的强震震中所对应的中间段(S10-2)进行了详细的分析，并在对应的正应力区域看到了库仑应力快速增加的现象(图7h).图8为S10-1和S10-3段对应的库仑应力变化曲线，可以看到S10-1段自1739年强震开始断层面上的应力有一个快速增加的趋势(图8c)，而对于S10-3段经历了1709年和1739年两次明显的库仑应力降低，在1888年强震作用下有轻微的库仑应力上升，整体呈现出应力演化后期显著降低的特点(图8d).

图8 1920年S10-1、S10-3断层面同震和震后库仑应力变化随时间演化过程

此外，就1920海原M81/2强震的整个破裂长度来看，从1219年强震开始到1920年海原强震发生前其断层上的库仑应力呈现出正负相间的分布模式(图4i和图9a)，这种相间出现的库仑应力变化如何形成两百多公里长的破裂范围？为了解释这一现象，我们根据1920年强震前断层面上的应力分布对断层进行了新的分段(图9a)，并计算了前序断层滑动导致的后续断层上库仑应力变化(图9b—9n).考虑到震源位置的不确定性，这里假定了三种破裂模式，分别从左依次向右破裂(图9b—9g)，从右向左依次破裂(图9h—9m)以及从中间向两侧开始破裂(图9n).结果表明，当正的库仑应力段落发生破裂后会导致后续破裂的断层面上产生极高的库仑应力变化，这种应力累积要高出前序1219—1888年强震造成累积库仑应力变化的数十倍，足以将原来负的应力极性改变，使原本相间分布的正应力区域贯通.在自左向右和自右向左破裂模式中，当破裂至S10-2-1时整个剩余断裂处于了极高的应力增加状态(图9e，k)，中间向两侧传递的模型中，中间S10-2段的破裂足以导致两侧剩余断裂全部处于极高的应力状态(图9n)，之后断层面上的应力会持续升高直至破裂.

在三种假设模型中，前序的断层滑动都会在端部引起后续破裂断层极高的库仑应力变化，这与各段断层走向变化较小的几何展布和均以左旋走滑为主的运动学性质具有紧密联系.海原强震长达200多公里的地表破裂以及跨越阶区的破裂很可能是由于这种特定的断层几何展布和断层力学性质所促成的.为了进一步验证这种独特的断层几何展布和走滑性质在海原M81/2强震中的作用，我们在弹性半空间介质模型中建立了垂直左旋走滑的源断层模型，采用均匀滑动模型，位移量大小为5 m，断层面长(L)、宽(W)比为1∶1和5∶1两种情况(该比值的变化范围涵盖了海原断裂不同分段的变化)，接收断层的走向为平行于源断层和顺时针变化30°两种情况(该走向的变化范围涵盖了海原断裂不同段的走向变化)(图9o—9r).结果表明破裂的两端会形成明显的端部应力升高区，其明显的升高区域可以达到源断层破裂长度的1～4倍.当接收断层的走向变化30°时，高应力花瓣的区域也会发生同向旋转(图9p,r)，这种旋转后的高应力分布模式与海原断裂的断层几何展布具有较好的一致性.

图9 海原强震断层面库仑应力累积破裂模式图

海原断裂独特的几何展布和运动学性质，造成了先前破裂部分对后续破裂断层面上产生较高的库仑应力积累，这种应力增加足以改变断层在1219—1888强震中形成的低应力变化区域，从而贯通原来正的库仑应力变化区域，促成了海原M81/2强震200多公里的破裂.此外，海原断裂带附近历史强震目录的完整性也会影响1920年发震断层面上震前库仑应力分布情况，这有待更多相关资料的考证和完善.

4.2 发震构造的不确定性及其影响

由于研究区存在多条断层相交汇和彼此邻近的情况，而历史地震发震构造的确定主要依赖于文献记载，并且历史地震破裂的调查受到野外工作条件等诸多因素影响，因此部分历史地震的发震构造尚难以被准确分辨出来.本研究涉及的强震序列中，关于1219年和1622年事件的发震构造尚存在争议，我们针对可能的发震构造构建了多个模型，以探讨这些发震构造的不确定性对库仑应力变化及迁移结果的影响.

袁道阳等(2008)认为六盘山东麓断裂与小关山断裂同为1219年强震的发震断层，但也有研究表明晚第四纪以来小关山断裂活动性较弱，在地形地貌和断层剖面上其晚第四纪以来活动的证据不足(Zhang et al., 1991)，该断层可能不是1219年强震的发震断层.上文的模型中，该两条断层均被设置为1219年强震的源断层，作为对比，此处另外引入了只包含六盘山东麓断裂的模型结果进行对比.如图10所示，不论模型中是否将小关山断裂当作1219年强震同震破裂的一部分，此次地震导致的1306年发震断层面上库仑应力增加的部分都主要集中在浅部8～9 km宽的断层面上.尽管库仑应力变化的强弱有明显变化，但小关山断裂的破裂对1306年强震断层面上库仑应力变化分布的总体特征没有本质影响(图10).此外，由于小关山断裂和1622年发震构造烟筒山断裂的走向近似，首尾相接并且变形方式相似，由于断层端部的应力效应，小关山断裂的破裂会有利于1622年发震断层面(S6-1)累积库仑应力的升高(图4e).

图10 1219年不同发震断层对1306年发震断层的库仑应力作用

关于1622年发震断层，有研究认为可能是云雾山断裂而不是烟筒山断裂(韩竹军等, 2008; 刘方斌等, 2014).计算显示1561年强震导致云雾山断裂(S6-2)上的同震以及同震和震后库仑应力增加十分明显(图11a,b)，说明云雾山断裂的确有可能受到应力触发作用而发生1622年强震.我们将云雾山断裂作为1622年强震发震断层，计算了1219—1920年强震发震断层面上库仑应力变化.结果显示，1622年地震将导致1709年强震断层面上的同震(图11c)及震后(图11d)库仑应力整体下降，但是1709年以前所有强震引起1709年发震断层面(S7-2)上的累积库仑应力依然增加(图11e,f).总体上，不论1622年强震发震构造是否是云雾山断裂，不改变1709年强震断层面上在前序强震作用下累积库仑应力变化明显升高的现象(图11e,f).且通过图4f和图11f可以看到虽然1622年发震构造不同，但1219—1622年强震的同震和震后库仑应力变化在后续强震断层面上的分布模式并没有显著不同.

图11 1622年发震断层的影响

4.3 模型介质参数不确定性及其影响

模型中介质参数的不同显然会影响模型的定量结果，但是否会改变模型结果的性质呢？我们通过计算优势深度断层面上库仑应力变化的中心值、平均值和达到阈值(≥0.01 MPa)的断层长度比例检验了黏度、等效摩擦因数和弹性层厚度对已有结果的影响(表4).

前文的黏度直接引用了孟令媛等(2016)对邻近祁连—海原断裂带上1920年M81/2强震后的6次M≥7强震的同震和震后库仑应力演化研究中采用的参数.其取值主要考虑了王庆良等(1997)利用1990年共和7.0级地震震后垂直形变资料反演震区岩石圈和软流圈均一化后的地球介质的有效黏度约为1018Pa·s量级，郝明等(郝明等，2010；郝明，2012)基于1990年共和地震，利用相邻水准点间的原始高差观测值约束连续介质位错模型，得到下地壳和上地幔内的黏度为1019～1020Pa·s量级，同时也考虑了邵志刚等(2008)利用震后短期内观测数据做约束反演给出的青藏高原东北部东昆仑断裂的最佳黏度为5.0×1017Pa·s.此外，许多学者用不同的方法对青藏高原及周缘的下地壳黏度进行了探讨，Royden(1996)利用解析方法研究青藏高原抬升高度及地壳厚度随时间变化问题时，认为其黏度范围很宽，介于1012～1021Pa·s 之间；Clark和Royden(2000)用地形梯度法对青藏高原周边地区下地壳黏滞性进行了估计，结果为1016～1021Pa·s；Beaumont等(2001)以湿黑山石英和马里兰辉绿岩为样品，根据热力学塑性流变定律计算得到的上地壳和中地壳等效黏度分别为1019Pa·s和1018Pa·s；Hilley等(2005)利用昆仑断裂附近的地质和大地测量数据，结合地震周期用贝叶斯方法得到的青藏高原下地壳黏度在1018～1021Pa·s；Ryder等(2011)使用 Burger 体流变学模型对2001年Mw7.8级昆仑山地震的震后余滑和黏弹性松弛等进行研究，得到下地壳震后初期的短期黏度和长期黏度分别为9×1017Pa·s和1×1019Pa·s，并与玛尼地震的流变学参数进行对比，得到流变学黏度有效约束为 5×1017Pa·s；魏聪敏等(2020)以下地壳管道流模型为基础，采用地貌解析方法估算了西秦岭—松潘构造节及邻区的下地壳黏度变化范围在1018～1020Pa·s.可见，青藏高原及其周边中下地壳黏度研究结果的范围介于1016～1021Pa·s，深部介质黏度不确定性较大.且在青藏高原及其周边强震序列的黏弹性库仑应力计算中(沈正康等，2003；万永革等，2007；雷兴林等，2013；Shan et al., 2013；Shao et al., 2016；徐晶等，2019)，其黏度的取值因研究区域以及流变模型等不同而有所差异.因此，作为简单对比，我们将原模型中的黏度统一降低一个量级进行了计算分析(表4).

表4 主要参数不确定性的影响

上文中等效摩擦因数为0.4，这里分别考虑了等效摩擦因数为0.2和0.6的情况.此外，在地球介质模型建立时，参考了宁夏地震工程研究院和宁夏回族自治区地震局(2016)给出的固原地区莫霍面埋藏深度大约为47～52 km这一结果，刘启民等(2014)利用接收函数方法研究表明在我们研究区附近大约以105°E为分界，西部比东部莫霍面平均深度至少深5 km.该部分，我们考虑了弹性层厚度由原来的25 km调整至30 km的情况.

整体上看黏度和弹性层厚度的变化，对模型中的计算结果影响较小，计算得到的高值比例变化基本维持在5%范围内(表4).不同等效摩擦因数下库仑应力变化结果表明，当摩擦因数较高时(μ′=0.6)，1622年和1888年库仑应力变化有所降低，且高值比例明显下降.但是，1739和1920年计算得到的中心值和平均值有所增加，且1920年强震断层面上的库仑应力变化高值比例增加5%(表4).可以看到，库仑应力变化随等效摩擦因数的改变不是简单的有效摩擦因数越大，断层就越稳定，库仑应力变化越小(朱守彪和缪淼, 2016; 徐杜远等, 2020).

4.4 震后黏弹性应力迁移的重要性

对比同震库仑应力变化与同震和震后库仑应力变化(图12)，可以看到仅考虑同震库仑应力时，应力演化的总体方向没有明显改变.1306和1352年发震断层面上同震累积库仑应力变化依旧处于负的应力值，1495年发震断层面上的同震累积库仑应力的平均值接近于0，最大约为0.0001 MPa，不足以触发1495年强震(图12a—12c).1561—1920年断层面上的库仑应力变化表明，除1561年事件外，其余同震库仑应力变化相较于考虑震后黏滞松弛效应有所降低.S5断层上仅考虑同震比考虑震后效应库仑应力变化增加了近0.002 MPa，占同震和震后累积库仑应力变化0.01 MPa的20%(图12d)，S6-1断层上仅考虑同震比考虑震后效应库仑应力变化降低了0.015 MPa，占考虑震后效应0.03 MPa的50%(图12e)，S7-2断层上仅考虑同震比考虑震后效应库仑应力变化降低了0.02 MPa，占考虑震后效应0.035 MPa的57%(图12f)，S8断层上仅考虑同震比考虑震后效应库仑应力变化降低了0.02 MPa，占考虑震后效应0.03 MPa的67%(图12g).1888年临震断层面上同震累积库仑应力变化仅为0.001 MPa，其库仑应力变化受1739年强震的影响很小，后期库仑应力的增加主要源于临近强震长期的黏滞松弛效应，并在1709年强震后这种震后黏滞松弛效应起到了主导作用，震后库仑应力增加明显加快(图12h).对S10-2断层上仅考虑同震效应得到的断层面上同震累积库仑应力平均值均为0.0013 MPa，其较低的平均库仑应力变化主要是由于临震断层东端存在较大的库仑应力变化减小区造成的，但是也可以看到1709—1888年震后效应极大的提升了断层面上的库仑应力(图12i).由此可见，本研究中的震后效应进一步抑制了1352和1495年强震的发生，并且对1622—1920年强震的触发作用具有重要贡献.

图12 震后黏滞松弛效应对库仑应力变化的影响

在我们的模型中如果仅考虑1561—1920年强震序列的同震库仑应力变化，1561—1709年强震发震断层面上平均库仑应力变化均达到或接近0.01 MPa，强震间存在明显的触发作用，这一结果同前人研究一致(韩竹军等, 2008; 刘方斌等, 2014).但是，对1739—1920年强震断层面上同震库仑应力变化的平均值并没有明显的库仑应力增加，这一结果可能是因为地球介质模型、强震序列、断层位错模型以及库仑应力变化的取样方式不同所导致的.

为了进一步解释库仑应力性质的演化，我们选择了图12i 中的1561—1622和1709—1739两个典型时间段，通过剖面应力的演化分析断层周边应力变化过程.剖面位置为图3中AA′和BB′两条剖线垂直下切50 km深的截面(图13a—13f, 13m—13r)，以及沿断层面S10-2延深至50 km深面上的库仑应力变化(图13g—13l, 13s—13x).可以看到1561年强震发生后在S10-2周边处于负的库仑应力变化区域(图13b)，随着1561年强震离逝时间增加，25 km以下及1561年发震构造靠近S10-2一侧的负应力持续迁移，负的库仑应力持续在断层S10-2附近累积(图13c—13f)，从而形成了图12i中1561—1622年应力持续降低.同时，S10-2断层面上的负的库仑应力在大小和分布范围上持续扩大(图13g—13l).1709年强震将震前S10-2断层附近的负应力极性转变(图13m、n)，且随着时间的演化，下地壳和上地幔中正的库仑应力向上迁移，使上地壳S10-2断层附近的应力持续增加(图13o—13r)，这种应力从下地壳和上地幔转移到上地壳的过程在断层面上(图13s—13x)同样可以观察到.这种同震阶段应力分布和黏弹性介质模型结构导致的后期库仑应力向相对高黏度地壳中的迁移，使上地壳接收断层周围的库仑应力变化的大小和范围发生改变，从而导致图12i中1561—1622和1709—1739时间段应力的演化.

图13 同震和震后库仑应力演化剖面

该部分仅选取了1920年强震断层S10-2段上两个典型应力升降事件相关剖面进行了库仑应力变化的分析，可以看到同震阶段的应力分布模式和黏弹性结构对震后应力的迁移具有重要影响.在其它的阶段中，我们也可以观察到受同震应力分布和黏弹性结构影响，下地壳和上地幔中不能长时间承受剪应力，会向黏度高的上地壳迁移的现象.此外，这种黏弹性库仑应力变化与断层的力学性质、几何展布以及地球介质模型参数的选取具有重要联系，相关工作仍需进一步研究.

4.5 强震间库仑应力作用

库仑应力变化只是断层面上应力状态的扰动，因此同震和震后库仑应力变化不会产生地震，但足够高的应力变化可以促进地震的发生，反之负的应力变化会抑制地震的发生(Harris, 1998).如果以0.01 MPa为触发阈值，我们的结果显示至临震时，1306年、1352年和1495年强震断层面上累积库仑应力变化平均值均为负值，而1561—1920年强震临震前断层面上累积库仑应力变化的平均值接近或超过了触发阈值，且具有较高的高值比例(图4,6，表3,4)，可以认为研究区内1219—1495年的强震之间没有应力触发作用，而1495—1920年的强震之间存在应力触发作用.

限于鄂尔多斯活动地块西缘历史地震记录时间区间长度，以及早期历史地震记录可能存在缺失，因此无法直接通过比较不同时期的地震活动性来判断库仑应力变化作用是否准确、有效地改变了强震的活动.该地区早在公元876年便记录有M≥61/2地震，震中位于宁夏雄州(今吴忠)附近.不妨假设公元1219年以后61/2级以上地震记录是完整的，那么可以看出大约在公元1495年以后，该地区历史强震应变能释放存在一个加速释放的过程(图14).对照1495年前后临震时各强震断层面上库仑应力变化所呈现的触发作用，累积库仑应力变化水平的持续以及加速升高(图6d—6i)很可能导致了鄂尔多斯活动地块西缘强震活动变得更加活跃，从而促成近400余年的强震丛集(韩竹军等, 2008).

图14 876—1920年M≥61/2强震应变能累积释放时间曲线

如果排除1495年以前的强震影响，只计算1495—1920年强震序列中各事件引起的累积库仑应力变化，结果显示这一阶段发震断层面上库仑应力变化还是明显增加(图15)，即强震间的触发作用依然存在，该阶段的应力作用的性质不受1219—1352年强震的影响.综合考虑1495年以前强震间库仑应力的抑制作用，进一步指示出强震间的库仑应力作用是强震活动是否进入活跃期的重要因素.

图15 1495年及以后强震引起的断层面上同震和震后累积库仑应力变化

需要指出的是，1495年以后强震发震时间普遍地比断层面上累积库仑应力变化超过阈值(0.01 MPa)的时间延迟数十年到数百年(图6,图7).以1622年强震为例，1219年地震导致该断层面上库仑应力变化超过0.02 MPa(图6e)，在断层面东南端库仑应力变化甚至达到0.04～0.1 MPa(图7c,g)，此后在1622年地震发生前的400多年里，累积库仑应力变化平均值一直在0.02 MPa附近或者更高(图6e).相似地，在1920年强震发震断层面上，编号7、8、9的区域在1219年以后库仑应力变化已超过阈值并持续了300多年(图7d,h)，但直到1561年都未发生类似1920年强震的事件.另外需要引起注意的是，1306年强震发生前，其发震断层面浅部7～8 km以上区域的库仑应力变化超过了阈值(图7a,e,图9).这些现象说明，强震间的库仑应力触发或抑制作用比强震对余震或中小地震活动的影响复杂得多，不宜简单以库仑应力变化的绝对大小来判别.如果1495年前后发震断层面上库仑应力变化的时程特征的差异是典型的，那么相较于应力变化大小，断层面上发生整体的平缓—快速的应力变化进程对判别区域强震活动的长期趋势更具可操作性.当然，目前无法完全排除该地区1495年以前61/2级以上历史地震记录缺失的可能性，这有可能会影响对强震间库仑应力作用的分析，更加细致和可靠的研究还有待于相关资料的更新与完善.

5 结论

黏弹性库仑应力模拟结果显示，1495—1920年发生于鄂尔多斯活动地块西缘的7个M≥61/2的历史强震序列中，前序强震活动均引起后序事件发震断层面上同震和震后累积库仑应力变化升高0.01～0.1 MPa，并且随着时间上临近后序事件的发生，同震和震后累积库仑应力变化趋势普遍经历了较高的增长速率，说明通过库仑应力作用，这些前序强震的活动对后序强震的发生有明显的促进作用.模型还显示，1219—1495年4个强震间累积库仑应力变化以负值为主，强震之间的触发作用不明显；另外，1219—1352年的3次强震对1495—1920年的7个强震断层的库仑应力作用较小，不影响模拟所得1495年后7次强震间同震和震后累积库仑应力变化总体分布和趋势.模型显示出鄂尔多斯活动地块西缘强震间库仑应力作用的特点在1495年前后发生了前后相反的变化，该时间节点与研究区强震的应变能开始快速释放的时间相符合，指示出强震间的库仑应力作用是强震活动是否进入活跃期的重要因素.因此，初步认为强震引起的库仑应力变化与迁移对鄂尔多斯活动地块西缘边界带的强震间的触发关系及强震活跃期发挥了重要作用.

致谢计算所用软件为PSGRN/PSCMP(Wang et al., 2006)和雷兴林研究员提供的Geotaos，部分图件采用GMT(Wessel and Smith, 1998)绘制，审稿专家对本文提出的修改建议，在此一并表示感谢！