继2008年5月12日汶川M_S8.0级地震后，时隔不到5年，2013年4月20日四川芦山县发生M_S7.0级地震.两次强震相距约90 km，以逆冲破裂为主，分别发生于龙门山断裂带中北段和南段，是主余型还是孤立型学界尚有争论，若为孤立型地震则重复地震周期与现有研究相差较大(张培震等，2008; 柳畅等，2012).这促使我们重新认识龙门山断裂带现今运动、受力状态及应力积累.在高原物质侧向挤出与华南地块阻挡作用的大背景下，川滇及周边形成巴颜喀拉、滇东、菱形块体等被断层切割、围限的地块体系(骆佳骥等，2012)，不同地块既区别又协调的相对运动共同构成区内主要断裂带与地震活动带——构造块体边界带——运动与变形的动力源.研究对比各断裂带运动、受力和应力积累及地震活动特征则有助于这一问题的认识.考虑到断裂带空间展布复杂多变，建立较为客观的模型和模拟计算是较好途径，而川滇地区近30年多期GPS结果为实现数值模拟提供了很好的约束条件.

本文构建川滇地区二维有限元接触模型，利用中国大陆网络工程项目1999—2007年GPS结果和“块体加载”有限元方法(刘峡等，2013)，模拟龙门山、鲜水河、小江、红河等川滇地区主要断裂带形变运动.对比、分析其变形特征、应力状态及与区域强震的关系，研究汶川地震、芦山地震的力学背景.

采用文献(李延兴等，2009)ITRF框架下中国大陆1999—2007年站速度结果，通过坐标旋转，得到区域无旋转坐标框架下GPS速度场，见图 1.该结果反映了川滇地区在汶川、芦山地震之前较长时期内的水平运动状态，可以清晰观察到断层两侧和各构造块体之间的GPS速度变化，方便于下文模拟结果与实测结果的对比和残差分析.

图 1 1999—2007年GPS站速度与有限元模拟结果(a)GPS速度和有限元模拟速度对比(图中蓝色框为有限元模型范围);(b)模拟计算与GPS观测结果的残差分布.Fig. 1 1999—2007 GPS speed and the simulated results of FEA(a)The GPS speed and the simulated results(blue frame denotes the FEA model area);(b)The residual between the simulated results and GPS speed.

在96°E—109°E、20°N—37°N范围内构建川滇地区地壳运动二维有限元接触模型，涵盖巴颜喀拉块体东部、滇东块体、滇西南块体和川滇菱形块体.见图 2.纳入模型的断层共43条，包括：龙门山中央、龙门山后山(龙门山中央断裂带与后山断裂带并为一条)、龙门山前山、岷江、虎牙、马尔康、文县、甘孜—玉树、东昆仑、鲜水河、小江(两个分支断裂合并为一条)、则木河、红河、小金河、大凉山、理塘—义敦—德巫、南华—楚雄—建水、楚雄—建水—通海、维西—巍山、怒江、马边—盐津、磨盘山—绿汁江、字嘎寺—德钦、镇远—普洱、大扁担山—景洪、大盈江、德钦—中缅、新津—德阳、龙泉山西缘断裂带等.采用ANSYS有限元软件，断裂带由二维接触单元 conta175、targe169组成，其他区域由面单元plane182 组成.采用静态弹性接触分析方法.非接触单元的杨氏模量和泊松比依据P波速度给定，分别为7.0×10¹⁰Pa、0.25(朱守彪等，2004)，四川盆地的杨氏模量高于设定值0.7～1.0倍(郑勇等，2006).断层接触面的摩擦系数统一设定为0.3(试算结果显示在本文约束条件下，在0.3～0.6范围内改变摩擦系数，应力分布变化低于5％，断层运动速率变化低于 15％，因此在设定断层摩擦系数时进行了简化处理).

图 2 川滇地区有限元二维接触模型(a)有限元模型构架；(b)有限元模型分网.Fig. 2 2D contact FEA model of Sichuan-Yunnan Region(a)The model region;(b)The mesh of the FEA model.

构造块体的边界带一般由大型断裂带组成，块体间的差异运动是其运动与变形的主要动力源.为模拟边界带的断层运动，作者提出“块体加载”有限元数值模拟方法(刘峡等，2013).该方法利用多个GPS结果提取块体运动参数并在有限元模型上进行加载计算，物理意义明确，能够避免实测数据的随 机变化引起的计算结果畸变，确保区域断裂带体系的形变运动相互协调和计算结果稳健可靠.按照该方法，本文第一步利用文献(李延兴等，2007)中公式(1)和本文GPS数据，得到川滇菱形块体、巴颜喀拉、滇东、滇西南块体的运动参数与形变参数见表 1(其中形变参数略)；第二步，根据上述参数计算得到有限元模型中相应块体区内的单元节点(位于图 2a阴影区内的节点)的速度值，形成约束条件；第三步，加载计算，模拟断裂带(位于图 2a空白区内的断层)形变运动.

表 1(Table 1)

表 1 川滇地区主要构造块体的运动参数 Table 1 The movement parameters of the main blocks in Sichuan-Yunnan Region 构造块体 ωx(10-9/a) ωy(10-9/a) ωz(10-9/a) 块体模型 rms(mm·a-1) 数值 rms 数值 rms 数值 rms 巴颜喀拉块体东部 1.1367 0.1888 -0.5549 0.7236 -0.4589 0.5009 1.3235 滇东块体 -1.444 0.3895 8.442 1.6162 4.5349 0.8539 0.9161 川滇菱形块体 -0.4523 0.2576 -2.9221 1.2714 -0.8339 0.6955 1.0412 滇西南块体 0.3965 0.1922 -6.2498 1.0268 -3.8133 0.4932 1.4933

表 1 川滇地区主要构造块体的运动参数 Table 1 The movement parameters of the main blocks in Sichuan-Yunnan Region

图 1a为模拟速度与GPS实测结果的对比，红色箭头显示的是模拟结果，其中大部分位于图 2a中的阴影区，为加载约束，位于断裂带附近的(图 2a的空白区)是实算结果.图 1b为模拟结果与GPS结果的残差分布.总体上看，模拟速度场的整体分布形态与实测GPS结果保持一致，表明加载计算客观反映了研究区断层运动的动力学条件.

依据单元接触对的位移可以得到断层错动速率，包括平行于断层的水平走滑分量(右旋/左旋，右旋为正)和垂直断层面的水平拉张分量(拉张/挤压，拉张为正)，见表 2和图 3.

图 3 川滇地壳主要断裂带运动速率模拟结果(a)走滑分量(右旋/左旋)；(b)拉张分量(拉张/挤压).Fig. 3 The simulated movement speed of the main faults in Sichuan-Yunnan Region(a)Slip component(right-lateral/left-lateral);(b)Tension component(tensile/compressive).

根据模拟结果，龙门山断裂带(模型将后山、中央和前山断裂带合并为两条断层)为右旋、挤压，右 旋合计不超过1 mm·a^-1，挤压不超过0.5 mm·a^-1. 右旋量沿断层分布不均匀，在汶川地震、芦山地震震源区附近右旋量很小，呈现一定程度闭锁，本文认为这主要受模型中断裂带在空间上蜿蜒展布的影响，有待进一步考证.地质研究显示(徐锡伟等，2003)该断裂带为右旋逆冲，水平速率约1 mm·a^-1.杜方等(2009)指出汶川地震前龙门山断裂带呈右旋、挤压但量值很小，均低于1 mm·a^-1；主要形变弥散于与其西侧巴颜喀拉块体之间约230 km范围内，水平压缩、右旋剪切和垂直隆升量分别达到1.3、2.6、 2～3 mm·a^-1.基于GPS数据和位错计算等，乔学军等(2004)给出龙门山水平运动速率约为1.7 mm·a^-1，王阎昭等(2008)给出 呈左旋挤压运动，速率分别为-0.5、-1.1 mm·a^-1. 本文结果与徐锡伟等(2003); 杜方等(2009)是一致的，反映出龙门山断裂带在汶川、芦山地震前较长时期内保持低变形.

表 2(Table 2)

表 2 川滇地区主要断裂带运动速率模拟结果及与地质数据比较Table 2 The simulated movement speed of the main faults in Sichuan-Yunnan Region compared with geologic data 断裂带 走滑分量(右旋/左旋)/mm·a-1 拉张分量(拉张/挤压)/mm·a-1 地质研究水平速率(徐锡伟，2013)/ mm·a-1 龙门山后山/中央 0~0.57 -0.10~-0.17 右旋逆，1 龙门山前山 0~0.38 -0.07~-0.19 鲜水河 -5.71~-8.31 0.58~4.64 左旋，14±2 则木河 -6.02~-8.09 0.20~3.70 左旋正，6.4±0.6 小江 -5.25~-6.62 0.04~4.41 左旋正，10±2 大凉山 -0.48~-2.31 -0.01~0 左旋逆，3.3±0.7 南华—楚雄—建水 2.18~4.22 0~1.33 右旋逆，3.5 红河北段 2.50~3.45 0~0.72 右旋正，3.3 红河中段 0.62~1.95 0~0.02 — 红河南段 2.12~3.93 0~1.61 右旋逆，3.5±1.5 维希—巍山 0.85-1.73 1.73~2.75 — 注：右旋、拉张为正.

表 2 川滇地区主要断裂带运动速率模拟结果及与地质数据比较Table 2 The simulated movement speed of the main faults in Sichuan-Yunnan Region compared with geologic data

鲜水河、则木河、小江和大凉山断裂带共同构成菱形块体东边界，本文结果显示鲜水河、则木河、小 江断裂带为左旋、拉张，大凉山断裂带为左旋挤压，后者速率明显低于前三条断裂带.前三条以鲜水河断裂带左旋量最大，向南渐小，变化比较均匀，变化范围是-5.25～-8.09 mm·a^-1，但拉张很不均匀，最大拉张4.64 mm·a^-1，这应是块体相对运动与块体边界特有的几何形状共同作用的结果.大凉山断裂带左旋速率约-2.31 mm·a^-1，挤压速率约-0.01 mm·a^-1.地质研究显示(徐锡伟等，2003)鲜水河、则木河、小江断裂带均为左旋和正断活动，左旋速率依次为14±2、6.4±0.6、10±2 mm·a^-1，鲜水河东南段正断速率3.2±0.7 mm·a^-1.大凉山为左旋逆冲，左旋速率3.3±0.7 mm·a^-1.模拟结果给出的断层活动性质与地质结果一致，但鲜水河、小江断裂带走滑速率偏低，以鲜水河断裂带最显著.实际上从现有GPS、地震矩计算结果看本文给出的速率也相对较低：文献中( Shen Z K et al.2005; 王阎昭等，2008; 乔学军等，2004; 申重阳等，2002)鲜水河走滑速率分别为8～10、8.9～17.1、约10.4、约30 mm·a^-1，地震矩反演结果(孙建中等，1994)为10.9 mm·a^-1.本文分析原因有二，其一，在块体加载计算中，没有考虑菱形块体内部川北、滇中两个次级块体，而这两块体的差异运动对鲜水河、则木河、小江断裂带的速度估值是有影响的；其二，基于GPS数据的位错反演仅考虑断层运动对速度场的影响，而本文的块体加载则兼顾了块体本身的变形，因此本文的断层错动量势必偏低.

模拟计算给出红河断裂带为右旋、拉张，部分区段略呈挤压，运动速率不均匀，南、北段高于中段，右旋、拉张速率变化范围0.62~3.93 mm·a^-1、0～1.61 mm·a^-1.与红河断裂带相比，西侧、北侧的维西—巍山、南华—楚雄—建水断裂带运动速率更高，显示其现今活动更突出.上述结果与地质结果比较接近(徐锡伟等，2003)，但值得注意的是，地质研究给出南华—楚雄—建水断裂带为逆冲，而本文为拉张，基于GPS数据的位错反演(王阎昭等，2008)也给出拉张变形约0.9 mm·a^-1.因此，该断裂带现今运动性质可能与历史情形相反.

应力场模拟结果既包含块体内部(即加载区)应力分布也包含断层两侧应力分布.模拟采用的是速度值加载，得到的应力结果是以年为单位的应力增量(张应力为正).结果显示，研究区主压应力方向自北而南由NEE转为NS、呈扇形分布，与现有研究(崔效锋，2006)一致.与连续变形的块体区相比，断层沿线的应力大小、方向变化突出.图 4为断层两侧应力分布以及整个区域第一、第二主应力和最大剪应力等值线分布.

图 4 川滇地区应力场模拟结果(a)沿断裂带应力分布；(b)区域主张应力等值线分布；(c)区域主压应力等值线分布；(d)区域最大剪切应力等值线分布.Fig. 4 The simulated stress distribution in Sichuan-Yunnan Region (a)The stress along the main faults(red narrows denote the principal tensile stress，black narrows denote the principal compressive stress);(b)The contour map of the principal tensile stress;(c)The contour map of the principal compressive stress;(d)The contour map of the maximum shear stress.

根据图 4a，龙门山断裂带主压应力大于主张应力，主压方向与断层走向近于垂直，形成一条很宽的挤压带；鲜水河、则木河、大凉山断裂带主压应力大于主张应力，主压方向EW-NWW，与断层呈40°～45°交角；小江断裂带北、中段挤压应力较大，南段拉张更突出，主压方向NW-NNW；维希—巍山、南华—楚雄—建水断裂带以拉张为主，主张方向与断裂带接近垂直；红河断裂带北段拉张稍强、中南段挤压稍强，但量值均小.钟继茂等(2006)给出的川滇地区中强震震源机制解的P轴分布与本文上述结果基本一致.为分析断层面受力，将断层两侧节点的应 力张量提取出来，剔除局部极大、极小点(这些极值 点是断层折拐或断层衔接处的应力集中引起的，占总数据量的7％以下)，取平均值投影至断层面上，得到断层面正应力、剪切应力见表 3.根据表 3，龙门山断裂带正压应力最大，为-333.74Pa，切向应力为0；南华—楚雄—建水的正张应力最大，为340.33 Pa，切向应力仅-0.40 Pa；鲜水河的切向应 力最大，为-251.79 Pa，其次是小江，为-197.84 Pa. 尽管鲜水河断裂带两侧主压应力均值达到-455.03 Pa，但由于应力方向与断层相交，致使断层面正压应力仅-155.23 Pa，不到龙门山断裂带的二分之一.

图 4(b、c、d)显示，研究区应力场呈非均匀分布，巴颜喀拉块体挤压、拉张均较高，其中龙门山断裂带附近挤压更为突出，其最大主张应力亦为负值；菱形块体东北部与西南部明显不同，分别呈现强烈挤压、强烈拉张；滇西南块体挤压和拉张量持平；四川盆地应力量值很小，说明该区域以刚性运动为主，变形很小.在块体衔接区，即龙门山、鲜水河、红河 等断裂带附近区域，应力变化剧烈，等值线沿断层 密集分布，形成梯度带，这一特征在图 4d中最为突出.在断层拐折区与交汇区出现局部应力集中区，应力值高于周边5～10倍，如小江断裂带中南部、玉树 与鲜水河断裂带交汇区、小江与红河断裂带交汇区等.

表 3(Table 3)

表 3 川滇地区断裂带应力分布的统计结果(张应力为正)Table 3 The pressure on the surfaces of the main faults in Sichuan-Yunnan Region(tensile pressure is positive) 断裂带 第一主应力平均值/Pa 第二主应力平均值/Pa 第二主应力与 断层平均夹角/(°) 断层面正应力/Pa 断层面剪切应力/Pa 龙门山 -26.43 -333.74 90 -333.74 0 鲜水河 56.23 -455.03 40 -155.23 -251.79 则木河 26.26 -218.52 45 -96.23 -122.39 小江 167.50 -234.22 40 1.34 -197.84 大凉山 -26.23 -226.23 38 -102.12 -97.052 南华—楚雄—建水 340.33 -157.88 0 340.33 -0.40 红河 132.627 -210.69 45 -39.17 -171.66 维希—巍山 310.84 -68.03 5 307.91 -33.18

表 3 川滇地区断裂带应力分布的统计结果(张应力为正)Table 3 The pressure on the surfaces of the main faults in Sichuan-Yunnan Region(tensile pressure is positive)

本文结果显示，以鲜水河、龙门山、小江、红河等断裂带为代表的川滇地区主要断裂带现今运动彼此差异很大.如鲜水河断裂带为左旋走滑，龙门山断裂带为右旋挤压，两者速度相差约10倍；红河断裂带北段、南段的运动速度高于中段.又如，位于川滇菱形块体内侧的断层运动和形变普遍高于外侧断层，如南华—楚雄—建水断裂带高于红河断裂带，则木河、小江断裂带高于大凉山断裂带.这些差异的总体特征符合现有研究结果.表明断裂带运动受块体差异运动所驱动，菱形块体的运动与变形较周边块体更为突出，导致其内侧断层形变运动较外侧更显著，由此形成复杂而彼此协调的断层运动体系.本文结果还显示，同一条断裂带不同区段的运动也不尽相同，如龙门山断裂带在汶川、芦山地震震源区附近的错动速率小于其他区段.表明在该断层运动体系之下，单一断层运动速率的不均匀性明显受其空间展布的影响，因此，在判断能否形成强震孕育区的过程中，需要将断层空间展布变化与动力学运动相结合，考察断层运动的不均匀特征和容易产生“运动亏损”的区段.

利用GPS数据获取断层现今运动速率并结合历史地震、地质数据，可以得到强震复发周期等具有重要价值的信息，然而不同研究给出的结果差异较大，从本文速率结果与已有结果(乔学军等，2004; 王阎昭等，2008;申重阳等，2002)的对比情况看亦是如此，显示出方法与模型的依赖性，需谨慎对待之.本文认为，作为高原物质侧向挤出通道，川滇地区地形、地壳厚度变化剧烈，地幔拖曳(朱守彪等，2004)可能造成挤出速度在不同深度上存在差异，导致其内部变形强烈、与刚性板块运动差异巨大.因此，考察研究区的运动与变形，不能仅考虑大型断裂带对区内次一级块体的控制作用，块体本身变形也不应忽视.本文在“块体加载”过程中，考虑了块体本身的形变，这是本文结果与位错反演结果(乔学军等，2004;王阎昭等，2008; 申重阳等，2002)相比偏低的重要原因，故可将本文结果视为断层速率估值的下限.

从本文应力结果反映的块体受力来看，以龙门山、小江断裂带为界，西侧各块体应力值高于东侧块体.巴颜喀拉块体东移在四川盆地的阻挡之下，受到强烈挤压，在红河断裂带以南，高原物质向南流动趋快，导致滇西南块体强烈挤压与拉张并存.在菱形块体内部，物质流动由东转南，受华南块体阻挡和菱形块体周边断裂带展布的围限，块体东北部挤压较强、西南部拉张较强.此外，与断层速率彼此相差较大迥然不同的是，本文给出的沿断裂带的主张或主压应力值较为接近，变化范围基本在200～450 Pa范围内(除断层交汇、拐折等部位)，反映了相对统一的力学背景.断层面的应力状态与沿断层发生的强震破裂方式基本一致，揭示地震孕育、发生与块体之间的相互作用密切相关，如鲜水河断裂带的切向应力在全区最高，且是法向应力的1.6倍，其强震主要为走滑型；而龙门山断裂带的法向应力在全区最高，但切向应力几乎是零，汶川地震、芦山地震均以逆冲破裂为主.

根据模拟结果，尽管龙门山断裂带运动速率远低于其他断裂带，仅约1 mm·a^-1，但断层面的应力与其他断裂带相当，甚至更高，表明在龙门山断裂带貌似平静的背后，其能量积累速度堪与鲜水河、小江等近期强震活跃带相匹敌.作为本文结果的印证，赵炜喆(2008)等利用地震资料发现1987年以来龙门山断裂带的累积Benioff应变与鲜水河、小江等断裂带相当.地震记录(石绍先等，2003)显示1900年以来鲜水河断裂带共发生3次7级以上地震，根据表 3，在每次强震的间隔期，断层面的切 向、法向应力积累平均约为9.45×10³ Pa、 约5.84×10³ Pa.若设龙门山断裂带的重复地震周期约3000年，其间的挤压应力积累约10⁶ Pa，远高于鲜水河断裂带的震间积累量，与岩石破裂强度相当.虽然这样的对比仅基于整条断裂带的平均量，没有考虑断层沿线的应力差异，也未顾及不同破裂类型所需的法向应力与切向应力的差别，但无论如何这样的差异还是惊人的.考虑到上述震间应力积累量是基于GPS结果得出的，并不适于千年时间尺度，因此最有可能的情况是目前龙门山断裂带的应力积累速度远远高于其千年时间尺度上的平均值.

研究表明川滇地区主要断裂带强震复发趋于随机、丛集行为，并不具有良好的周期性(易桂喜等，2012).结合上文并考虑到龙门山断裂带在汶川地震前的低形变状态，本文推测目前巴颜喀拉块体东部和龙门山断裂带的构造挤压应力呈快速增长，并已有相当长的时间积累，导致其实际应力处于高水平，并进入地震活跃时期.在此期间强震呈丛集性发生，故而有汶川、芦山地震在很短时间内重复发生的现象.

本文基于1999—2007年GPS数据和有限元“块体加载”方法实现川滇地区断层运动数值模拟，认识如下：

第一，作为衔接不同地块的构造边界带，研究区主要断裂带——龙门山、鲜水河、红河、小江断裂带等——形变运动及受力特征主要由块体间的差异运动所决定，形成复杂而彼此协调的断层运动体系.1999—2007年期间，龙门山断裂带右旋错动不超过1 mm·a^-1，挤压不超过0.5 mm·a^-1，但断层面法向挤压应力积累为-333.74 Pa·a^-1，显著高于其他断裂带，并且主压应力轴垂直断裂带走向，形成一条很宽的挤压带.这一低变形和高应力增长现象提示我们，现今龙门山断裂带和巴颜喀拉块体东部的实际挤压应力可能已处于较高水平，并诱发强震活跃，此期间强震呈丛集性发生，出现汶川、芦山地震接连发生的现象.

其二，模拟计算显示在汶川、芦山震源区附近，断层运动速率明显低于其他区段，尽管有待进一步证实，但提示我们，断层速率的不均匀性明显受其空间展布的影响，因此，在判断能否形成强震孕育区的过程中，需要将断层的空间变化与动力学运动相结合，考察断层运动的不均匀特征和容易产生“运动亏损”的区段.另外，研究区地壳变形不仅受大型断裂带所控制，块体本身变形也不容忽视，故而利用GPS数据估算断裂带运动速率需兼顾这两方面因素，本文的速率结果可视为现今断裂带形变运动的下限.

需要指出的是，本文的模拟计算未考虑各断层地下延伸状态和断层部分锁定的形变影响.汶川地震的主、余震研究显示破裂并不完全与龙门山断裂带走向一致，在理县NW向分支也有分布，揭示震源区地下结构复杂(易桂喜等，2012; Zheng et al., 2010; Zheng et al., 2009).震源区局部特殊构造可能加快应力积累速度并诱发地震.此外，川滇地区地形差异很大，如龙门山断裂带东西两侧重力侧滑以及岩石层横向不均匀性引起的蠕变“流”也可对断层形成逐年变化的影响.在今后工作中，通过建立3D数值模型以及采用分布更为密集的GPS测站结果，相信可以取得更好的结果.