0 引言

巴颜喀拉块体东北缘位于南北地震带中段，区域内分布有龙门山断裂、东昆仑断裂、西秦岭北缘断裂等多条大型活动断裂带，为现今地壳变形显著、地震活跃的区域。仅20世纪以来，该区域就发生了1933年叠溪7.5级、1976年松潘平武7.2级、2008年汶川8.0级、2017年九寨沟7.0级等多次7级以上强震。对该区域的地震危险性跟踪一直是地震科技工作者关注的焦点，M7专项工作组根据历史地震活动、大地测量及地球物理场数据对该区域给出了甘青川交界危险区、西秦岭中西段危险区以及六盘山南-西秦岭东危险区等强震高风险区域。汶川8.0级地震后，一些学者(Xiong et al，2010；Toda et al，2008；Wan et al，2010；单斌等，2009)基于历史强震演化的断层库仑应力变化的研究结果表明，东昆仑断裂东段以及西秦岭北缘断裂西段的地震危险性值得关注。此外，该区域还存在历史地震破裂空段、低b值、小震活动图像、地震矩加速释放等异常(陈为涛等，2013；M7专项工作组，2012)以及跨断层水准等地球物理场观测数据异常^①。

① 中国地震局第二监测中心, 2016, 2017年度地震趋势会商报告。

库仑应力研究作为研究地震轮回中构造应力场时空变化的重要组成部分，已成为研究断层应力状态及区域地震危险性的重要手段之一。King等(1994)研究指出，1992年美国Landers 7.3级地震在Big Bear产生的0.2MPa左右的应力加载量加速了Big Bear 6.5级地震的发生。1999年Izmit 7.4级地震发生前，Stein等(1997)等基于库仑应力计算指出该区域地震危险性较高。而2004年Parkfield 6.0级地震的发生则受1983年Coalinga 6.5级地震及Nuñez多个中等地震的应力卸载的影响，出现了约11年的时间延迟(Bakun et al，1985)。此次九寨沟地震就发生在M7专项工作组(2012)给出的甘青川交界危险区内，基于库仑应力的强震作用研究表明，区域内1976年松潘平武7.2级地震、2008年汶川8.0级地震等都对此次九寨沟地震发震断层有较强的库仑应力加载作用(徐晶等，2017；单斌等，2017)。

有限元数值模拟技术作为探索构造变形等的虚拟实验手段，在地震危险性分析中得到了广泛应用(李玉江等，2013；Moreno et al，2010；Hergert et al，2010；陈连旺等，2008；祝爱玉等，2015)。为进一步分析九寨沟地震对巴颜喀拉块体东北缘主要活动断裂的加卸载效应以及该区域未来可能的强震风险，本文构建了巴颜喀拉块体东北缘三维粘弹性有限元模型，模拟计算了九寨沟地震对周边活动断裂的加卸载效应，分析了震后地震活动及其与九寨沟地震间的可能关系以及未来可能的强震风险。

1 模型构建 1.1 几何模型及物性参数

本文在详细调研地震地质资料的基础上，通过对所调研的地质资料进行假设构建了巴颜喀拉块体东北缘有限元模型：模型纵向深100km，地壳分层、速度结构以及密度结构根据Crust1.0速度结构给出(Laske et al，2013)；模型中包含礼县-罗家堡断裂、虎牙断裂、雪山断裂、东昆仑断裂、龙门山断裂、青川断裂、青川-平武断裂、临潭-宕昌断裂、白龙江断裂、塔藏断裂、龙日坝断裂、岷江断裂、西秦岭北缘断裂以及基于有限断层反演给出的此次九寨沟地震发震断裂等等；区域孕震层深度(段星北，1997；张国民等，2002)，除九寨沟地震发震断层外，假定其他断层面下边界深度为20km；除九寨沟地震发震断层施加位错荷载外，其余断层均为摩擦接触断层，断层空间几何结构简化为由断层地表出露沿着倾角向下延伸的空间斜面；下地壳及地幔粘滞性结构根据邵志刚等(2008)、徐晶等(2017)的研究给出。地壳单元尺寸控制在10km内，地幔单元尺寸控制在20km内，断层面上单元尺寸控制在3km内，模型共包含359676个单元、67796个节点，构建的有限元模型及对应的物性参数分别见图 1和表 1。

1.2 边界约束及模拟方案

地震是构造应力作用下断层应力状态达到一定程度而断层破坏形成的自然灾害。在数值模拟中，由于我们并不知道断层当前绝对应力状态以及断层自身的承载极限，因此，很难做到对地震自发的模拟。为此，本文首先对模型施加以GPS约束的边界荷载来模拟初始构造应力场；然后，再施加地震同震位错，模拟计算给出九寨沟地震后区域的模拟构造应力场；最后，从震后模拟应力场中扣除模拟初始应力场，给出九寨沟地震的发生对区域主要活动断裂的加卸载效应。其中，九寨沟地震位错模型参考波形反演结果(http://www.cea-igp.ac.cntpxw275883.html)给出。所有模拟计算过程均在有限元开源代码Pylith(https://geodynamics.org/cig/software/pylith/)平台下进行。

1.3 库仑应力计算

三维有限元模拟软件并不能直接给出断层面上的库仑应力，因此，需要通过计算给出的应力张量来计算断层面上的正应力及剪应力，进而计算断层面库仑应力变化。图 2为断层面主方向($\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over 1} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over 2} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over 3} $)与模型坐标系($\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over x} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over y} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over z} $)相对关系，其中，$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over 1} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over 2} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over 3} $分别为断层走向(面向走向，右手在上盘)、倾向(逆冲为正)以及张拉方向(张为正)；$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over x} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over y} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over z} $分别为地理E、N向以及垂直地表向上的方向。α为断层走向；β为断层倾角。将断层面应力分量在模型坐标系($\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over x} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over y} $、$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over z} $)下表示为

$ s = {\left[ {{\sigma _{xx}}{\sigma _{yy}}{\sigma _{zz}}{\sigma _{xy}}{\sigma _{yz}}{\sigma _{xz}}} \right]^T} $

(1)

其中，T为矩阵转置符号。

由图 2(b)可知，将模型坐标系下的应力张量绕z轴逆时针旋转90°-α，然后，再绕主轴1逆时针旋转β，便可得到主方向下的应力张量。则

$ \left[ \begin{array}{l} {\sigma _{11}}\;\;{\sigma _{12}}\;\;{\sigma _{13}}\\ {\sigma _{12}}\;\;{\sigma _{22}}\;\;{\sigma _{23}}\\ {\sigma _{13}}\;\;{\sigma _{23}}\;\;{\sigma _{33}} \end{array} \right] = \mathit{\boldsymbol{T}} * \left[ \begin{array}{l} {\sigma _{xx}}\;\;{\sigma _{xy}}\;\;{\sigma _{xz}}\\ {\sigma _{xy}}\;\;{\sigma _{yy}}\;\;{\sigma _{yz}}\\ {\sigma _{xz}}\;\;{\sigma _{yz}}\;\;{\sigma _{zz}} \end{array} \right] * \mathit{\boldsymbol{T'}} $

(2)

其中，T = A*B* $\left[ \begin{array}{l} 1\;\;0\;\;0\\ 0\;\;1\;\;0\\ 0\;\;0\;\;1 \end{array} \right] = \left[ \begin{array}{l} {\sigma _{11}}\;\;{\sigma _{12}}\;\;{\sigma _{13}}\\ {\sigma _{12}}\;\;{\sigma _{22}}\;\;{\sigma _{23}}\\ {\sigma _{13}}\;\;{\sigma _{23}}\;\;{\sigma _{33}} \end{array} \right]$，且

$ \mathit{\boldsymbol{A = }}\left[ \begin{array}{l} 1\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;0\\ 0\;\;\;\cos \beta \;\;\sin \beta \\ 0\;\; - \sin \beta \;\;\cos \beta \end{array} \right], \mathit{\boldsymbol{B = }}\left[ \begin{array}{l} \;\cos \left({90^\circ - \alpha } \right)\;\;\;\;\;\cos \left({90^\circ - \alpha } \right)\;\;\;\;0\\ - \sin \left({90^\circ - \alpha } \right)\;\;\;\;\cos \left({90^\circ - \alpha } \right)\;\;\;\;0\\ \;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;1 \end{array} \right] $

由式(1)、(2)给出的新的应力分量σ₃₃即为断层面上的正应力变化Δσ_n，σ₂₃、σ₁₃分别为断层面上沿倾向和走向的剪应力分量。令断层面优势破裂方向的滑动角为θ，则断层面优势破裂滑动方向的剪应力变化为

$ \Delta \tau = {\sigma _{23}}\;\sin \theta + {\sigma _{13}}\;\cos \theta $

(3)

根据断层面库仑应力变化的定义(King et al，1994)，给出断层面优势破裂方向的库仑应力变化为

$ \Delta {\sigma _{\rm{f}}} = \Delta \tau + \mu \Delta {\sigma _{\rm{n}}} = \mu \Delta {\sigma _{33}} + {\sigma _{{\rm{23}}}}\sin \theta + {\sigma _{{\rm{13}}}}\cos \theta $

(4)

其中，μ为断层面有效摩擦系数。

2 模拟结果分析 2.1 模拟初始应力场

以1999~2007年GPS观测速度为约束，通过本文构建的有限元动力学模型计算了巴颜喀拉块体东北缘的主应力场(图 3)。模拟结果与基于小震震源机制反演的区域构造应力场(王晓山等，2015)以及基于GPS推算的区域应变率场(陈长云等，2013)基本一致。区域主要处于NW-SE拉张、NE-SW西向挤压的动力环境。九寨沟地震震中南侧至龙门山断裂，主张应力逐渐减小，主压方向逐渐偏转并呈现出近EW向挤压态势。九寨沟地震震中附近主应力方向与USGS给出的震源机制解(https://earthquake.usgs.gov/earthquakes/eventpage/us2000a5x1#moment-tensor)的P轴(方位角104°)和T轴(方位角204°)基本一致。

2.2 九寨沟地震及其对周边断裂的加卸载效应

以同震位错为荷载模拟给出的九寨沟地震震中附近的同震位移场如图 4所示。由图 4可见，震中最大同震位移约为9cm，九寨沟地震的发生对东昆仑断裂、塔藏断裂、西秦岭北缘断裂东段以及龙门山断裂均为左旋作用，同震位移场的左旋作用有利于左旋活动为主的东昆仑断裂、塔藏断裂以及西秦岭北缘断裂东段断层面的应力积累；而对于具有右旋逆冲活动的龙门山断裂则为卸载作用。同理，对西秦岭北缘断裂西段应为卸载影响。为进一步定量分析九寨沟地震对各断裂的加卸载效应，我们基于本文模拟结果，根据式(2)、(4)计算了九寨沟地震在这些断裂断层面上引起的库仑应力变化(图 5)。

由图 5可见，九寨沟地震对塔藏断裂西段、虎牙断裂加载效应显著，达7056Pa；东昆仑断裂东段、龙日坝断裂、迭部-白龙江断裂西端和东端、临潭-宕昌断裂东段以及西秦岭北缘断裂东段也均受到不同程度的加载作用；塔藏断裂东段以及岷江断裂受九寨沟地震的影响，断层面库仑应力变化表现为卸载，最大卸载量达73984Pa。位于此次九寨沟地震破裂面正前方的雪山断裂断层面库仑应力变化较为复杂：由图 4模拟同震位移场可知，九寨沟地震对雪山断裂西段的影响表现为挤压，对东段的影响表现为拉张，根据库仑应力计算公式可知，当断层面应力变化主要表现为垂直于断层面的挤压/拉张时，有效摩擦系数的取值对库仑应力变化影响显著，甚至可以影响其极性变化。Wan等(2010)在分析汶川地震对周边断裂影响时也给出同样的结论。从本文模拟结果看，当有效摩擦系数较小(μ=0.1)时，九寨沟地震对雪山断裂西段表现为库仑应力加载，对东段表现为卸载；当有效摩擦系数较大(μ=0.7)时，对西段表现为卸载作用，而对东段表现为加载(图 5)。

此外，此次九寨沟地震还对汶川地震破裂区东北端的青川断裂以及青川-平武断裂的西南段有一定的库仑应力加载，最大加载量可达1139Pa。邵志刚等(2010)研究汶川地震的震后影响发现，青川断裂库仑应力加载非常显著，实际上青川地区也是汶川地震余震活动非常活跃的区域(石军，2008)，但从2013年开始，该区域地震活动水平明显减弱。此次九寨沟地震后不到2个月，9月30日青川(32.27°N，105°E)再次发生5.4级地震。分析认为，此次青川地震可能是由九寨沟地震的触发作用所致。

3 结论与讨论

基于构建的三维有限元模型，以同震位错为荷载模拟了九寨沟地震的发生及其对巴颜喀拉块体东北缘主要活动断裂的加卸载效应。数值模拟结果显示，九寨沟地震的发生对处于甘青川交界危险区内的东昆仑断裂东段至塔藏断裂西段库仑应力加载效应显著，最大可达7056Pa；对六盘山南-西秦岭东危险区的西秦岭北缘断裂东段也有一定程度的加载效应，约为346Pa。此外，计算结果还显示，九寨沟地震对龙日坝断裂、虎牙断裂、青川-平武断裂西段、迭部-白龙江断裂西段和东段、临潭-宕昌断裂东段等也均有不同程度的库仑应力加载效应；对岷江断裂、塔藏断裂东段库仑应力卸载效应显著，最大卸载量可达73984Pa。

地震是区域构造应力积累到一定程度且超越断层面承载极限而引起断层剪切破坏的动力过程。但实际上我们很难获取断层面的绝对应力状态和断层面的承载极限。库仑应力作为由地震发生而引起的构造应力场变化在断层面优势破裂方向的体现，在地震危险性定性分析判断中有着广泛的应用。汶川地震后，许多学者通过库仑应力计算分析了汶川地震对区域活动断层的加卸载效应及地震危险性(Toda et al，2008；单斌等，2009；邵志刚等，2010；李玉江等，2013)。随后的2013年芦山7.0级地震、岷县漳县6.6级地震以及此次九寨沟7.0级地震均发生于文献报道的库仑应力加载区域。这些强地震的时空丛集活动也表明，巴颜喀拉地块东段在未来一段时间内可能继续保持这种活跃态势。从库仑应力变化的角度来看，此次九寨沟地震对甘青川交界危险区内的主要活动断裂(东昆仑断裂东段和塔藏断裂)以及位于六盘山南-西秦岭东危险区的主要活动断裂(西秦岭北缘断裂东段)，均具有不同程度的库仑应力加载效应。