2017年8月8日，四川省阿坝州九寨沟县发生M7.0地震，是继2008年汶川M8.0地震和2013年芦山M7.0地震后发生在巴颜喀拉块体东部的又一强震，引起了社会及学术界的广泛关注.震中位于巴颜喀拉块体东部岷江、虎牙和塔藏断裂带之间，见图 1.该地震震源机制解显示发震断层为走向北西-南东的左旋走滑断层，余震分布和地震破裂反演结果指示地震破裂深度约为5~20 km (http://www.csi.ac.cn/manage/eqDown/05LargeEQ/201708082119M7.0/zonghe.html).

图 1

Fig. 1

图 1 巴颜喀拉地块东部构造背景图 黑色粗线条为研究区主要活动断裂带，地震震源球为区内1976年以来M6.0以上地震分布，黑色走滑型和逆冲型震源球分别代表九寨沟地震和汶川地震，其中北西-南东向节面为断层面(http://www.globalcmt.org)；九寨沟地震震中附近深灰色圆点为九寨沟地震M4.0以上余震分布(https://earthquake.usgs.gov/earthquakes/).灰色箭头为汶川地震之前1999—2007年之间的平均GPS速度场；黑色箭头为汶川地震之后2011—2016年间平均GPS速度场. Fig. 1 Geological settings in Eastern Bayan Har block The thick black lines represents the actve faults in the research area. The focal spheres are the earthquakes greater than M6.0 since 1976, among which the black strike-slip and thrust ones represent Jiuzhaigou earthquake and Wenchuan earthquake, with the former striking in NW-SE direction(http://www.globalcmt.org). The dark gray dots near Jiuzhaigou earthquake are aftershocks greater than M4.0(https://earthquake.usgs.gov/earthquakes/). The gray arrows and black arrows are the average GPS velocities during 1999—2007 before Wenchuan earthquake, and 2011—2016 after Wenchuan earthquake, respectively.

巴颜喀拉块体位于青藏高原中部，是青藏高原物质向东“逃逸”的重要构造单元之一.其运动模式受控于南、北边界大型走滑断裂带(邓起东等，2010).块体北边界的左旋走滑昆仑断裂带，是调整青藏高原变形的重要构造单元(Kirby et al., 2007；Meyer et al., 1998; Van der Woerd et al., 2002; 张军龙等，2012; Ren et al., 2013; 胡朝忠等, 2017)，东边界为龙门山逆冲断裂，兼有右旋走滑分量(邓起东，2010；徐锡伟等，2005；Densmore et al., 2007).块体东部现今地震活动强烈，是南北地震带中部的重要组成单元(丁国瑜，1990)，区内伴随多次历史强震活动(徐锡伟等，2008；邓起东等，2010；闻学泽等，2011).本次M7.0九寨沟地震发生在东昆仑断裂带的分支左旋走滑断裂上，断层走向呈北西-南东向.

在印度—欧亚板块碰撞及高原侧向扩展作用下，巴颜喀拉块体现今整体呈北东向运动(Avouac and Tapponnier, 1993; 张培震，2003).东部受到华南地块的阻挡，地壳呈近东-西向挤压缩短变形，造成块体东缘活动构造分布密集(徐锡伟等，2005).现今地表GPS速度观测数据显示巴颜喀拉块体整体呈NE方向运动，在邻近龙门山断裂的东部向东南方向运动(张培震，2003; Gan et al., 2007; Wu et al, 2011；Liang et al., 2013).基于中国大陆构造环境监测网络近20年GPS观测数据，解算得到了1999—2007和2011—2016期GPS速度场结果(图 1)(Wu et al, 2015；占伟等，2015；梁洪宝等，2016)；结果显示，2008年汶川地震前后巴颜喀拉块体东部地表速度场存在明显差异：岷江断裂带以南，2011—2016年平均GPS速度场东-西向分量较1999—2007年明显增加，且速度场方向往北东方向偏转；岷江断裂以北，汶川地震之后区域速度减小，速度场方向呈现小幅的北东向偏转；2011—2016年间龙门山断裂带东侧的四川盆地水平速度场较汶川地震之前明显减小，速度场方向无明显变化，见图 1.汶川地震前后区域应变和应力场如何变化？2008年汶川地震与区域变形时间差异及九寨沟地震是否存在关联？

地震的发生是地壳应力释放的过程，往往会对周围区域的形变场、应力场及断层上的应力积累和区域地震活动性有着较大的影响(Stein, 2003；Lin et al., 2004；黄禄渊等，2017).汶川地震之后，许多学者从地震静态应力触发的角度计算分析了2008年汶川地震的同震效应.例如，Parsons等(2008)和Toda等(2008)利用均匀半无限空间地壳模型解析解计算得到汶川地震使左旋走滑型的东昆仑断裂库仑应力增加，增加了其发生破裂的危险性.单斌等(2009)基于分层地壳模型的数值方法，计算得到汶川地震降低了逆冲型岷江断裂和虎牙断裂的库仑应力，但会因破裂模型的不同而有所差异.以上研究均发现汶川地震对巴颜喀拉块体东部活动断层有不同程度的影响.

GPS观测数据有效反映现今地壳形变状态(Zhu and Shi, 2011)，利用不同时段GPS速度场数据可以得到区域应变场的动态演化趋势.同时，分析地震孕育过程中地壳形变和应力变化对认识强震机理具有理论和实际意义.据此，本文基于近20年的GPS速度场观测数据，尝试从构造加载角度分析巴颜喀拉块体东部地壳形变、应力及应变能积累过程.通过建立三维高分辨率有限元地质模型，计算现今巴颜喀拉块体东部及九寨沟地震周缘地壳形变和应变能积累的趋势.进一步对比汶川地震前后两个时间段(1999—2007年和2011—2016年)的GPS速度场约束下区域形变特征，探讨汶川地震对巴颜喀拉地块东部形变及九寨沟地震发生的影响.

基于巴颜喀拉块体东缘的区域活动构造特征和九寨沟地震的震中位置，本文选取(102—104.5°E, 32.5—34°N)为主要研究区域，覆盖巴颜喀拉块体东部的龙日坝断裂、岷江断裂、虎牙断裂和塔藏断裂等大型活动断裂带.其中，北西-南东向的塔藏断裂为以左旋滑动为主；近南北向的虎牙断裂和岷江断裂逆冲滑动为主，兼具有走滑运动分量(邓起东等，1994；马保起等，2005；周荣军等，2006)；龙日坝断裂带近期呈现右旋走滑的运动模式(徐锡伟等，2008).本次研究的深度范围为上地壳脆性孕震层，深度为25 km，暂未考虑地壳垂向精细分层结构，图 2为研究区域的三维有限元模型.模型采用四面体网格剖分，平均网格大小约为4 km，总节点数为61, 918，单元数为327, 290，三维模型网格划分见图 2b.由于万年以内的短时间尺度地球介质变形主要表现为弹性特征(Turcotte and Schubert, 2014)，因此，分析现今几十年内的地壳变形时可将地壳视为弹性介质，其动力学控制方程采用线弹性本构关系.

图 2

Fig. 2

图 2 三维有限元计算模型和计算边界条件 (a)区域实测GPS速度场时空分布，白色虚线框为研究区位置；(b)三维有限元网格模型；(c) 1999—2007年平均GPS速度场约束得到的模型边界条件；(d) 2011—2016年间GPS速度场约束得到的模型边界条件. Fig. 2 The 3-D numerical model and boundary conditions (a) GPS velocities in the research area, the white dotted rectangle is the location of the 3-D model; (b) The 3-D model with meshes; (c—d) The boundary conditions contrained by GPS velocities during 1999—2007, and 2011—2016 respectively.

根据已有的巴颜喀拉块体东部地壳S波速度结构及地壳厚度分布(王椿镛等，2008；Zhang et al., 2011)及Crust1.0数据(http://igppweb.ucsd.edu/~gabi/rem.html)，获得了巴颜喀拉块体上地壳的弹性力学性质作为模型的力学参数，见表 1.相对周围岩石，断层力学性质较弱(Ben-zion and Sammis, 2003)，其弹性力学参数低于周围基岩两个数量级.模型中断层的走向参考实际野外测量，视为具有2 km宽的软弱带.模型分别采用中国大陆构造环境监测网络1999—2007年和2011—2016年两个时期的GPS平均速度场作为边界约束.利用最小二乘配置方法分别将两个时间段内离散的GPS速度场插值到模型边界上(Wu et al., 2017)，两个模型边界速度约束在方向和大小上均有明显的差异，见图 2(c—d).由于GPS反映地球表层水平运动特征，模型边界仅为水平方向的速度约束，边界速度不随深度变化.

表 1 Table 1

表 1 模型力学参数 Table 1 Mechanical parameters of the model 弹性模量/ Pa 泊松比 深(厚)度 /km 密度 /(kg·m-3) 上地壳 9.0×1010 0.26 25 2750 断裂带 9.0×109 0.26 2 2750

表 1 模型力学参数 Table 1 Mechanical parameters of the model

图 3为不同时段GPS观测值约束下计算得到的区域地壳速度场分布，计算得到的速度场与实际GPS观测速度场方向和大小基本一致, 受区域活动断裂带间断错动等因素的影响，研究区个别观测点的计算结果与实际观测稍有差异.1999—2007年GPS速度约束下，区域地表速度场计算值与观测值的均方根误差东-西向分量为0.38 mm·a^-1，南-北向分量为0.41 mm·a^-1；2011—2016年GPS速度约束下的速度场计算值与观测值的均方根误差，东-西向分量为0.577 mm·a^-1，南-北向分量为0.32 mm·a^-1；地表速度场的计算值与观测值吻合较好.2008年汶川地震发生之前的1999—2007年期间，巴颜喀拉块体东部的地表速度南北向分量显著，由西往东速度的东-西向分量呈递减趋势.而汶川地震之后的2011—2016期间，区域地表速度场相对之前发生明显调整，速度场向南运动分量整体减小，自西往东东-西向速度分量差异增加.综上，汶川地震前后两个时段研究区速度场时空差异明显.

图 3

Fig. 3

图 3 研究区域水平速度场分布 (a) 1999—2007年间GPS约束下区域速度场分布；(b) 2011—2016年间GPS约束下区域速度场分布. Fig. 3 Distribution of horizontal velocities in the research area (a) Horizontal velocities derived from GPS velocities during 1999—2007; (b) Horizontal velocities derived from GPS velocities during 2011—2016.

通过计算地壳应变率可以更直观地分析地壳的变形状态，图 4给出了计算得到的巴颜喀拉块体东部地壳水平主应变率分布.其中，汶川地震之前的1999—2007年间区域地壳应变率计算结果与Wu等(2015)和陈长云等(2013)利用GPS观测数据推算的水平主应变率分布基本一致.由于断裂带具有较弱的力学性质，在构造驱动下较周围基岩容易变形，导致模型中断裂带上存在较大的速度和应变率，例如模型中塔藏断裂带内的计算节点上明显的应变率“奇点”异常现象.汶川地震之前，在现今构造驱动作用下，巴颜喀拉块体东缘主要处于北西-南东向拉张和北东-南西向挤压变形状态.九寨沟地震南侧区域，水平主应变率方向稍有偏离，特别是岷江断裂和虎牙断裂之间的区域，水平主应变率调整为北西-南东向挤压和北东-南西向拉张；在虎牙断裂东侧，北东-南西向挤压变形占主导，见图 4a.汶川地震之后，区域地壳同样以北西-南东的拉张和北东-南西的挤压变形为主，水平应变率相对于汶川地震之前明显增加；九寨沟地震以南的区域，应变调整为北西-南东的拉张和北东-南西的挤压变形；白龙江断裂周缘北西-南东向的挤压变形显著.综上，汶川地震前后两个时间段内地壳水平应变率主轴方向整体近似，说明汶川地震前后巴颜喀拉块体东部变形趋势是一致的.汶川地震前，本次九寨沟地震震中位置水平主压应变率约为3.0×10^-8·a^-1，水平主张速率约为2.0×10^-8·a^-1；汶川地震后，对应区域的水平主压应变率增至8×10^-8·a^-1，水平主张应变率为5×10^-8·a^-1，其水平主应变率增至汶川地震前约1.5~3.0倍.龙日坝断裂带以东，汶川地震前后区域水平主应变率方向调整较为明显.由此可知，汶川地震之后该区域地壳变形加剧，水平应变率主轴方向也有所调整.由于模型不考虑地壳分层耦合作用，且边界约束未考虑随深度变化，计算得到的地壳应变率不随深度变化.

图 4

Fig. 4

图 4 研究区域地壳孕震层水平主应变率场分布 (a) 1999—2007年间GPS约束下区域10 km深处水平应变率分布；(b) 2011—2016年间GPS约束下区域10 km深处水平应变率分布. Fig. 4 Distribution of horizontal strain in the crust of the research area (a) Horizontal strain at the depth of 10 km in the crust contrained by GPS velocities during 1999—2007; (b) Horizontal strain at the depth of 10 km in the crust contrained by GPS velocities during 2011—2016.Red and black lines represent extension and compression, respectively.

地震孕育过程中构造应力不断积累，同时地壳应变能也持续增加.因此，地壳应变能积累的定量计算也是分析地震危险性的重要手段.GPS观测速度场约束下的三维地球动力学模拟，可以计算得到现今地壳应变和应力的变化率，在此基础上计算区域地壳应变能密度变化率中应力和应变变化引起的部分贡献δN，在未知地壳绝对应力的情况下δN可以表征部分应变能密度变化率，计算见公式(1).图 5为汶川地震前后两个时期GPS观测速度场约束下，区域地壳应变能密度变化率相关量δN的分布.由图可知，断裂带具有较低的弹性模量，存在较大的应变量，总的效应使断裂带的应变能密度变化率相关量高于周围地块，这也与大地震主要沿断裂带发生的观测事实是吻合的.但是，特别应该注意的是，汶川地震前后区域应变能密度变化率相关量δN的变化.根据计算结果可知，汶川地震之后区域δN整体呈现明显增加，达到原来的3.0倍左右，汶川地震显著改变了区域构造运动状态，大大增加了应变能密度变化率相关量δN，这可能不仅有助于九寨沟地震的发生，而且对研究区域其他部位的地震发生的可能性也有影响，应引起关注.

图 5

Fig. 5

图 5 数值计算得出的地壳孕震层弹性应变能密度变化率相关量δN分布 (a)汶川地震之前区域10km深处应变能密度变化率相关量δN分布; (b)汶川地震之后区域10km深处应变能密度变化率相关量δN分布. Fig. 5 The related rate of crustal strain energy density δN in the research area derived from numerical modeling (a) Distribution of δN before the Wenchuan earthquake at the depth of 10 km, (b) Distribution of δN after the Wenchuan earthquake at the depth of 10 km.

(1)

根据区域水平应变率计算结果可知，巴颜喀拉块体东部地壳水平北东-南西向主压应变与北西-南东向的主拉应变大小相近，无明显占主导的拉张或挤压变形.结合区域构造背景和应力积累结果可以更直观地分析地壳构造变形趋势和活动断层上地震潜在危险性.图 6为研究区GPS速度场约束下得到的反映构造应力积累状态的震源球及断层上库仑应力变化率结果，其中各个震源球的P、T和N轴为对应区域地壳应力增量主轴的最大、最小和中间主应力方向.汶川地震前后两个时期内，反映应力增量的震源球都以走滑型地震形式为主，说明现今区域应力积累整体上促进了巴颜喀拉块体东缘北东-南西方向右旋走滑断层(如龙日坝断裂)和北西-南东左旋走滑断层(如塔藏断裂和白龙江断裂)上地震的孕育过程.而汶川地震之前，白龙江断裂和塔藏断裂之间，受东侧边界显著的东向运动影响，近北西-南东向拉应力变化占主导，该阶段更有利于两断层间正断层型地震的孕育；汶川地震之后，对应区域地壳形变调整为利于走滑型地震的孕育.虎牙断裂带东侧区域北东-南西向挤压变形占主导，利于逆冲断层型地震的孕育.

图 6

Fig. 6

图 6 区域应力累计趋势分析 (a)汶川地震之前区域应力变化率趋势分布；(b)汶川地震之后区域应力变化率趋势分布；图中红色震源球为区域实际地震震源机制，黑色震源球代表计算得到的地壳主应力年变化率趋势，反映地壳应力孕育地震的趋势；断层上彩色条带代表GPS速度约束下断层上库仑应力变化率. Fig. 6 The analysis on stress accumulation in the research area (a) The trend of crustal stress accumulation before Wenchuan earthquake; (b) The trend of crustal stress accumulation after Wenchuan earthquake; The red focal spheres represent the real earthquake; Black focal spheres are derived from crustal stress changes in our numerical models, representing the trend for earthquae nucleation process; The colored lines represent the rate of Coulomb stress changes on the active faults contrained by GPS velocities.

库仑应力变化是有效分析地壳应力变化与地震危险性关系的重要手段.利用计算得到的地壳应力变化，参考研究区活动断层产状和历史地震滑动信息，见表 2，本文计算了GPS速度场约束下活动断层上库仑应力变化率，见图 6.其中，九寨沟地震的发震断层产状等参数主要参考本次地震的破裂反演及余震分布结果.根据结果可知，巴颜喀拉块体东部现今构造变形导致主要走滑断层上库仑应力整体增加.1999—2007年间，左旋走滑型白龙江断裂库仑应力增长率约为1.0 kPa·a^-1，左旋走滑型塔藏断裂和九寨沟地震发震断裂库仑应力增长率约为0.7 kPa·a^-1，右旋走滑型龙日坝断裂上平均库仑应力增长率约为0.5 kPa·a^-1；逆冲型的岷江断裂北段库仑应力降低，中南端呈小幅增加趋势；逆冲型虎牙断裂带北段库仑应力增加约为0.4 kPa·a^-1，其南段库仑应力呈小幅降低.2011—2016年，可能受汶川地震震后效应的影响，研究区库仑应力变化幅度增加，龙日坝断裂和塔藏断裂带上，库仑应力增长率约为汶川地震之前的1.5倍；而虎牙断裂和岷江断裂上库仑应力整体降低；九寨沟地震发震断裂上库仑应力增长率平均约为1.0 kPa·a^-1，最大值可达1.2 kPa·a^-1.因此，从库仑应力变化角度进一步验证，巴颜喀拉地块东部构造驱动作用有利于孕育北西-南东向的左旋走滑和北东-南西向的右旋走滑型地震.通过对比汶川地震前后两期地壳应力及库仑应力变化结果，我们认为汶川地震的发生在一定程度上增强了区域走滑型断裂带上地震的孕育过程.

表 2 Table 2

表 2 研究区域断层形态及计算参数 Table 2 Geometry and slip parameters of active faults 断层名称 走向 倾向 滑动角 断层摩擦系数 参考文献 塔藏断裂 100~145 89 0 0.6 (He and Chéry, 2008) 龙日坝断裂 205~229 60 135 0.6 (Kirby et al., 2000; 邓起东等，2002; 张培震等，2003) 白龙江断裂 290 68~85 45 0.6 岷江断裂 180 45~56 45 0.6 虎牙断裂 150 75 45 0.6 九寨沟发震断裂 326 62 -15 0.6

表 2 研究区域断层形态及计算参数 Table 2 Geometry and slip parameters of active faults

本文基于巴颜喀拉块体东部强震活跃区的构造背景，利用中国大陆构造环境监测网络1999—2007年和2011—2016年两个时段GPS速度场观测分布，通过建立三维有限元地球动力学模型，计算了汶川地震前后两个时期GPS约束下九寨沟地震周缘地壳形变特征、应变能和应力积累趋势.

计算结果显示在印度-欧亚板块碰撞驱动青藏高原物质向东逃逸的前提下，现今巴颜喀拉块体东部呈东南向运动特征，地壳主要为北西-南东向水平拉张和北东-南西向的水平挤压变形.现今构造驱动下地壳应力积累过程有助于区内北西-南东向左旋走滑或北东-南西向右旋走滑型断裂上地震的孕育，而对逆冲型岷江和虎牙断裂的地震孕育的驱动作用相对较小.由此可知，巴颜喀拉块体东部现今的变形模式依然以左旋走滑构造占主导.而虎牙断裂以东的区域，地壳水平挤压变形占主导，有利于逆冲型断裂构造发育.

汶川地震之后，巴颜喀拉块体东部地壳速度场、形变场等发生明显的调整，主要表现为: (1) 汶川地震之后，区域地壳速度场方向和大小均较之前发生调整, 速度向南运动分量减小；(2) 从形变角度分析，汶川地震之后区域水平应变率增加，局部地区主应变率方向发生变化，应变能密度变化率相关量也明显增加，说明区域应变能积累增强；(3) 汶川地震之后，地壳应力增长率明显增加，走滑型断裂上库仑应力变化率大于汶川地震之前.

前人围绕汶川地震静态应力触发相关研究表明，汶川地震同震和震后效应会引起九寨沟地震震中及周缘地区库仑应力的明显变化，调整区域地壳形变及应力积累过程(Parsons et al., 2008；Toda et al., 2008；单斌等，2009；程佳等，2011).根据近20年该区域的GPS速度场观测数据，利用有限元模拟方法，从库仑应力角度分析得出汶川地震之后九寨沟地震发震断裂上库仑应力增长率约为1.2 kPa·a^-1，而汶川地震之前构造驱动该断层库仑应力增加约为0.7 kPa·a^-1.综合对比汶川地震前后区域地壳形变、应变能、应力积累结果，本文初步得出现今构造环境有利于九寨沟地震发震断层的孕震过程，汶川地震在一定程度上加速了本次地震的孕育过程.

两位审稿专家对本文提出的建设性修改意见，中国地震局第一监测中心武艳强研究员在文章撰写过程中提供的帮助，中国大陆构造环境监测网络提供GPS观测数据、中国地震局第一监测中心GPS工作组提供1999—2007年和2011—2016年两个时段GPS速度场结果，在此一并致谢.