基于潘家铮法的滑坡涌浪参数敏感性分析

李太清，曾树元，王 瑞，张 斌，刘 羿

（1. 贵州乌江水电开发有限责任公司索风营发电厂，贵州 贵阳 550215；2. 中国电建集团贵阳勘测设计研究院有限公司，贵州 贵阳 550081）

0 引 言

降雨、水位变动等作用通常会影响到河道或水库岸坡的稳定性，一旦失稳将与水体发生强烈的相互作用及能量传递形成涌浪，并沿着河道传播，一方面在对岸形成较大的涌浪爬坡，越岸及冲刷，另一方面，涌浪向上、下游传播，对航道航行构成巨大威胁，对沿岸的水工建筑物、居民区造成破坏，对下游大坝产生巨大冲击，甚至发生漫顶、冲毁等灾害，进而造成水库下游产生洪水灾害。滑坡涌浪问题作为近水滑坡灾害链中的重要问题，已经被国内外学者广泛关注和研究。滑坡运动是涌浪预测计算的基础，滑坡体形状、滑坡体前缘形态、滑坡体体积、滑面摩擦角及水面宽度等因素均对涌浪计算产生影响［1］。

目前滑坡涌浪灾害的主要研究手段有物理模型试验，解析求解和数值模拟等。其中物理模型实验受实验条件限制，实验结论从小尺度推广到实际大尺度时受相似准则制约，误差较大；数值模拟计算准确，模拟效果好，黄宇云、王志超［4，9］等使用CFD、SPH 等数值方法对涌浪过程进行了数值模拟，但存在建模复杂，在实际应用中计算规模大、计算效率低的问题，在实际工程应用中成本较高；基于解析求解的经验估算方法在实际工程中具有计算简便、快捷，计算成本低的优势，也得到了广泛应用。其中常用的涌浪计算方法有美国土木工程学会建议法、Node 法［2］、水科院经验公式法、潘家铮法，黄锦林、滕光亮等对各方法进行了对比［3，5］，万志勇等进行了潘家铮法计算参数的敏感性分析［10］。

潘家铮法［6-8］作为目前我国现行相关规范中规定的滑坡涌浪分析方法之一，为了分析采用该法进行库区滑坡涌浪分析时计算参数取值对计算结果的影响，本文以贵州省索风营库区卞家寨滑坡为例，对其涌浪灾害进行了评估和分析，并对滑坡强度折减参数、涌浪高度计算中级数项数选取、涌浪传播计算中壁面反射系数等指标进行了敏感性分析，为潘家铮法计算滑坡涌浪中计算参数的选取提供依据。

1 潘家铮法滑坡涌浪计算原理概述

采用潘家铮法进行滑坡涌浪计算分析。计算分为滑坡变形过程计算，涌浪传播及爬坡高度计算两个步骤。

1.1 潘家铮法滑坡计算原理

根据《建筑物的抗滑稳定和滑坡分析》［6］中潘家铮法，滑坡滑动过程采用条分法进行估算，得到滑坡体运动速度。 滑坡滑动过程计算中基于条分法的基本思想，将滑坡体垂直划分为有限条块。对单个条块进行受力分析如图1，假设条块在发生滑动时按照刚体运动规律滑动，不同条块之间在滑动过程中发生变形错动。滑动面上的反力有，法向力Ni，切向力fiNi+Ci，其中i表示条块编号，fi= tanφ，φ为滑动面处摩擦角，Ci为黏聚力大小。垂直界面上的反力有，条块间垂直作用面的法向力H，H+ ΔHi，沿作用面的切向力Q，Q+ΔQi，考虑水压力项Ui，条块自重Wi（地下水位以下使用饱和容重计算），αi为条块底边线坡脚，图1中ax，ay分别表示水平和竖直方向加速度方向。

图1 条块i受力分析Fig.1 Forces analysis of block i

使用简单条分法计算，求解过程中假定条块间完全无摩擦及其他阻力，忽略各条块间剪切作用，则剪切力Q为0。基于两个方向平衡方程，可求出滑坡整体滑动的水平加速度ax，即此时刻所有条块的水平加速度。基于运动学方程，各条块的水平速度为：

其中ΔL为条块宽度，移动过该水平距离ΔL所需时间：

基于以上过程循环迭代计算，直至阻力大于下滑力，或求解得到加速度为负值，产生减速作用，最终滑动停止。

1.2 潘家铮法涌浪计算原理

潘家铮法中涌浪计算包括滑坡前缘初始涌浪高度、涌浪在滑坡对岸的爬坡高度、及涌浪传播至下游不同位置处的涌浪高度计算。

滑坡变形滑动方式分为水平运动和垂直运动两种模式，且各自对初始水波产生不同程度的影响。《建筑物的抗滑稳定和滑坡分析》［6］给出了滑坡不同运动模式示意图。

1.2.1 初始涌浪高度计算

水平变形运动产生的初始涌浪按公式（3）计算：

式中：h为当量水深，1.17为经验系数；vx为滑坡入水的水平滑动速度；ζ0x为初始涌浪高度水平运动项分量。

垂直变形运动初始涌浪按式（4）、式（5）计算，其中v为垂直运动模式下滑坡体沿滑坡床滑动的速度：

影响涌浪初始高度的主要因素是滑坡速度，且相对滑动速度在水平运动模式中与初始涌浪高度呈线性关系，在垂直运动模式中与初始涌浪高度呈非线性关系。

1.2.2 涌浪传播计算

基于初始涌浪高度，涌浪传播计算遵循孤立波的变形传播及反射规律。图2为卞家寨滑坡涌浪传播计算示意图。

图2 卞家寨滑坡涌浪计算示意图Fig.2 Calculation diagram of surge propagation of Bianjiazhai landslide

其中A点为滑坡断面对岸的爬坡点，A'点为滑坡下游任意位置处爬坡点，xi为下游某位置到滑坡距离，L为滑坡体宽度，B为涌浪传播的水面宽度。k为壁面反射系数，对于直立、光滑、不透水的理想岸壁，考虑水波到达后完全反射，此时反射系数为1.00，在实际应用中，反射系数取值约在0.85 与0.95 之间［6］。入射波与岸坡法线夹角为θ，对滑坡下游任意位置，反射系数常取kcosθ。

滑坡对岸A点最大涌浪高度ζAmax按公式（6）计算：

其中ζ0为初始涌浪高度，级数项数n取决于滑坡历时T及涌浪从本岸传至对岸所需时间，c为预估波速：c=，根据表1 选择合适的级数项数取值。

表1 涌浪高度计算级数求和中项数取值Tab.1 Number of items in summation of surge height calculation series

按公式（7）计算得到滑坡下游对岸A'位置最大涌浪高度ζ：

其中θn为传到A'点的第n次入射线与岸坡法线的夹角，计算滑坡上下游其他位置处最高涌浪时，反射系数k=0.85～0.95。

2 卞家寨滑坡涌浪计算方案

2.1 工程概况

卞家寨滑坡位于索风营水电站上游水库库区左岸，滑坡体位置距坝址直线距离约1.50 km，隶属黔西县素朴乡马路村沿河一组，其地质成因为早期后缘陡壁大规模崩滑堆积，堆积体分布高程760～1 100 m，堆积体平行河床平均长800 m，垂直河床方向平均长600 m，推测厚度50～130 m，体积约3 300 万m3。

根据实际开裂情况，为便于涌浪计算，创建模型滑坡体所在水库水位为837.00 m，滑坡体部分位于水面以下，根据滑坡体形状，滑动趋势方向及地质勘探得到的目前滑坡体的变形破坏范围，滑坡体前缘不稳定部位的体积约为4×106m³（图3）。本文将针对该不稳定滑体失稳可能产生的涌浪灾害进行分析和评估。

图3 条分法计算断面位置及库区坝址位置示意图Fig.3 Schematic diagram of section position calculated by slice method and dam position

2.2 滑坡变形滑动计算参数

根据图3 所示滑坡体形态，及现场滑坡体不同区域滑动方向，综合考虑滑坡体整体运动趋势，选取将滑坡体两侧体积均分的断面Ⅰ-Ⅰ作为滑坡运动计算剖面（图4），两侧滑体体积均为2×106m³，滑动方向与对面岸坡法线夹角约18°。根据断面形状，将整个滑坡体分为14个条块，每个条块宽度为ΔL=17.50 m，条块形状、条块编号、条块侧边线编号（图5）。

图4 卞家寨滑坡I-I断面Fig.4 Section I-I of Bianjiazhai landslide

图5 I-I断面分块图Fig.5 Block diagram of section I-I

计算选取滑坡体密度为2 100 kg∕m³，内摩擦角φ= 30°，黏聚力为C=20 kPa。同时为了分析滑坡体强度参数对计算得到的运动速度的影响，计算中考虑对滑坡体进行强度折减：Fos为折减系数，且（Fos＞1.0）。取不同折减系数，对应滑坡体滑动不同强度参数，从而得到滑坡体不同滑动运动过程：

2.3 涌浪传播计算参数

根据图4 断面Ⅰ-Ⅰ，计算得到滑坡体断面面积为9 667 m2，同时根据断面两侧滑体体积，计算得到断面两侧滑坡影响宽度均为316.40 m，故滑坡体宽度L=632.83 m。库水位为837.0 m，滑坡体平均厚度为λ=39.39 m（最大厚度56.94 m），断面当量水深为h=52.73 m。根据断面处水面宽度，及至坝址段水面形状，根据模型，坝址距滑坡距离为1 676.19 m，水面宽度为200.35 m，涌浪最大高度计算取宽度平均值B=220.60 m（滑坡下游至坝址处间隔100.00 m断面宽度平均值）。

根据图4、5 滑坡断面与库水位的关系，卞家寨滑坡入水产生涌浪需同时考虑水平运动与垂直运动的影响，按公式ahv+bλv'分别计算水平与垂直两部分的影响，其中为滑坡变形率，a为水平运动所占比例系数，b为垂直运动所占比例系数，根据滑坡变形计算中的最大水平速度与垂直速度占比，确定比例系数a，b。

根据初始涌浪高度选取不同级数项数按公式（6）计算得到滑坡对岸最大涌浪高度ζAmax，并分析级数项数取值对滑坡对岸涌浪高度的影响。在涌浪传播计算中，不同位置处河岸壁面粗糙程度不同，根据《建筑物的抗滑稳定分和滑坡分析》［4］中壁面反射系数的经验取值范围0.85～0.95，计算分析壁面反射系数取值对涌浪传播计算的影响。

3 卞家寨滑坡涌浪计算及敏感性分析

3.1 滑坡计算敏感性分析

滑坡运动过程计算，主要考虑对滑坡强度的不同程度折减，根据条分法计算，当滑坡体强度折减系数Fos（Fos＞1.0）为1.2时滑坡开始滑动，选取Fos从1.2 至3.0，步长取0.1 进行取值，计算得到不同折减系数下滑坡滑动结果（图6）。

图6 不同折减系数下滑坡滑动速度曲线Fig.6 Sliding velocity curve of landslide under different reduction coefficients

滑坡滑动入水的最大冲击速度及滑动时间T，将用于后续涌浪计算。从图6 可以看出，当对滑坡体强度折减系数为1.2时，滑坡体开始滑动，随折减系数增加，滑坡体滑动入水速度增加，符合实际规律。当折减系数为1.5 时，滑坡滑动入水时水平速度为6.14 m∕s，垂直速度为0.35 m∕s，合速度为6.15 m∕s。结合实际滑坡滑动过程中滑坡体前缘入水速度，滑坡体滑动阶段条分法计算结果与实际误差最小，并选择该计算工况作为参考工况，用于涌浪计算的敏感性计算分析。

3.2 涌浪计算敏感性分析

根据2.1 节计算结果，选取Fos=1.5 工况下滑坡滑动计算结果为参考工况，进行涌浪计算。根据滑坡变形计算中的最大水平速度与垂直速度占比，分别取a=0.95，b=0.05。按公式（3）～（5）计算得到涌浪初始高度为ζ0=16.34 m。

3.2.1 级数项数对涌浪计算影响

分析级数项数选取对涌浪计算的影响，参考工况下，滑坡滑动时间计算为9.76 s，波速=27.98 m∕s，根据表1级数项数取值为1，计算按公式（6）计算，滑坡对岸涌浪高度为13.87 m。根据表 1取值分析级数项数影响，级数项数选取从1～10，取值步长为1，对岸A点最大涌浪高度变化曲线如图7，随级数项数增加，滑坡对岸涌浪高度增加，增长速率降低；考虑对坝址处涌浪影响，根据壁面反射系数经验取值范围0.85～0.95，计算中取值选择中值k=0.90，计算得到下游不同断面涌浪高度，选择下游坝址处断面为典型断面，坝址断面距离滑坡发生位置1 676.00 m。不同级数项数下坝址涌浪高度变化曲线见图7，随级数项数增加，坝址处涌浪高度增加，增长速率降低，与对岸涌浪高度变化规律相同。

图7 涌浪高度随级数项数变化曲线Fig.7 Variation curve of surge height with the number of series terms

当项数取值大于5时，滑坡对岸涌浪高度级数第i项相对第i-1 项相对误差小于10%；当项数取值大于7 时，坝址处涌浪高度级数第i项相对第i-1 项相对误差小于10%，并趋于相对稳定，且计算涌浪高度较按照表1 取值均偏大较多。故在实际取值时，应根据实际滑坡滑动速度、时间以及水面宽度按照表1中公式T∕Δt计算，对级数项数进行取值。同时根据实际情况，只有对体量较大、滑动过程持续时间长、河面较窄的滑坡，才有可能选择较高的级数项数；对体量较小，滑动快速完成，且河面宽阔的情况，级数项数取值较小，取值为1 较为常见。从能量角度，体量越大，滑动时间越长，滑坡体动能传递至水体的能量越多，且河面越窄，涌浪传播中的能量耗散越少，涌浪高度越大。本文滑坡持续时间较短，滑动完成后，河面较宽，根据公式计算级数项数取1时为合理值。

3.2.2 壁面反射系数k对涌浪传播计算影响

分析壁面反射系数k，控制单一分析指标，基于参考工况，级数项数选取为1，壁面反射系数取0.85～0.95，取值步长为0.01，不同k值坝址涌浪高度变化曲线见图8。结果显示坝址处涌浪高度随壁面反射系数k增加而增大，且呈现线性相关的规律，k取中值0.90 时，坝址涌浪高度为1.094 m，在推荐取值范围k=0.85～0.95 之间，计算坝址涌浪高度为1.088～1.101 m，相较中值变动幅度在0.6%以内。综上，下游涌浪高度计算对壁面反射系数不敏感，在计算建议范围内取值时，对计算断面处壁面情况不明确时可取中值0.90；当明确计算位置壁面情况时，如遇地形变化平缓或光滑岩壁时可在0.85～0.90范围取值，如遇地形变化剧烈岸壁参差情况可在0.90～0.95 范围取值，由k值引起的误差较小。

图8 坝址处涌浪爬坡高度随壁面反射系数k变化曲线Fig.8 Variation curve of surge climbing height with wall reflection coefficient k at dam site

综上，基于潘家铮法计算涌浪高度时，对级数项数选取敏感，且不同取值差别较大，在取值时需根据表1中公式T∕Δt计算并根据表1中范围取值，同时考虑滑坡特点及河面宽度等特征，选择合适的项数；对壁面反射系数k不敏感，在建议范围0.85～0.95 之间取值时，引入误差减小。本文中根据滑坡滑动时间计算结果，选取级数项数为1，壁面反射系数取中值k=0.90 时，初始涌浪高度为16.34 m，对岸涌浪高度为13.87 m，坝址处涌浪高度为1.094 m。

4 结 论

基于潘家铮法对卞家寨滑坡滑动过程、涌浪传播过程进行了计算分析，并对级数项数、壁面反射系数k对涌浪传播计算的影响进行了敏感性分析，主要结论如下。

（1）基于潘家铮法估算卞家寨滑坡滑动过程，采用条分发思想计算，当滑坡体强度折减系数为1.2 时开始滑动，当为1.5时，滑坡体滑动入水水平速度为6.14 m∕s，垂直速度为0.35 m∕s，合速度为6.15 m∕s，滑动时间为9.76 s；

（2）涌浪计算中，根据滑动断面形状，同时考虑滑坡体的水平运动与垂直运动，且比例系数分别为0.95 与0.05，初始涌浪高度为16.34 m，基于表1，同时考虑涌浪传播的能量耗散，级数项数取值为1。滑坡对岸涌浪高度为16.12 m；坝址处壁面反射系数k取中值0.90，坝址处涌浪高度为1.094 m；

（3）基于潘家铮法的涌浪传播计算对级数项数的选取敏感，随着保留的级数项数增大，估算的涌浪高度增大，变化率减小，对计算影响较大。在取值时应根据滑坡滑动时间、入水速度、河面宽度计算T∕Δt，并参照表1 取值，同时需考虑滑坡及附近河面特点，只有对体量较大、滑动过程持续时间长、河面较窄的滑坡，才有可能选择较高的级数项数；对体量较小，滑动快速完成，且河面宽阔的情况，级数项数取值较小，取值为1 较为常见；

（4）基于潘家铮法的涌浪传播计算对壁面反射系数的选取不敏感，随k值增大涌浪高度呈现线性增加的规律，在建议范围0.85～0.95 内变化时与中值计算高度误差在0.6%以内，对计算结果影响较小。故当对壁面情况不明确时可取中值0.90；当、遇地形变化平缓或光滑岩壁时可在0.85～0.90范围取值，遇地形变化剧烈岸壁参差情况可在0.90～0.95范围取值。