基于Scoops 3D模型的库岸边坡稳定性研究

王志浩，党进谦，郭 钊

（西北农林科技大学水利与建筑工程学院，陕西杨凌 712100）

水库在运行期间，为了满足防洪、灌溉、发电等需要，库水位需要进行周期性的涨落，在这过程中岸坡的稳定性会出现不同程度的下降，进而发生库岸滑坡灾害。国内外库岸滑坡灾害屡见不鲜，如湖南省资水河柘溪水库的塘岩光滑坡［1］、三峡库区的千江坪滑坡［2］、意大利Vajont 水库的左岸滑坡［3］等均对人民生命财产安全造成重大损失。

目前，国内外针对库岸边坡有着大量的研究成果。郑颖人［4］、吴琼［5］通过建立岸坡饱和与非饱和非稳定渗流模型，探讨了库水位波动情况下岸坡内地下水渗流场的变化情况；江强强［6］、肖诗荣［7］采用物理模型试验研究岸坡在不同水位波动、降雨作用下的稳定性演化规律，揭示滑坡的变形过程；仉文岗［8］等使用FLAC 进行数值模拟，表明库水位下降对岸坡稳定性影响较大，且水位下降过程存在最危险水位；杨敏［9］等采用离散元软件FINAL 探索了锚索加固后岸坡稳定性分析方法；Xia［10］等通过对三峡区域石榴树包滑坡多年的监测研究，表明岸坡浅层滑坡主要受降雨影响，而库水变动是造成岸坡深层滑坡的主要原因，D R Montgomery 和W E Dietrich［11］建立了可考虑降雨的SHALSTAB 稳定渗流模型；J W Godt 等［12］建立了TRIGRS 模型。史笑凡［13］等提出了不确定性边坡稳定性Bootstrap-BCa估计方法；王禹玺［14］等基于Copula 理论和Bootstrap 方法提出了边坡系统可靠度区间的高效估计方法；但是，对三维真实岸坡整体区域稳定性进行的动态研究较少。

所以，本文以山美水库为研究对象，首先采用解析法即包辛涅斯克潜水非稳定性方程计算库水位升降过程中岸坡地下水位，然后采用三维斜坡稳定性计算模型Scoops 3D 进行库岸边坡稳定性的计算，并就岩土力力学参数对水库岸坡的稳定性影响程度进行分析。此论文主要目的是研究库水位的升降速率和岩土体力学参数对岸坡稳定性的影响，为水库运行和库岸滑坡评价提供参考。

1 山美水库区域概况

1.1 地质、水文概况

山美水库位于福建省泉州市西北部，地处晋江上游的东溪，是晋江流域唯一一座综合利用的大型水利枢纽工程，水库集雨面积1 023 km2，总库容6.55 亿m3。研究区雨量年内分布极不均匀，多集中在4-9月，占全年雨量的70%～80%，而10-12月仅占6%～10%，且夏季多暴雨、大雨。研究区域地层发育不完全，出露的地层主要为上三叠统-侏罗系、第四纪。其流域地势由西北向东南显著倾斜，西部山岭纵横，地形多以低山与丘陵为主，东部呈阶梯状，散布着山间小盆地。

1.2 库水位变动特征

近年来，山美水库库水位在70～95 m 之间呈现规律性地波动，可以将库水位的波动过程分为4 个阶段（见图1），具体划分如表1所示。本论文主要针对库水位蓄水、下降阶段进行库岸边坡的稳定性研究。

表1 山美水库运行阶段划分表Tab.1 Shanmei Reservoir Operation Stage Division Table

2 库岸边坡稳定性研究方法

2.1 库岸边坡地下水浸润线计算

地下水是影响库岸边坡稳定性的重要因素，所以，在库岸边坡稳定性研究中准确、合理地确定库水位升降作用下的浸润线分布具有重要的意义。在进行库水位变动引起的岸坡地下水位变化计算时，解析法具有方便快捷的优势。目前常用的解析解法有线性解、近似解、简化解，此论文计算库岸边坡浸润线变化时使用Boussinesq 一维非稳定渗流方程。计算中的基本假设为：采用裘布依假定、潜水层为均质含水层。由这些假设条件下的潜水非稳定性运动微分方程可由包辛涅斯克（Boussinesq）方程得到，即：

式中：h是指岸坡内部浸润线高程，m；K表示岩土体渗透系数，m/d；μ为给水度；t表示水位变动时间，d；H表示岩土体含水层厚度。

计算简图如图2所示。依据基本假定可以得到水位升降的半无限含水层中岸坡地下水非稳定渗流计算渗流模型：

式中，u表示土体中水位变化高度；用水位升降始末高程的平均值hm代替岩土体含水层厚度，简化方法将其看为常量，即取hm等于H。

将上述的渗流计算模型利用拉普拉斯积分变换和逆变换即可得到微分方程如下：

为方便计算对η-M(η)函数曲线进行拟合，得到的多项式拟合公式为：

于是就可以得到库水位升降时浸润线的简化计算公式，其表达式为：

式中：h0，0为库水升降前的水位，d；v0为库水位升降速率，m/d，水位下降时取负值，水位上升时取为正值。

2.2 Scoops 3D模型

Scoops 3D 是在2015 由美国地质调查局的M E Reid 等［15］开发的一个三维斜坡稳定性计算模型，此模型以数字高程模型为基础，采用球形滑面，可系统全面的搜索大量的潜在滑坡并计算其稳定性。在进行稳定性计算时，此模型可以耦合多种不同环境、力学等条件参数进行分析，常见的有：斜坡体的孔隙水压力、坡体的地下水位高度、岩土体力学参数、地震荷载、真实复杂地形等。

Scoops 3D 模型基于一个椭球形的潜在滑面利用三维柱体极限平衡分析法计算潜在滑坡的稳定系数［16，17］，进行稳定性计算的基本单元为三维柱体，主要过程为，首先在设定的数字高程模型范围上方形成大量的球心，然后生成若干的搜索球体和潜在滑坡三维柱体。球面和数字高程模型相交形成的不规则坡体内所有三维柱体构成潜在的滑动体，模型中潜在三维滑坡体如图3所示。

模型使用力矩平衡法即简化的Bishop 法进行三维斜坡稳定性系数计算，Scoops 3D中基本计算公式为［18］：

式中：Am，n为岸坡潜在滑面的表面积；αm，n为岸坡潜在滑面的视倾角；βm，n为岸坡潜在滑面的倾角；Wm，n为潜在滑动三维柱体的重量；um，n为孔隙水压力，作用在三维柱体潜在滑面上；km，n为三维柱体中心受到的水平的振动荷载；cm，n和φm，n滑面处的岩土体抗剪强度参数；Rm，n为转动轴到潜在三维柱体滑面的距离，在三维中随着柱体的位置等改变Rm，n也不断发生变化。

3 研究区数据参数

3.1 DEM数据

此次研究中，山美水库的表面地形DEM数据由地理信息数据库得到，使用ArcGIS 软件进行处理后得到分辨率为5 m 的高精度DEM，将其作为模型的基础地形高程数据。库岸边坡高程以及坡度数据分别如图4、5所示。研究区库岸边坡的面积约为16.1 km2。

3.2 岩土体参数

模型计算所需的主要物理力学参数为：库岸边坡岩土体的土体重度、饱和抗剪强度以及土体的渗透系数。在Scoops 3D模型计算时，根据山美水库地质勘察资料选取岩土体物理力学参数，其中黏聚力为20 kPa，内摩擦角为16°，渗透系数为0.05 m/d，土体部分饱和重度为18 kN/m3，饱和重度为20 kN/m3。

3.3 计算工况的确定

3.3.1 不同库水位升降速率下库岸边破稳定性

近年来，山美水库水位在90～75 m之间涨落，水位下降期的平均水位降速为0.8 m/d。为了研究库水位升降速率对库岸边坡稳定性的影响，设置水位变动范围为90 m 至75 m，选取库水位升降速率分别为0.5、0.8、2 和3 m/d，计算不同水位升降速率下水位变动最后时刻瞬时状态下山美水库库岸边坡稳定性的变化情况，据此评价水位升降速率对库岸边坡稳定性的影响。

3.3.2 库岸边坡稳定性对岩土参数的敏感性

模拟计算过程中使用毕肖普强度理论进行岸坡稳定性计算，朱锐［19］等通过室内直剪试验，揭示了干湿循环作用下膨胀土边坡土体黏聚力出现明显降低、内摩擦角呈小幅度减小的衰减规律，据此确定黏聚力c、内摩擦角φ、土体重度γ为岸坡体稳定性主要影响因素。对库岸边坡稳定性进行大量的数值模拟，来获得库岸边坡在主要影响因素变化时稳定性的变化规律，进而进行敏感性分析。模拟中水位波动条件设定为由90 m 以0.8 m/d 的速度降至75 m。通过控制其中两个影响因素，另外一个因素作为变量进行计算，具体模拟方案见表2。

表2 模型岩土力学参数变化表Tab.2 Model rock and soil mechanics parameter change table

3.4 模型参数

为提高模型计算效率和精度，需预先设定搜索范围即三维搜索矩阵，在此范围内Scoops 3D 模型将生成大量的球形滑面球心。研究的区域分辨率为5 m，岸坡最高点高程为261 m，共有1 464×1 956 个栅格。据此定义三维搜索矩阵的起点为1×1号栅格，搜索矩阵的最低高程设定为300 m，最高高程为500 m，搜索半径增量为0.5 m。库岸区域模型参数见表3。

表3 山美水库库岸区域模型参数Tab.3 Model parameters of Shanmei reservoir bank area

4 结果分析

根据节2 的岸坡稳定性研究方法和节3 的参数，并按照3.3节中设定工况进行山美水库库岸边坡稳定性分析。结果分析中参考相关标准［20］依据稳定性系数对岸破稳定性进行分级。分级标准见表4。

表4 岸坡稳定性划分标准Tab.4 Bank slope stability classification standard

4.1 不同水位升降速率下岸坡稳定性计算

4.1.1 不同水位降速下岸坡稳定性计算结果

将不同水位降速条件下，库水位降至75 m时岸坡体瞬时浸润线分布状态作为地下水位参数，并将岩土体力学等参数一并输入到Scoops 3D 模型中进行计算，得到山美水库库岸边坡在90 m 常水位和不同水位下降速率下库水位为75 m 时的岸坡安全系数，并进行稳定性分级，结果如图6所示。其中水位下降速率为0.5 和0.8 m/d 时岸坡稳定性变化极小，稳定性分级图基本一致，故在图6中只绘制了v=0.8 m/d 降速下的稳定性分级图。在各工况间的趋势对比如图7所示。

由图6、7 可以得到：①库水位在90 m 高水位运行期间，库岸边坡大部分具有较高的稳定性，稳定性等级在较高以上的岸坡占全库岸段的68.77%，低稳定性等级的岸坡占全库岸段的14.37%。②库水位从90 m 下降至75 m 过程中，当下降速度为0.8 m/d 时，低稳定性岸坡占全库岸段的19.44%，较90 m 常水位时增加了5.07%。岸坡出现了新的零星分布的低稳定性区域，右侧库岸段低稳定和较低稳定性区域较之前明显扩大，部分低稳定性区域有点状分布连通成片。③库水位下降速度由0.5 m/d 增加到0.8 m/d 时，库岸边坡稳定性基本不变。④库岸边坡的稳定性随着水位下降速度的增大逐渐较小，下降速度增加到2、3 m/d 时，低稳定性岸坡占全库岸段的百分比分别为21.82%、24.84%，相对0.8 m/d面积分别增加了0.38、0.87 km2。因此山美水库在泄洪等排水过程中，库水位下降速率应控制在0.8 m/d之内，且在此范围内根据需要可适当增加减速，但应减少2 m/d、3 m/d高降速的情况。

4.1.2 不同水位升速下岸坡稳定性计算结果

将不同水位上升速率条件下库水位升至90 m 瞬时地下水位状态等条件输入到Scoops 3D 模型中进行库岸稳定性分析，计算得到库水位上升速率分别为0.8、2、3 m/d 时库库边坡稳定性系数并进行分级，各工况间的变化趋势对比如图8所示。

由统计可知：①库水位上升工况下，库岸边坡整体稳定性系数较库水位下降工况下大，处于高稳定区域占全库岸段的百分比增加了5%～10%，低稳定性区域也较水位下降时有所降低，降低约3%～6%。②库水位上升过程中，随着库水位上升速度的增加岸坡稳定性逐渐增大，上升速度为2、3 m/d 相对0.8 m/d 低稳定性区域分别降低了1.17%、1.57%。所以在水库运行期间，水位上升速度在3 m/d范围内可适当增大。

据此可得到，研究区库岸滑坡稳定性对库水位升降速很敏感。库水位下降高度相同条件下，随着库水位下降速率的增大库岸稳定性显著降低；而库水位上升高度保持不变时，库水位上升速度的增大对库岸边坡的稳定性有利。

4.2 岩土体参数敏感性分析

采用单因素敏感性分析方法按照3.3.2 节设计工况对库岸边坡滑坡稳定性进行研究，库岸滑坡稳定性分级计算结果如图9所示。

由图9可知，岩土体力学参数的变化对库岸边坡稳定性影响很显著，将这些指标按照标准归一化处理进而进行分析。采用敏感度计算公式［21，22］：

式中：L(xk)为岩土体参数xk的敏感度；F*为xk取力学参数标准值xk*时对应的稳定性区域占比；Δxk、ΔF为岩土体力学参数x*k和稳定性占比F的变动值。

计算过程中，对每个力学参数，选取库岸边坡低稳定性区域占岸坡总库岸段的百分比进行敏感度计算。计算得到，改变黏聚力c、内摩擦角φ、土体重度γ，敏感度分别为1.14、3.77、1.26。由此可知，各个岩土体力学参数对岸坡稳定性的影响程度大小为：内摩擦角φ＞土体重度γ＞黏聚力c。

5 结 论

（1）运用Scoops 3D 模型计算不同工况下岸坡稳定性的结果表明，库水位下降速率为0.5、0.8、2、3 m/d 时，低稳定性库岸区域占去全库岸段的百分比为：18.95%、19.44%、21.82%、24.84%；库水位上升速率为0.8、2、3 m/d 时，低稳定性库岸区域占去全库岸段的百分比为：18.59%、17.42%、17.02%。即库岸边坡稳定性随库水位下降速率的增大而降低，而较大的水位上升速率对岸坡稳定性有利。

（2）采用单因素敏感性分析方法研究了库岸边坡稳定性对岩土体力学参数的敏感性，结果表明，敏感度由大到小依次是内摩擦角φ、土体重度γ、黏聚力c。□