基于运动波壤中流理论与无限边坡稳定分析理论的滑坡预测模型

金保明，林 鹏，李光敦

(1.福州大学 土木工程学院，福建 福州 350108；2.台湾海洋大学 河海工程学系，台湾 基隆 20224)

降雨型滑坡是山区流域常见的地质灾害之一[1]，往往导致严重的人员伤亡与经济损失[2]。滑坡还会诱发泥石流与堰塞坝决堤等次生灾害[3]。降雨作为滑坡的触发机制入渗至坡体[4]，将坡面表层土壤软化、泥化，使其抗剪强度降低，从而诱发滑坡。考虑水文变化的坡面稳定分析模型较少，因此，构建考虑降雨变化的滑坡预测模型，为开展滑坡研究提供参考。滑坡灾害与地形特征、降雨强度、土壤特性等因子有着密切关系[5]。坡面表层土壤特性指标是滑坡预测模型的重要参数[6]，但由于土壤厚度、水文地质条件分布不均匀，饱和水位随时间与空间变化特性机理复杂[7]，因此土壤特性与坡面稳定性之间具有复杂关系。

对于一般浅层边坡，其滑坡需考虑径流影响[8]。预测浅层滑坡的方法可以分为两个方面：一方面，传统多采用临界雨量进行滑坡预测[9]，但由于不同流域的地形地貌以及气候条件，所能触发滑坡的降雨条件也有所不同，因此模型应用具有局限性[10]。另一方面，针对滑坡物理机制构建预测模型。预测因降雨引起的浅层滑坡灾害较为常用方法为无限边坡稳定分析。Skempton等[11]以摩尔-库仑破坏准则为基础构建无限边坡稳定分析模型，评估土体边坡稳定。Montgomery等[12]以地形模式为基础，结合无限边坡理论进行坡面稳定性评估。Iverson[13]利用Richards方程的简单解析解，计算斜坡入渗量，并结合无限边坡稳定分析，模拟饱和水位上升造成边坡破坏情形。Lee等[14]利用修正地形指数模式计算流域内饱和水位的时间与空间变化，配合无限边坡稳定分析理论，计算流域内不稳定格点与平均安全系数；根据评估指标分析预报结果，均满足滑坡预测要求[15]。Talebi等[16]将地下水运动理论与坡面稳定分析理论相结合，进一步将其应用于9种不同山坡[17]。马世国等[18]推导强降雨入渗条件下边坡的安全系数随湿润锋深度的变化关系。然而，目前，尚缺乏土壤蓄水与坡面稳定性的相关研究。

浅层滑坡中饱和水位随时间变化是影响坡面稳定的重要因素之一[19]，研究水流运动可以计算出饱和水位变化值[20]。Boussinesq模型描述了一个均匀坡面中无承压含水层运动，但由于其非线性形式，求解过程往往很困难[21]。Fan等[22]将土壤蓄水层简化为1维水流运动，基于特征线法求解无承压含水层运动方程，提出了运动波壤中流方程的解析解。Troch等[23]拟合山坡地形为双变量函数，推导运动波壤中流方程，并应用此模型分析山坡地形对于土壤蓄水量、地下径流量影响。Liu等[24]基于地形因子预测地下蓄水与流量关系，并较好地拟合洪水时期的退水曲线。Afshar Ardekani等[25]基于蓄满机制研究所设计的降雨条件下不同理想山坡土壤含水量变化，结果表明土壤厚度对于壤中流也有着较大影响。金保明等[26]分析了山区流域表层土壤蓄水特性对壤中流运动的影响。然而，鲜有论文讨论土壤蓄水与坡面稳定性之间的关系。

综合以上分析，壤中流运动对于滑坡致灾机理至关重要，滑坡体内孔隙水压力和地下水作用是滑坡的最主要致灾因子，但目前尚缺乏滑坡与水文模型相结合的应用创新。本文将运动波壤中流理论与无限边坡理论相结合进行坡面稳定性分析，并依据土壤蓄水机制，建立滑坡预测模型，并将其应用于实际流域。此外，设计两场双峰对称雨型的降雨过程，分析降雨期间坡面出口处饱和水位变化过程；研究两种雨型和不同土壤厚度条件下坡面流量过程线与安全系数变化过程，为滑坡预测模型提供理论依据。

1 结合水文模型构建坡面稳定性分析模型

以无限边坡理论进行坡面稳定性分析，饱和水位是无限边坡稳定分析理论中需要考虑的重要因素。本文采用运动波壤中流理论估算饱和水位的时空变化。

1.1 模型理论

1.1.1 无限边坡稳定分析理论

通常因降雨入渗而导致的滑坡深度远小于山坡的长度与宽度。模型采用Skempton等[11]提出的无限边坡稳定分析理论：

式中：WFS(x,t) 为坡面上x处t时刻的安全系数；ct为有效黏聚力，N/m2； γm为 天然土重度，N/m3； γs为饱和土重度，N/m3； γw为水重度，N/m3； γb为浮重度，N/m3，γb=γs-γw； β 为斜坡倾角，(°)；φ为土壤内摩擦角，(°)；g为重力加速度，m/s2；d(x) 为土壤厚度，m；h(x,t) 为x处t时刻饱和水位，m。

1.1.2 运动波壤中流理论

由式（1）可发现，安全系数会随着饱和水位的升降而变化。在陡峭山坡区域，可假设地下水流动方向与基岩面平行。应用Fan等[22]所提出1维简化的山坡运动波壤中流理论描述复杂山坡水文过程，并计算随时间变化的饱和水位。引入山坡蓄水能力Sc(x)的概念：

式中：Sc(x)为 单宽土壤蓄水能力，m2；ε为土壤孔隙率；w(x)为坡面宽度，m。土壤蓄水能力由山坡地形决定，土壤中的水流满足达西公式：

式中：Q1(x,t) 为t时刻在x处壤中流流量，m3/s；S1(x,t)为t时刻在x处 单宽土壤蓄水量，m2；K为饱和渗透系数，m/s； dz/dx为基岩面坡度。假定坡面上饱和渗透系数、孔隙率固定不变。山坡地表面形状采用二次函数拟合，且土壤厚度固定不变，因此地表面与基岩面表达式相同。

土壤颗粒之间存在孔隙，雨水受到重力作用沿着孔隙流动。壤中流符合水流连续方程：

式中，r(t)为降雨强度，m/s。将达西公式（3）与水流连续方程（4）结合，基于特征线法解微分方程，将1维山坡剖面地形方程代入特征线方程可以得到土壤蓄水方程式：

式中，ξ为初始位置，A(x)为 0至x处 坡面面积，α、 β、γ为拟合基岩面是二次多项式z=α+βx+γx2的形状系数。通过式（5）计算t时刻在x处单宽土壤蓄水量，联立式（3）与（5）可以得到坡面出口处壤中流流量方程式（6）。进而，计算出降雨时期坡面出口（x=L）处的流量过程线Q1(L,t)：

式（1）中饱和水位可由土壤蓄水方程（5）求解，可得出口处饱和水位h(L,t)随时间变化：

将饱和水位式（7）代入无限边坡稳定分析理论式（1），并对整个坡面积分，得到坡面平均安全系数：

1.2 模型构建

式（1）中饱和水位是触发滑坡的关键因素。本文结合运动波壤中流理论与无限边坡稳定分析理论，构建一个结合浅层滑坡预测模型。具体步骤为：首先，根据式（5）对坡面进行产汇流计算得到土壤蓄水量随时间与空间变化；然后，用式（7）研究坡面出口处饱和水位随时间的变化特性；最后，根据式（6）计算降雨径流过程线，式（8）计算坡面平均安全系数。

2 模型应用

收集研究区域的降雨与流量记录、地形、土壤等自然地理要素资料，根据植被覆盖指数NDVI判别何时何地发生滑坡灾害。图1显示在2009年第8号莫拉克（MORAKOT）台风期间台湾高雄市那玛夏区南沙鲁里发生滑坡灾害的区域，其总面积为1.052 km2，滑坡坡面面积为0.266 km2。表1为无限边坡稳定分析所需的土壤参数，可经由现场土壤采样数据进行土壤颗粒比重与三轴压缩实验、卫星影像识别等方式获取，以此模型参数进行径流与滑坡模拟。

表1 土壤参数Tab. 1 Soil parameters

图1 滑坡区域示意图Fig.1 Schematic diagram of landslide areas

莫拉克台风带来的特大暴雨袭击台湾中南部、东部及东北部地区，并发生多处重大滑坡灾情。2009年8月6日上午8时30分气象局发布第8号台风莫拉克的海上陆上台风警报，8月7日5时台风暴风圈开始接触台湾至脱离（8月9日18时）一共经历61 h，其中心登陆至出海经历14 h。临近雨量站有达喀尔努瓦里雨量站，海拔高度为1 040 m，东经120°42′19″，北纬23°16′47″，台风期间累计降雨量高达2 580 mm，最大降雨强度为76 mm/h。

在本文模型应用之前需对坡面地形特征参数进行提取。基于GIS提取基本单位的地形特征，拟合坡面形状为二次函数，坡面地形参数如表2所示。根据相对数值高程理论将研究区域划分为12个坡面（图1）。

表2 坡面产汇流参数Tab. 2 Slope runoff yield and confluence parameters

采用滑坡预测模型对研究坡面进行模拟，图2为特征稳定坡面NL1、NL5、NL6的降雨径流过程线，图3为同等降雨情况下特征滑坡坡面L1、L2、L3的径流过程线。

图2 坡面NL1、NL5、NL6降雨径流过程线Fig.2 Rainfall runoff hydrograph of NL1、NL5、NL6 slope

图3 坡面L1、L2、L3径流过程线Fig.3 Runoff hydrograph of L1、L2、L3 slope

根据第1.1节理论计算坡面的安全系数，在莫拉克暴雨期间，有8个坡面安全系数均大于1，即为稳定状态（NL1、NL2、NL3、NL4、NL5、NL6、NL7、NL8）；其余4个坡面出现安全系数小于1的时段，即为不稳定状态（L1、L2、L3、L4）。同等降雨情况（图2（a））下，图4为特征稳定坡面NL1、NL5、NL6的安全系数变化，图5为滑坡坡面L1、L2、L3的安全系数变化。

图4 坡面NL1、NL5、NL6安全系数变化过程线Fig.4 Temporal variation of safety factor of NL1、NL5、NL6 slope

图5 坡面L1、L2、L3安全系数变化过程线Fig.5 Temporal variation of safety factor of L1、L2、L3 slope

由图2（a）、4、5可知，安全系数随降雨强度大小而产生变化，其原因是由于降雨量增加，导致更多的雨水入渗到土壤中，造成饱和含水层的水位逐渐抬升，即距滑动面以上的饱和水位升高，发生浅层滑坡的可能性亦随之增加。研究发现基于降雨型浅层滑坡预测模型所得到的坡面分析结果，与滑坡发生事件实测记录相符，说明本文模型具有一定的可靠性，可用于实际流域的滑坡预测。

3 模型验证与讨论

3.1 模型参数设置

为了进一步探讨降雨条件下土壤厚度对于坡面稳定性的影响，设计了一种坡面，对该坡面进行降雨径流模拟与稳定性分析。设计坡面为平直型山坡（图6），坡面长度、宽度、高度分别为100 m、100 m、60 m。

图6 设计坡面示意图Fig.6 Schematic diagram of design slope

设计降雨延时24 h、降雨量相同，但雨型对称的两场降雨。图7（a）为所设计的先大后小的雨型Ⅰ，其中：第1个降雨尖峰发生在6 h，最大降雨强度为30 mm/h；随后12 h发生第2个降雨尖峰，降雨强度为20 mm/h；图7（b）为所设计的先小后大的降雨型态Ⅱ，第1个降雨尖峰为20 mm/h，第2个降雨尖峰为30 mm/h。

图7 两种设计雨型的降雨过程线Fig.7 Rainfall hydrograph of two designed rainfall types

土壤参数采用表1中的参数，饱和渗透系数0.000 5 m/s、孔隙率0.5。设置3组不同土壤厚度，分别模拟设计降雨过程在设计坡面所产生的坡面出口处饱和水位变化过程、流量过程线以及坡面平均安全系数随时间变化过程。

3.2 降雨时期坡面出口处饱和水位变化

山坡表层松散土壤往往是浅层滑坡的来源，而土壤厚度决定山坡表层蓄水能力大小，土壤厚度越大蓄水能力越强。讨论土壤厚度对坡面出口处饱和水位影响，图8（a）、（b）分别为雨型Ⅰ、Ⅱ降雨条件下不同土壤厚度饱和水位变化。

由图8（a）可知：当d=0.1 m时，随着坡面出口汇集的水量增加，饱和水位上涨速率逐渐加快，直到4.50 h土壤饱和，饱和水位上涨至0.1 m；接着土壤排水，饱和水位缓慢下降。当d=0.5 m时，饱和水位上涨速率随着降雨强度减低而减缓，直到19.00 h饱和水位上涨至0.5 m。当d=1.0 m时，饱和水位上涨速率逐渐减缓，直到23.67 h饱和水位上涨至0.57 m；接着土壤排水，饱和水位缓慢下降。

由图8（b）可知：当d=0.1 m时，随着坡面出口汇集的水量增加，饱和水位上涨，直到5.50 h饱和水位上涨至0.1 m。当d=0.5 m时，随着坡面出口汇集的水量增加，饱和水位上涨，直到19.83 h饱和水位上涨至0.5 m。当d=1.0 m时，随着坡面出口汇集的水量增加，饱和水位上涨，直到23.67 h饱和水位上涨至0.57 m；接着土壤排水，饱和水位缓慢下降。

图8 两种雨型条件下不同土壤厚度饱和水位变化Fig.8 Temporal variation of saturated water level using different soil thickness assumptions under two rainstorm types

在整个降雨过程中，初始时刻土壤未饱和，由于饱和渗透系数均为0.000 5 m/s，饱和水位上涨速率相同，雨水下渗至土壤，向坡面出口方向汇流；随着降雨的持续，饱和水位上升至地表面；而后降雨量减小，土壤排水速率大于降雨强度，饱和水位缓慢下降。

在土壤蓄水阶段，蓄水前期汇集雨水速率较快，而后逐渐缓慢，降雨时期饱和水位呈现快涨慢退的性质。土壤厚度主要决定蓄水量，土壤厚度越大，土壤达到饱和状态所需的时间越长。虽然两场降雨雨型不同，但是随着土壤厚度的增大，饱和点的移动过程越相似。研究降雨时期饱和水位变化，有助于了解滑坡物理机理。

3.3 降雨时期流量过程线与安全系数变化

图9为在两种雨型和不同土壤厚度条件下坡面流量过程线与安全系数变化过程。

若山坡表层有0.1 m厚度的土壤，如图9（a）所示：雨型Ⅰ降雨输入同样产生先大后小双峰流量过程线；流量过程线第1个尖峰发生在6.17 h，尖峰量为0.082 m3/s；第2个尖峰发生在18.17 h，尖峰量为0.055 m3/s。初始时刻，坡面安全系数为最大值3.522；随着降雨量的增大，安全系数快速减小；安全系数最小值发生在6.50 h，两个降雨尖峰之间，此时安全系数降低至最小值2.760；随着土壤饱和，安全系数并没有很大变化；降雨停止时，安全系数缓慢上升至初始值。同等土壤厚度下，如图9（b）所示：雨型Ⅱ降雨输入同样产生先小后大双峰流量过程线；流量过程线第1个尖峰发生在6.17 h，尖峰量为0.055 m3/s；第2个尖峰发生在18.17 h，尖峰量为0.082 m3/s。安全系数最大值3.522；随着降雨量的增大，安全系数快速减小；安全系数最小值发生在18.50 h，安全系数降低至最小值2.760。

然后，假设表层土壤厚度0.5 m。如图9（c）所示：雨型Ⅰ降雨输入模型中产生单峰的流量过程线，尖峰发生在19.33 h，尖峰量为0.047 m3/s。安全系数最大值1.391；随着降雨量的增大，安全系数快速减小；安全系数最小值发生在19.67 h，在第2个降雨尖峰之后，此时安全系数降低至最小值1.011。同等土壤厚度下，如图9（d）所示：雨型Ⅱ降雨输入模型中产生单峰的流量过程线，尖峰发生在20.17 h，尖峰量为0.054 m3/s。安全系数最大值1.391；随着降雨量的增大，安全系数快速减小；安全系数最小值发生在20.83 h，此时安全系数降低至最小值0.990。

最后，假设土壤厚度1.0 m足够大时，较大的土壤厚度能蓄积更多的水量。如图9（e）所示：雨型Ⅰ降雨输入模型中产生平滑的单峰流量过程线，尖峰发生在23.33 h，尖峰量为0.019 m3/s。安全系数最大值1.125；随着降雨量的增大，安全系数快速减小；安全系数最小值发生在21.67 h，安全系数降低至最小值0.926。同等土壤厚度下，如图9（f）所示：雨型Ⅱ降雨输入模型中产生单峰的流量过程线，尖峰发生在23.67 h，尖峰量为0.020 m3/s。安全系数最大值1.125；随着降雨量的增大，安全系数快速减小；安全系数最小值发生在21.83 h，安全系数降低至最小值0.913。

为了直观比较，根据图9中的结果，表3列出了两种雨型和不同土壤厚度条件下坡面出口流量峰值与坡面最小安全系数统计。

由图9和表3可知：初始时刻，降雨未进入表层土壤，含水量低，土壤抗剪强度大，山坡不易发生滑坡，此时安全系数为最大值。随着降雨量增大，含水量增加，抗剪强度降低，安全系数快速减小。当降雨量减少，土壤排水速率大于降雨强度，土壤含水量减少，抗剪强度升高，安全系数增大，逐渐恢复至天然状态下的安全系数。坡面平均安全系数与壤中流运动密切相关，在降雨条件下，表层土壤厚度增大，含水量增加，壤中流流量增大，最小安全系数减小且延后，对称雨型的流量过程线越相似，发生最小安全系数的时间也越接近。在土壤厚度d=0.5 m条件下，先大后小的降雨输入土壤汇集水量较少，最小坡面安全系数为1.011；而先小后大的降雨输入土壤汇集水量较多，最小安全系数为0.990。这是由于前期少量降雨输入，土壤来不及排水，而后较大的降雨强度使得雨水蓄积在土壤中不易排出，土壤含水量增大，导致抗剪强度降低，坡面安全系数减小。

表3 坡面出口流量峰值与坡面最小安全系数统计Tab. 3 Maximum discharge and the minimum factor of safety of different soil thickness

图9 不同土壤厚度条件下流量过程线与安全系数变化过程Fig.9 Discharge hydrographs and temporal variation of safety factor using different soil thickness assumptions

4 结 论

将无限边坡模型与水文模型相结合，构建了坡面稳定性分析模型，探究了降雨过程中土壤厚度对于坡面稳定性影响，主要结论如下：

1）应用山坡运动波壤中流理论与无限边坡理论分析坡面稳定性，构建了结合水文模型的降雨型浅层滑坡的预测模型。与传统滑坡预测模型相比，模型参数较少。

2）验证了滑坡预测模型的可靠性。模型应用于台湾高雄市那玛夏区达喀尔努瓦里小流域内，根据实测雨量，推求出降雨径流过程线和坡面平均安全系数，进而分析的坡面稳定性结果与实际卫星影像图所得到的结果一致，率定的土壤物理参数符合实际，可以认为所构建滑坡预测模型，在实际应用中具有一定的可靠性。

3）研究了双峰连续降雨过程中坡面出口处饱和水位运动过程，以及坡面安全系数变化过程。随着降雨的进行，饱和水位与壤中流呈现快涨慢退的性质，而坡面安全系数的变化规律为快速减小缓慢增大。在壤中流汇流过程中，土壤含水量逐渐增大，安全系数最小值往往发生在最大降雨强度之后。

4）探讨了双峰对称降雨条件下土壤厚度对于坡面出口流量与稳定性影响。降雨时坡面稳定性与壤中流相关，土壤厚度增大，壤中流流量增大，安全系数减小，降雨对于流量与坡面安全系数的影响降低，在对称的两种雨型条件下的流量过程线越相似，可能发生滑坡的时间越接近。在相同降雨量条件下，先小后大的降雨型态所汇集的水量较多，更易发生滑坡灾害。

以小流域为例，所建立的降雨型浅层滑坡预测模型，运行效率快，可应用于不同地形地貌条件的区域，满足滑坡预警要求。将来有条件时，可进一步以中等流域或大流域为研究对象，进行深入分析，结合实时水文资料与天气预报降雨资料建立更为实用的滑坡预测模型，为可能发生的滑坡灾害提供预警依据。