引言

汶川地震是迄今为止记录到的在一次地震中产生滑坡、崩塌灾害分布范围最广、数量最多的事件，其中，滑坡崩塌灾害在断裂带两侧的不对称分布引起了相关领域科学家的广泛关注。

地震诱发滑坡是自然边坡的原始平衡状态遭到破坏后寻求新平衡状态的过程。已有研究表明活动断裂的几何形状、运动方式对破裂速度、地表变形和地震动参数等起着重要的影响（Abrahamson，Somerville，1996；Oglesby et al，2000a，b；Oglesby，Mai，2012）。地震过程中，地震断层面的错动会改变震区周围的应力状态，一方面是通过降低坡体滑动面的摩擦力，破坏坡体的平衡条件从而导致崩滑发生；另一方面，地震波作用于坡体，使岩土内部产生附加应力，造成了应力集中，从而致使区域内的岩石破裂，岩石块体与斜坡分离形成滑坡。

目前对边坡稳定性的计算研究方法大致可分为极限平衡法、数值分析法和概率法（詹军，于清杨，2002；谭儒蛟等，2004）。相较于其它两种方法，数值分析法可通过滑坡体内部应力应变来分析滑动面的变形及稳定性问题，并可计算出滑坡内应力、变形规律及其对稳定性的影响，因此受到广泛关注（Lam，Fredlund，1993；张伯艳，陈厚群，2001；兰恒星等，2002）。在不同的数值分析方法中，有限元法有较好的灵活性，且求解精度较高，虽然计算量较大，但随着计算机硬件与大型并行机的发展，越来越多的被应用到与地震过程有关的动力学分析和边坡稳定性的计算中（焦玉勇等，2000；刘礼领，殷坤龙，2003；邓琴等，2010）。

大量研究实例和部分动力学模拟成果证实，相较于走滑断层和正断层，逆冲型断层能够在更大的范围内触发更多的滑坡和崩塌。张冬丽等（2009）对汶川地震强震记录的研究认为断层距和上盘效应控制了地震动强度的变化，在断层上盘地震动的峰值加速度明显大于断层下盘等距离台站的相应值，距断层较近的区域内地震动峰值加速度较大，随断层距的增加峰值加速度迅速衰减，且断层下盘衰减明显快于上盘。黄润秋和李为乐（2008）等在汶川地震断层效应的研究中认为，逆冲型断层的上下盘效应是导致滑坡崩塌灾害在上盘分布明显多于下盘的主要原因。陈晓利等（2011a）通过对北川地区滑坡崩塌与岩性、坡度、高程等的影响因素相互关系的分析表明，断层错动方式是影响该区滑坡在断裂带两侧不对称分布的主要原因。Oglesby 等（2000b）针对1999年台湾集集地震中断裂的几何形态对该地震的动力学影响进行了三维数值模拟研究，认为集集地震中逆冲型断层的几何形态不仅是引起断层上盘地震动强度较下盘强烈的主要原因，同时也是较高的地震动加速度分布在上盘近断层区域的主要原因。陈晓利等（2011b）建立了地震动载响应的有限元计算模型，对汶川地震滑坡在地震波传播过程中的触发进行模拟分析，讨论了地震波传播过程的滑坡危险区分布范围。

本文拟基于静态和动态理想模型，使用有限元方法，建立了地震动载效应的动态计算模型，对汶川地震滑坡在地震波传播过程中的触发进行动力学计算模拟分析，以期得到地震波传播过程对边坡稳定性的影响。

1.   地震动载对铲形断层边坡稳定性影响的有限元计算模拟

地震动态载荷即地震波对滑坡体的波动震荡。波动震荡在岩土体内部产生了相应的附加应力，与活动断裂构造和滑坡地形相互作用，造成了应力集中，进而模拟点脉冲源地震波对边坡稳定性的影响。实现地震波的模拟，对震源（深度为14 km）处节点输入指定加速度，得到相应的位移分布。

1.1   研究区区域概况

龙门山地处西部青藏高原腹地与东部扬子地台的结合部位，两侧构造复杂，地貌相差较大。龙门山前陆盆地地壳厚度约为44 km，向西部高原过渡过程中地壳厚度急剧增厚，厚度可达63 km，形成一个莫霍面陡变带（徐锡伟等，2010；陶玮等，2011）。爆破地震剖面结果显示松潘—甘孜地块及龙门山推覆体在深度20—30 km处存在一厚度约为10 km的低速层（刘启元等，2009）。

龙门山断裂带可以分为以走滑断裂为主的东北段和以逆冲断裂为主的西南段。汶川地震破裂带起始于西南段，初始破裂段断裂性质以逆冲为主，破裂的最大值发生在西南段，该段破裂主要由垂直于断层走向的水平向推挤造成，在研究中可以近似为二维平面模型（陶玮等，2011）。前人研究表明映秀—北川断裂随着深度增加而倾角变缓，在15—20 km深度上逐渐趋于水平，并入低速层（张培震等，2008；Shen et al，2009；徐锡伟等，2010）。

汶川地震中，单薄山脊、孤立山体、地形陡缓变化急剧部位、山体上部高程较高部位等对地震震动有明显的放大效应，通常是滑坡、崩塌的易发部位。大部分的地质灾害发生在高程相对较高部位，形成大量的高位滑坡与崩塌（黄润秋，李为乐，2008）。而汶川地震地质灾害的极强发育区是在断层上盘的10 km范围以内；上盘10—20 km范围是地质灾害的强发育区；上盘20—30 km和下盘0—10 km范围是地质灾害中等发育区；此外为弱发育区。绝大多数的大型滑坡都分布在距断层5 km范围内，大于10 km发生大型滑坡的概率很小（黄润秋，李为乐，2008）。

本文根据以上地质背景及汶川动力学背景，参考陶玮等（2011）的铲形逆冲断层模型及材料参数，制定本文研究模型。

1.2   计算方法

本文以FEPG为有限元平台，构建了含有边坡结构的活动断裂带的二维动载模型，加载脉冲震源，模拟地震振动对边坡稳定性的影响。

弹性波动方程考虑介质阻尼，经过空间离散后方程为

式中，M为质量矩阵，A为加速度向量，V为速度向量，U为位移向量，C为阻尼矩阵，S为刚度矩阵，F为载荷向量。

计算采用罚函数法在已知节点上指定加速度，加载单位脉冲函数，即式（1）等号两侧分别加上一个“大质量”，如要实现第i个分量指定节点的加速度使$\ddot{u}_i＝a_0 $，则需将第i个方程

用 $ M_{ii}＋M_0 $代替$ M_{ii} $，其中 $ M_0\gg M_{ij} $，相应地在方程右边加上$ M_0a_0 $，即：

则可得到：

即实现了指定节点加速度。

之后采用wave速度法对方程进行时间离散，经过时间离散的方程，可以采用先求位移再求速度方式求解，也可以采用先求速度再求位移的方式求解，本研究中采用先求位移的方式。

1.3   几何模型的建立

弹性波理论认为，结构对地震滑坡的产生有重要的影响。设计模型如图1所示，采用宽500 km，深100 km的二维简化模型，大尺度的模型可以减少有限元计算的边界影响。断层宽度约为200 m；在0—15 km深度处断层倾角为45°，15—20 km深度处断层倾角逐渐变缓并在20 km以下连接中地壳的低速层；地壳厚度从西向东由70 km递减为40 km；最下层为上地幔。在地表设计了边坡作为地貌的简化，为了说明滑坡的上下盘影响，选取的边坡在断层出露点左右对称分布，用L和R分别表示边坡位于断层左侧和右侧，边坡a，b，c分别距断层1 km，30 km，100 km；每个边坡高5 km，宽10 km。

图  1  二维计算模型

Figure  1.  2-D computational model

下载: 全尺寸图片 幻灯片

1.4   初始条件、边界条件及模型参数

汶川地震滑坡由灰岩、白云岩、玄武岩等较坚硬岩石以及页岩、千枚岩等软岩所组成，灾区中部多为软岩（陈晓利等，2011a）。根据地质构造情况，模型介质分为五类：地壳、低速层、地幔、断层以及边坡。

本文计算基于FEPG有限元生成程序。计算中采用四边形网格，共划分了45 540个单元，4万6 245个节点，每个单元平均约为1 km²，断层处加密。

汶川地震震源深度约为14 km （徐锡伟等，2008），本文模型的震源设计在断层上14 km深度处，以双力偶δ源模拟。模型在地震波传播的时间尺度内（共计算100 s），忽略青藏高原的挤压作用引起的地表位移，右边界受四川盆地的阻挡，底部受地幔的支撑作用，将三个边界面均设定为滚桶约束，即面的法向固定，切向自由；因地表无法积累应力，为与现实情况相同，将其设为自由表面。为消除边界效应，模型边界采用吸收介质。模型材料参数列于表1。

表  1  模型材料参数（修改自陶玮等，2011）

Table  1.  Material parameters of the model （modified from Tao et al ，2011）

介质	密度/（kg·m−3）	泊松比	弹性模量/1010 Pa		介质	密度/（kg·m−3）	泊松比	弹性模量/1010 Pa
地壳	2.7	0.25	8.37		断层	2.7	0.25	5.45
低速层	2.7	0.31	3.35		边坡	2.7	0.25	8.37
地幔	3.3	0.30	12.10

下载: 导出CSV 

| 显示表格

2.   计算结果分析

为方便比较和分析，二维简化模型结果图中字母和数字标识如图2所示，边坡上编号1， 2， 3以远离断层为序的三个点位，其中1和3号点位置同高度。La2即代表断层左侧第一个边坡2号点。

图  2  边坡点位分布示意图

R与L分别代表点在断层的右侧与左侧；a为边坡编号；1号和3号点等高

Figure  2.  Sketch map of the results analysis

R and L represent points on the right and left sides of the fault，respectively；a is the slope numbers of the fault；points 1 and 3 at the same height

下载: 全尺寸图片 幻灯片

计算模型中各点的x方向即时位移振动，如图3—5所示。图中，曲线La2即代表左侧第一个边坡2号点位置位移随时间变化。

图  3  水平方向断层左（a）右（b）两侧不同边坡相同位置位移变化比较

L和R表示断层的左右两侧，a，b，c代表边坡，1，2，3为边坡的不同位置，La1即代表断层左侧第一个边坡1号点位振动图，下同

Figure  3.  Displacement variation comparison of different slopes in the same location on both sides of the fault in horizontal direction

Different color represent the vibration figure of different location；L，R represent the left and right sides of the fault；a，b，c represent the side slope from the near to the distant；the number represents different location on the side slope；La1 represents the vibration figure of the No.1 point of first side slope on the left of the fault，the same below

下载: 全尺寸图片 幻灯片

图  5  水平方向断层左（a）右（b）两侧相同边坡不同位置位移变化比较

Figure  5.  Displacement variation comparison of the same side slope in different location on left side （a） and right side （b） of the fault in horizontal direction

下载: 全尺寸图片 幻灯片

断层同侧不同边坡的相同点号位置处振动如图3所示。断层上盘：距离断层1 km的边坡La最大振幅为1.5 m， 1和2号点位最大振幅均为1.5 m，3号点位最大振幅约为0.8 m；距离断层30 km的边坡Lb最大振幅为0.7 m，1，2和3号点位最大振幅分别为0.5，0.6和0.7 m；距离断层100 km的边坡Lc振幅值相对微弱。断层下盘：距离断层1 km的边坡Ra最大振幅为1.2 m， 1，2，3号点位最大振幅分别为1，1.2和1 m；距离断层30 km的边坡Rb最大振幅约0.2 m， 3个点位最大振幅相近；距离断层100 km的边坡Rc振幅值相对微弱。

图4为断层两侧相同距离的边坡上相同点号位置的振动图。断层上盘（左侧），距离断层1 km的边坡La，1号点位最大振幅为1.5 m，3号点位最大振幅约为0.8 m；距离断层30 km的边坡Lb，1号点位最大振幅为0.5 m，3号点位最大振幅约为0.75 m。断层下盘（右侧），距离断层1 km的边坡Ra，1号和3号点位最大振幅相近，均为1 m。距离断层30 km的边坡Rb，1号点位最大振幅均为0.22 m，3号点位最大振幅约为0.18 m

图  4  水平方向距断层两侧边坡a，b在1号 （a）和3号 （b）点位移变化比较

Figure  4.  Displacement variation comparison of point 1 （a），3 （b） of slope a and b from both sides of the fault in horizontal direction

下载: 全尺寸图片 幻灯片

断层两侧同一边坡上不同点号位置的振动如图5所示，最大振幅值列于表2.

表  2  x方向不同位置最大振幅及时刻

Table  2.  Maximum amplitude and its time in different places of along x direction

边坡位置	边坡a振动方向						边坡b振动方向						边坡c振动方向
	正方向			负方向			正方向			负方向			正方向			负方向
	时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m
上盘1号点	6.25	1.50		4.40	−0.35		5.60	0.50		-	-		22.1	0.024 0		-	-
上盘2号点	6.70	1.50		5.60	−0.60		6.00	0.60		-	-		27.90	0.029 0		-	-
上盘3号点	5.60	0.82		3.75	−0.30		6.25	0.65		-	-		31.80	0.025 0		-	-
下盘1号点	10.00	1.00		10.06	−0.35		13.75	0.22		-	-		45.0	0.012 5		-	-
下盘2号点	8.00	1.20		6.90	−0.30		13.00	0.20		-	-		45.0	0.012 5		-	-
下盘3号点	7.50	1.06		9.40	−0.24		10.06	0.18		-	-		45.0	0.012 5		-	-

下载: 导出CSV 

| 显示表格

上盘La边坡中，1，2号点振幅及衰减相似，3号点振幅相对变小；上盘Lb边坡中，三点的振幅及衰减速度相近。下盘Ra，Rb边坡中的三点振幅与衰减速度都比较相近，Rb边坡的衰减时间更长；Rc边坡的最大振幅均比较小，上盘比下盘略大。

由图3—5可以看出，离断层最近的两个边坡La和Ra最大振幅较大，上下盘分别达到了1.5 m与1 m；振幅向远离断层两侧衰减，幅值降低，A_c＜A_b＜A_a ，且最远的边坡振动幅值基本为0；相同距离的点位，上盘比下盘振幅大（a边坡3号点除外）；相同边坡不同位置的最大振幅没有统一规律；振动图中，振动到达时间各不相同，尤其是断层右侧的考察点的振动图可以明显看到，边坡a与边坡b的到时差达到近5 s。

参考点的z方向振动图如图6—8所示。图6为断层同侧不同边坡上相同点号位置的振动图。断层上盘：距离断层1 km的边坡La最大振幅为0.9 m，1，3号点位最大振幅分别为0.9和−1 m，2号点位最大振幅约为0.5 m；距离断层30 km的边坡Lb最大振幅绝对值为0.4 m，1，2和3号点位最大振幅分别为−0.35、−0.4和−0.3 m；距离断层100 km的边坡Lc振幅值相对微弱。断层下盘：距离断层10 km的边坡Ra最大振幅为0.8 m，1，2，3号点位最大振幅分别为0.8，0.75和0.6 m；距离断层30 km的边坡Rb最大振幅约0.1 m，3个点位最大振幅相近；距离断层100 km的边坡Rc振幅值相对微弱。

图  6  垂直方向断层左（a）右（b）两侧不同边坡相同位置位移变化比较

Figure  6.  Displacement variation comparison of different side slope in the same location on left side （a） and right side （b） of the fault in vertical direction

下载: 全尺寸图片 幻灯片

图  8  垂直方向断层左（a）右（b）两侧相同边坡不同位置位移变化比较

Figure  8.  Displacement variation comparison of the same side slope in different location on left side （a） and right side （b） of the fault in vertical direction

下载: 全尺寸图片 幻灯片

断层两侧相同距离的边坡上相同点号位置的振动如图7所示。断层上盘：距离断层1 km的边坡La，1和3号点位最大振幅为0.9和−1 m；距离断层30 km的边坡Lb，1号点位最大振幅为−0.35 m，3号点位最大振幅约为−0.3 m。断层下盘：距离断层10 km的边坡Ra，1号点和3号点位最大振幅分别为0.8和0.6 m；距离断层30 km的边坡Rb，1，3号点位最大振幅均为0.1 m，图8中，每个子图为断层两侧，同一边坡不同点号位置的振动图。最大振幅值列于表3。

表  3  z方向不同位置最大振幅及时刻

Table  3.  Maximum amplitude and its time in different place of z direction

边坡位置	边坡a振动方向						边坡b振动方向						边坡c振动方向
	正方向			负方向			正方向			负方向			正方向			负方向
	时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m		时刻 /s	最大振幅 /m
上盘1号点	3.50	0.93		7.35	−0.74		8.8	0.14		5.59	−0.36		21.11	0.007 4		17.22	−0.014
上盘2号点	7.78	0.57		6.67	−0.82		9.44	0.29		6.11	−0.46		40.55	0.006 1		17.78	−0.013
上盘3号点	7.77	0.91		5.55	−1.04		14.11	0.14		6.11	−0.29		21.11	0.007 8		21.94	−0.013
下盘1号点	4.70	0.81		6.47	−0.59		9.41	0.15		12.78	−0.09		32.22	0.002 3		-	-
下盘2号点	5.55	0.71		10.55	−0.64		17.78	0.14		16.67	0.11		33.33	0.003 4		-	-
下盘3号点	6.11	0.59		8.89	−0.38		19.41	0.13		20.00	−0.11		31.67	0.002 1		-	-

下载: 导出CSV 

| 显示表格

图  7  断层两侧边坡a，b垂直方向上1号（a）和3号（b）点位移变化比较

Figure  7.  Displacement variation comparison of point 1 （a），3 （b） of side slopes a，b from both sides of the fault in vertical direction

下载: 全尺寸图片 幻灯片

上盘La边坡中，1，3号点振幅相似，2号点振幅相对变小；上盘Lb边坡中，三点的振幅相近；下盘Ra，Rb边坡中的三点振幅相差不大，3号点振幅相对小一些；Rc边坡的最大振幅均比较小，上盘比下盘略大。

图6—8中，离断层最近的两个La和Ra边坡最大振幅较大，上下盘分别达到了1 m与0.8 m；相同边坡中点1和点3处振幅相近；类似x方向，振动振幅向远离断层两侧衰减，幅值降低，A_c＜A_b＜A_a ，且最远的边坡振动幅值基本为0；相同距离的位置，上盘比下盘振幅绝对值大（a边坡2号点除外）；振动图中，振动到达时间各不相同，类似x方向，下盘边坡Ra与边坡Rb的到时差也近5 s。

3.   讨论与结论

本文基于FEPG平台，建立了二维动态铲形断层模型，并在地表断层两侧构建对称边坡，模拟滑坡地形，通过对震源（深14 km）处节点指定加速度的方法，加载脉冲震源，模拟计算了地震发生后100 s内地震波传播对离断层分别1 km，30 km及100 km的不同边坡的边坡震动产生的影响，得出以下结论

1） 水平方向，上盘最大振幅为1.5 m，下盘最大振幅为1.2 m；垂直方向，上盘最大振幅为0.9 m，下盘最大振幅为0.8 m。

2） 距断层相等距离处，上盘比下盘振动幅度大。

3） 相同边坡、垂直方向、相同高程位置振幅相近，高程大的位置，振幅相对变小；水平方向，三点振幅无明显规律。

4） 随着与断层距离的增加，边坡的水平和垂直方向振幅都逐渐变小。

黄润秋和李为乐（2008）指出汶川地震的地质灾害具有显著的断层效应，即表现出了明显的“上/下盘效应”，发震断层上盘较下盘地质灾害分布密度大、分布范围更广、地质灾害规模也较大。 陈晓利等（2011a）指出汶川地震灾区的滑坡分布，87% 的滑坡发生在断裂上盘，且滑坡的空间分布与发震断层的空间位置关系密切。随着与断裂距离的增大，滑坡的数量在减少，80%的滑坡发生在距离断层30 km的区域内。本文的结果与已有的研究结果具有很好的一致性。

在三维静态模型中，最大主应力降均出现在断层附近，地震滑坡危险区域主要分布于断层附近的上盘区域，且呈弧形分布，断层距为100 km的区域，滑坡危险性大大降低。与三维静态模型相比，本研究地震各点的振动情况可以直接反映该点的边坡稳定性。可以看出，在远离震源的地方振幅比近震源振幅小很多，也说明在近震的地方振动强度大；近断层的边坡，高程相对小的振幅大，考虑重力作用，高程大的点滑坡可能性变大。计算结果显示，在与断层相等距离的地方，上盘比下盘的点振动强度大，衰减慢，即如果此类断层发震，其他条件相似的时候，近震源上盘的地方发生滑坡的可能性更高，与三维静态模型结果一致。

地震断层两侧边坡遭遇的地震动载，还要充分考虑震源破裂过程、行波效应等的影响，得到相对完整的地震波数据（包含振幅、频谱和持时），本文并未涉及，下一步将更加细致深入地揭示地震发生对同震滑坡的控制性影响。