2019, Vol. 62
2. 中国科学院计算地球动力学重点实验室, 中国科学院大学, 北京 100049;
3. 同济大学, 海洋地质国家重点实验室, 上海 200092
2. Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences, Beijing 100049, China;
3. State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China
地震的发生通常认为是地下岩体所承受的应力超出其极限而发生破裂和能量释放的过程.构造地震学研究表明区域构造应力场往往对地震活动性具有控制作用，地震沿着已有的活动断裂带发生，当活动断裂带所承受的应力环境达到极限强度，断层失稳破裂滑动(Scholz, 1998).近几十年来，不少学者对构造应力场展开测量和分析研究，以期探讨地震与应力的相互作用等.例如，通过钻孔测量方法(陈群策等, 2012; 孟文等, 2017)和水压致裂法(郭啟良等, 2009; 王成虎等, 2014)等获取不同震区的浅部应力状态来分析区域构造应力场，进而探讨邻近区域的地震活动性.Liao等(2003)对昆仑山地震破裂区附近开展了浅部数十米的地应力测量，发现地震发生后测点上最大水平主压应力显著降低.2008年汶川地震后，许多学者(郭啟良等, 2009; 陈群策等, 2012; 杜建军等, 2013; Meng et al., 2015)在龙门山断裂带及其周边区域的200~500 m深度处开展了一系列原地应力测量工作，获得了区域最大水平主应力方向和地应力积累数据，分析认为2013年4月20日芦山MW6.7地震的发生与西南段地应力高度集中有关.孟文等(2017)在拉萨块体东部地区开展地表浅层百米地应力测量，认为该区地壳处于临界破裂状态，2017年11月18日在米林地区发生了MS6.9地震.
虽然，地应力测量能够揭示浅部地壳应力水平，可为地震危险性评估提供一定的定量化参考.但是，陆壳大震通常发生在十几公里深度处，甚至部分地震达到几十公里(Prieto et al., 2012)，仅仅使用有限的浅部地应力测量数据明显不足以描述深部孕震应力环境.一方面，地应力测量深度通常在数百米范围内且容易受局部岩体结构影响.例如，孟文等(2017)在拉萨块体东部地应力测量结果显示朗县测点应力值约为相距百余公里的林芝测点的1.5~3倍.另一方面，地应力随深度变化复杂，中国大陆浅层地应力测量数据揭示百米深度内地壳水平主应力偏高，垂向应力偏低，应力结构易以逆冲为主；深度约200 m下，垂向应力不断增加，超过最小水平主应力，易形成走滑型应力结构；当达到一定深度(约1 km)以下，垂向应力超过最大水平主应力，则易转变为拉张型应力结构为主(杨树新等, 2012b).同时，地应力测点相对较少，因此仅根据个别测点观测结果难以对区域应力状况做出准确判断.另外，虽然也可通过震源机制解信息推测深部震源区域的最大主应力方向(崔效峰等，2005)，但却很难确定应力的大小.
随着高性能数值计算的快速发展，对于孕震环境以及地震活动性的数值模拟研究日渐丰富(杨树新等, 2012a; Liu et al., 2016; 庞亚瑾等, 2017; Pang et al., 2018)，地震数值预测概念不断被提出(刘启元和吴建春, 2003; 石耀霖等, 2013, 2018).数值计算中对于区域初始构造应力场的选取尤为重要，然而，目前有限的浅部地应力数据和深部仅能描述地应力方向的资料不足以支撑数值计算的需求.据此，不少学者对初始构造应力场的选择采用了均匀估计值(胡才博等, 2009; 朱守彪等, 2017)，或者视初始应力值为零(程惠红等, 2014; 柳畅等, 2014; Liu et al., 2016)，一定程度上可以应用于应力场演化过程计算和解释孕震环境等，但是仍有其局限性，给出的往往是应力的变化量.
本文在前人研究基础上，以巴颜喀拉块体为例，提出利用Monte Carlo方法，结合区域历史大震信息、地质背景及GPS等观测数据，以库仑-摩尔破裂准则为计算依据来反演该块体的初始构造应力场.同时，将反演过程中不确定部分限定在一个合理的上下限范围内进行独立的重复性随机试验，运用统计学方法给出1997玛尼地震前区域初始构造应力场.
1 数值计算模型近20年来巴颜喀拉块体及其边界断裂带陆续发生了7次中强震(见图 1)，是青藏高原内地震最活跃的区域之一.不少学者(闻学泽等, 2011; 邓起东等, 2014)研究认为巴颜喀拉块体仍是青藏高原地震活动的主体区域，特别是其东边界和北边界仍存在较高地震危险性.据此，本文以巴颜喀拉块体及其周边区域为研究区域，主要包括：昆仑山断裂带、甘孜—玉树断裂带(包括鲜水河断裂带西北段)和龙门山断裂带，如图 2中线框所示.据前人研究(王辉等, 2006; 曹建玲等, 2009)，采用线弹性本构方程，构建了有限元数值计算地质模型(断裂带宽度~10 km)，给出了各个块体和断层参数(见表 1，断裂带弹性模量较周边区域弱).以现今的GPS观测速率(Gan et al., 2007)为位移边界的约束条件.
|
图 1 巴颜喀拉块体及邻区的地质构造背景 B1—B4分别为：柴达木块体、巴颜喀拉块体(图中粉色区域)、羌塘块体、拉萨块体；F1—F3分别为：昆仑山断裂带、甘孜—玉树断裂带、龙门山断裂带；红色震源球为1997—2014年间MW6.5以上的强震震源机制；深蓝色箭头为GPS观测速度场(引自Gan et al., 2007). Fig. 1 Tectonic sketch map of Bayan Har block and surrounding regions B1—B4 : Qaidam block, Bayan Har block (pink region), Qiangtang block, Lhasa block; F1—F3: Kunlun fault, Garzê-Yushu fault, Longmenshan fault; Red beach balls denote the focal mechanism of earthquakes larger than MW6.5 occurred between 1997—2014; Dark blue arrows represent the GPS velocity (from Gan et al., 2007). |
|
图 2 研究区域的数值计算模型及剖分网格(改自董培育等, 2016) Fig. 2 Numerical model and mesh of the study area (modified from Dong et al., 2016) |
|
表 1 介质参数表 Table 1 Material parameter |
现有研究(Royden et al., 1997; 曹建玲等, 2009; 朱守彪和石耀霖, 2005)指出在印度—欧亚板块的长期碰撞挤压作用下，高原物质整体东向—东南向运动，现今高原运动形变场主要是由水平运动速率不均匀的柔性下地壳拖动脆性上地壳而形成的.据此，我们通过加载体力来实现柔性下地壳对脆性上地壳的拖曳力作用，以期将计算得出的现今地表形变与GPS实测结果很好地吻合，进而获得区域年应变增长速率(董培育等, 2016).
2 计算方法地震视为断层在长期构造活动作用下应力持续积累直至超出其承受极限而发生错动并释放应力的过程.那么，历史地震序列在一定程度上可以揭示震源破裂区在地震前后的应力状态变化.据此，基于库仑-摩尔破裂准则，不同的构造应力场会引起不同类型的断层错动，反之，不同的震源特征可大致反推构造应力场.由此可为推测震前应力场提供一些参考.
2.1 库仑-摩尔破裂准则库仑-摩尔(Coulomb-Mohr)破裂准则可很好地描述壳内岩石应力状态与破裂关系，如式(1).
|
(1) |
其中，τ为断层面上剪应力，μ为摩擦系数(也为岩石摩擦强度)，τ0为岩石内聚力，σn为面上正应力.由于破裂往往发生在已有断层上，根据拜尔利定律(Byerlee′s Law)，当σn < 200 MPa时，τ0=0，μ=0.85；当σn>200 MPa时，τ0=50 MPa，μ=0.6，断层发生滑动，产生破裂.因内聚力与最大主应力相比可忽略(Goodman, 1989)，在此作为初步估算暂不考虑内聚力τ0.则式(1)可变为式(2)
|
(2) |
判断准则为：当断层上应力状态为最大主应力σ′1，最小主应力σ′3，应力摩尔圆未碰触破裂线；当差应力(σ1-σ3)增大，即最大主应力增大到σ1，最小主应力减小到σ3时，摩尔圆与破裂线相切，岩石可能发生破裂，见图 3.其中，φ为断层摩擦角，摩擦系数μ等于tanφ.
|
图 3 应力摩尔圆示意图 Fig. 3 Sketch map of Mohr′s stress circle |
据区域历史地震记录，我们对初始应力的估计分为以下两种情况：(a)历史地震震源区的初始应力估计；(b)历史无震区的初始应力估计.
如图 4a所示为历史上发生中强震的地区：在某个地震发生前，该地震的震源破裂区应力水平至少达到了断层临近破裂的状态或者半临界状态，据地震的应力降、震源区断层摩擦强度以及区域构造应力场的加载速率，我们可估计出震前初始状态的应力值.其中，在本文初步模型中的地震应力降为主震静态应力降，可以采用有限元横向各向同性“杀伤单元”法(董培育和石耀霖，2013)来计算；震源区断层摩擦强度，则根据不同区域构造活动性质来获得，活动频繁区域对应较低摩擦强度；区域构造应力增长速率可由区域应变率(董培育等, 2016)通过弹性力学胡克定律(如式(3)所示)计算得到.
|
图 4
初始应力估计原则示意图
(a)历史地震震源区的初始应力估计; (b)历史无震区的初始应力估计.其中，σs：断层破裂强度；σr：地震释放应力后的状态；σ0：初始应力； ：应力加载速率；σd：地震释放的应力即应力降；σu：历史无震区的应力上限；σl：历史无震区的应力下限.
Fig. 4
Sketch map of initial tectonic stress estimation
(a) Stress estimation in the rupture zone of historical earthquakes; (b) Stress estimation in the region without earthquake.σs: strength of fault; σr: the state after the earthquake releases stress; σ0: initial stress; : stress increase rate; σd: stress drop; σu: upper limit of stress in the region without earthquake; σl: lower limit of stress in the region without earthquake.
|
|
(3) |
其中，σij为应力分量，εij为应变分量，E为杨氏模量，ν为泊松比(材料参数见表 1).
图 4b所示为历史上没有发生地震的较稳定地区：该区域的应力水平未曾达到破裂临界值，可估计其上限值；同时参考该区域已有观测数据及其他相关资料即可限定于某一个应力下限值.即历史无震区初始应力场仅能大致限定其上下限范围.
2.3 巴颜喀拉块体及邻区初始应力场估算选取巴颜喀拉块体及邻区1997年至2014年的7个MW≥6.5破坏性地震：1997年11月8日玛尼MW7.5大震、2001年11月14日昆仑山MW7.8大地震、2008年3月20日于田MW6.9地震、2008年5月12日汶川MW7.8大地震、2010年4月13日玉树MW6.9地震、2013年4月20日芦山MW6.7地震和2014年2月12日于田MW6.8地震.具体的地震震源参数等信息，见表 2.
|
|
表 2 地震序列目录 Table 2 Earthquakes catalog |
研究区域震源深度范围约为0~20 km，暂且选择10 km深度作为计算深度面，该深度处的静岩压力(Pv=ρrock×g×h，其中ρrock=2700 kg·m-3)约为270 MPa(按照岩石力学的定义，压应力为正值，拉应力为负值).若将二维平面模型应力值叠加上10 km深度处静岩压力，可视为准三维应力.假设计算得出的二维水平主压应力分别为最大值P10，最小值P20，那么准三维的应力值为：P1=P10+Pv、P2=P20+Pv、Pv.根据Anderson(1951)断层理论，由P1，P2，Pv三个主应力状态即可判断断层类型，即：P1>P2>Pv为逆冲型；P1>Pv>P2为走滑型；Pv>P1>P2为拉张型.若假设最大主压应力为Pmax，最小为Pmin.以库仑-摩尔破裂准则为约束，对于断层面上的正应力σn和剪应力τ满足条件τ=μσn，则可转换为最大主应力和最小主应力之间的关系：
|
(4) |
当实际主应力比值r≥r0时，岩石发生破裂，且不同类型的地震对应不同的主应力比值.由此可推测历次地震震前应力状态满足条件：(1)历次地震破裂区，震前临界主应力比值r≥r0；(2)历史无震区震前应力虽无法唯一限定，但可以确定主应力比值r上限小于r0，下限可根据静水压大小、构造活动性质及浅部地应力观测资料等大致推测.对于不同地区给出不同的主应力比值下限值.最终取值可在上下限范围内(如图 4b中的σu~σl之间)给出.
考虑现今构造应力场必然是由多重因素(如重力、板块边界力、地壳介质的非均匀力等)共同作用形成的，可假设线弹性有限元模型中每个单元的应力通过一定时间(t)的积累(以此来代表构造应力场形成过程中的多种因素作用)，达到震前不同区域的不同状态(破裂区的临界或半临界状态，历史无震区的稳定状态).在给定一个粗糙的现今构造应力场的基础上，往前倒推，减去震间应力积累和区域地震应力扰动作用，即可得到某一时刻的初始构造应力场.具体计算如下所示：
(1) 二维平面模型中每个单元计算得到的应力增长速率为
|
(5) |
根据主应力计算公式(6)，求得最大/小水平主应力P1, 2，
|
(6) |
根据线弹性计算特征，可以直接根据应力增长速率
|
(7) |
(2) 依据实际断层特征及区域构造特征等信息，确定模型中不同单元的滑动类型来反推积累时间t，计算如下：
1) 逆冲型：
|
(8) |
2) 走滑型：
|
(9) |
3) 拉张型：
|
(10) |
根据不同区域单元的不同主应力比值r(震源破裂区单元r≥r0，稳定地区r < r0)、主应力增长速率
(3) 计算过程中，首先可获得地震序列中最近一个地震(2014年于田地震)的震前应力场(假设于田地震震源区的应力处于临界破裂状态)；然后倒推计算，该应力场减掉与前一个地震(2013年芦山地震)的震间应力积累量，以及同震应力变化量(利用有限元“杀伤单元”法计算所得)，并利用步骤(2)计算芦山地震震源区应力状态，可得到芦山地震震前应力场(芦山地震震源区应力处于临界破裂状态).依次倒推推算到1997年玛尼地震前的初始应力场.
须强调，这个估计出的初始应力场只是根据已有的历史地震数据推测出的一种可能值，并非是一定真实的状态，只满足在历史地震区域应力超出其强度，其他地区相对的应力较低的假设.
3 计算结果及讨论 3.1 年应力增长速率图 5给出了研究区水平年应力增长速率，可以看出，约以96°E为分界，东西向正应力Sxx年增量呈现出西部为拉张、东部为压缩的特征；南北向应力Syy年增量以94°E分界，呈现出西部为南北向压缩、东部为拉张的特征；剪应力Sxy年增量在甘孜—玉树断裂和阿尔金断裂交汇处、昆仑山断裂带中部和东部龙门山与鲜水河断裂交汇处有较高值；水平主应力年增量则呈现出西部以东西向拉张为主、南北向压缩为辅，东部则反之，东西向压缩、南北向拉伸的特征，且主应力年增量在84°E—92°E之间有较高值，在东缘龙门山断裂带附近有最小值.最大主压应力增量方向自西向东由近NS逐渐转变为NNE、NE及近EW向，与块体运动方向保持一致，呈现顺时针旋转趋势.块体东部计算结果同Liu等(2016)计算得出的上地壳应力增长速率(最大值约3 kPa·a-1)量级相当.
|
图 5 计算构造应力年增长量 (a) Sxx年增量；(b) Syy年增量；(c) Sxy年增量；(d)水平主应力年增量.其中x为东西向，y为南北向，张应力为正值；红色箭头表示拉应力，黑色表示压应力；其中块体和断层名称同图 1. Fig. 5 Calculated stress accumulation rate of each year (a) Sxx rate; (b) Syy rate; (c) Sxy rate; (d) Horizontal principal stress rate. Red arrow show the tensile stress; Black ones show the compressive stress. The name of the blocks and faults are the same as in Fig. 1. |
在应力场反演过程中对于模型中历史无震稳定区部分的应力场估计仅能确定上下限范围，得到的初始应力场只是一种可能的模型.因此，采用Monte Carlo方法进行随机试验，即通过对随机变量的统计试验来获取近似解的数值方法(朱守彪和石耀霖，2007).利用随机数发生器，在初始应力场不确定部分的上下限范围内随机取值，通过多次大量独立试验(受计算限制，本文重复10000次计算)，得到10000组可能的初始应力场.图 6展示研究区域某一点(87.93°E，35.57°N，见图 7中的粉色五角星点)处10000种随机的水平主应力分布的统计直方图，呈现正态分布.在该点位处，最大、最小水平主应力平均值分别为：422 MPa，273 MPa.
|
图 6 某一点位处的水平主应力分布直方图 (a)最大水平主应力; (b)最小水平主应力. Fig. 6 Histograms of the computed horizontal principal stress at one point (a) Maximum principal stress; (b) Minimum principal stress. |
|
图 7 计算地应力与实测值对比 红色圆圈数字代表实际地应力测量的点，红色的箭头表示实际测量浅部的主应力方向；其中块体和断层名称同图 1. Fig. 7 Comparison between the simulation stress value and the observed value The numbers in the red circle represent the actual stress measurement point, and the red arrows indicate the actual measured principal stress direction; The name of the blocks and faults are the same as in Fig. 1. |
本研究中，我们通过对所有计算节点进行10000组统计分析，得到空间上各点的主应力分量的平均值作为初始应力场值，如图 7所示.可以看出，巴颜喀拉块体内最大水平主应力方向自西向东由NWW向逐渐转变为近NS向、NEE向，在块体东缘呈现出顺时针旋转趋势.其最大水平主应力值和最小水平主应力值分别为~400 MPa和~200 MPa.这与中国大陆地壳应力环境基础数据库中收录的数据基本吻合(谢富仁等, 2004；孟文等, 2017).进一步说明，由于印度板块对欧亚板块的碰撞俯冲及板块间相互作用，巴颜喀拉块体与南北相邻块体存在东向“逃逸”速度差，且其东缘与四川盆地存在着挤压作用.
另外，将计算结果进一步同青藏高原及龙门山地区开展的一些点位的地应力测量结果(吴满路等, 2005; 陈群策等, 2012; 杜建军等, 2013; Meng et al., 2015)相对比，如表 3及图 7所示.可以看出，模拟计算最大主压应力方向与实测浅部主应力方向大致相同，在1号点位和4号点位相差仅在5°左右；2号和3号点位略大，分别为24°和14°，原因在于：一方面浅部实测应力方向受岩体结构影响，测得主应力方向在浅部不同深度本身存在差异；另一方面，浅部主应力方向不一定与深部一致，且模拟计算深部应力存在误差.
|
|
表 3 实测地应力与计算值对比 Table 3 Comparison between the observed stress value and the simulation value |
孟文等(2017)给出拉萨地块浅部450 m范围内最大水平主压应力(PH)，最小水平主压应力(Ph)的梯度值约为54 MPa/km，33.2 MPa/km.假设巴颜喀拉块体深部10 km处主应力梯度值与拉萨地块浅部类同，则最大与最小水平主压应力量级分别约500 MPa和300 MPa.德国深钻测得陆壳9 km深度处差应力(P1- P3)可达230 MPa(Zoback and Harjes, 1997).本文计算1号点位于2001年可可西里大地震附近，是典型走滑型应力场结构(PH>Pv>Ph)，模拟计算PH=424 MPa，Ph=237 MPa，PH及Ph量级与前述推测量级相当；若给定Pv=270 MPa，则计算应力场与走滑型结构吻合；且计算差应力为PH-Ph=187 MPa，量级亦与德国深钻实测值相当.巴颜喀拉块体区域内整体以剪切形变为主，发震构造多为走滑型，最大主应力和最小主应力均在水平面上.
图 8a—8c分别为最大水平主应力、最小水平主应力和二者差值分布图，结果显示，块体自西向东，三者的值均呈现逐渐增加的趋势.西部阿尔金断裂与昆仑断裂交汇处为低值区，最大水平主应力约380 MPa，最小水平主应力约220 MPa，差应力约150 MPa；东部昆仑山断裂带东段和甘孜—玉树断裂带东南段部分区段为高值区，最大水平主应力约450 MPa，最小水平主应力约250 MPa，差应力约200 MPa.根据库仑-摩尔破裂准则，主应力差值的大小可以反映出区域的稳定性，差值越大越易发生破裂，该结果表明在昆仑山断裂带与阿尔金断裂带交汇处以及甘孜—玉树断裂带中西段地震活动性较低，而昆仑山断裂带东段及甘孜—玉树断裂带东南段部分区段地震活动性较高.野外地质调查、遥感影像解译及探槽研究等认为甘孜—玉树断裂带东南段全新世以来活动较强，有强震可能性(石峰等, 2013; 吴中海等, 2014).地震地质资料分析、古地磁资料研究以及库仑应力理论计算研究结果认为块体北边界的东昆仑断裂东段仍有较高强震可能性(李正芳等, 2012; 闻学泽等, 2011; 万永革等, 2007)，本文研究与前人研究结果一致.
|
图 8 (a) 最大水平主应力值，(b)最小水平主应力值，(c)最大与最小水平主应力差值 其中块体和断层名称同图 1. Fig. 8 (a) Maximum horizontal principal stress; (b) Minimum horizontal principal stress; (c) Principal stress differences The name of the blocks and faults are the same as in Fig. 1. |
另外，巴颜喀拉块体整体东向—东南向运移过程中，在块体运动前缘(东部龙门山地区)伴随有地壳缩短的逆冲推覆构造，块体运动后缘(块体西部于田地区)伴随有局部伸展拉张构造.发生在龙门山断裂带上的2008年汶川地震和2013年芦山地震均为逆冲兼具走滑型地震(陈运泰等，2013)，在于田地区2008年、2014年相距约100 km的两次地震分别为拉张型和走滑型地震(程惠红等，2014).其发震构造反映出区域构造运动特征，揭示区域构造应力场的复杂性.如前文所示，逆冲型地震和拉张型地震的最大和最小主应力分别为水平向和垂直向.但是由于受二维平面模型的限制，本文尚且仅能计算水平向应力场，对于垂向应力仅给出估计均值，无法给出局部逆冲型或拉张型地震应力场状况.因此，对于图 8c中块体东缘龙门山地区以及西缘于田地区的水平向差应力不能直接反映这两个区域的稳定性.
但本文计算中仍存在问题值得深入研究：(1)初始应力场的估算较为粗略.主要根据库仑-摩尔破裂准则，以及历史地震释放的应力降，来进行反演估算.对于历史地震的选取，采用了近期7个较典型的大震，仍有部分历史地震没有参与计算，但不能排除这些地震的影响.(2)对于历史地震应力降的计算较为简单，仅计算同震应力降.实际发生地震后，岩石的弹性瞬态响应会产生同震应力变化，同时，地壳岩层的黏弹性响应会在长时间内产生应力松弛效应，需在进一步研究中考虑该效应.(3)本文采用二维平面弹性模型和静水压力准三维应力场进行分析，结果虽大致可行，但巴颜喀拉块体的物质结构的横向和纵向不均匀性柔性下地壳的黏弹性在应力场演化中的动力学作用，以及对于巴颜喀拉块体挤压的东边界(龙门山断裂)及局部拉伸的西边界于田地区尚无法真实反映.
4 结论震前初始构造应力场是决定孕震环境的重要因素，是数值计算分析应力场演化及地震序列演化的重要初始条件，但是目前暂且不能通过直接观测手段获取深部区域构造应力场数据.本文基于已有的历史地震资料和观测数据，初步反演计算得到震前初始构造应力场.以近二十年来巴颜喀拉块体及其边界断裂带发生的7次中强震为例，基于库仑-摩尔破裂准则和历史地震信息、地质背景及GPS观测数据，建立有限元数值模型.并提出利用Monte Carlo方法，对初始应力场反演中的不确定部分限定在一个合理的上下限范围内进行独立的重复性随机试验，运用统计学方法给出其均值，以此得到巴颜喀拉块体于1997年玛尼MW7.5地震前区域初始应力场.其结果显示，巴颜喀拉块体区域的最大水平主压应力方向，自西向东逐渐由近NS向转变为近EW向，呈现出顺时针旋转趋势.最大与最小水平主应力量级分别为~400 MPa和~200 MPa，其中最大水平主应力、最小水平主应力以及二者应力差值均在块体西部有较低值，分别为：380 MPa，220 MPa，150 MPa；块体东部有较高值，分别为：450 MPa，250 MPa，200 MPa.在昆仑山断裂带与阿尔金断裂带交汇处以及甘孜—玉树断裂带中西段区域差应力值较低，地震活动性较低；在昆仑山断裂带东段及甘孜—玉树断裂带东南段有较高差应力值，地震活动性较高.该结果可与浅部地应力观测值以及历史地震活动类型有较好一致性.
初始应力场归根结底是需要测出来的，但是以目前的科技水平，尚不能短期内实现.本文在已有的观测条件基础上，提出利用数值模拟方法来推测初始构造应力场，以期对构造应力场与地震活动性相关研究提供参考.
致谢 感谢两位匿名审稿人对本文提出的宝贵指导意见！
Anderson E M. 1951. The Dynamics of Faulting and Dyke Formation with Applications to Britain. Edinburgh, Oliver and Boyd, U.K.: Hafner Pub. Co..
|
Cao J L, Shi Y L, Zhang H, et al. 2009. Numerical simulation of GPS observed clockwise rotation around the eastern Himalayan syntax in the Tibetan Plateau. Chinese Science Bulletin, 54(8): 1398-1410. DOI:10.1007/s11434-008-0588-7 |
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 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 |
Cheng H H, Pang Y J, Dong P Y, et al. 2014. Analysis of the stress environment of the 2008 and 2014 Yutian MS7.3 earthquakes. Chinese Journal of Geophysics (in Chinese), 57(10): 3238-3246. DOI:10.6038/cjg20141012 |
Cui X F, Xie F R, Zhao J T. 2005. The regional characteristics of focal mechanism solutions in China and its adjacent areas. Seismology and Geology (in Chinese), 27(2): 298-307. |
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, 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 J J, Chen Q C, Ma Y S, et al. 2013. Faults activity and stress state in the northeast segment of Longmenshan faults zone. Progress in Geophysics (in Chinese), 28(3): 1161-1170. DOI:10.6038/pg20130307 |
Funning G J, Parsons B, Wright T J. 2007. Fault slip in the 1997 Manyi, Tibet earthquake from linear elastic modelling of InSAR displacements. Geophysical Journal International, 169(3): 988-1008. DOI:10.1111/j.1365-246X.2006.03318.x |
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 |
Goodman R E. 1989. Introduction to Rock Mechanics. 2nd ed. New York: Wiley.
|
Guo Q L, Wang C H, Ma H S, et al. 2009. In-situ hydro-fracture stress measurement before and after the Wenchuan MS8.0 earthquake of China. Chinese Journal of Geophysics (in Chinese), 52(5): 1395-1401. DOI:10.3969/j.issn.0001-5733.2009.05.029 |
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 |
Li Z F, Zhou B G, Ran H L. 2012. Strong earthquake risk assessment of eastern segment on the East Kunlun fault in the next 100 years based on paleo-earthquake data. Chinese Journal of Geophysics (in Chinese), 55(9): 3051-3065. DOI:10.6038/j.issn.0001-5733.2012.09.023 |
Liao C T, Zhang C S, Wu M L, et al. 2003. Stress change near the Kunlun fault before and after the MS8.1 Kunlun earthquake. Geophysical Research Letters, 30(20): 2027. DOI:10.1029/2003gl018106 |
Liu C, Shi Y L, Zhu B J, et al. 2014. Crustal rheology control on the mechanism of the earthquake generation at the Longmen Shan fault. Chinese Journal of Geophysics (in Chinese), 57(2): 404-418. DOI:10.6038/cjg20140207 |
Liu C, Zhu B J, Yang X L, et al. 2016. Geodynamic background of the 2008 Wenchuan earthquake based on 3D visco-elastic numerical modelling. Physics of the Earth and Planetary Interiors, 252: 23-36. DOI:10.1016/j.pepi.2016.01.003 |
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 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 |
Meng W, Guo C B, Zhang C Y, et al. 2017. In situ stress measurements and implications in the Lhasa Terrane, Tibetan Plateau. Chinese Journal of Geophysics (in Chinese), 60(6): 2159-2171. DOI:10.6038/cjg20170611 |
Meng W, Chen Q C, Zhao Z, et al. 2015. Characteristics and implications of the stress state in the Longmen Shan fault zone, eastern margin of the Tibetan Plateau. Tectonophysics, 656: 1-19. DOI:10.1016/j.tecto.2015.04.010 |
Pang Y J, Cheng H H, Zhang H, et al. 2017. Numerical modeling of crustal deformation in the eastern margin of the Bayan Har block and analysis of seismogenic environment of the 2017 Jiuzhaigou earthquake. Chinese Journal of Geophysics (in Chinese), 60(10): 4046-4055. DOI:10.6038/cjg20171030 |
Pang Y J, Zhang H, Gerya T V, et al. 2018. The Mechanism and dynamics of N-S rifting in southern Tibet:insight from 3-D thermomechanical modeling. Journal of Geophysical Research:Solid Earth, 123(1): 859-877. DOI:10.1002/2017JB014011 |
Prieto G A, Beroza G C, Barrett S A, et al. 2012. Earthquake nests as natural laboratories for the study of intermediate-depth earthquake mechanics. Tectonophysics, 570-571: 42-56. DOI:10.1016/j.tecto.2012.07.019 |
Royden L H, Burchfiel B C, King R W, et al. 1997. Surface deformation and lower crustal flow in eastern Tibet. Science, 276(5313): 788-790. DOI:10.1126/science.276.5313.788 |
Scholz C H. 1998. Earthquakes and friction laws. Nature, 391(391): 37-42. DOI:10.1038/34097 |
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. Physics (in Chinese), 42(4): 237-255. DOI:10.7693/wl20130402 |
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 |
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. DOI:10.3321/j.issn:0253-3782.2007.02.001 |
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 H, Zhang G M, Shi Y L, et al. 2006. Numerical simulation of movement and deformation of Qinghai-Tibet Plateau. Journal of Geodesy and Geodynamics (in Chinese), 26(2): 15-23. DOI:10.3969/j.issn.1671-5942.2006.02.004 |
Wang Q, Qiao X J, Lan Q G, et al. 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan. Nature Geoscience, 4(9): 634-640. DOI:10.1038/ngeo1210 |
Wen X Z, Du F, Zhang P Z, et al. 2011. Correlation of major earthquake sequences on the northern and eastern boundaries of the Bayan Har block, and its relation to the 2008 Wenchuan earthquake. Chinese Journal of Geophysics (in Chinese), 54(3): 706-716. DOI:10.3969/j.issn.0001-5733.2011.03.010 |
Wu M L, Zhang C S, Liao C T, et al. 2005. The recent state of stress in the central Qinghai-tibet plateau according to in-situ stress measurements. Chinese Journal of Geophysics (in Chinese), 48(2): 327-332. |
Wu Z H, Zhou C J, 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. DOI:10.3969/j.issn.1671-2552.2014.04.003 |
Xie F R, Cui X F, Zhao J T, et al. 2004. Regional division of the recent tectonic stress field in China and adjacent areas. Chinese Journal of Geophysics (in Chinese), 47(4): 654-662. |
Xu L S, Chen Y T. 2004. The inversion of temporal and spatial rupture process of the 14 Nov. 2001 Kunlun earthquake based on the global long-period waveform data. Science in China Series D:Earth Sciences (in Chinese), 34(3): 256-264. |
Yang S X, Lu Y Z, Chen L W, et al. 2012a. 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 |
Yang S X, Yao R, Cui X F, et al. 2012b. Analysis of the characteristics of measured stress in Chinese mainland and its active blocks and North-South seismic belt. Chinese Journal of Geophysics (in Chinese), 55(12): 4207-4217. DOI:10.6038/j.issn.0001-5733.2012.12.032 |
Zhang G H, Qu C Y, Shan X J, et al. 2011. The coseismic InSAR measurements of 2008 Yutian earthquake and its inversion for source parameters. Chinese Journal of Geophysics (in Chinese), 54(11): 2753-2760. DOI:10.3969/j.issn.0001-5733.2011.11.005 |
Zhang Y, Xu L, Chen Y T. 2010. Fast inversion of rupture process for 14 April 2010 Yushu, Qinghai, earthquake. Acta Seismologica Sinica (in Chinese), 32(3): 361-365. DOI:10.3969/j.issn.0253-3782.2010.03.011 |
Zhou Y, Wang W M, Xiong L, et al. 2015. Rupture process of 12 February 2014, Yutian MW6.9 earthquake and stress change on nearby faults. Chinese Journal of Geophysics (in Chinese), 58(1): 184-193. DOI:10.6038/cjg20150116 |
Zhu S B, Shi Y L. 2005. Genetic algorithm-finite element inversion of topographic spreading forces and drag forces of lower crust to upper crust in Tibetan plateau. Acta Scicentiarum Naturalum Universitis Pekinesis (in Chinese), 41(2): 225-234. DOI:10.3321/j.issn:0479-8023.2005.02.008 |
Zhu S B, Shi Y L. 2007. Error analysis of strain rates from GPS measurements based on Monte Carlo method. Chinese Journal of Geophysics (in Chinese), 50(3): 806-811. |
Zhu S B, Yuan J, Miao M. 2017. Dynamic mechanisms for supershear rupture processes of the Yushu earthquake (MS=7.1). Chinese Journal of Geophysics (in Chinese), 60(10): 3832-3843. DOI:10.6038/cjg20171013 |
Zoback M D, Harjes H P. 1997. Injection-induced earthquakes and crustal stress at 9 km depth at the KTB deep drilling site, Germany. Journal of Geophysical Research:Solid Earth, 102(B8): 18477-18491. DOI:10.1029/96jb02814 |
曹建玲, 石耀霖, 张怀, 等. 2009. 青藏高原GPS位移绕喜马拉雅东构造结顺时针旋转成因的数值模拟. 科学通报, 54(2): 224-234. |
陈群策, 丰成君, 孟文, 等. 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 |
程惠红, 庞亚瑾, 董培育, 等. 2014. 于田2008年和2014年两次MS7.3地震孕育的应力环境. 地球物理学报, 57(10): 3238-3246. DOI:10.6038/cjg20141012 |
崔效峰, 谢富仁, 赵建涛. 2005. 中国及邻区震源机制解的分区特征. 地震地质, 27(2): 298-307. DOI:10.3969/j.issn.0253-4967.2005.02.012 |
邓起东, 程绍平, 马冀, 等. 2014. 青藏高原地震活动特征及当前地震活动形势. 地球物理学报, 57(7): 2025-2042. DOI:10.6038/cjg20140701 |
董培育, 石耀霖. 2013. 关于"用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移"的讨论——横向各向同性"杀伤单元"才是更好的途径. 地球物理学报, 56(6): 2133-2139. DOI:10.6038/cjg20130633 |
董培育, 胡才博, 石耀霖. 2016. 青藏高原及周边区域地表长期变形数值模拟. 地震地质, 38(2): 410-422. DOI:10.3969/j.issn.0253-4967.2016.02.014 |
杜建军, 陈群策, 马寅生, 等. 2013. 龙门山断裂带东北段地应力状态与断裂活动性研究. 地球物理学进展, 28(3): 1161-1170. DOI:10.6038/pg20130307 |
郭啟良, 王成虎, 马洪生, 等. 2009. 汶川MS8.0级大震前后的水压致裂原地应力测量. 地球物理学报, 52(5): 1395-1401. DOI:10.3969/j.issn.0001-5733.2009.05.029 |
胡才博, 周一杰, 蔡永恩. 2009. 如何用有限元新模型研究地震触发和应力场连续演化. 中国科学D辑:地球科学, 39(5): 546-555. |
李正芳, 周本刚, 冉洪流. 2012. 运用古地震数据评价东昆仑断裂带东段未来百年的强震危险性. 地球物理学报, 55(9): 3051-3065. DOI:10.6038/j.issn.0001-5733.2012.09.023 |
柳畅, 石耀霖, 朱伯靖, 等. 2014. 地壳流变结构控制作用下的龙门山断裂带地震发生机理. 地球物理学报, 57(2): 404-418. DOI:10.6038/cjg20140207 |
刘启元, 吴建春. 2003. 论地震数值预报——关于我国地震预报研究发展战略的思考. 地学前缘, 10(S1): 217-224. DOI:10.3321/j.issn:1005-2321.2003.z1.030 |
孟文, 郭长宝, 张重远, 等. 2017. 青藏高原拉萨块体地应力测量及其意义. 地球物理学报, 60(6): 2159-2171. DOI:10.6038/cjg20170611 |
庞亚瑾, 程惠红, 张怀, 等. 2017. 巴颜喀拉块体东缘形变及九寨沟地震孕震环境数值分析. 地球物理学报, 60(10): 4046-4055. DOI:10.6038/cjg20171030 |
石峰, 李安, 杨晓平, 等. 2013. 甘孜-玉树断裂带东南段晚第四纪活动性研究. 地震地质, 35(1): 50-63. DOI:10.3969/j.issn.0253-4967.2013.01.004 |
石耀霖, 孙云强, 罗纲, 等. 2018. 关于我国地震数值预报路线图的设想——汶川地震十周年反思. 科学通报, 63(19): 1865-1881. DOI:10.1360/N972018-00335 |
石耀霖, 张贝, 张斯奇, 等. 2013. 地震数值预报. 物理, 42(4): 237-255. DOI:10.7693/wl20130402 |
万永革, 沈正康, 曾跃华, 等. 2007. 青藏高原东北部的库仑应力积累演化对大地震发生的影响. 地震学报, 29(2): 115-129. DOI:10.3321/j.issn:0253-3782.2007.02.001 |
王成虎, 宋成科, 郭启良, 等. 2014. 利用原地应力实测资料分析芦山地震震前浅部地壳应力积累. 地球物理学报, 57(1): 102-114. DOI:10.6038/cjg20140110 |
王辉, 张国民, 石耀霖, 等. 2006. 青藏活动地块区运动与变形特征的数值模拟. 大地测量与地球动力学, 26(2): 15-23. DOI:10.3969/j.issn.1671-5942.2006.02.004 |
闻学泽, 杜方, 张培震. 2011. 巴颜喀拉块体北和东边界大地震序列的关联性与2008年汶川地震. 地球物理学报, 54(3): 706-716. DOI:10.3969/j.issn.0001-5733.2011.03.010 |
吴满路, 张春山, 廖椿庭, 等. 2005. 青藏高原腹地现今地应力测量与应力状态研究. 地球物理学报, 48(2): 327-332. DOI:10.3321/j.issn:0001-5733.2005.02.014 |
吴中海, 周春景, 冯卉, 等. 2014. 青海玉树地区活动断裂与地震. 地质通报, 33(4): 419-469. DOI:10.3969/j.issn.1671-2552.2014.04.003 |
谢富仁, 崔效锋, 赵建涛, 等. 2004. 中国大陆及邻区现代构造应力场分区. 地球物理学报, 47(4): 654-662. DOI:10.3321/j.issn:0001-5733.2004.04.016 |
许力生, 陈运泰. 2004. 从全球长周期波形资料反演2001年11月14日昆仑山口地震时空破裂过程. 中国科学D辑:地球科学, 34(3): 256-264. |
杨树新, 陆远忠, 陈连旺, 等. 2012a. 用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移. 地球物理学报, 55(1): 105-116. DOI:10.6038/j.issn.0001-5733.2012.01.010 |
杨树新, 姚瑞, 崔效锋, 等. 2012b. 中国大陆与各活动地块、南北地震带实测应力特征分析. 地球物理学报, 55(12): 4207-4217. DOI:10.6038/j.issn.0001-5733.2012.12.032 |
张国宏, 屈春燕, 单新建, 等. 2011. 2008年MS7.1于田地震InSAR同震形变场及其震源滑动反演. 地球物理学报, 54(11): 2753-2760. DOI:10.3969/j.issn.0001-5733.2011.11.005 |
张勇, 许力生, 陈运泰. 2010. 2010年4月14日青海玉树地震破裂过程快速反演. 地震学报, 32(3): 361-365. DOI:10.3969/j.issn.0253-3782.2010.03.011 |
周云, 王卫民, 熊林, 等. 2015. 2014年2月12日MW6.9于田地震震源破裂过程及对周围断层的应力影响. 地球物理学报, 58(1): 184-193. DOI:10.6038/cjg20150116 |
朱守彪, 石耀霖. 2005. 青藏高原地形扩展力以及下地壳对上地壳的拖曳力的遗传有限单元法反演. 北京大学学报(自然科学版), 41(2): 225-234. DOI:10.3321/j.issn:0479-8023.2005.02.008 |
朱守彪, 石耀霖. 2007. 基于Monte Carlo方法的由GPS观测计算地应变率的误差分析. 地球物理学报, 50(3): 806-811. DOI:10.3321/j.issn:0001-5733.2007.03.020 |
朱守彪, 袁杰, 缪淼. 2017. 青海玉树地震(MS=7.1)产生超剪切破裂过程的动力学机制研究. 地球物理学报, 60(10): 3832-3843. DOI:10.6038/cjg20171013 |