2020, Vol. 63
2. 中国科学院大学, 中国科学院计算地球动力学重点实验室, 北京 100049
2. Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
印度—欧亚板块的碰撞挤压，形成了构造环境复杂的青藏高原，在高原南部发育了喜马拉雅巨型俯冲带；高原内部物质被侧向挤出，发育了数条巨型走滑断裂带，主要有昆仑山、甘孜—玉树—鲜水河、阿尔金、祁连山—海原断裂带等；高原东缘受古老而坚硬的四川盆地阻挡，发育出逆冲推覆的龙门山断裂带(Molnar and Tapponnier, 1975，England and McKenzie, 1982; 邓起东等, 2002), 如图 1所示.该区域是中国大陆地震发生最频繁、强度最大的地区，近百年间有300多个震级大于MW6.0的地震发生于此(张培震等，2013).该地区的孕震环境、应力场演化及未来的地震活动性趋势等一直是地学界研究的热点及难点问题.
|
图 1 青藏高原及邻区的构造简介 F1—阿尔金断裂带；F2—祁连—海原断裂带；F3—昆仑山断裂带；F4—甘孜—玉树断裂带；F5—嘉黎断裂带；F6—喀喇昆仑—雅鲁藏布江断裂带；F7—喜马拉雅俯冲带：F8—理塘断裂带；F9—鲜水河断裂带；F10—龙门山断裂带；F11—岷江断裂带；R1—塔里木盆地；R2—柴达木盆地；R3—鄂尔多斯地块；R4—巴颜喀拉块体；R5—羌塘块体；R6—拉萨地块；R7—雅鲁藏布江块体；R8—四川盆地.图中沙滩球为1904—2014年间MW6.5以上的强震震源机制. Fig. 1 Tectonic setting of the Tibetan Plateau and its surrounding area F1—Altyn Fault; F2—Qilian-Haiyuan Fault; F3—Kunlunshan Fault; F4—Garzê-Yushu Fault; F5—Jiali Fault; F6—Karakoram-Yarlung Zangbo River Fault; F7—Himalaya Subduction Zone; F8—Litang Fault; F9—Xianshuihe Fault; F10—Longmenshan Fault; F11—Minjiang Fault; R1—Tarim Basin; R2—Qaidam Basin; R3—Ordos Block; R4—Bayan Har Block; R5—Qiangtang Block; R6—Lhasa Block; R7—Yarlung Zangbo River Block; R8—Sichuan Basin. The beach balls in the figure are the focal mechanisms of MW > 6.5 strong earthquakes occurred from 1904 to 2014. |
青藏高原及其邻区的地震活动性在受板间(印度—欧亚板块)作用力远程驱动的同时，还受到板内构造运动的控制，即地壳岩石层在板内构造作用力下发生变形，内部能量/应力非均匀积累，在某处积累的应力超过其承受极限，则发生地震释放应力.且地震释放的应力并没有凭空消失，而是转移到周边区域，对其应力场产生扰动.因此，对于区域地震活动性的数值研究，需要同时考虑构造活动产生的应力场及地震产生的应力扰动的共同作用.在我国，王仁等(1980)最早展开数值模拟应力场演化问题的研究工作，考虑背景构造应力场的影响，重现了华北地区从1966年到1976年之间的五次大震的演化过程，但背景应力场的选取非常粗略单一.随后大量学者对于该问题开展了一系列有益的工作(陆远忠等, 2007; 杨树新等, 2012; Hu et al., 2009).杨树新等(2012)利用有限元数值计算中的“降低单元刚度法”，模拟分析了中国大陆地区历史强震发生位置长距离跳迁的问题，但其应力场计算简单并未能重复真实历史地震序列的发震顺序.Hu等(2009, 2012, 2013)同时考虑构造应力场、断层相互作用以及下地壳和上地幔黏弹性松弛等因素对地震序列演化的影响，数值计算分析了1976年唐山大地震、1992年美国Landers地震以及2008年汶川大地震在数年尺度范围内的余震触发，但未能计算长时间尺度和大空间范围内的地震远距离迁移和应力场演化问题.
近年来，对于区域地震活动性趋势的研判，地震科学界认为在发展已有的经验预报的同时，需要逐步展开地震数值预报的探索性工作(刘启元和吴建春，2003；石耀霖等, 2013, 2018).其中石耀霖等(2018)提出了我国地震数值预测技术线路图的设想，其中指出对于中长期时间相关的地震数值预报，可以从一个特定时间点，给定初始应力状态，开始模拟重现历史大震序列，然后再继续计算未来可能的地震危险性，其关键问题为如何选取可靠的初始应力场.董培育等(2019)根据已有的历史地震信息及其他相关资料，提出了一种反演估算初始构造应力场的方法，估算了巴颜喀拉块体内的应力场，其计算结果可与实测应力场对比.本文将采用该方法，选取青藏高原及其邻区为研究对象，建立二维平面弹性有限元地质模型，并考虑静水压力作用，得到准三维应力场.在反演估算的震前初始应力的基础上，数值模拟区域应力场演化过程，借此重现历史大震发震过程，并尝试给出该区域未来的地震危险性概率预估.
1 数值计算模型及计算方法 1.1 数值计算模型本文计算模型范围为75°E—110°E，25°N—41°N，包括青藏高原及邻区的主要大型断裂带，但去除西南角的印度板块区域.边界条件采用现今GPS观测速率(Gan et al., 2007)在模型边界上的插值, 如图 2所示.已有研究表明高原内物质在印度板块的推挤作用下被侧向挤出，其中柔性下地壳运动速度高于脆性上地壳，对脆性上地壳有拖曳力作用，本文采用已有研究成果(董培育等, 2016, 2019)中的模型及材料参数进行计算，材料参数如表 1所示.
|
图 2 研究区域数值模型及网格剖分，修改自董培育等(2016, 2019) Fig. 2 Numerical model with mesh of the study region, modified from Dong et al.(2016, 2019) |
|
表 1 介质参数 Table 1 Material parameters |
应力场演化过程受多重因素影响，本文将此复杂问题简单化处理，简化为图 3所示：区域各处从某一初始状态开始，在构造活动加载作用下(假设在百年尺度范围内构造活动加载速率保持恒定)，其应力不断积累至强度极限，发生破裂并释放应力，随后继续加载构造应力，直至发生下一次地震，依次循环.在此过程中，可能受到周边区域发生地震或者其他地质事件产生的应力扰动增加或者减小.
|
图 3
应力场演化示意简图
其中，σs：断层破裂强度；σr：地震释放应力后的状态；σ0：初始应力； ：应力加载速率；σd：地震释放的应力即应力降；Δσ1：应力扰动减小量；Δσ2：应力扰动增加量.
Fig. 3
Sketch map of tectonic stress evolution
σs: strength of fault; σr: the state after the earthquake releases stress; σ0: initial stress; : stress increase rate; σd: stress drop; Δσ1: stress disturbance reduction; Δσ2: stress disturbance increase.
|
该过程涉及到的构造应力积累速率，在模型中施加GPS速率边界条件，通过数值模拟计算得到，本文直接利用已有的计算结果(董培育等, 2016, 2019)；地震造成的应力扰动(暂且仅考虑同震变化)，在地震破裂滑动模型已知的情况下，有较多成熟的计算方法(Wang et al., 2003; 张贝等，2015)，但本文为了提高计算效率，在二维平面模型中采用横向各向同性“杀伤”单元法，计算更加快速便捷(董培育等，2013).本文最关键问题是，如何选取初始应力，原则上应该通过观测手段获取，但以目前的科学水平和技术手段，几乎不可能直接观测到广大深部区域的应力(陈群策等，2012；王成虎等，2014)，对构造应力场认识仍存在很大的局限.
1.3 初始应力场估算方法如图 3所示，依据应力场演化过程可知，地壳内岩石层应力持续积累，直至达到其强度极限，发生破裂释放应力，为一次地震事件.那么历史地震序列信息至少揭示了不同区域上的历史应力状态，可以为初始应力场反演提供一些关键性参考.根据库仑-摩尔破裂准则和研究区内不同地区的岩石强度参数，可以推测，历史地震破裂区的震前应力至少达到了临界或半临界破裂状态，即岩石强度极限值.那么根据该次地震释放的应力降以及岩石强度值，则可以大致估计震源区震前应力状态.然而历史无震区，其应力水平必然未曾达到岩石强度极限值，不可直接估计，但是可参考其他相关研究资料和观测数据等，确定其初始应力下限值及上限值.因此，历史震源区的初始应力可限定某一个值，但历史无震区初始应力场尽管未知，但能大致给出其上下限范围.且历史地震破裂总是在断层的某个区间发生，至少说明在该区间上是高应力区或低强度区，而区间的两端是低应力区或者高强度区，才导致破裂不能持续下去.因此在反演估算初始应力场时，在强度相当的条件下，给出历史破裂区较高初始应力，而破裂区两端初始应力较低.另外，二维平面模型计算的平面应力场，考虑一定深度上的静岩压力作用，即可视为准三维应力场.本文选择发生在该研究区内，自1904年至2014年之间的较典型的26个M≥6.5破坏性强震进行计算，具体地震信息详见表 2.
|
|
表 2 26个历史大震信息 Table 2 Information of the 26 historical earthquakes |
基于库仑-摩尔破裂准则，可利用区域最大主应力Pmax与最小主应力Pmin比值之间的关系来判断岩石是否达到破裂条件(董培育等，2019)，如公式(1)所示：
|
(1) |
其中μ为岩石内摩擦系数，r0为根据内摩擦系数μ计算得到的破裂临界主应力比值，r为实际主应力比值，当r < r0时，安全；当r≥r0时，破裂.另外，本文提出危险度(risk degree, 简称rd)的概念，令rd=r/r0-1.0，则rd < 0，安全；rd≥0，危险.依据不同地区的地震类型以及岩石内摩擦系数μ，推断其对应的震源区在震前应该处于临界破裂状态，即rd=0或略大于0，则可确定震源区震前临界主应力值；无震区在震前长时间内rd＜0，并结合静水压大小、地应力测量数据以及构造活动资料等，可给出其rd上限值和下限值，即主应力的上下限值.具体应力场反演估算方法，详见董培育等(2019).
2 初始应力场估算由于本文所选历史大震序列的震源深度在10~30 km范围内，暂且选取10 km深度面为计算深度.根据前述方法粗略地反演估算出一种可能的1904年道孚地震震前的初始应力场，如图 4所示.10 km深部的静岩压力约为270 MPa(Pv=ρrock×g×h，其中ρrock=2700 kg·m-3)，反演得到的最小和最大水平主应力值在200~600 MPa之间.在高原南缘，受印度板块挤压作用强烈，区域主压应力水平较高，主要方向为NNE-SN向；高原北部受坚硬的塔里木盆地的阻挡作用，高原内部水平主压应力方向以NNE-NE向为主，但主应力值自南向北逐渐减小；在高原东缘，高原内侧向挤出的物质受到四川盆地的阻挡，最大水平主压应力方向在该地区发生偏转至近EW向，东南缘转为NW向.最大水平主压应力的方向基本与地表GPS观测到的形变方向一致.
|
图 4 反演估算出的一种可能1904年震前初始水平主应力场 Fig. 4 A estimated possible initial horizontal principal stress before the 1904 earthquake occurred |
需指出，这个估计出的初始应力场只是根据仅有的地震序列信息推测出的一种可能的模型，并非是真实的状态，只满足在历史地震震源区应力较高，其他地区应力相对较低的假设.
3 应力场及历史大震序列演化过程图 5给出了本文采用的应力场演化过程计算流程图，从震前初始应力场Sij_before(n)开始计算，判断某处应力是否超出其破裂临界值.若未达到破裂条件，则持续加载构造应力；若达到破裂条件，则计算同震应力变化量ΔSij(n)，则震后应力场为Sij_after(n)=Sij_before(n)-ΔSij(n).在该次地震震后应力场Sij_after(n)基础上，震间继续加载构造应力，以及可能的应力扰动，直至下次地震.依次循环计算，重复大震序列中的每次地震.
|
图 5 应力场演化过程计算流程图 Fig. 5 Calculation flow chart of the stress field evolution |
其中，历次地震释放的应力量，在本文弹性模型中，仅考虑同震变化，采用已经过验证的更加简洁高效的有限元横向各向同性“杀伤”单元法(董培育和石耀霖，2013)，数值计算结果可以与实际地震同震变化对比.
由于篇幅有限，我们不一一列举所有历史强震的同震变化，只以1920年的M8.0海原地震为例.利用有限元横向各向同性“杀伤”单元法，计算海原地震的同震位移和应力变化结果如图 6a和b所示，其主震震源区应力降ΔSxy约为1.0 MPa左右，位错4 m左右.实际上海原地震为多段破裂，前人研究给出其加权平均水平位移约为4 m(见表 2)，最大同震应力变化在兆帕量级，本文计算结果与此基本吻合(刘桂萍和傅征祥，2001).同时，我们也计算了海原地震主震破裂方向上的同震库仑破裂应力变化，如图 6c所示.图 7a、b和c分别为海原地震前后的危险度分布图以及地震前后的危险度变化图，由图可见，地震前震源区的危险度高于0，其他地区危险度均低于0，尤其块体区域非常安全，断层带不同区段危险度不同.海原地震之后，震源区应力释放，危险度大幅下降，而破裂区两端受应力扰动影响最大(图 6b)，库仑应力变化也为正值，且最大可达到兆帕量级(图 6c)，危险度明显上升(图 7c).但由于破裂区两端自身应力水平较低，即便应力扰动增加，仍未达到破裂极限，危险度仍低于0(图 7b).地震发生前后的危险度变化(图 7c)更进一步揭示了地震破裂区的危险度大幅下降，破裂区两端的值则大幅升高，其分布云图与库仑应力变化图类似，同为相对量.由此可见，利用库仑应力变化理论仅能够揭示哪些区域相对危险或者安全，而本文考虑背景构造应力场的计算，可以给出更加科学的值得参考的绝对的危险状态.
|
图 6 (a) 1920年海原地震产生的位错; (b) 1920年海原地震产生的应力降; (c)库仑破裂应力变化 Fig. 6 (a) Displacement caused by the 1920 Haiyuan earthquake; (b) Stress drop caused by the 1920 Haiyuan earthquake; (c) Coulomb failure stress change |
|
图 7 (a) 1920年海原地震震前的危险度分布图；(b) 1920年海原地震震后的危险度分布图；(c)地震前后的危险度变化图(图中红色表示危险区，蓝色表示安全区) Fig. 7 (a) Earthquake risk degree before the 1920 Haiyuan earthquake occurred; (b) Earthquake risk degree after the 1920 Haiyuan earthquake occurred; (c) Earthquake risk degree change before and after Haiyuan earthquake occurred (The red area indicates danger and the blue area indicates safe in the figure) |
|
图 8 不同历史地震震源区平均危险度值随时间变化曲线 Fig. 8 Time-varying curve of average earthquake risk value in different historical earthquake source regions |
图 8为不同历史地震震源区平均危险度值随时间变化曲线.整体而言所有历史地震，在其各自地震发生前，震源区平均危险度基本处于零线以下，但维持较高水平，且随时间推移，有逐步升高趋势，直至地震发生临界时间，危险度高出零线.震后应力释放，危险度亦大幅下降.其中在震间期，还可能受到周边地震的应力扰动影响以及构造应力加载的影响，致使其危险度有所调整，如1904年道孚地震震源区和1923年炉霍地震震源区可能受邻近的1973年道孚地震的影响，在1973年其各自的平均危险度均有明显变化；受2008年汶川大地震的影响，在2008年的时间线上，1933年叠溪震源区和1976年松潘震源区的平均危险度有明显下降，这两个区域更加安全，而2013年芦山震源区的平均危险度明显升高.另外1904年道孚地震、1923年炉霍地震、1955年康定地震以及1973年道孚地震，都发生在鲜水河断裂带上，其震源区位置非常接近甚至有重叠部分(见图 1)，导致在历史非震时间这些震源区的平均危险度也高出了零值.主要原因可能是在地震发生前后，震源区尤其是断层上的强度即介质属性实际上会发生变化，同震或震后短期内弱化，震后逐渐恢复强度(Taira et al., 2009)，然而本计算尚未考虑该因素，因此重复发生地震的地区或邻近地区，计算误差可能较大，平均危险度计算结果显示与实际情况有差异.
图 9为本文所选强震序列中最后一个地震——2014年于田地震后的危险度分布，可以看出受构造应力加载和历史地震共同作用的控制，断裂带的部分区段以及块体内部一些地区的危险度明显升高，另外一些地区的危险度明显下降，形成“安全岛”区域，如格尔木和玉树以东、道孚以北地区等.但这只是其中的一种可能的初始应力场算例计算所得结果，并非完全可靠.
|
图 9 2014年于田地震后的危险度分布 Fig. 9 Earthquake risk degree after the 2014 Yutian earthquake occurred |
完成历史地震序列的演化计算之后，可得到现今应力场状况，在此基础上，继续加载构造应力，可观察后续的地震活动性趋势.由于反演估算初始应力场时，对于模型中的历史无震区，其应力难以准确约束，仅能给定其上下限范围，在此范围内给出一个随机值.因此得到的初始应力场并非真实状态，只是一种在一定假设条件下的可能的算例，进而得到的现今应力场状况必然存在很大的不确定性，仅据此随机的单一算例来分析未来地震活动性趋势，可参考性较低.
因此本文选择采用Monte Carlo方法，进行随机试验，在合理范围内，随机产生千种可能的初始应力场算例，且保证每种算例都能够重复历史地震序列(董培育等，2019).如果有一些算例不能重复历史地震序列，则会微调初始模型，直至其可完全重复历史地震序列发震过程为止，每个算例的初始应力条件不同，对后续应力场演化也不尽相同.有些地区在多个算例中的危险度rd均大于0，说明在该地区的地震危险性可能较大；反之，有些地区在多个算例甚至全部算例中的危险度rd均小于0，则说明该地区可能为安全区.对于成功重复历史地震的千种算例，统计其不同地区在不同时间段的地震危险度rd＞0出现的频次，最终得到其危险性概率(频次/总数)分布，以此反映未来的地震活动性.计算结果如图 10所示.
|
图 10 不同时间点的危险性概率分布图 (a) 2014年；(b) 2114年.蓝色圆点表示从2014年3月至2019年3月以来发生的M>5.0地震的位置. Fig. 10 Risk probability distribution map at different time (a) 2014; (b) 2114. The blue dots show the location of the M>5.0 earthquake occurred from Mar. 2014 to Mar. 2019. |
图 10a和b分别为2014年于田地震后以及百年后——2114年的地震活动性概率分布图.在2014年于田地震之后，阿尔金断裂带西端与昆仑山断裂带交叉区域(于田地区)地震危险性概率较小，百年后阿尔金断裂带西南段和东北段局部地区危险性概率有所上升.且2014年之后在地震危险性概率较高地区发生了两次M>5.0地震，一定程度上验证了本文计算的可靠性.在祁连—海原断裂带上，2014年及百年后，该断裂带上除了历次地震(分别是1920海原、1927古浪及1932昌马大地震)造成的各自的破裂区两端有较低概率的危险性，其他地区仍然较为安全，说明该断裂带目前处于较稳定状态.
昆仑山断裂带历史地震较多，断裂带上未破裂的地区受到历史地震的影响较大.本文计算结果显示，百年后该断裂带上多数地区均有较高概率的地震危险性，昆仑山断裂带的平均滑动速率也较高，约为12.0 mm·a-1(Fu and Awata, 2007).地震地质资料以及GPS观测资料分析认为东昆仑断裂带东端存在较高地震危险性(邓起东等，2014；徐锡伟等, 2014, 2017).需要说明的是2017年8月8日，在东昆仑断裂带东端的隐伏分支断层——虎牙断裂带北段, 发生了MS7.0九寨沟地震.在该地震发生前，地震地质勘察尚未确定该断层的存在，本文模型中未包含此区段，仅有虎牙断裂带南段.虽然本模型中未考虑九寨沟地震发震断层，但计算结果显示在岷江断裂带区域附近存在约15%的地震危险性概率，一定程度上表明了该区域的危险性.
甘孜—玉树断裂带，发生过两次玛尼地震(分别是1973年和1997年)以及2010年玉树地震，这三次地震对整个断裂带影响较大.计算结果显示，现今，仅在玛尼地震破裂区东西两侧及玉树地区西侧小区域有最高约40%的危险性概率.百年后该断裂带上玛尼东西两侧附近区域均有较高概率危险性(最高可达60%)，在断裂带东南段出现约10%左右的危险性概率.显现出构造加载在该断裂带上的重要作用.野外地质调查、遥感影像解译及探槽研究等认为玉树断裂带东南段全新世以来活动较强，有强震可能性(石峰等，2013；吴中海等，2014).
鲜水河断裂带是活跃的左旋走滑断裂带，自1725—2008年间发生14次M>6.5以上的强震(Allen et al., 1991).本文计算结果显示，2014年该断裂带上，仅在道孚地区及其南侧与龙门山断裂带交叉地区，有小概率(约5%)的危险性，大部分区域为安全区.百年后，受构造活动加载作用，断裂带中段危险性概率明显升高，约至20%.该断裂带大部分区段发生过强震，其应力释放的同时，断层内岩石承载力下降，摩擦系数降低.且地震地质学及重磁资料表明鲜水河断裂带东南段存在高异常，地震活动性值得关注(刘冠中等，2013；李大虎等，2015).2014年至今在断裂带中南段发生了两次中强地震.
嘉黎—雅鲁藏布江断裂带，本文计算结果显示，在2014年于田地震后，仅在那曲附近区域有小概率(约10%左右)的地震危险性.百年后，在断裂带的西侧大部分区段危险性概率增加至最高约30%左右.同样，在现今一段时间内，喀喇昆仑断裂带处于安全期，百年后，仅在中段有小概率(约5%)的地震危险性.实际上这两条断裂带及其内部拉萨地块内，发育较多正断层.由于正断层活动需要的应力水平相比于走滑型和逆掩型断层小很多，据地质学资料分析认为该区域相较于其他地块更易发生地震(傅征祥等，2010).受限于本文的二维平面模型计算拉张或逆冲型地震误差较大，因此该区域的计算结果可能存在误差.
龙门山断裂带，百年内发生过两次强震，2008汶川MS8.0地震和2013芦山MS7.0地震.本文结果显示，现今地震危险性概率较高地区主要集中在龙门山断裂带的东北段和西南段部分地区，尤其与鲜水河断裂带交叉区域，且在百年后，这些区域地震危险性概率值继续升高.汶川地震和芦山地震破裂带中间存在约50 km的空区，陈运泰等(2013)认为汶川地震之后，龙门山断裂带的东北段和西南段很有可能发生强余震，芦山地震的发生初步证实了其西南段危险性估计，但应力并未完全释放.宽频带数字地震波形资料计算得到的视应力结果表明汶川地震后，龙门山断裂带东北段视应力下降，然而西南部仍维持较高水平(易桂喜等, 2013).杜方等(2015)认为龙门山断裂带的中北段历史地震释放能量不充分，仍有危险性.计算结果显示，龙门山西北侧的岷江断裂带、虎牙断裂带，以及鲜水河西侧的理塘断裂带，这些地区的地震危险性在现今一段时间内维持较低水平，百年后，受构造应力加载影响，危险性概率有所升高.这些断裂带均处于构造活跃、地震频发的川滇地区，未来的地震活动性不容忽视.
喜马拉雅断裂带是构造活动最为复杂的板间逆冲巨型断裂带，计算结果显示，在2014年于田地震之后，该断裂带均处于安全状态.百年后，仅在1905年和1934年大地震震源区的地震危险性极高.2015年4月25日，在喜马拉雅俯冲带上发生了MS8.1低角度逆冲型强震，震源位置28.1°N，84.7°E.由于本文计算过程中，此次地震尚未发生，因此该地震尚未参与历史地震序列的演化计算过程.此外，该断裂带位于本文计算模型的南部边界区，受有限元数值计算的边界效应影响较大.同时，在二维平面模型中，一方面对于逆掩型大地震的同震变化计算误差较大；另一方面，在该复杂的逆冲环境区域中，垂向应力占主导作用，而本文中的垂向应力采用的是均匀的静水压力，误差更为突出.因此该断裂带上的计算结果误差较大，对未来的强震危险性评估可参考性较小.
综上，本文计算结果表明未来地震活动性概率较高区域主要集中在昆仑山断裂带中东段区域、玛尼西北地区，青藏高原东缘(龙门山断裂带南段与鲜水河断裂带交叉区域)，以及嘉黎断裂带那曲段等.
5 结论与讨论本文充分利用了历史大震序列信息，基于岩石库仑-摩尔破裂准则，反演估算了区域初始构造应力场，在此基础上重现了历史大震发震过程.其中对于初始应力场的反演做了大量的试验工作，在推测的合理范围内随机生成数千种可能的有差异性的初始应力场算例，所有算例都能够满足历史地震序列的有序发生，但最终计算结果有所差异，对这些结果进行数理统计，仅能给出中长期尺度的大致估计，时间窗口难以确定.虽不能确切计算哪里是危险的，但是缩小了范围，给出危险性可能较高的地区，对未来的地震发展趋势有一定程度的估计.本文通过综合考虑背景初始构造应力场、构造应力增长速率以及地震断层错动产生的同震应力变化，可以更加科学地判断地震危险/安全区，给出未来百年范围内的地震危险性概率图, 为中长期地震危险性评估提供科学定量的参考依据.根据本文结果，从中长期尺度而言，未来地震活动危险性概率较高的地区主要集中在玛尼、格尔木以及二者的中间地带，玉树西南部，龙门山断裂带的南部附近区域，鲜水河断裂带周围地区，这些地区的未来强震活动性值得关注.
但本文计算中仍存在问题值得深入研究：(1)初始应力场的估算较为粗略.主要根据库仑-摩尔破裂准则，以及历史地震释放的应力降以及其他地应力观测资料等来进行反演估算.对于历史地震的选取，由于本文研究区域较大，受计算量的限制，仅采用了26个较典型的历史大震，仍有很多历史地震没有参与计算，例如，最近的2015年尼泊尔大地震，以及2017年九寨沟地震等，以及MW < 6.5的历史强震以及历史强震的大量余震等，虽然小震释放的能量小，对周边区域产生的应力扰动也较小，但不能排除这些地震对计算结果的影响.在开展局部区域精细研究中，可考虑更多的小震对计算结果的影响.(2)本文采用简捷的二维模型计算，给定均匀的垂向应力(静水压)作为“假”三维应力场进行分析，虽然可行，但终究未能考虑真实区域物质结构的横向和纵向不均匀性，以及黏性下地壳在应力场演化中的动力学作用.尤其对于拉张和逆冲型地震，受垂向应力影响较大，因此给定的均匀的垂向应力值会导致计算误差较大.(3)计算走滑型地震产生的同震变化，二维平面模型基本可行，然而逆冲型或拉张型地震，计算误差较大.因此本文的计算结果可能在喜马拉雅俯冲带和龙门山断裂带上可靠性较低.未来仍需要建立三维模型，才能够更好的计算拉张或逆冲型地震的同震变化，以及考虑黏弹性模型，计算震后长期的应力调整，以期更接近实际情况.
综上所述，该工作中的任何一个计算环节都可能存在误差，如模型材料参数、边界条件、断层强度等.尤其对于初始应力场的估计是比较粗略的，但是在目前仅有的观测数据基础上，该工作是一种尝试和探索, 比常规的零初始应力场或均匀初始应力场的计算更进一步.总言之，本文是对中长期时间相关地震数值预报的探索性试验工作，所得结果有一定的参考价值.
致谢 感谢审稿专家对本文提出的宝贵建议.
Allen C R, Luo Z L, Qian H, et al. 1991. Field study of a highly active fault zone:The Xianshuihe fault of southwestern China. Geological Society of America Bulletin, 103(9): 1178-1199. DOI:10.1130/0016-7606(1991)103<1178:FSOAHA>2.3.CO;2 |
Chen Q C, Feng C J, Meng W, et al. 2012. Analysis of in situ stress measurements at the northeastern section of the Longmenshan fault zone after the 5.12 Wenchuan earthquake. Chinese Journal of Geophysics (in Chinese), 55(12): 3923-3932. DOI:10.6038/j.issn.0001-5733.2012.12.005 |
Chen S F, Wilson C J L, Deng Q D, et al. 1994. Active faulting and block movement associated with large earthquakes in the Min Shan and Longmen Mountains, northeastern Tibetan Plateau. Journal of Geophysical Research:Solid Earth, 99(B12): 24025-24038. DOI:10.1029/94JB02132 |
Chen Y T, Yang Z X, Zhang Y, et al. 2013. From 2008 Wenchuan earthquake to 2013 Lushan earthquake. Scientia Sinica Terrae (in Chinese), 43(6): 1064-1072. DOI:10.1360/zd-2013-43-6-1064 |
Deng Q D, Zhang P Z, Ran Y K, et al. 2003. Basic characteristics of active tectonics of China. Science in China Series D:Earth Sciences, 46(4): 356-372. DOI:10.1360/03yd9032 |
Deng Q D, Cheng S P, Ma J, et al. 2014. Seismic activities and earthquake potential in the Tibetan Plateau. Chinese Journal of Geophysics (in Chinese), 57(7): 2025-2042. DOI:10.6038/cjg20140701 |
Dong P Y, Cheng H H, Shi Y L, et al. 2019. Numerical inversion of regional initial tectonic stress based on Monte Carlo method-A case study of Bayan Har block. Chinese Journal of Geophysics (in Chinese), 62(8): 2858-2870. DOI:10.6038/cjg2019M0393 |
Dong P Y, Hu C B, Shi Y L. 2016. Numerical simulation of long-term deformation of Tibetan Plateau and surrounding area. Seismology and Geology (in Chinese), 38(2): 410-422. DOI:10.3969/j.issn.0253-4967.2016.02.014 |
Dong P Y, Shi Y L. 2013. A discussion on "The mechanism of long-distance jumping and the migration of main active areas for strong earthquakes occurred in the Chinese continent"-Transverse isotropic "wounded element" is a better method. Chinese Journal of Geophysics (in Chinese), 56(6): 2133-2139. DOI:10.6038/cjg20130633 |
Du F, Su Y J, Long F, et al. 2015. The M7.0 Lushan, the M8.0 Wenchuan earthquakes and the history earthquakes in the eastern border area of the Bayan Har block. Progress in Geophysics (in Chinese), 30(4): 1569-1579. DOI:10.6038/pg20150411 |
England P, McKenzie D. 1982. A thin viscous sheet model for continental deformation. Geophysical Journal International, 70(2): 295-321. DOI:10.1111/j.1365-246X.1982.tb04969.x |
Fu B, Awata Y. 2007. Displacement and timing of left-lateral faulting in the Kunlun Fault Zone, northern Tibet, inferred from geologic and geomorphic features. Journal of Asian Earth Sciences, 29(2-3): 253-265. DOI:10.1016/j.jseaes.2006.03.004 |
Fu Z X, Liu G P, Chen Q F. 2001. Dynamic analysis on interaction between the Haiyuan-Gulang-Changma great earthquake in the north boundary of the Tibetan Plateau. Seismology and Geology (in Chinese), 23(1): 35-42. |
Fu Z X, Lu X J, Zhang Y S, et al. 2010. Seismicity and precursory heat anormaly in the Lasa and Qiangtang Blocks. Earthquake (in Chinese), 30(1): 20-27. |
Gan W J, Zhang P Z, Shen Z K, et al. 2007. Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements. Journal of Geophysical Research:Solid Earth, 112(B8): B08416. DOI:10.1029/2005JB004120 |
Hu C B, Zhou Y J, Cai Y E. 2009. A new finite element model in studying earthquake triggering and continuous evolution of stress field. Science in China Series D:Earth Sciences, 52(7): 994-1004. DOI:10.1007/s11430-009-0082-3 |
Hu C B, Cai Y E, Wang Z M. 2012. Effects of large historical earthquakes, viscous relaxation, and tectonic loading on the 2008 Wenchuan earthquake. Journal of Geophysical Research:Solid Earth, 117(B6): B06410. DOI:10.1029/2011JB009046 |
Hu C B, Cai Y E, Liu M, et al. 2013. Aftershocks due to property variations in the fault zone:A mechanical model. Tectonophysics, 588: 179-188. DOI:10.1016/j.tecto.2012.12.013 |
Jones L M, Han W B, Hauksson E, et al. 1984. Focal mechanisms and aftershock locations of the Songpan earthquakes of August 1976 in Sichuan, China. Journal of Geophysical Research:Solid Earth, 89(B9): 7697-7707. DOI:10.1029/JB089iB09p07697 |
Li D H, Ding Z F, Wu P P, et al. 2015. The deep seismogenic environment of the southeastern section of the Xianshuihe fault zone and the 2014 Kangding MS6.3 earthquake. Chinese Journal of Geophysics (in Chinese), 58(6): 1941-1953. DOI:10.6038/cjg20150610 |
Liu C L, Zheng Y, Ge C, et al. 2013. Rupture process of the MS7.0 Lushan earthquake, 2013. Science China Earth Sciences, 56(7): 1187-1192. DOI:10.1007/s11430-013-4639-9 |
Liu G P, Fu Z X. 2001. A study on the great Gulang earthquake triggered probably by static stress change resulting from the great Haiyuan earthquake. Chinese Journal of Geophysics (in Chinese), 44(Suppl.): 107-115. DOI:10.3321/j.issn:0001-5733.2001.z1.014 |
Liu G Z, Ma J, Zhang H X, et al. 2013. Study on activity features of Xianshuihe fault zone with fault creep and short baseline observation for the last 20 years. Chinese Journal of Geophysics (in Chinese), 56(3): 878-891. DOI:10.6038/cjg20130317 |
Liu Q Y, Wu J C. 2003. On numerical forecast of earthquakes-thinking about the strategy for promoting earthquake prediction. Earth Science Frontiers (in Chinese), 10(S1): 217-224. DOI:10.3321/j.issn:1005-2321.2003.z1.030 |
Lu Y Z, Ye J D, Jiang C, et al. 2007. 3D numerical simulation on the mechanism of precursory seismicity pattern before strong earthquake in China. Chinese Journal of Geophysics (in Chinese), 50(2): 499-508. |
Luo G, Liu M. 2010. Stress evolution and fault interactions before and after the 2008 Great Wenchuan earthquake. Tectonophysics, 491(1-4): 127-140. DOI:10.1016/j.tecto.2009.12.019 |
Molnar P, Tapponnier P. 1975. Cenozoic tectonics of Asia:Effects of a continental collision. Science, 189(4201): 419-426. DOI:10.1126/science.189.4201.419 |
Molnar P, Deng Q D. 1984. Faulting associated with large earthquakes and the average rate of deformation in central and eastern Asia. Journal of Geophysical Research:Solid Earth, 89(B7): 6203-6227. DOI:10.1029/JB089iB07p06203 |
Papadimitriou E, Wen X Z, Karakostas V, et al. 2004. Earthquake Triggering along the Xianshuihe Fault Zone of Western Sichuan, China. Pure and Applied Geophysics, 161(8): 1683-1707. DOI:10.1007/s00024-003-2471-4 |
Shen Z K, Wan Y G, Gan W J, et al. 2003. Viscoelastic triggering among large earthquakes along the east Kunlun fault system. Chinese Journal of Geophysics (in Chinese), 46(6): 786-795. |
Shi F, Li A, Yang X P, et al. 2013. Research on late quaternary activity of the southeastern segment of Ganzi-Yushu fault zone. Seismology and Geology (in Chinese), 35(1): 50-63. DOI:10.3969/j.issn.0253-4967.2013.01.004 |
Shi Y L, Zhang B, Zhang S Q, et al. 2013. Numerical earthquake prediction. Physical (in Chinese), 42(4): 237-255. |
Shi Y L, Sun Y Q, Luo G, et al. 2018. Roadmap for earthquake numerical forecasting in China-reflection on the tenth anniversary of Wenchuan earthquake. Chinese Science Bulletin (in Chinese), 63(19): 1865-1881. DOI:10.1360/N972018-00335 |
Taira T, Silver P G, Niu F L, et al. 2009. Remote triggering of fault-strength changes on the San Andreas fault at Parkfield. Nature, 461(7264): 636-639. DOI:10.1038/nature08395 |
Wan Y G, Shen Z K, Zeng Y H, et al. 2007. Evolution of cumulative coulomb failure stress in northeastern Qinghai-Xizang (Tibetan) Plateau and its effect on large earthquake occurrence. Acta Seismologica Sinica (in Chinese), 29(2): 115-129. |
Wang C H, Song C K, Guo Q L, et al. 2014. Stress build-up in the shallow crust before the Lushan Earthquake based on the in-situ stress measurements. Chinese Journal of Geophysics (in Chinese), 57(1): 102-114. DOI:10.6038/cjg20140110 |
Wang R, He G Q, Yin Y Q, et al. 1980. A mathematical simulation for the pattern of seismic transference in north China. Acta Seismological Sinica (in Chinese), 2(1): 32-42. |
Wang R J, Martín F L, Roth F. 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust-FORTRAN programs EDGRN/EDCMP. Computers & Geosciences, 29(2): 195-207. |
Wu Z H, Zhou J C, Feng H, et al. 2014. Active faults and earthquake around Yushu in eastern Tibetan Plateau. Geological Bulletin of China (in Chinese), 33(4): 419-469. |
Xu L S, Chen Y T. 2004. Temporal and spatial rupture process of the great Kunlun Mountain Pass earthquake of November 14, 2001 from the GDSN long period waveform data. Science in China Series D:Earth Sciences, 48(1): 112-122. |
Xu X W, Cheng W, Xu C, et al. 2014. Discussion on block kinematic model and future themed areas for earthquake occurrence in the Tibetan Plateau:inspiration from the Ludian and Jinggu earthquakes. Seismology and Geology (in Chinese), 36(4): 1116-1134. DOI:10.3969/j.issn.0253-4967.2014.04.015 |
Xu X W, Chen G H, Wang Q X, et al. 2017. Discussion on seismogenic structure of Jiuzhaigou earthquake and its implication for current strain state in the southeastern Qinghai-Tibet Plateau. Chinese Journal of Geophysics (in Chinese), 60(10): 4018-4026. DOI:10.6038/cjg20171028 |
Yang S X, Lu Y Z, Chen L W, et al. 2012. The mechanism of long-distance jumping and the migration of main active areas for strong earthquakes occurred in the Chinese continent. Chinese Journal of Geophysics (in Chinese), 55(1): 105-116. DOI:10.6038/j.issn.0001-5733.2012.01.010 |
Yi G X, Wen X Z, Xin H, et al. 2013. Stress state and major-earthquake risk on the southern segment of the Longmen Shan fault zone. Chinese Journal of Geophysics (in Chinese), 56(4): 1112-1120. DOI:10.6038/cjg20130407 |
Yue H, Zhang Z Q, Chen J. 2008. Interaction between adjacent left-lateral strike-slip faults and thrust faults:the 1976 Songpan earthquake. Chinese Science Bulletin, 53(16): 2520-2526. DOI:10.1007/s11434-008-0210-z |
Zhang B, Zhang H, Shi Y L. 2015. Equivalent-bodyforce approach on modeling elastic dislocation problem using finite element method. Chinese Journal of Geophysics (in Chinese), 58(5): 1666-1674. DOI:10.6038/cjg20150518 |
Zhang P Z, Deng Q D, Zhang Z Q, et al. 2013. Active faults, earthquake hazards and associated geodynamic processes in continental China. Science China Earth Sciences (in Chinese), 43(10): 1607-1620. |
Zhang Y, Xu L S, Chen Y T, et al. 2014. Fast inversion for the rupture process of the 12 February 2014 Yutian MW6.9 earthquake:Discussion on the impacts of focal mechanisms on rupture process inversions. Acta Seismologica Sinica (in Chinese), 36(2): 159-164. DOI:10.3969/j.issn.0253-3782.2014.02.001 |
陈群策, 丰成君, 孟文, 等. 2012. 5.12汶川地震后龙门山断裂带东北段现今地应力测量结果分析. 地球物理学报, 55(12): 3923-3932. DOI:10.6038/j.issn.0001-5733.2012.12.005 |
陈运泰, 杨智娴, 张勇, 等. 2013. 从汶川地震到芦山地震. 中国科学:地球科学, 43(6): 1064-1072. DOI:10.1360/zd-2013-43-6-1064 |
邓起东, 张培震, 冉勇康, 等. 2002. 中国活动构造基本特征. 中国科学:D辑, 32(12): 1020-1030. DOI:10.3969/j.issn.1674-7240.2002.12.007 |
邓起东, 程绍平, 马冀, 等. 2014. 青藏高原地震活动特征及当前地震活动形势. 地球物理学报, 57(7): 2025-2042. DOI:10.6038/cjg20140701 |
董培育, 程惠红, 石耀霖, 等. 2019. 基于Monte Carlo方法数值反演区域初始构造应力场——以巴颜喀拉块体为例. 地球物理学报, 62(8): 2858-2870. DOI:10.6038/cjg2019M0393 |
董培育, 胡才博, 石耀霖. 2016. 青藏高原及周边区域地表长期变形数值模拟. 地震地质, 38(2): 410-422. DOI:10.3969/j.issn.0253-4967.2016.02.014 |
董培育, 石耀霖. 2013. 关于"用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移"的讨论——横向各向同性"杀伤单元"才是更好的途径. 地球物理学报, 56(6): 2133-2139. DOI:10.6038/cjg20130633 |
杜方, 苏有锦, 龙锋, 等. 2015. 芦山7.0级, 汶川8.0级地震与巴颜喀拉地块东缘区域历史地震. 地球物理学进展, 30(4): 1569-1579. DOI:10.6038/pg20150411 |
傅征祥, 刘桂萍, 陈棋福. 2001. 青藏高原北缘海原、古浪、昌马大地震间相互作用的动力学分析. 地震地质, 23(1): 35-42. DOI:10.3969/j.issn.0253-4967.2001.01.004 |
傅征祥, 吕晓健, 张元生, 等. 2010. 拉萨和羌塘地块地震活动和热前兆现象及其机理研究. 地震, 30(1): 20-27. DOI:10.3969/j.issn.1000-3274.2010.01.003 |
李大虎, 丁志峰, 吴萍萍, 等. 2015. 鲜水河断裂带南东段的深部孕震环境与2014年康定MS6.3地震. 地球物理学报, 58(6): 1941-1953. DOI:10.6038/cjg20150610 |
刘成利, 郑勇, 葛粲, 等. 2013. 2013年芦山7.0级地震的动态破裂过程. 中国科学:地球科学, 43(6): 1020-1026. DOI:10.1360/zd-2013-43-6-1020 |
刘桂萍, 傅征祥. 2001. 海原大地震对古浪大地震的静应力触发研究. 地球物理学报, 44. DOI:10.3321/j.issn:0001-5733.2001.z1.014 |
刘冠中, 马瑾, 张鸿旭, 等. 2013. 二十年来蠕变和短基线观测反映的鲜水河断裂带活动特征. 地球物理学报, 56(3): 878-891. DOI:10.6038/cjg20130317 |
刘启元, 吴建春. 2003. 论地震数值预报——关于我国地震预报研究发展战略的思考. 地学前缘, 10(S1): 217-224. DOI:10.3321/j.issn:1005-2321.2003.z1.030 |
陆远忠, 叶金铎, 蒋淳, 等. 2007. 中国强震前兆地震活动图像机理的三维数值模拟研究. 地球物理学报, 50(2): 499-508. DOI:10.3321/j.issn:0001-5733.2007.02.022 |
沈正康, 万永革, 甘卫军, 等. 2003. 东昆仑活动断裂带大地震之间的黏弹性应力触发研究. 地球物理学报, 46(6): 786-795. DOI:10.3321/j.issn:0001-5733.2003.06.010 |
石峰, 李安, 杨晓平, 等. 2013. 甘孜-玉树断裂带东南段晚第四纪活动性研究. 地震地质, 35(1): 50-63. DOI:10.3969/j.issn.0253-4967.2013.01.004 |
石耀霖, 张贝, 张斯奇, 等. 2013. 地震数值预报. 物理, 42(4): 237-255. |
石耀霖, 孙云强, 罗纲, 等. 2018. 关于我国地震数值预报路线图的设想——汶川地震十周年反思. 科学通报, 63(19): 1865-1881. |
万永革, 沈正康, 曾跃华, 等. 2007. 青藏高原东北部的库仑应力积累演化对大地震发生的影响. 地震学报, 29(2): 115-129. DOI:10.3321/j.issn:0253-3782.2007.02.001 |
王成虎, 宋成科, 郭启良, 等. 2014. 利用原地应力实测资料分析芦山地震震前浅部地壳应力积累. 地球物理学报, 57(1): 102-114. DOI:10.6038/cjg20140110 |
王仁, 何国琦, 殷有泉, 等. 1980. 华北地区地震迁移规律的数学模拟. 地震学报, 2(1): 32-42. |
吴中海, 周春景, 冯卉, 等. 2014. 青海玉树地区活动断裂与地震. 地质通报, 33(4): 419-469. DOI:10.3969/j.issn.1671-2552.2014.04.003 |
许力生, 陈运泰. 2004. 从全球长周期波形资料反演2001年11月14日昆仑山口地震时空破裂过程. 中国科学D辑:地球科学, 34(3): 256-264. |
徐锡伟, 程伟, 许冲, 等. 2014. 青藏高原块体运动模型与地震活动主体地区讨论:鲁甸和景谷地震的启示. 地震地质, 36(4): 1116-1134. DOI:10.3969/j.issn.0253-4967.2014.04.015 |
徐锡伟, 陈桂华, 王启欣, 等. 2017. 九寨沟地震发震断层属性及青藏高原东南缘现今应变状态讨论. 地球物理学报, 60(10): 4018-4026. DOI:10.6038/cjg20171028 |
杨树新, 陆远忠, 陈连旺, 等. 2012. 用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移. 地球物理学报, 55(1): 105-116. DOI:10.6038/j.issn.0001-5733.2012.01.010 |
易桂喜, 闻学泽, 辛华, 等. 2013. 龙门山断裂带南段应力状态与强震危险性研究. 地球物理学报, 56(4): 1112-1120. DOI:10.6038/cjg20130407 |
岳汉, 张竹琪, 陈永顺. 2008. 相邻左旋走滑和逆冲断层之间的相互作用:1976年松潘震群. 科学通报, 53(13): 1582-1588. DOI:10.1360/csb2008-53-13-1582 |
张贝, 张怀, 石耀霖. 2015. 有限元模拟弹性位错的等效体力方法. 地球物理学报, 58(5): 1666-1674. DOI:10.6038/cjg20150518 |
张培震, 邓起东, 张竹琪, 等. 2013. 中国大陆的活动断裂、地震灾害及其动力过程. 中国科学:地球科学, 43(10): 1607-1620. |
张勇, 许力生, 陈运泰, 等. 2014. 2014年2月12日于田MW6.9地震破裂过程初步反演:兼论震源机制对地震破裂过程反演的影响. 地震学报, 36(2): 159-164. DOI:10.3969/j.issn.0253-3782.2014.02.001 |