2017-08-08四川九寨沟地区发生M_S7.0地震，震中位于103.89°E、33.21°N，震源深度约为20 km，发震断层为NW走向、左旋走滑型的虎牙断裂北段^[1]，余震呈NWW向条带状分布，长约42 km，主震位于余震带中央^[2]，如图 1所示。震区位于青藏高原巴颜喀拉块体东北边界，是高原内物质侧向挤出与四川盆地碰撞挤压的边缘地带，具有孕育强震的构造背景。2017年九寨沟地震是该区继2008年汶川M_S8.0地震和2013年芦山M_S7.0地震之后发生的又一次强震，其中，汶川地震和芦山地震均发生在逆冲推覆的龙门山断裂带上。由于芦山地震震级相对较小，且与九寨沟地震震中距离较远，对九寨沟地震几乎无影响^[3-6]。然而，2008年汶川地震是罕见的逆冲兼走滑型特大地震，主震及余震的破裂自震源处沿龙门山断裂带向NE方向延伸，长度达360 km，地表破裂带长达240 km^[7]，破裂带北端距九寨沟地震震中约150 km。GPS观测资料显示，受2008年汶川地震影响，巴颜喀拉块体东缘的水平运动速率明显比震前增强1.5~3倍^[8]。

地震在释放震源区地壳积累应力的同时，会对周边区域产生应力扰动，即介质瞬态弹性形变的同震应力变化和长期粘弹性效应的震后应力调整。根据弹性位错理论，以地震位错为扰动源可以计算同震及震后应力变化量^[9]，并采用库仑应力变化理论^[10]定量计算库仑应力变化量(以下简称ΔCFS)，据此判断对周边断层是否产生促进或抑制滑动的作用。该方法在研究主余震关系及区域地震触发关系中得到广泛应用^[11-13]。

前人根据该理论定量计算并探讨了汶川地震与九寨沟地震的关系，但不同研究者之间的结论各异。徐晶等^[4]的研究表明，若采用王卫民^[14]的汶川地震破裂模型计算九寨沟地区的ΔCFS，结果为负值，若采用Ji^[15]的破裂模型计算则结果为正值，不能确定汶川地震对九寨沟地震的影响效应。Jia等^[6]则认为汶川地震延迟了九寨沟地震的发生。

本文利用2008年汶川地震破裂模型、2017年九寨沟地震接收断层参数和断层面有效摩擦系数等不同模型参数，计算汶川地震造成的库仑应力变化，分析前人研究结论产生差异的原因，并探讨汶川地震对九寨沟地震的触发作用。

1 计算方法

本文采用Wang等^[9]根据半无限空间位错理论开发的计算程序PSGRN/PSCMP，以2008年汶川地震滑动分布为扰动源，计算其产生的同震及震后应力变化，将其投影到接收断层面上，得到断层面上的剪应力变化量Δτ及正应力变化量Δσ_n，并根据库仑应力变化理论计算断层面上的ΔCFS：

$ \Delta {\rm{CFS}} = \Delta \tau - \mu ^{\prime} \times \Delta {\sigma _n} $

(1)

式中，μ′为考虑介质孔隙压等因素的有效摩擦系数。可以看出，不仅作为扰动源的地震滑动分布模型能够影响计算结果，接收断层的几何参数和有效摩擦系数等同样会对结果产生影响。

地震破裂模型是由不同的学者和研究机构根据不同类型的大地测量数据(地震波、GPS等)通过反演计算得到。汶川地震发生之后，不同学者和机构发布了大量的滑动分布模型^[14-17]，本文选取3个基于不同数据资料的模型以示对比：1)模型1，来源于USGS网站，基于宽频带远震P波、S_H波及长周期面波数据反演；2)模型2，由Wang等^[16]根据高分辨率GPS观测数据反演；3)模型3，Fielding等^[17]联合远场地震波、GPS和InSAR数据反演。

根据研究区P波、S波速度信息及岩石圈粘滞系数等^[18]，本文计算采用的地壳分层结构模型主要分为4层：上地壳、中地壳、下地壳和岩石圈地幔，具体信息见表 1。

2 计算结果 2.1 不同的接收断层参数对计算结果的影响

选取不同来源的九寨沟地震震源机制解节面之一作为接收断层参数，分析其对计算结果的影响。其来源主要有3种：1)GCMT(http://www.globalcmt.org/CMTsearch.html)，走向为151°，倾角为79°，滑动角为－8°；2)CEA-IGP(http://www.cea-igp.ac.cn/tpxw/275881.html)，走向为151°，倾角为83°，滑动角为0°；3)USGS(https://earthquake.usgs.gov/earthquakes/eventpage/us2000a5x1#moment-tensor)，走向为153°，倾角为84°，滑动角为－33°。

采用汶川地震破裂模型2^[16]，有效摩擦系数μ′取值均为0.4，计算汶川地震后虎牙断裂带北部(图 1中AA′剖面)的同震及震后(到2017年)ΔCFS总量，计算深度范围0~24 km，包含余震精定位深度(4~20 km)和震源深度(20 km)。同时计算九寨沟地震震源处的ΔCFS，结果见图 2和表 2(单位kPa)，其中，图 2(a)和2(b)分别为采用GCMT和CEA-IGP接收断层参数时计算的ΔCFS，二者走向一致，倾角相差4°，滑动角相差8°。结果显示，在这两种接收断层参数条件下，ΔCFS的分布形态基本一致，计算深度10 km处的ΔCFS分别为2.13 kPa和2.08 kPa，16 km处的ΔCFS分别为2.67 kPa和2.57 kPa，差异较小。

图 2(c)为采用USGS接收断层参数时计算的ΔCFS结果，由于参数中的走向和倾角与GCMT和CEA-IGP的差异较小，但滑动角差异较大，分别相差25°和33°，导致该条件下计算的ΔCFS分布形态明显不同，且计算深度10 km处的ΔCFS为1.77 kPa，16 km处为2.39 kPa，同样与前两个参数的结果相差较大。由此可知，与单斌等^[3]的结论一致，差异较大的滑动角对计算结果的影响非常明显。

2.2 不同汶川地震破裂模型对计算结果的影响

采用3种不同的汶川地震破裂模型，利用USGS给出的九寨沟地震震源机制解节面之一(走向为153°，倾角为84°，滑动角为－33°)作为接收断层参数，在相同的有效摩擦系数(μ′取0.4)下分别计算汶川地震后虎牙断裂带北部(图 1中AA′剖面)的同震及震后(到2017年)ΔCFS，结果见图 3。由破裂模型1(图 3(a)和3(b))的计算结果可知，九寨沟地震主震及余震均落在ΔCFS增加的区域——断裂带北(右)段，同震值约为2~3 kPa，震后受粘弹性松弛效应影响，ΔCFS持续增加，变化总量约为3~8 kPa；而在断裂带南(左)段的非震源区，同震ΔCFS为负值，震后继续减小。破裂模型2(图 3(c)和3(d))的结果显示，九寨沟地震震源处的同震ΔCFS约为3 kPa，75%的余震落在ΔCFS增加区，仅断裂带南段有少量余震落在ΔCFS减少(约－1 kPa)区；断裂带北段震后ΔCFS略有增加，但增幅较小，断裂带南段ΔCFS则小幅减小。破裂模型3(图 3(e)和3(f))的结果显示，主震和余震均在同震ΔCFS为2~3 kPa的区间。在断裂带南段深度18 km以下，同震ΔCFS为正值，约为2 kPa，而另两个模型的结果均为负值。另外，模型3断裂带南段同震及震后ΔCFS的总量减小变为负值，北段则为正值且有所增加，达到3~6 kPa。

综上所述，3种不同模型的计算结果存在差异性，但也存在共性。受震后粘弹性松弛效应的影响，虎牙断裂带剖面上ΔCFS的震后变化总量与同震相比，在断裂带北段均有所增加，如在九寨沟地震震源(图 3中黑色五角星)处，ΔCFS可达7 kPa，而同震ΔCFS仅为3 kPa左右；断裂带南段ΔCFS则与之相反，有继续减小的趋势。由于本文选取的地壳结构模型在深度22 km处为分界面，上地壳和中地壳有明显的粘滞系数变化，因此同震及震后总量值显示，在深度22 km处有ΔCFS阶变。

另外，在有效摩擦系数μ′取0.4、USGS为接收断层参数的条件下，分别计算3种破裂模型在九寨沟地震震源处的ΔCFS(同震及震后)，结果见表 3(单位kPa)。由表 3可知，在深度10 km处计算结果分别为3.28 kPa、1.77 kPa和4.27 kPa，16 km处计算结果分别为3.53 kPa、2.39 kPa和4.75 kPa，采用不同破裂模型计算的结果差异明显。由此可知，作为扰动源的汶川地震破裂模型的选取很大程度上影响了ΔCFS的计算结果。

2.3 不同的有效摩擦系数μ′对九寨沟地震的影响

为研究在库仑应力计算中有效摩擦系数的取值对计算结果的影响，利用汶川地震破裂模型3，选取USGS为接收断层参数，μ′分别取0.2、0.4、0.6、0.8，计算九寨沟地震震源处ΔCFS，结果见表 4。由表 4可知，μ′取0.2~0.8，在计算深度分别为10 km和16 km处，ΔCFS分别为5.12~2.57 kPa和5.58~3.08 kPa。从式(1)可以看出，在断层几何形态一致的情况下，断层面上的剪应力变化量和正应力变化量相同，μ′越大则ΔCFS越小，符合摩擦系数越大物体越不易滑动的物理规律，与徐晶等^[4]的研究成果相同。

3 讨论

为分析汶川地震对九寨沟地震的应力扰动作用，分别选取3个不同的汶川地震破裂模型并展开计算。结果表明，汶川地震对九寨沟地震震源处产生的同震ΔCFS约为1~5 kPa，同震及震后ΔCFS总量约为3~7 kPa。九寨沟地震主震及大部分余震均落在发震的虎牙断裂带北段ΔCFS增加(1~7 kPa)区。在震源深度约10 km处的同震计算结果与前人对比结果见表 5，其中，徐晶等^[4]选取CRUST1.0地壳模型，在2个不同接收断层参数(GCMT和USGS)条件下分别展开计算，采用王卫民等^[14]的汶川地震破裂模型计算的结果分别为－15.9 kPa和－13.2 kPa，而采用Ji等^[15]模型计算结果则完全相反，分别为15.7 kPa和11.0 kPa，均超出触发阈值(10 kPa)。汪建军等^[5]选取巴颜喀拉块体的CRUST1.0平均地壳模型，破裂模型采用USGS网站更新模型(模型1)，相同计算条件下结果分别为17.4 kPa和9.6 kPa；采用Wang等^[16]的模型(模型2)，计算结果均比模型1的值大，分别为26.3 kPa和18.3 kPa。本文选取区域地壳模型，在相同计算条件下，采用Ji等^[15]破裂模型计算的震源处ΔCFS分别为8.4 kPa和5.59 kPa，采用模型1、模型2、模型3计算的结果分别为1.6 kPa和0.94 kPa、5.14 kPa和2.94 kPa、2.78 kPa和1.74 kPa，得到的结果值最小。

通过对比分析发现，在相同的地壳模型和接收断层参数条件下，采用不同的破裂模型，甚至可以得到完全相反的结果。Ji等^[15]和王卫民等^[14]的汶川地震破裂模型是震后根据地震波快速反演产出的结果，仅有一个断层面，虽然能够迅速为震后应急救援等工作提供重要参考，但由于受精度限制，对于进一步深入分析研究略显不足。模型1、模型2及模型3则是在地震发生较长时间后，根据大量精细数据(如宽频带远场及近场地震波、高频GPS和InSAR等)通过复杂的计算方法反演得到的，断层面离散程度较高。其中，模型1划分为3个子断层，断层呈现斜立状，模型2和模型3均划分为5个子断层，断层呈叠瓦状分布的弯曲平滑状，更加接近真实情况。因此，在研究强震对后续地震的触发及相关作用时，选取精度较高且更为可靠的破裂模型尤为重要。其次，在采用相同的地震破裂模型条件下，选取不同的地壳分层结构也会得到有差异的计算结果。CRUST1.0是全球三维地壳模型，但由于巴颜喀拉块体东缘地形高度的变化和地壳厚度及结构的复杂性^[19]，CRUST1.0模型在该区域不一定适用，需要参考更多的区域研究成果，采用相对精确的局部区域地壳模型。另外，选取不同的接收断层参数及有效摩擦系数也会对计算结果产生影响。

在本文的计算条件下，汶川地震对九寨沟地震震源处产生的同震及震后库仑应力总量仅为3~7 kPa，尚未达到触发阈值(10 kPa)。GPS观测数据揭示，2008年汶川地震之后该区域的剪应变速率一直较高^[1]，左旋走滑的九寨沟地震发生在虎牙断裂带北段，而与其相隔约50 km的虎牙断裂带南段1976年发生了3次逆冲型强震，揭示该断裂带北段以走滑为主、南段以逆冲为主的复杂构造应力场。地震地质调查资料显示，青藏高原内大型左旋走滑的东昆仑断裂带东端向东扩展，形成岷江断裂、虎牙断裂、塔藏断裂等分支断裂，这些断裂吸收了东昆仑断裂带大部分的左旋走滑分量，承担着区域构造的转换形变，造成地壳缩短和岷山隆起，控制着区域的强震孕育。大地电磁资料分析认为，岷江断块中上地壳存在低阻的粘性流体，导致高原内软流圈物质侧向流动，驱动巴颜喀拉东缘上地壳沿着中上地壳低阻层向东运移，导致虎牙断裂带北部走滑、南部逆冲^[19]。由此认为，2017年九寨沟地震是由区域地壳结构和构造环境共同主导，断层应力受控于构造活动加载，汶川地震对九寨沟地区的影响作用非常有限，尚未达到触发地震的阈值。

4 结语

通过计算汶川地震对九寨沟地震震源处产生的同震及震后ΔCFS，分析其对九寨沟地震的影响，得出以下结论：

1) 在本文选取的汶川地震破裂模型和地壳模型基础上，采用不同接收断层参数下的同震ΔCFS值为1~5 kPa，受粘弹性效应影响的同震及震后共同库仑应力增量约为3~7 kPa，尚未达到地震触发阈值(10 kPa)。汶川地震对九寨沟地震的触发作用有限，区域地壳流变结构和构造活动的加载可能是诱发九寨沟地震的主要因素。

2) 不同研究机构给出的九寨沟地震震源机制解的走向和倾角基本一致，但滑动角差异较大，对计算结果影响较大。另外，有效摩擦系数μ′也影响计算结果，μ′值越大，ΔCFS越小，断层越不易滑动。

3) 不同的汶川地震破裂模型选取对计算结果的影响最大，甚至会得到相反结果，这也可能是造成前人得到不同研究结论的重要原因之一。在大震对后续地震触发作用的相关研究中，可靠的破裂模型选取十分重要，不可忽视。