淤地坝减轻坡沟系统滑坡侵蚀的数值模拟

邹兵华，袁 洁，李占斌，李 鹏

（1.中国水电顾问集团成都勘测设计研究院，四川 成都610072；2.西安理工大学 西北水资源与环境生态教育部重点实验室，陕西 西安710048；3.西安理工大学 水利水电学院，陕西 西安710048）

黄土高原有大小沟壑27万多条，沟壑纵横，支离破碎，沟壑密度3～5km／km2，沟壑面积占土地总面积的20%～40%，其支离破碎的地貌主要是由沟蚀造成的，沟蚀的发展使沟壑面积日益扩大，耕地面积日趋缩小［1］，是世界上水土流失最严重的地区之一。黄土高原的各类沟壑中以沟头前进、沟底下切、沟岸扩张3种形式的沟蚀危害最为严重［1］。沟蚀加剧了面蚀的发展，造成了更多的陡峭临空面，加剧了重力侵蚀［2－3］。重力侵蚀作为黄土区主要侵蚀类型，其研究一直处于滞后状态［4］。理论研究方面，重力侵蚀的微观及宏观力学机理尚不明晰，计算模型的构建刚起步；试验及观测方面，至今仍未形成重力侵蚀系统观测和试验研究体系［4］。根据黄河水利委员会西峰、天水、绥德3个水土保持站在典型小流域的调查，黄土高原沟壑区的西峰南小河沟，重力侵蚀面积占流失面积的9.1%，重力侵蚀量占总流失量的57.5%；黄土丘陵沟壑区第Ⅲ副区的天水吕二沟，重力侵蚀面积占流失面积的30.7%，重力侵蚀量占总流失量的68.0%；黄土丘陵沟壑区第Ⅰ副区的绥德韭园沟，重力侵蚀面积占总流失面积的12.9%，重力侵蚀量占总流失量的20.2%。在黄土塬区和丘陵区，沟头前进多以土体崩塌形式进行，沟岸扩张是崩塌与滑坡共同作用的结果［5－7］。

淤地坝系作为治理水土流失的主要工程措施，在蓄水拦沙、防 洪 保 收 等 方 面 起 到 了 重 要 作 用［1，5，8－10］，有机统一了当地致富和治河的关系，深受黄土高原地区群众的喜爱，同时又为治河部门所关注。打坝淤地、蓄水拦沙是流域水土保持综合治理的一项重要措施，目前国内外有关淤地坝的研究主要集中在坝系相对稳定、减水减沙和生态环境效应等方面［10－11］，而对淤地坝减轻坡沟系统重力侵蚀特别是滑坡侵蚀方面的研究则少有人提及。淤地坝通过抬高沟底侵蚀基准面，提高了沟坡的稳定性，减少了沟坡重力滑坡侵蚀发生的可能性。本研究采用数值模拟方法对随着坝地逐渐淤高沟道坡沟系统的稳定性、滑塌概率以及坡沟系统滑坡侵蚀破坏的部位等方面进行研究，以期为坡沟系统水土保持工程措施及生物措施的配置提供有益的参考，并为评价坡沟系统的稳定性提供一定可靠度的依据。

1 研究区概况

本研究选取的关地沟是王茂沟流域左岸的一条典型支沟。王茂沟是陕北绥德县韭园沟中游左岸的一条试验性治理支沟，是无定河的一条2级支沟，地理坐标为东经110°20′26″—110°22′46″，北纬37°34′13″—37°36′03″，海 拔 高 度 940～1 188m，流 域 面 积5.97km2，主沟长3.75km，沟道平均比降为2.7%，沟谷地面积2.97km2，占流域总面积的46.7%，沟壑密度4.3km／km2。流域内共有沟道46条，其中＜0.1km2的沟道24条，0.1～0.5km2的沟道18条，0.5～1.0km2的沟道2条，＞1.0km2的沟道2条。地面坡度一般在20°以上，坡度＜10°的面积仅占2.2%，坡度＞45°以上的面积占34.7%。地形地貌主要由梁、峁以及分割梁峁的沟谷组成［12］，地形破碎，沟壑纵横，坡陡沟深。梁顶和峁顶面积不大，均略成穹形，坡度在8°～10°，梁、峁坡的坡度一般在20°～35°。梁、峁坡以下为沟谷，沟谷横断面在上游及支沟均呈“V”字型，在下游略呈“U”字型。沟谷坡极陡，一般都在35°以上。按土壤侵蚀区划，属黄土高原丘陵沟壑区第Ⅰ副区。

流域内地质构造比较单纯，表层多被质地匀细、组织疏松的黄绵土覆盖，厚度20～30m，在梁、峁顶，梁、峁坡上均有分布。其下为红色黄土，厚度50～100m，多出露于沟谷。再下为基岩，基岩主要为三迭纪的砂页岩，岩层倾角甚小，接近水平，其露头只在干沟中下游沟床及其两侧有所见。

流域属北温带干旱大陆性气候，年平均气温10.2℃，无霜期160d左右，夏季多东南风，春秋多西北风。流域多年平均降水量513.1mm，降雨量的年际变化率大，年最大降雨量是年最小降水量的3.5倍；降雨的年内分配极不均匀，年内降雨量主要集中在汛期7—9这3个月，汛期降水量占年降水量的73.1%，且多以暴雨形式出现，一次暴雨产沙量为全年产沙量的60%以上，土壤侵蚀以水蚀及重力侵蚀为主，治理前流域年平均侵蚀模数为18 000t／（km2·a），属剧烈侵蚀区，在黄土丘陵沟壑区具有一定的代表性。

2 坡沟系统滑坡侵蚀的数值分析原理与方法

2.1 崩塌滑坡的蒙特卡洛概率计算

对于土质边坡，根据其土体结构、破坏机理和受力状况，可以建立坡体地质条件和环境因素的状态函数：Z＝g（X1，X2，…，Xm）其中 X1，X2，…，Xm为m个具有一定分布、独立统计的随机变量，假定它们的统计量已知。如果把状态函数定义为安全系数，且随机地从诸随机变量Xi的全体中抽取同分布的变量X′1，X′2，…，X′m，则由上述状态函数可求得安全系数的一个随机样本。如此重复，直到达到预期精度的充分次数N，则可得到N个相对独立的安全系数Z1，Z2，…，Zn。安全系数所表征的极限状态为Z＝1，可构造一个随机变量R：当Z≤1时，R＝1；当Z＞1时，R＝0。设在N次随机抽样的试验中，出现R＝1即Z≤1的次数为M，则坡体破坏的概率pf为［13－14］：

式中：pf——坡体破坏的概率；M——出现安全系数Z≤1的次数；N ——随机抽样试验的次数。

此式即为直接蒙特卡洛法计算坡体破坏概率的公式。显然当N足够大时，由安全系数的统计样本Z1，Z2，…，Zn可比较精确地近似安全系数的分布函数G（z），并估计其分布参数，其均值和标准差分别为：

式中：zi——安全系数的样本值；μz——安全系数的均值；σz——安全系数的标准差。

进而可根据G（z）拟合的理论分布，通过积分方法求得破坏概率。本文中N取值为10万次，滑坡计算选用Spencer法。

2.2 有限元强度系数折减法的基本原理

有限元强度系数折减法的基本原理是将坡体强度参数黏聚力c和内摩擦角φ，同时除以一个折减系数ω，得到一组新的c′，φ′，其中：

式中：c——岩土体原本的黏聚力；φ——岩土体原本的内摩擦角；ω——岩土体强度折减系数；c′——岩土体折减后的黏聚力；φ′——岩土体折减后的内摩擦角。

然后作为新的资料参数输入，再进行试算，利用相应的稳定判断准则，程序可以自动根据弹塑性计算结果得到破坏滑动面，确定相应的ω值为坡体的最小稳定安全系数，此时坡体达到极限状态，发生剪切破坏，同时又可得到坡体的破坏滑动面［15－17］。有关研究［15－19］表明，有限元强度折减法的安全系数在本质上与传统方法是一致的。郑颖人等［15］通过多种比较计算说明有限元折系数法用于分析土坡稳定问题是可行的，但必须合理地选用屈服条件以及严格地控制有限元法的计精度。

2.3 强度系数折减法失稳判断准则与屈服准则

采用有限元强度折减法分析边坡稳定性的一个关键问题，是如何根据有限元计算结果来判别边坡是否处于破坏状态。目前的失稳判据主要有两类：第一类是在有限元计算过程中采用力和位移的不收敛作为边坡失稳的标志［18］；第二类以广义塑性应变或等效塑性应变从坡脚到坡顶贯通作为坝坡破坏的标志［19］，以上两种判据得到的安全系数相差不大［15］。本文采用节点不平衡力的不收敛作为判据，便于计算和可视化，同时根据塑性区的范围及其连通状态确定相应的安全系数，以此评价坡体的稳定性。当前流行的有限元软件中采用的均是广义米赛斯准则［15］，由于强度折减有限元法计算的安全系数与选用的屈服准则密切相关，郑颖人等［15，18］等研究发现不等角六边形外接圆（DP1）屈服准则计算的安全系数较传统的莫尔—库仑准则会偏大许多，而郑颖人［15］研究提出的莫尔—库仑等面积圆准则的计算精度与Spencer的平均误差在5%左右，可满足工程要求，并给出了不等角六边形外接圆（DP1）屈服准则求得的安全系数与莫尔—库仑等面积圆准则求得的安全系数之间的转换关系式（5）。本文把在外接圆（DP1）屈服准则下求得的安全系数转换为莫尔—库仑等面积圆屈服准则下的安全系数：

式中：η——外接圆（DP1）屈服准则下安全系数转换为莫尔—库仑等面积圆屈服准则下安全系数的转换系数；φ——岩土体的内摩擦角。

3 坡沟系统滑坡侵蚀的有限元数值模拟与分析

3.1 坡沟系统概化模型及有限元计算模型

根据相关地质资料，关地沟地层构造主要是马兰黄土（Qeol3），梁、峁顶、峁坡均有分布，厚度20～30m，其下为离石黄土（Qeol2），厚度50～100m，多出露于沟坡上，再下主要是三叠纪砂页岩层，基本接近水平，多出露于干沟、支沟的下游及其两侧。通过实测数据统计分析得知，坡沟系统峁坡坡度在21°～25°所占比例最大，达到50%，而实测坡度在21°～25°的平均值为23.1°，取峁坡坡度24°；沟坡坡度在36°～45°所占比例最大，达到43.6%，而实测坡度在36°～45°范围内的平均值为40.8°，取沟坡坡度41°。通过对关地沟基本坡沟地貌特征的综合分析，建立了坡沟系统概化模型。该重力滑坡侵蚀模型取峁坡坡度24°，沟坡坡度为41°。沟坡概化模型土层从上到下分别为马兰黄土（Qeol3）和离石黄土（Qeol2），厚度分别取为土层厚度的平均值25m和75m。

关地沟流域典型坡沟系统的有限元网格见图1。坡沟系统的计算采用平面应变条件下摩尔—库仑不等角六边形外接圆屈服准则（DP1）和非关联流动的理想弹塑性模型［20－21］。坡沟系统二维模型基岩底边界为固定约束，左边界为水平约束，其他自由。弹塑性有限元计算采用大变形静态分析选项，最大迭代次数为1 000。

图1 坡沟系统三维有限元网格

3.2 计算参数取值

根据边坡工程经验［13－14］、现场资料分析、现场及室内岩土物理力学试验和有限元计算的需要，沟坡概化模型各土层材料特征取值如表1所示。

表1 坡沟系统的计算参数

3.3 沟坡系统稳定性随坝地淤高度的变化

由于关地沟4号坝从1959年修建到1987年被洪水冲毁的28a运行过程中，累计於高4.86m，年均淤高0.18m；垮坝后由于群众对坝体的零星修复填筑及2005年的修复加高，坝地被填高，坝地已不再是原来垮坝前的自然淤高，截至2007年坝地累计高10.8m。由于关地沟4号淤地坝跨坝前坝地的年均淤高仅为0.18m，在跨坝前的淤高模拟计算中以5a为1个时间跨度，计算坝地每淤高0.9m时的坡沟系统稳定特征参数；垮坝后为了与垮坝前保持一致性，坝地以0.9m递增。在有限元程序中，模拟淤地坝坝地按0～0.9，0.9～1.8，1.8～2.7，2.7～3.6，3.6～4.5，4.5～4.86，4.86～5.4，5.4～6.3，6.3～6.48，6.48～7.2，7.2～8.1，8.1～9.0，9.0～9.9，9.9～10.8m的淤积过程每淤积（增高）0.9m其相应的稳定系数、滑塌概率和滑塌量的散点图如图2—4所示。

图2 随坝地淤高坡沟系统稳定系数的变化规律

图3 随坝地淤高坡沟系统滑坡侵蚀滑塌概率的变化规律

图4 随坝地淤高坡沟系统滑坡侵蚀滑塌量的变化规律

由图2—4分别拟合出相关系数最高的随坝地淤高坡沟系统的稳定系数、滑塌概率和滑塌量变化规律的相关方程如表2所示。

表2 随坝地淤高坡沟系统的稳定系数、滑塌概率和滑塌量的变化规律

当坝地分别淤高到11.7，12.6，3.5，14.4和15.3m时，由程序计算得到坡沟系统的稳定系数、滑塌概率和滑塌量与通过上述拟合出的相关方程计算得到的对应值的比较如表3所示。

由表2—3可以得出，随坝地淤高坡沟系统的稳定系数和滑塌概率满足二项式分布规律，而滑塌量则满足线性减少分布规律。稳定系数和滑塌概率的相关系数最高，达到0.980；滑塌量的相关系数为0.977，也相当的高。通过坝地累计淤高11.7，12.6，13.5，14.4和15.3m时坡沟系统的稳定系数、滑塌概率和滑塌量的程序计算值和拟合方程计算值的比较，稳定系数的最大相对误差为1.61%，可见稳定系数的变化是很有规律的，而且稳定性好；滑塌概率的相对误差较大，最大达到15.67%。这是由于滑塌概率本身就是另一个概率计算的结果，还有别的影响因素。滑塌量的最小误差7.92%，最大误差20.09%，这可能是由于随着坝地的淤高，坡沟系统形态的改变，滑坡侵蚀发生过程中最危险和最可能滑动面变化，导致滑坡侵蚀中量的波动比较大，故而滑塌量的相对误差较大的原因。总体来说，三者的相对误差均在满足的范围内，拟合出的方程精度较高，可应用于该沟道重力滑坡侵蚀的定量分析和计算。

表3 坡沟系统的稳定系数、滑塌概率和滑塌量软件计算值和相关方程计算值的比较

3.4 稳定系数最小、坡沟系统濒临破坏时的应力和位移分布

由图2可得，在坝地相对原始沟道累计淤高0m的时候，其稳定系数最小，滑塌的概率最大。为了优化配置坡沟系统的水土保持工程措施和生物措施，达到水土保持措施的最佳防治效果，从力学机理上对坡沟系统的应力场和位移场进行分析，得出坡沟系统应力场和位移场的最危险区域，有针对性地进行措施强化，保证坡沟系统的土壤侵蚀减到最小。由于坡沟系统的近似对称性，这里只考虑沟道的1／2来进行有限元分析。关地沟4号坝上游典型坡沟系统的有限元网格见图1。由有限元强度折减法得摩尔—库仑不等角六边形外接圆屈服准则（DP1）下求得此时的安全系数为1.507，将其按式（5）转换为莫尔—库仑等面积圆屈服准则下的安全系数为1.224，其与滑坡分析采用Spencer法求得此时的安全系数1.255的相对误差为2.42%，在郑颖人等［16，18］研究指出莫尔—库仑等面积圆屈服准则的安全系数与Spencer法的平均误差5%范围内，精度满足要求。

由计算可得，坡沟系统中X方向的位移均为正，在沟坡中下部位垂直向坡体内约15m的范围为X方向位移的最大值区域，X方向的其他位移从右向左、从沟坡向峁坡递减，直至减小为最小值0。Y方向均为下沉位移，最大值在从峁顶向左向右各延伸10m左右的扇形区域，其他下沉位移从峁顶向沟坡坡脚递减，直至减小为最小值0。

坡沟系统中除峁坡和沟坡坡缘线向坡体内一定范围内出现了局部拉应力外第1主应力和第3主应力几乎全为压应力，且第1主应力的拉应力区域几乎包括在第3主应力拉应力区域内。由于土体是抗压不抗拉的，出现拉应力就会使土体发生破坏，进而发生侵蚀，这就从力学机理上说明了坡沟兼治的正确性和峁坡和沟坡进行生物措施和工程措施治理的必要性。坡沟系统中第1主应力的最大拉应力值区域位于从峁顶到沟坡坡脚的整个坡缘线在峁坡向坡体内延伸6m左右、在沟坡向坡体延伸3m左右的带形区域内；第3主应力的最大拉应力值在从峁顶经峁坡向沟坡坡缘线向坡体内延伸9m左右的带形区域内。

4 结论

基于边坡分析软件和有限元技术，通过建立关地沟典型坡沟系统的概化模型，得到并验证了随淤地坝坝地淤高坡沟系统的稳定系数和滑塌概率满足二项式分布规律，而滑塌量则满足线性减少分布规律，三者的精度均较高，可用于该沟道重力滑坡侵蚀的预报、滑坡侵蚀的定量计算和指导淤地坝的分期加高。

应用有限元强度折减法对坡沟系统的应力场和位移场进行了分析，指出了稳定系数最小、滑坡概率最大、坡沟系统濒临滑坡侵蚀时的最大应力和位移分布区域。（1）X方向的最大位移：沟坡中下部位垂直向坡体内约15m的范围；（2）Y方向的最大位移：从峁顶向左向右各延伸10m左右的扇形区域；（3）拉应力最大值区域：从峁顶经峁坡向沟坡坡缘线向坡体内延伸9m左右的带形区域。这些区域是坡沟系统最容易发生侵蚀的危险区域，在配置坡面水土保持工程措施（水平沟和鱼鳞坑等）和生物措施的时候要有针对性地进行重点防护，以使坡沟系统的土壤侵蚀降低到最低限度。

［1］唐克丽.中国水土保持［M］.北京：科学出版社，2004.

［2］Williams J R，Berndt H D.Sediment yield prediction based on watershed hydrology［J］.Transaction of the ASAE，1977，20（6）：1100－1104.

［3］Nachtergaele J，Poesen J，Vandekerck H L，et al.Testing the ephemeral gully erosion model（CEGEM）for two mediter ranean environment［J］.Earth Surface Process and Landforms，2001，26（1）：17－30.

［4］薛海，王光谦，李铁键.黄河中游区重力侵蚀研究综述［J］.水科学进展，2009，7，20（4）：599－606.

［5］崔灵周，李占斌，朱永清，等.流域侵蚀强度空间分异及动态变化模拟研究［J］.农业工程学报，2006，22（12）：17－22.

［6］高鹏，刘作新，邹桂霞.丘陵半干旱区小流域土地资源定量化评价研究［J］.农业工程学报，2003，19（6）：298－301.

［7］于国强，李占斌，李鹏，等.黄土高原小流域重力侵蚀数值模拟［J］.农业工程学报，2009，25（12）：74－79.

［8］陈浩.黄土丘陵沟壑区流域系统侵蚀与产沙关系［J］.地理学报，2000，55（3）：354－363.

［9］郭文元，贾志军.晋西淤地坝试验研究文集［M］.郑州：黄河水利出版社，2001：102－165.

［10］郑宝明.黄土高原丘陵沟壑区淤地坝建设效益与存在问题［J］.水土保持通报，2003，23（6）：32－35.

［11］黄自强.黄土高原地区淤地坝建设的地位及发展思路［J］.中国水利，2003（17）：8－11.

［12］郑宝明，王晓，田永宏，等.淤地坝实验研究与实践［M］.郑州：黄河水利出版社，2003：152－367.

［13］王家臣.边坡工程随机分析原理［M］.北京：煤炭工业出版社，1996：68－139.

［14］祝玉学.边坡可靠性分析［M］.北京：冶金工业出版社，1993：45－96.

［15］郑颖人，赵尚毅.有限元强度折减法在土坡与岩坡中的应用［J］.岩石力学与工程学报，2004，23（19）：3381－3388.

［16］Zienkiewicz O C，Taylor R L.The finite element method［M］.New York：McGraw－Hill，1989.

［17］Zienkiewicz O C，Humpeson C，Lewis R W.Associated and nonassociated visco－plasticity in soil mechanics［J］.Geotechnique，1975，25（4）：671－689.

［18］赵尚毅，郑颖人，时卫民，等.用有限元强度折减法求边坡稳定安全系数［J］.岩土工程学报，2002，24（3）：343－346.

［19］吴海真，顾冲时.有限元强度折减法在土石坝边坡稳定分析中的应用［J］.水电能源科学，2006，8（4）：54－58.

［20］石长.赵新铭.用有限元强度折减法分析边坡稳定［J］.隧道建设，2006，26（3）：5－8.

［21］张彩双，李俊杰，胡军.有限元强度折减法的边坡稳定分析［J］.中国农村水利水电，2006（5）：72－74.