0 引言

2008年5月12日在龙门山地区发生汶川M_S 8.0地震，2013年4月20日该地区发生芦山M_S 7.0地震，2次地震分别发生在龙门山断裂带的中段和南段。龙门山断裂带位于巴颜喀拉块体和华南地块交界处，自北西到南东由龙门山后山断裂、中央断裂、前山断裂和山前断裂呈叠瓦状组成(徐锡伟等，2008；高翔等，2013)。龙门山断裂带在汶川地震前未发生过7级以上地震，GPS观测和地质研究结果表明，该断裂自晚更新世以来活动性较弱(Gan W J et al，2007)，横跨整个断裂带的滑动速率不超过2 mm/a，单条断裂的滑动速率不到1 mm/a(邓起东等，1994；张培震等，2008)，因此，汶川地震和芦山地震的发生引起了地震学家的关注。张培震等(2008)认为，汶川地震是一次发生在大陆内部的高角度逆冲型地震，可见青藏高原和华南地块之间的相对运动是2次地震的主要动力来源。张培震等(2009, 2012)进一步提出汶川地震的多单元组合模式，Zhu S B等(2010，2013)使用数值模拟方法验证了该结论。

区域构造应力场控制着地震的孕育、发展和发生，因此研究构造应力场的分布和演化，有助于了解地震的孕育和发生过程(李玉江等，2013)。数值模拟方法因其灵活性、低成本和短周期等特点，成为定量研究构造应力场的重要方法(陈连旺等，2008；王辉等，2008)。汶川地震和芦山地震发生后，国内外诸多学者采用数值模拟方法，对龙门山地区构造应力场进行研究。如：Parsons T等(2008)采用弹性有限断层模型，分析汶川地震对龙门山地区应力分布的影响，认为本次地震使龙门山断裂带南段应力增加超过0.1 MPa；Toda S等(2008)和万永革等(2009)计算汶川地震后周围断层的库仑破裂应力变化，认为龙门山断裂带北部和最南端、鲜水河断裂带南段、岷江断裂的地震危险性增加；单斌等(2009, 2013)和董培育等(2013)采用弹性位错理论和岩石圈分层模型，计算汶川地震和芦山地震后周边断层应力变化，结果认为，汶川地震可能有效促进了芦山地震的发生，Wang Y Z等(2014)采用三维粘弹性分层位错模型分析，也得到相同结论。这些工作对于认识龙门山地区构造应力场的分布和演化具有重要意义。近年来，随着对龙门山断裂带及地震孕育机制研究的进一步深入(张培震等, 2008, 2009, 2012；邓起东等, 2010, 2011)，以及在地球深部结构探测与反演、GPS观测等方面取得的重要进展，结合不断成熟的有限元建模技术，可以建立更加合理的模型来研究关注的问题。

本文在前人研究基础上，考虑龙门山地区主要活动断层分布特征，使用最新的深部反演结果和GPS观测资料，建立反映地表起伏和岩石圈分层的龙门山地区三维粘弹性有限元模型。采用Maxwell粘弹性本构模型，模拟岩石圈介质在长期地质演化中表现出的流变性质及发生地震时的弹性效应。使用GPS观测值对模型进行约束，重新构建龙门山地区现今构造应力场。在此基础上，采用准静态方法，±次模拟1900年以来龙门山地区发生的7级以上历史地震序列，分析应力场演化与强震之间的关系。

1 龙门山地区三维粘弹性有限元模型 1.1 模型建立

根据龙门山地区地质构造特征及地震分布特点，选取101.1°E—106.85°E，28.85°N—33.7°N的范围作为研究区域，建立三维地质模型，模型中包括龙门山断裂带(后山、中央和前山断裂)、岷江断裂、虎牙断裂、龙日坝断裂、鲜水河断裂带南段等该区域的主要活动断裂，为了更加接近实际地质构造，在模型中将龙门山断裂带的后山、中央和前山断裂设置为倾角70°的高角度倾斜断层，其他断裂设置成直立断层，见图 1。模型中断层分布在地表至地下20 km处，断层宽度设为5 km。由于选取的研究区域不大，可以用笛卡尔坐标系代替球坐标系进行计算(曹建玲等，2009)，模型中设东向为X轴正向，北向为Y轴正向，垂直向上为Z轴正向，采用Guess-Kruger投影算法，将经纬度坐标投影计算成XYZ直角坐标。根据深部反演结果，所选区域的青藏高原东部和四川盆地的地壳厚度、介质参数等均存在较大差异，为此，在模型中引入地球深部反演结果，将模型在横向分为东侧的青藏高原东缘和西侧的四川盆地2地质部分，在纵向(Z轴)沿地表至地下100 km处分为4层，分别为上地壳、下地壳、岩石圈上地幔和软流圈上地幔。模型中使用分辨率为1弧度分的ETOPO1高程数据获得地表的Z坐标，Moho面深度采用Li Y H等(2014)的结果。使用4节点的四面体单元对三维模型进行剖分，得到83 758个节点和403 930个四面体单元。本文三维模型的建立、剖分和计算均在通用有限元计算程序ADINA(Automatic Dynamic Incremental Nonlinear Analysis)(ADINA R&D，2013)中完成。

1.2 模型材料参数及本构关系

岩石圈在漫长地质演化过程中主要表现为流变性，而短时间变化可视为弹性过程，该现象可以用粘弹性本构关系来描述。本文采用Maxwell粘弹性本构模型，其本构关系为

$K\left( t \right) = K\quad \quad G\left( t \right) = \mu {e^{ - \mu t/\eta }}$

(1)

其中，K(t)和G(t)分别为Maxwell模型的体积模量和剪切模量，η为介质的粘滞系数，K为弹性体积模量，μ为弹性剪切模量。弹性参数可由线弹性本构模型的杨氏模量E和泊松比ν计算得到。杨氏模量E和泊松比ν可以根据研究区深部反演得到的三维速度模型和密度模型(Li Y H et al，2013)计算得到。然后将弹性参数在每个单元的质心处插值，得到每个单元的弹性参数。为了在模型中体现断层单元，根据断块介质比非断块介质软的特点，在模型中将断层单元的杨氏模量取为其两侧介质杨氏模量平均值的1/3，泊松比取为比两侧介质的平均值高0.02(王辉等，2006；胡幸平等，2012；李玉江等，2013)。粘滞系数η参考石耀霖等(2008)、孙玉军等(2013)的结果。该方法充分利用深部反演成果，使模型材料的分辨率与反演的速度结构分辨率一致。

1.3 模型加载条件和边界条件

为了使模型可以计算，需要对模型施加加载条件和边界条件。使用地表观测的GPS值(张培震等，2008)，在模型的4个侧边界插值得到位移载荷，见图 2。由于地球深部水平位移与地表GPS的关系尚无定论，因此本文假设加载条件从地表到地下为同一值(李玉江等，2013)。另外，由于GPS观测值表示地壳水平位移的年变化量，因此将加载时间函数设置为每年增加1个载荷单位，使加载条件在物理意义上是合理的。模型底边界条件为沿垂直(Z轴)方向固定，沿XY水平方向自由，上表面为自由边界。

1.4 模型可靠性检验

设置好材料参数和约束条件后，对模型加载10万年，使用ADINA的重启算法解决因引入重力而导致的模型塌陷问题，最终得到稳定的应力状态，将其作为研究区现今构造应力场。为了检验模型的可靠性，将模拟结果与GPS观察值进行比较，见图 3。

由图 3可以看出，在大部分区域，地壳水平运动速率的模拟值与观测值的速率和方向基本吻合，只在少数观测点上存在突出差异，模拟结果基本能够反映龙门山地区现今构造变形状态，可以作为研究区背景应力场，用来模拟历史强震序列。

2 历史地震序列模拟 2.1 地震模拟方法

地震是在构造应力场下发生的，与漫长的地质演化历史相比，地震的发生过程是一个短暂的瞬态，可视为力学上的一次弹性过程。如果忽略地震发生时断层的动态破裂过程，只关注地震发生前后的应力状态，可以采用准静态的方法来模拟地震的发生。为此，本文使用ADINA的重启动算法引入研究区现今构造应力场，作为地震发生前的背景应力场，并采用降低发震断层材料参数的方法模拟实际地震的发生。根据弹性力学理论(王敏中等，2011)，对于弹性材料，泊松比ν、体积模量K、剪切模量G的关系为

$\nu = \frac{K}{{2\left( {K + G} \right)}}$

(2)

一般认为材料的泊松比是一个不变量，因此本文在模拟地震时，对发震断层单元的体积模量K和剪切模量G按照相同的比例系数降低，这样既模拟了地震，又使泊松比不发生变化。为了在模拟地震时使断层错动尺度与实际地震相同，选取发生地震时的断层水平滑移量作为约束条件。

2.2 地震序列模拟

选取龙门山地区1900年以来发生的5次M_S 7.0以上地震进行模拟，表 1给出5次地震的相关参数，图 4为所模拟的5次地震的空间分布。模型中，根据震后现场调查得到的地表破裂长度预先在模型中设置发震断层，对于发生时间较早、缺乏现场调查资料的地震，由震级与破裂长度关系的统计公式(Wells D L et al，1994)，计算得到断层的破裂长度。

对于表 1给出的数据，可知：1933年叠溪地震震源参数参考阚荣举等(1977)、万永革等(2007)的结果；1955年康定地震震源参数参考Molnar P等(1984)、阚荣举等(1997)的结果；1976年松潘地震震源参数、破裂长度、水平错动来自岳汉等(2008)；2008年汶川地震震源参数、破裂长度、水平错动数据来自中国地震台网中心、Global CMT、徐锡伟等(2008)；2013年芦山地震震源参数、破裂长度、水平错动数据来自USGS、曾祥方等(2013)、张勇等(2013)。对于无明确来源的破裂长度、水平错动数据，使用Wells D L等(1994)的统计关系给出。

为了方便分析芦山地震对其后续强震的影响，需用到2014年康定M_S 6.3地震的震源参数。康定M_S 6.3地震发生在鲜水河断裂带磨西段，根据Fang L H等(2015)的研究结果，本次地震的震中位置为(101.7°E，30.3°N)，走向142°，倾角82°，滑动角-8°。

3 计算结果及分析 3.1 库仑破裂应力计算方法

根据库仑破裂准则，当岩石趋近破裂时库仑破裂应力(CFS)为

${\rm{CFS = }}\tau {\rm{ + }}\mu \left( {{\sigma _n} - p} \right)$

(3)

其中，μ为摩擦系数，τ为主震产生的应力张量在断层面上的剪切应力大小，以沿断层滑动方向为正，σ_n为断层面上正应力大小，以张应力为正，P为孔隙压力。通常更关心的是库仑破裂应力的变化，当μ不随时间变化时，由式(2)可得库仑破裂应力的变化ΔCFS为

$\Delta {\rm{CFS = }}\Delta \tau {\rm{ + }}\mu \left( {\Delta {\sigma _n} - \Delta p} \right)$

(4)

孔隙压力变化控制着断层面上的有效正应力，对于各向同性均匀介质，孔隙压力对摩擦系数的影响可用等效摩擦系数μ'=μ(1-B)表示，其中B为Skempton系数，范围为0—1(万永革等，2000)，则式(3)变为

$\Delta {\rm{CFS = }}\Delta \tau + \mu '\Delta {\sigma _n}$

(5)

等效摩擦系数μ取值通常在0—0.75，平均值约0.4(Freed A M，2006)。一般认为，库仑破裂应力变化大于0.01 MPa是触发后续地震的阈值(Harris R A，1998)。

3.2 计算结果分析

根据数值模拟结果，由库仑破裂应力计算公式(5)计算20 km处的同震库仑应力变化，其中以现今构造应力场作为地震发生前的背景构造应力场，计算时将应力张量投影到后续地震的发震断层上。计算结果见图 5(色标柱为ΔCFS数值，单位MPa)，图中红色五角星表示后续地震震源位置，黑色虚线和红色虚线分别表示库仑破裂应力变化为0和0.01 MPa的等值线。

由图 5可知，在模拟的5次地震中，有3次对后续地震的同震库仑破裂应力变化大于0，分别为：1933年叠溪地震对1955年康定地震的同震库仑应力变化1.66×10^-4 MPa；1976年松潘地震对2008年汶川地震的同震库仑应力变化2.536 2×10^-4 MPa；2008年汶川地震对2013年芦山地震的同震库仑应力变化1.136 1×10^-2 MPa。表明3次地震加快了其后续地震的发生，特别是汶川地震对芦山地震的库仑破裂应力达到触发阈值0.01 MPa，说明汶川地震可能触发了芦山地震。

3.3 等效应力计算结果

除了从库仑破裂应力角度外，还可以从等效应力角度分析应力演化与地震分布、强震间的相互关系。计算的等效应力由于不受震源机制参数的影响，得到的结果可能更加准确。等效应力是将物体受到的6个应力分量等效成为一个单向拉伸应力。在一般应力状态下，等效应力σ_e表示(杨兴悦等，2013)为

${\sigma _{\rm{e}}} = \frac{1}{{\sqrt 2 }}\sqrt {{{\left( {{\sigma _x} - {\sigma _y}} \right)}^2} + {{\left( {{\sigma _y} - {\sigma _z}} \right)}^2} + {{\left( {{\sigma _z} - {\sigma _x}} \right)}^2} + 6\left( {\sigma _{xy}^2 + \sigma _{yz}^2 + \sigma _{zx}^2} \right)} $

(6)

其中，σ_x、σ_y、σ_zx、σ_xy、σ_yz、σ_zx为应力张量的6个分量。

等效应力反映介质趋于屈服程度的强度，其值的增大会加速介质失稳，使地震提前发生。使用公式(6)分别计算5次模拟地震产生的等效应力变化，计算深度为20 km。图 6给出每次地震对其后续地震的等效应力变化，图中红色五角星表示后续地震震源位置，色标柱数值单位MPa。

由图 6可知，对于5次模拟地震，除2013年芦山地震对2014年康定地震的等效应力变化小于0外，其余4次地震对其后续地震的等效应力变化均为正，表明对后续地震的发生具有加速作用，与采用库仑破裂应力计算得到的结果较为一致。

从等效应力角度还可以研究强震间的相互关系。研究区历史强震震源处等效应力随时间变化曲线见图 7。图中横轴为每次地震前的地震序列，纵轴为所有前续地震在该次地震震源处的等效应力之和。

由图 7可知：1933年叠溪地震对1955年康定地震有促进作用；1933年叠溪地震对1976年松潘地震具有抑制作用，1955年康定地震对1976年松潘地震具有促进作用；1933年叠溪地震和1976年松潘地震对2008年汶川地震具有促进作用，1955年康定地震对2008年汶川地震具有抑制作用；1976年松潘地震和2008年汶川地震对2013年芦山地震具有促进作用，1933年叠溪地震和1955年康定地震对2013年芦山地震具有抑制作用；除叠溪地震和汶川地震对2014年康定地震影响较小外，其余3次地震对2014年康定地震具有抑制作用。

4 结论与讨论

本文基于地质构造、最新地球深部反演结果，建立龙门山地区三维粘弹性有限元模型，以GPS观测值为约束模拟研究区现今构造应力场，将所得应力场与观测值对比，两者吻合较好。以现今构造应力场为背景，模拟龙门山地区1900年以来发生的5次7级以上强震，分别从库仑破裂应力和等效应力角度分析应力场演化对强震的影响及强震之间的相互作用关系，得到以下结果。

(1) 从库仑破裂应力角度，1933年叠溪地震对1955年康定地震的库仑应力变化为1.66×10^-4 MPa；1976年松潘地震对2008年汶川地震的库仑应力变化为2.536 2×10^-4 MPa；2008年汶川地震对2013年芦山地震的库仑应力变化为1.1361×10^-2 MPa，表明模拟的5次地震中，有3次地震对其后续地震的库仑破裂应力变化大于0，特别是汶川地震对芦山地震的库仑破裂应力变化大于0.01 MPa，达到触发阈值。

(2) 从等效应力角度，除2013年芦山地震在2014年康定地震震源处的等效应力小于0外，其余4次地震在其后续地震震源处的等效应力变化均为正，表明4次地震加速了其后续地震的发生。从等效应力角度还能得出地震间的相互关系：1933年叠溪地震对1955年康定地震、2008年汶川地震有较明显的促进作用；1955年康定地震对1976年松潘地震有明显的促进作用；1976年松潘地震对2008年汶川地震、2013年芦山地震有较明显的促进作用；2008年汶川地震对2013年芦山地震有明显促进作用。

对比库仑应力和等效应力2种分析指标，二者计算结果具有较高的一致性，主要差异是，因使用参数不同，计算库仑应力需要用到震源位置和震源机制解，而计算等效应力只需震源位置。因此，在分析不易获得震源机制解的地震时，使用等效应力更加方便。