汶川地区区域构造应力场特征及汶川MS8.0地震对周围主要断层面的影响＊

何 金 许 鑫 吴 彪 刘晓磊 杨东辉 路 彤

1) 河北省地震局易县地震台，河北保定 074200

2) 防灾科技学院，河北三河 065201

3) 河北省地震动力学重点实验室，河北三河 065201

4) 中国地质大学（北京）地球物理与信息技术学院，北京 100083

5) 河北省地震局承德地震监测中心站，河北承德 067000

6) 河北省地震局唐山地震监测中心站，河北唐山 063000

引言

2008 年5 月12 日14 时28 分，四川省汶川县（31.00°N，103.47°E）发生了MS8.0 地震。汶川地震是新中国成立以来，中国大陆破坏性最强、影响范围最广的地震，它给我国人民的生命财产造成了巨大损失。

汶川地震发生于青藏高原内部的巴颜喀拉块体东南缘。青藏高原是我国现代构造活动与地震活动最强烈的地区[1]，参考邓起东等[2]根据时空分布提出的地震系列，本次地震与同样发生在巴颜喀拉块体的2010 年玉树MS7.1 地震和2017 年九寨沟MS7.0 地震为同一系列地震。所以计算出本次地震的同震库仑应力变化对邻域地震事件及周围主要断层的影响，对了解地震事件间的触发与否等关联性具有一定的参考意义。大量的震例研究表明，地震引起的静态库仑应力变化，会对周围地区后续地震的发生起到提前或延迟的作用[3-6]。McCloskey 等[3]通过2004年苏门答腊MW9.0 地震对周围断层的静态应力触发研究后发现，Nias-Simeulue 段是未来可能发生地震的危险区域，2005 年发生在Nias 的MW8.7 地震证实了他们的预测。万永革等[4]分别给出了2008 年汶川MS8.0 地震和于田MS7.3 地震在周围断层的不同断层段上产生的库仑应力变化，并指出部分断层上的大震的发震被提前或滞后，2013 年4 月20 日发生在汶川地震南部的芦山M7.0 地震和2014 年2 月12 日发生在于田地震震中东北部的MS7.3 地震就验证了他们的研究成果。万永革等[5]计算了2015 年尼泊尔强震序列在中国大陆地区造成的静态库仑应力变化，结果表明该次强震序列产生的应力加载区主要集中在青藏断块西南地区的逆冲断裂和大部分拉张断裂。靳志同等[6]研究了2017 年九寨沟MS7.0 地震对周围断层的应力影响，发现该次地震的应力加载区集中在虎牙断裂中段和南段、岷江断裂北段南部和塔藏断裂西段。

由于应力场的状态与强震的发震密切相关，并且震后很多学者或机构利用不同方法给出此次地震序列的震源机制解，故而研究区域构造应力场成为汶川地震的主要研究方向之一。程方正等[7]利用2.5＜M＜6.0 地震震源机制解研究了震前区域应力场动态。陈学忠等[8]基于丰富的数字地震波记录，对区域构造应力场及其变化特征进行研究，发现汶川地震前后的应力水平也有明显差异。张勇等[9]将利用远场地震波形资料获取大震震源机制的方法应用于汶川地震，发现地震断层的南端为逆冲，随着破裂向东北方向延伸，断层走滑分量逐渐加大，从以逆冲为主到以走滑为主的转折点在震中东北大约190 km 的位置。杜义等[10]通过对龙门山断裂带震后断层擦痕测量，得到311 条断层擦痕数据，利用断层资料反演构造应力张量的计算方法，得到研究区8 个测点的构造应力张量数据，并获得研究区构造应力场特征。郑勇等[11]利用“裁剪—粘贴”（cut and paste，即CAP）法，对MS4.0 以上余震震源机制解进行反演给出该地区现今应力场。基于前人给出的震源机制解结果，进一步对2008 年汶川地震的区域构造应力场进行分段研究。

研究区域大震前后的应力变化对评估地震危险性是非常重要的，但由于其受限于区域构造背景、源场模型和接受断层面参数等因素。通常，利用应力传输理论研究应力场变化主要从大震影响本地断层应力入手，从而会忽略地震带来的构造应力变化，并且单一的库仑应力研究也只能提供某个特定断层上的应力变化，所以将构造应力场及库仑应力变化结合起来研究，能为地震危险性的评估提供更全面的参考。

本文搜集了相对全面的汶川地震序列的震源机制解，使用中心震源机制解剔除重复地震对应力场结果的影响，并使用网格搜索法反演该区域构造应力场；然后结合中心震源机制解法得到的主震的接受断层面参数，结合USGS 给出的地震破裂模型，计算了本次地震对周围断层的影响和水平面应力变化；最后综合区域构造应力场结果和静态库仑应力变化，计算出此次地震对背景构造应力场的影响。

1 汶川地震区域构造应力场

1.1 研究数据

自2008 年5 月12 日汶川大地震发生之后，越来越多的学者参与到汶川地震的研究中，由于各个组织和个人利用不同的方法，得出的数据结果也大不相同，为了更好地研究汶川地震，本文在知网以及地震局相关网站等查找出汶川地震序列的震源机制解。本文所用的震源机制来源于胡幸平等[12]、崔效锋等[13]、易桂喜等[14]、杨宜海等[15]、郭祥云等[16]、蒋海昆等[17]、郑勇等[11]、Tian 等[18]和中国地震台网网站发布的震源机制解，共搜集到2008 年5 月12 日—2013年4 月19 日，（27.39°—33.15°N，102.00°—106.19°E）范围内的1 660 条地震的震源机制数据。将这些数据按照时间顺序排序，根据发震时间、经度、纬度、震级区别出同一次地震，并将同一次地震的震源机制整理到一起，共找到230 个同一地震的多个震源机制。利用万永革[19]提出的利用同一地震的多个震源机制解确定其中心解的方法，求解汶川地震序列的中心震源机制解，用中心震源机制解替换同一地震的多个震源机制解。将自己收集到的震源机制解资料用GMT 绘图，根据Zoback[20]给出的对震源机制解的划分标准（表1），将这些数据划分成正断型、走滑型和逆断型，用于分析研究。

表1 震源机制解划分标准[20]Table 1 Division criteria of focal mechanism solutions[20]

通过震源机制的分析，我们可以发现汶川地震的余震沿着NE 方向贯通了整个龙门山断裂带的北段和中段。利用本文搜集到的震源机制解资料作Mt图（图1）发现，在汶川MS8.0 地震发生后的4 个月内，每月有超过20 个ML＞3.5 的余震，说明这个地区的构造应力场发生了比较大的改变。而在2008 年9 月之后的每个月则少于20 个ML＞3.5 的余震，这说明这个地区的构造应力场开始变得平稳了[21]。

图1 2008 年5 月12 日—2013 年4 月19 日汶川地震后ML＞3.0 地震的M-t 图Fig.1 M-t diagram with magnitude greater than 3.0 after Wenchuan earthquake from May 12，2008 to April 19，2013

1.2 研究方法

确定中心震源机制解主要分为两个步骤[19]：①求解两震源机制间的最小空间旋转角。引入两个震源机制a和b，以PBT轴为基：BPBT=Rab APBT，可得矩阵Rab以及Rab的主对角线元素之和tr(Rab)；同理以xyz轴为基：。两种基之间一定存在可逆线性变换关系，由相似矩阵定义知，1+2cosθab，可用最小值表示两者的差别，即最小空间旋转角(θab)min。②Levenberg-Marquardt 方法迭代求解：所求的中心解实际上就是与机构和学者给出的所有的震源机制的最小空间旋转角平方和最小的那个解。故将目标函数转化为线性问题，其中，(θci)min为中心震源机制与测量震源机制的最小空间旋转角；再使用Levenberg-Marquardt 方法迭代求解，为使其快速收敛加入了一个可变阻尼项；最后计算各震源机制解作为初始解与全部震源机制解的标准差其中，s采用类似于万永革表述应力场误差的方法[22]），取标准差S最小作为最终的结果，还可得到与其他震源机制解之间差别最小的空间旋转角。

用于应力场反演的方法很多。许忠淮等[23]在假定区域构造应力场一致的前提下，利用中小地震的综合节面解求得区域构造应力场；Gephart 等[24]利用网格搜索法通过区域内若干震源机制解约束应力张量；Michael[25]把断层剪切面上的应力分量归一化，变为线性方程求解应力场。相比前几种方法，Wan 等[26]采用全局网格搜索法求最优解，并引入了震源机制数据精度的权重因子，可以反演出更精确的应力场，所以本研究采用Wan 等[26]的全局网格搜索法。

计算震源机制解权重的公式为：

式中，W为权重；D为衰减系数；r为震级的相对大小。

应力场中描述3 个主轴应力间相对大小的应力形因子R值的计算公式为：

式中，σ1、σ2、σ3分别为最大、中间、最小主压应力。

1.3 汶川地区震源机制分区及应力场研究

本节将前文计算得到的汶川地区多个震源机制中心解数据（详见附表I）联同去除同一次地震的其余震源机制整理在一起，共计911 个震源机制解数据，结合地质构造背景[10,27]以及汶川地震同震破裂过程（图2）可将汶川地区分为3 区，分别是西南区、中区和东北区。

按照上述分区，第一个部分是东北区，共有231条数据；第二个部分是中区，共有150 条数据；第三部分是西南区，共有530 条数据。根据这些分区对所在区域内的震源机制解数据，使用Wan 等[26]提出的基于震源机制解的网格搜索法进行应力场反演。

从反演结果（表2）可看出：在置信度为90%时，东北区的最大主压应力轴走向为281.84°—282.90°，最大主张应力轴走向为157.79°—159.29°，主要受WNW-ESE 向的挤压；中区的最大主压应力轴走向为255.31°—264.31°，最大主张应力轴走向为-2.42°—6.58°，主要受WSW-ENE 向的挤压；西南区最大主压应力轴走向为267.50°—270.00°，最大主张应力轴走向357.50°—360.00°，主要受近W-E 向的挤压。由反演得到的应力场示意图（图3）可知，东北区、中区和西南区的R值分别为0.4、0.2 和0.5。中间应力轴在东北区和中区呈现出压应力轴特征，在西南区不表现出任何应力轴特征。最大主压应力轴倾伏角东北区为1.50°—3.00°，中区为0.99°—2.48°，西南区为-0.50°—0.80°，由此看出主压应力轴的倾伏角都比较小。主张应力轴倾伏角东北区为85.89°—86.89°，中区为79.78°—82.78°，西南区为52.50°—54.00°，由此可见，最大主张应力轴倾伏角在东北区和中区都比西南区大。

图3 分区的应力场及辐射花样图Fig.3 Stress field and radiation pattern of zoning

2 汶川地震的同震库仑应力变化对周围主要断层面的影响

研究区域大震后的应力变化对评估地震危险性是非常重要的，但由于其受限于区域构造背景、源场模型和接收断层面参数等因素。通常，利用应力传输理论研究应力场变化主要从大震影响本地断层应力入手，从而会忽略地震带来的构造应力变化，并且单一的库仑应力研究也只能提供某个特定断层上的应力变化，所以将构造应力场和库仑应力变化结合起来研究，能更全面地为评估地震危险性提供参考。

2.1 研究方法

地震的发生是由地下岩石的破裂造成的。Jaeger[28]认为，岩石的抗剪强度与作用于节面上的剪应力和正应力有关。地震专家和学者基于地震活动性观测和应力变化，提出一次地震的发生会对后续的地震产生影响（应力触发）[29-31]。在计算静态库仑应力时，主要采用基于弹性半空间的解析表达式[32]。如果已知发震断层面的几何模型和滑动量，则可求出弹性体产生的应力变化张量[29,33]，从而把应力变化张量投影到接收断层面的滑动方向和法向方向，求得库仑破裂应力变化：

式中，Δτ为沿接收断层面滑动方向的剪切应力变化（为正时，Δτ与滑动方向一致）；Δσ为沿接收断层面法向的正应力变化（拉张为正）；μ为视摩擦系数，由于孔隙流体也会引起应力变化，同时考虑到断层面的介质特性，一般μ在0.2—0.8 之间[34-36]。本研究延续前人的经验[37-39]，μ取0.4。当Δσf为正时，说明该次地震促进了研究区域断层的破裂，未来地震发生在研究区域的时间将会提前；当Δσf为负时，说明该次地震抑制了研究区域断层的破裂（也称应力影区），未来地震发生在研究区域的时间将滞后。同时研究表明，库仑应力变化大于0.01 MPa（静态应力触发阈值）时，能有效地对后续地震的时空分布造成影响[40-41]。

在计算库仑破裂应力时，需要知道接收断层的精确参数（走向、倾角和滑动角)。活动断层的走向可以根据其在地表的几何展布获取，结果较为精确；但是活动断层的倾角通常只给出一个大概的范围，精度无法满足计算的要求；对于活动断层的滑动角，在地质调查中只给出基础的滑动性质描述（正（逆）断层或左（右）旋走滑断层），导致无法获取一个精确的值。由于以上参数的精确度问题，在计算沿接收断层面滑动方向的剪应力变化和法向方向的正应力变化会产生误差，从而导致库仑破裂应力变化也存在误差。另外，在研究地震对周围地区的影响时，水平应力作用（应力变化张量水平分量）也是一个重要因素，这样可以直观、简捷地看出该地震在地图上的动力学作用特征。

在只考虑应力变化张量水平分量时，应力变化张量可改写成如下形式：

水平面应力则可以表示为：

式中，σHmax表示水平面最大主应力；σHmin表示水平面最小主应力；σHc＞0 表示处于膨胀区，反之，表示处于压缩区。

最大主应力与x轴的夹角：

2.2 汶川地震在震中附近主要断层面上产生的同震库仑应力变化

根据计算方法可知，无论是研究主震对周围断层的静态库仑破裂应力影响还是对周围地区产生的水平应力影响，首先需要确定汶川主震的破裂模型。本研究基于Okada[32]给出的弹性半空间的位错理论，采用USGS 给出的发震断层面几何模型和滑动量，接收断层面参数根据邓起东等[1]给出的中国大陆断层的几何形状和运动特征（表3）得到，以此来计算该次地震对周围断层造成的影响。

由图4 可见，汶川地震对2013 年芦山MS7.0 地震及2021 年玛多MS7.4 地震触发作用明显，对2010年玉树MS7.1 地震、2017 年九寨沟MS7.0 地震、2021 年漾濞MS6.4 地震及2021 年泸州MS6.0 地震均存在抑制作用。除此之外，本次地震造成龙门山断裂南北两端、秦岭南缘断裂、鲜水河断裂东南端、东昆仑断裂、白玉断裂的库仑破裂应力增加。其中龙门山断裂西南端增加最为明显，最大达258 700 Pa（约0.259 MPa），远比静态应力触发阈值（0.01 MPa）大，秦岭南缘断裂次之（0.047 MPa）（表3）。在白玉断裂、东昆仑断裂的库仑应力值均在触发阈值附近（图4）。

图4 2008 年汶川地震14 km 深度处在周围主要断裂上产生的库仑破裂应力变化Fig.4 Variation of Coulomb fracture stress on the surrounding main faults at the depth of 14 km of Wenchuan earthquake in 2008

表3 本研究所用的断层面性质和库仑破裂应力增加段Table 3 Fault plane properties and Coulomb fracture stress increasing section used in this study

2.3 汶川地震在震中附近产生的水平面应力变化

本节基于2.1 节给出的方法，计算出汶川地震产生的水平面应力变化对周围地区的影响（图5）。由图5 可见，此次地震是一次单侧破裂事件，地震发生时，震源处会释放大量的应变能，震中的面应力值下降，导致震中处于面挤压应力区，同样可知，此次破裂沿NE 向分为3 段。

图5 汶川地震产生的应力场在深度14 km 处的水平面投影Fig.5 Horizontal projection of stress field generated by Wenchuan earthquake at the depth of 14 km

从水平面最大主应力和最小主应力方向来看，震中附近水平面最小主应力多为近SN 向，最大主应力多为近EW 向，这与上文计算出的区域构造应力场方向基本相反，说明本次地震产生的应力场在一定程度上抵消了该区域构造应力场。

3 结论和讨论

本文首先搜集了2008 年5 月12 日—2013 年4 月19 日汶川地震及其邻区的1 660 条震源机制解，并采用同一地震多个震源机制解的中心解的方法筛去重复地震事件的震源机制解，之后通过网格搜索法分段反演出区域构造应力场。然后，基于USGS 给出的汶川MS8.0 地震的破裂模型，计算出该地震对附近强震的触发关系，以及汶川MS8.0 地震对周围断层的同震库仑应力变化。最后，结合构造应力场结果和同震库仑应力结果，计算本次地震对背景构造应力场产生的影响。

（1）东北区、中区和西南区最大主压应力轴走向分别为281.84°—282.90°，255.31°—264.31°和267.50°—270.00°。杜义等[10]通过擦痕资料反演应力场的最大主应力走向为256°—301°；盛书中等[42]分南北段反演应力场结果显示，最大主压应力轴走向北端为249°，南端为255°，与本文反演结果在误差允许范围内。对于倾伏角而言，东北区、中区和西南区3 个区域的P轴倾伏角都较小，T轴倾伏角自西南向东北倾伏角由52.50°增加到85.89°；张勇等[9]利用远场地震波资料反演的破裂模型显示，从西南至东北延伸，在震中约190 km 处，由逆冲为主转为走滑为主，这与T轴倾伏角的变化相吻合。东北区主要受WNWESE 向的挤压，西南区受W-E 向的挤压，中区受WSW-ENE 向的挤压。由此可以观察到，自西南到东北，主压应力轴方向发生了小的方向变化。汶川地区位于青藏高原东南缘，其内多发育断层，地质构造活动剧烈，受到来自印度板块北北东方向的俯冲推挤、四川盆地的阻挡和巴颜喀拉块体东南向挤压的联合作用。由应力场反演得到的结果与汶川地区的构造应力场背景有关，也与前人的研究结果相似。

（2）汶川地震对2013 年芦山MS7.0 地震[43]及2021 年玛多MS7.4 地震触发作用明显，对2010 年玉树MS7.1 地震、2017 年九寨沟MS7.0 地震、2021 年漾濞MS6.4 地震及2021 年泸州MS6.0 地震均存在抑制作用。汶川地震产生的库仑应力变化对周围断层的影响结果与Parsons 等[44]、王连捷等[45]、万永革等[4]的结论相对应，具体表现为：造成龙门山断裂南北两端、秦岭南缘断裂、鲜水河断裂东南端、东昆仑断裂、白玉断裂的库仑破裂应力增加。其中龙门山断裂南北端增加最为明显，最大达258 700 Pa（约0.259 MPa），远比静态应力触发阈值（0.01 MPa）大，秦岭南缘断裂次之（0.047 MPa），白玉断裂、东昆仑断裂的库仑应力值均在触发阈值附近，因此，这4 条断裂的地震危险性值得重视。另外，水平面应力场结果显示，此次地震为一次单侧破裂事件，这与盛书中等[43]利用构造应力场均匀性和杜义等[10]测量擦痕的地质方法的研究结果较一致。主应力方向与背景构造应力场几近相反，一定程度上抵消了构造应力场。

本文在计算汶川地震产生的水平面应力变化时基于半无限均匀弹性空间模型，但在计算中未考虑介质体存在背景构造应力场。所以本研究根据背景构造应力场的主张应力轴、中间应力轴、主压应力轴的方向和R值，假设主张应力大小为20 MPa，利用Matlab 中的interp2 函数进行二维数据内插值（插值方法：双三次插值bicubic），计算出研究区域的构造应力场（图6），并以此讨论本次地震产生的水平面应力变化对构造应力场的影响。

从图6 可以看出，在考虑本次地震产生的应力变化后，构造应力场的水平面应力仅在震中附近（31.0°—31.6°N，103.5°—104.5°E）产生了明显的变化，而在远离震中处引起的变化几乎可以忽略不计。从构造应力场水平分量的最大主应力和最小主应力的方向来看，最大偏转量达19°左右，对背景构造应力场有一定影响。同时需要注意的是，在计算研究区域构造应力时，主张应力假设为20 MPa，随着主张应力值的逐渐扩大，本次地震对构造应力场的影响会逐渐减弱。这从另一方面说明背景构造应力场是相对稳定的，即使是大地震，也只能改变地震震源附近的小范围内的背景构造应力场，不能改变构造应力场的大趋势分布。

图6 2008 年汶川地震震前（a） 及震后（b） 构造应力场的水平分量Fig.6 Horizontal components of tectonic stress field before （a） and after （b） the 2008 Wenchuan earthquake

本文在计算本次地震对周围断层的影响时为一级近似估算，因此，采用USGS 提供的相对简单的弹性半空间破裂模型，实际的不均匀地球介质一定程度上会影响计算结果的准确性；震后地球的粘弹性松弛效应会导致应变扩散[46-47]，但是这种粘弹性松弛效应对发震间隔小于30 年的地震的影响可以忽略不计[3]。

附录

附表I 见电子版（10.19987/j.dzkxjz.2021-078）。

致谢

本研究使用万永革研究员提供的程序，采用GMT 绘图软件[48]，在此一并表示感谢！

表I 同一地震多个震源机制的中心解结果Table I Central solution results of multiple focal mechanisms of the same earthquake

续表I

续表I

续表I

续表I

续表I

续表I