重庆小南海滑坡形成机制离散元模拟分析

申 通，王运生，吴龙科

（成都理工大学 地质灾害防治与地质环境保护国家重点实验室，成都 610059）

1 引言

重庆黔江小南海是1856年6月10 日在湖北咸丰县大路坝与重庆黔江区小南海之间发生的6.25级地震作用下诱发形成的大型滑坡，长期以来一些学者对该滑坡开展了大量调查研究工作，如1968年初四川和湖北两省地质大队的研究人员先后深入到小南海实地考察，得出小南海由地震形成的结论；李庆海[1]通过调查和搜集实际资料，对小南海地震发生的地质背景、成因机制进行了论述。崔芳鹏等[2]在对小南海地震崩滑堆积体进行野外地质调查和工程地质勘查等相关工作的基础上，阐述了该崩滑堆积体的基本特征和形成机制。2000年，中国科学院及国家地震局有关专家考察小南海[3]，把震级定为6.25 级，极震区烈度定为Ⅷ度。然而，地震是怎样形成小南海滑坡的，至今没有合理的解释。在震级和烈度较低的情况下能够形成如此大规模的滑坡值得深入研究。

汶川地震以来，许多学者对地震滑坡开展了大量研究，特别是汶川地震余震及芦山地震主震与余震斜坡地震动监测数据揭示地震中在高位斜坡地震突出部位，地震动峰值加速度相对于背景值会放大2～3 倍[4－7]，由此联想到小南海滑坡是否也是由于高位斜坡地震过程中地形放大效应所致？为了证实这一推断，本文进行了数值模拟。以往研究的数值方法大都基于有限元和有限差分法[8]，这些程序都具有“界面单元”或“节理单元”，使程序能够模拟问题中的不连续面，扩大了程序的应用范围，但公式在一个或多个方面通常受到限制。首先，考虑很多相互切割的节理就可能打乱系统的逻辑关系；其次，不可能自动识别新的接触面进行自动考虑；第三，计算公式可能有小位移和无转动条件限制，通常更适用连续介质的程序。离散单元法允许离散块体发生有限的位移和转动（包括完全脱离），能够更好地解决这种含有众多不连续面的岩质边坡大变形动力分析问题。本文以小南海顺层岩质边坡为例，运用离散元软件UDEC对地震动力条件下斜坡变形破坏及运动机制进行初步分析研究。

2 小南海滑坡基本特征

崩滑体位于小南海北岸，距黔江约32 km。崩滑体周围地形以中低山为主，海拔400～1 900 m，属浅、中切割的中、低山地形见图1，滑坡平面图见图2。崩滑体包括大垮岩崩滑体和小垮岩崩滑体两处，大小垮岩处山体顶面原始高程均约1 300～1 500 m，其中大垮岩堆积面积约为79×104m2，小垮岩堆积面积约为71×104m2，堆积总面积约150×104m2，均厚约40 m，总体积约为6 000×104m3。大垮岩崩滑体后缘位于大垮岩陡崖下部，高程约900 m，见岩体崩滑形成的崩滑壁；前缘位于小南海岸边，高程约670 m。崩滑体纵向上呈凸形，整体北高南低，地形坡度为10°～35°，中上部相对较缓，前缘陡坡临空见图3。

图1 小南海滑坡体全貌[2]Fig.1 Panoram a of collapse-slidemass induced by Xiaonanhai earthquake in Chongqing[2]

图2 小南海滑坡体平面[2]（单位：m）Fig.2 Plan view of collapse-slidemass induced by Xiaonanhai earthquake[2](unit:m)

图3 地震前大跨岩地质剖面复原图Fig.3 Recovered geological section of Dakuayan before earthquake

崩滑体整体由志留系页岩碎块石组成，局部崩坡积物中还含有二叠系灰岩块石。崩滑体下基岩为志留系中统罗惹坪组（S2lr）页岩,完整性较好，岩层产状130°～140°∠5°～8°。坡体中节理发育，主要有3 组构造节理：（1）N45°～55°E∠70°～85°SE；（2）N50°～60°W∠80°～85°NE；（3）N55°70°W∠80°～90°SW[2]。

3 小南海崩滑体地震触发下失稳过程分析

3.1 小南海陡崖危岩块体地震抗倾覆稳定性分析模型

小南海陡崖岩体在地震力作用下发生了崩滑，通过对陡崖危岩块体的抗倾覆稳定性进行力学分析，定量分析出小南海陡崖危岩块体在多大的地震力作用下失稳。

建立如图4 所示的简化块体模型，模型是一个宽度为b，高度为h 的截面为矩形的块体，块体的自重为W。

图4 倾覆破坏分析块体模型Fig.4 Block model for block toppling failure

考虑块体向坡外的倾覆，取块体受到的各个作用力相对于块体右下角点O 的力矩，顺时针力矩为倾覆力矩，逆时针力矩为抗倾覆力矩。计算并比较抗倾覆力矩与倾覆力矩的相对大小，即可得出对块体抗倾覆稳定性的认识。

为简化分析，设陡崖块体与其他部分岩体的接触面无连接强度，块体仅受自重W和指向坡外的水平外力H 的作用，对于陡崖松动岩体的倾覆稳定性而言是相对不利的情况。这时，水平力H 产生的力矩MH为倾覆力矩，块体自重W 产生的力矩MW是抗倾覆力矩。定义抗倾覆力矩MW与倾覆力矩MH的比值为块体的抗倾覆稳定系数KM，即

当MW＞MH时，KM＞ 1，块体处于稳定状态；当MW＜MH时，KM＜ 1，块体有倾覆的危险；当MW=MH时KM=1，块体处于临界状态。

参照图所标注的块体尺度，取h=500 m，b=300 m。考虑水平地震力的作用，则水平力H 可按式（2）取值。

式中：β为水平地震力系数，等于场地水平地震动峰值加速度α 与重力加速度g 的比值，即β=α/g 。

3.2 分析结果

利用式（1）、（2），代入相关数据，计算得出了分析模型可能发生倾覆崩落的临界水平地震加速度为0.6g。2000年底，中国科学院及国家地震局有关专家考察小南海，把当年发生的那场黔江地震极震区烈度定为8 度[3]，通过查阅中国地震烈度表（1980年）规定，烈度为Ⅷ时相对应的峰值加速度为0.25g，明显小于临界水平地震加速度0.6g，似乎不足以对岩体产生破坏作用。但是，汶川、芦山大量震害调查发现，有些局部场地地形对地震作用下的次生灾害发育具有明显的影响或控制作用，这种地形效应会引起实际结构在地震时的严重破坏[9]。地震监测发现，高台上的地震动参数与平地及高台陡坎下部的有很大差别，陡峭坡的地形放大效应大于平缓坡的地形放大效应[10]。已有研究表明，三面临空的山体、孤立山脊、沟－脊－沟等地形部位地震加速度得到放大，坡顶的水平峰值加速度可以增大到谷底的2～3 倍[11]。因此，考虑到地形放大的因素，烈度为Ⅷ的黔江地震能够使小南海陡崖岩体发生倾覆崩落。为了进一步验证此结论，本文通过离散元数值模拟的方法，在模型底部施加相应的地震波来模拟地震力作用下小南海陡崖岩体的变形破坏过程。

4 数值模型和边界条件的建立

离散单元法属于块体理论数值方法，该方法允许离散块体发生有限的位移和转动（包括完全脱离），弥补了有限元法或边界元法对介质连续及小变形的限制。特别适合具有节理结构的块裂介质的不连续变形问题分析。UDEC[12]软件正是基于离散元的二维数值程序，能够模拟不连续块体集合体系在静力或动力荷载条件下的响应。所以，UDEC 能够很好地模拟地震力作用下边坡节理化岩体变形破坏后的运动过程。

小南海是1856年（清咸丰六年）黔江地震形成的，而有关斜坡地震动力响应的研究表明[13－15]。对于层状结构岩质边坡来说，水平地震力是影响其稳定的主要因素，本文主要考虑水平地震力对坡体的影响，以峰值加速度平均值为0.25g 的地震波（见图5）作为输入地震波，该地震波采用5·12 汶川地震文县台网地面加速度实测数据，近似模拟小南海地震滑坡变形破坏及运动过程。

图5 水平向加速度时程曲线Fig.5 Time-history curves of horizontal acceleration

4.1 离散元数值模型的建立

本文根据小南海地质剖面建立概化模型，由于篇幅所限，只对具有代表性的大垮岩做以研究。如图6 所示，模型底边宽2 138 m，左侧边界高1 446 m，右侧边界高735 m。在模型中设置了一些监测点，用以监测小南海滑坡在地震作用下位移、速度随时间的变化过程。

从小南海崩滑体物质组成调查来看，崩滑体整体由志留系页岩碎块石组成，局部崩坡积物中还含有灰岩块石，岩性较简单。为计算方便，将坡体岩性进行概化。把滑体岩性视为强风化岩，由页岩、灰岩组成；把滑床岩体岩性视为弱风化岩，由页岩夹粉砂岩和灰岩组成；结构面分为节理面和岩层面。坡体及结构面的物理参数（根据崔芳鹏等[16]）见表1，2。

4.2 模型边界条件的选取

斜坡失稳后，由于滑体中的单元变形很小，相对其运动位移可以忽略不计，因此滑体可以选取刚性本构模型（见图7），滑床可以选取塑性本构模型以及Mohr-Coulomb 屈服准则，结构面则选取库仑滑动模型。

图6 小南海滑坡数值模型和监测点分布Fig.6 Numerical model and distribution of monitoring points of Xiaonanhai landslide

表1 坡体的力学参数Table 1 Mechanical parameters of slope rock mass

表2 结构面的物理力学参数Table 2 Physico-mechanical parameters of structural plancs

图7 小南海坡体离散元数值模型网格剖分Fig.7 Meshing generation of numerical model for Xiaonanhai landslide

本文的边界条件：（1）两侧边界水平方向约束、底部边界竖直方向约束，在自重应力下使其达到平衡。（2）将两侧位移边界条件改为左右侧自由场边界，底部边界竖直向位移约束并设置黏滞边界（见图6），能够有效减少因反射地震波而损失的能量，在坡体底部施加地震力。

图8 最大不平衡力时程曲线Fig.8 Time-history curve of the maximum unbalanced force

4.3 输入地震力的条件

加速度代表的地震力荷载不能直接施加于黏滞边界上，需利用数值积分的方法先转化为速度时程曲线，再转化为应力曲线加在底部边界上[17]。由加速度曲线可知，施加的地震能量主要分布于90 s 以内，本文把前90 s 作为地震动力条件输入。

5 数值模拟结果的分析

5.1 最大不平衡力

随着地震动力荷载输入，最大不平衡力随时间的变化曲线如图8 所示。由图可见，约32 s 时达到3 800 MPa 的最大值，70 s 后不平衡力已经变得很小，90 s 后不平衡力趋近于0，表明基本达到稳定，计算收敛。

5.2 加速度放大效应的研究

为研究在地震作用下坡体中质点加速度的响应规律，沿坡体竖向设置了监测点A1～A5，其中监测点A5 位于结构面上，其他4 点则位于基岩上。通过模拟可知，加速度沿斜坡竖向存在显著的放大效应，如图9 所示。通过对比实测加速度时程曲线，可以得出A1～A5 的加速度放大系数见表3。从图表中可知，沿坡体竖向地震水平向加速度的放大效应很明显，监测点中A5 的水平加速度放大系数较A1～A4 显著增大。

图9 监测点A1、A5 水平向加速度时程曲线Fig.9 Time-history curves of horizontal acceleration of monitoring points A1 and A5

表3 不同监测点加速度峰值及放大系数Table.3 Peak values and amplification factors of acceleration of different monitoring points

地震波在坡体内发生多向反射及散射，从而产生了坡体波动震荡加速效应；此外坡体锁固段变形体中贮存的弹性应变能得到释放后也产生了加速效应，最终在崩滑体失稳启动时具有一个较大的加速度。当加速度达到最大值时，坡体失稳破坏而开始运动，因此在动力响应分析时，结构面的加速度放大系数相比稳定的坡体要高出许多，也是监测点A5加速度放大系数如此巨大的原因所在。从地震动响应的分析中得出，该坡体具有较大的启动加速度，也表明具平行坡面陡倾控制性结构面的高陡突出地形对地震波有明显的放大作用，加速度放大2 倍以上，这一结果与汶川震后监测结果吻合。

5.3 滑坡运动过程的模拟

变形累积达到质变后斜坡就进入到运动阶段，见图10。运行10 s 时，坡体中上部开始发生错动和位移，坡顶表层开始有少量岩块崩落；运行20 s 时，坡体获得了一个较大的启动加速度，块体发生大范围的错动，坡顶尤为明显，在地震波作用下岩块发生振荡抛射现象，岩体崩滑形成的崩滑壁已经明显可见；运行30 s 时，坡体进入整体崩滑阶段，上部岩块发生了显著垮落；运行45 s 时，整体崩塌基本结束，堆积在斜坡上的岩块沿滑带迅速下滑；运行75 s 时，整体滑动基本结束，滑体沿滑带一直滑动到小南海湖底，使堆积体前缘边界淹没于小南海湖区水位之下；运行75～90 s，滑坡整体不再滑动，表层有部分岩块被震荡弹落、翻转；运行90～130 s为滑体自稳过程；运行130 s 时，达到稳定状态，滑坡运动结束。通过与现场地质调查资料进行对比可以得出，该数值模拟获得的滑坡堆积情况与现状具有较好的一致性。

图10 各个时间模型运行所对应的不同运动状态Fig.10 Movement states of landslide element at different running states of numerical model

通过以上数值分析可知，小南海的失稳模式为振动拉裂-震荡抛射-崩滑堆积。山高坡陡，在地震动作用下地形放大效应显著，坡体振动拉裂、破裂解体，上部岩土体首先失稳破坏，崩塌滑落后给下部岩土体施加推力，导致滑体整体破坏[18]。地形放大效应随坡度比的增大而增大（坡度比是指坡高与地形坡宽的比值）[19]，所以通常这种高台陡坎坡体的山顶地震动力响应最为强烈，缓坡及山麓破坏相对不明显。

为了研究在地震作用过程中滑体的运动特征，在滑体上设置了监测点B1～B8（见图6）。通过比较滑体内监测点B1～B3和滑体表层监测点B4～B8 的位移时程曲线（见图11），滑体内岩体位移时程曲线较为平滑，而岩体表层位移时程曲线则呈波浪形，由此可以说明，在下滑过程中，由于受到地震纵波竖向的震荡作用，表层岩块体存在垂直抛落现象。也与许强及黄润秋的调查结果相符合[20]。

图11 滑体中监测点的位移随时间变化曲线Fig.11 Time-history curves of displacement of monitoring points in landslide mass

6 结论

（1）斜坡地形放大效应在强震区是客观存在的。小南海滑坡不是瞬间地、均衡地失稳破坏，由于地形放大效应的影响，坡顶首先形成应力集中，随着持续的地震动力输入，应力不断向中前部锁固段集中，坡体振动拉裂、破裂解体，上部岩土体崩塌滑落后给下部岩土体施加推力，导致滑体整体破坏。

（2）斜坡体中质点加速度在地震过程中具有地形放大效应，加速度放大在2 倍以上。在动力响应分析时，结构面监测点的加速度放大系数相比稳定的坡体要大的多，该坡体失稳破坏时具有较大的启动加速度，也表明具平行坡面陡倾控制性结构面的高陡突出地形对地震波有明显的放大作用。

（3）地震作用下的运动过程表明：滑坡失稳模式为峰值加速度放大-增加的振幅迫使岩体顺平行坡面陡倾控制性结构面迅速拉裂-沿缓倾层面滑移-高速脱离滑源区-巨大的势能和动能驱动坡地做长距离运动。

[1]李庆海.1856年黔江小南海地震初探[J].四川地质学报,2001,21(3):184－185.LI Qing-hai.A preliminary study of the Xiaonanhai earthquake in Qianjiang in 1856[J].Acta Geologica Sichuan,2001,21(3):184－185.

[2]崔芳鹏,胡瑞林,张明,等.重庆小南海地震崩滑体的基本特征及形成机制研究[J].工程地质学报,2009,17(2):167－174.CUI Fang-peng,HU Rui-lin,ZhANG Ming,et al.Basic characteristics and formation mechanism of collapsed and slided debris mass induced by Xiaonanhai earthquake in Chongqing city[J].Journal of Engineering Geology,2009,17(2):167－174.

[3]杜海忠,刘坤梅,丰志美.灾变景观旅游地开发研究—以小南海地震遗址景区开发为例[J].桂林旅游高等专科学校学报,2006,17(4):430－434.DU Hai-zhong,LIU Kun-mei,FENG Zhi-mei.Research on the development of tourist destination with cataclysmic sight—A case study of the tourism planning of Xiaonanhai[J].Journal of Guilin Institute of Tourism,2006,17(4):430－434.

[4]王运生,徐鸿彪,罗永红,等.地震高位滑坡形成条件及抛射运动程式研究[J].岩石力学与工程学报,2009,28(11):2360－2369.WANG Yun-sheng,XU Hong-biao,LUO Yong-hung,et al.Study of formation conditions and toss motion program of high landslides induced by earthquake[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(11):2360－2369.

[5]CELEBI M.Topographical and geological amplifications determined from strong—motion and aftershock records of the 3 march 1985 Chile earthquake[J].Bulletin of the Seismological Society of America,1987,77:1147－1107

[6]ASHTORD S A,SITAR N.Analysis of topographic amplification of inclined shear waves in a steep coastal bluff[J].Bulletin of the Seismological Society of America,1997,87(3):692－700

[7]何蕴龙,陆述远.岩石边坡地震作用近似计算方法[J].岩土工程学报,I998,20(2):66－68.HE Yun-long,LU Shu-yuan.A method for calculating the seismic action in rock slope[J].Chinese Journal of Geotechnical Engineering,1998,20(2):66－68.

[8]曹琰波,戴福初,许冲,等.唐家山滑坡变形运动机制的离散元模拟[J].岩石力学与工程学报,2011,30(增刊1):2878－2887.CAO Yan-bo,DAI Fu-chu,XU Chong,et al.Discrete element simulation of deformation and movement mechanism for Tangjiashan landslide[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(Supp.1):2878－2887.

[9]张建毅,薄景山,王振宇,等.汶川地震局部地形对地震动的影响[J].自然灾害学报,2012,21(3):164－169.ZhANG Jian-yi,BO Jing-shan,WANG Zhen-yu,et al.Influence of local topography on seismic ground motion in Wenchuan earthquake[J].Journal of Natural Disasters,2012,21(3):164－169.

[10]刘琨,卢育霞,李倩.复杂地形条件对强地震动特征的影响[J].世界地震工程,2012,28(4):80－88.LIU Kun,LU Yu-xia,LI Qian.Influence of complicated topography conditions on strong ground motion characteristics[J].World Earthquake Engineering,2012,28(4):80－88.

[11]罗永红,王运生.汶川地震诱发山地斜坡地震动的地形放大效应[J].山地学报,2013,31(2):200－210.LUO Yong-hong,WANG Yun-sheng.Mountain slope ground motion topography amplification effect induced by Wenchuan earthquake[J].Journal of Mountain Science,2013,31(2):200－210.

[12]Itasca Consulting Group Inc.Universal distinct element code(version 4.0)[M]// User′s guide.Minneapolis:Itasca Consulting Group,Inc.2004.

[13]王存玉,王思敬.边坡模型振动实验研究[C]//岩体工程地质力学问题(七).北京:科学出版社,1987.

[14]祁生文,伍法权.边坡动力响规律研究[J].中国科学(E辑),2003,33(增刊1):28－40.QI Sheng-wen,WU Fa-quan.Study on dynamic responses of slope[J].Science in China(Series E),2003,33(Supp.1):28－40.

[15]门玉明,彭建兵,李寻昌,等.层状结构岩质边坡动力稳定性试验研究［J］.世界地震工程,2004,20（4）：131－136.MEN Yu-ming,PENG Jian-bing,LI Xun-chang,et al.Research on vibration testing of models for dynamic stability of rock slope with layered structures[J].World Earthquake Engineering,2004,20（4）：131－136.

[16]崔芳鹏,胡瑞林,殷跃平,等.纵横波时差耦合作用的斜坡崩滑效应离散元分析—以北唐家山滑坡为例[J].岩石力学与工程学报,2010,29(2):319－327.CUI Fang-peng,HU Rui-lin,YIN Yue-ping,et al.Discrete element analysis of collapsing and sliding response of slope triggered by time difference coupling effects of P and S seismic waves—Taking Tangjiashan landslide in Beichuan County for example[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(2):319－327.

[17]KUHLEMEYER R L,LYSMER J.Finite element method accuracy for wave propagation problems[J].Journal of the Soil Mechanics and Foundations Division,ASCE,1973,99(5):421－417.

[18]裴向军,黄润秋.“4·20”芦山地震地质灾害特征分析[J].成都理工大学学报(自然科学版),2013,40(3):257－263.PEI Xiang-jun,HUANG Run-qi.Analysis of characteristics of geological hazards by“4·20”Lushan earthquake in Sichuan,China[J].Journal of Chengdu University of Technology(Science &Technology Edition),2013,40(3):257－263.

[19]刘洪兵,朱晞.地震中地形放大效应的观测和研究进展[J].世界地震工程,1999,15(3):20－25.LIU Hong-bing,ZHU Xi.Advance on topographic amplification effects of seismic response[J].World Earthquake Engineering,1999,15(3):20－25.

[20]许强,黄润秋.5·12 汶川大地震诱发大型崩塌灾害动力特征初探[J].工程地质学报,2008,16(6):721－729.XU Qiang,HUANG Run-qiu.Kinetics characteristics of large landslides triggered by May 12th Wenchuan earthquake[J].Journal of Engineering Geology,2008,16(6):721－729.