蓄水后典型土质岸坡降雨与地震稳定性区域分析

汪功茂，石安池，陈建平，周燕国，汪玉冰*

（1.浙江大学 软弱土与环境土工教育部重点实验室，浙江 杭州 310058；2.浙江大学 超重力研究中心，浙江 杭州 310058；3.中国电建集团华东勘测设计研究院有限公司，浙江 杭州 311122；4.浙江华东建设工程有限公司，浙江 杭州 310014）

0 引 言

我国水资源丰富，西南山区河流落差大，开发利用河流资源可以实现当地社会经济高质量发展，同时推动国家碳达峰、碳中和目标的实现。但是水电的大力开发建设，附带对人居环境的影响也越来越大[1]。以金沙江白鹤滩水电站为例，其位于金沙江下游的攀枝花至宜宾河段，坝址距上游巧家县45 km。该水电站蓄水后，淹没巧家县沿江一带岸坡，库水会对县城沿江地带土质岸坡产生较大的影响。同时巧家县地属亚热带与温带共存的高原立体气候，近年降雨量呈明显上升的趋势[2]。该地区地质构造背景复杂，地震活动频繁，先后发生了2014 年8 月鲁甸MS 6.5 及2020 年5 月巧家MS 5.0 破坏性地震[3]。这些因素将影响岸坡的稳定性，进而影响到岸坡附近移民安置工程的安全。

国内外学者借助地理信息系统（GIS）进行不同工况下的区域边坡危险性研究，技术相对成熟，应用较为广泛。王萌[4]运用GIS 软件计算三峡库区水位上升及下降到某一高度时的水库岸坡稳定性系数，对区域尺度的水库滑坡危险性进行了预测研究。BAUM 等[5]基于数字高程模型（DEM）耦合降雨入渗导致的瞬态孔隙压力变化模型和无限边坡模型，提出了定量评价降雨致浅层滑坡失稳的TRIGRS 模型。MARIN 等[6]则运用该模型研究哥伦比亚热带山区盆地降雨引起的滑坡失稳，成功估算降雨强度和持续时间阈值。刘甲美等[7]基于25 m 分辨率DEM 数据，考虑潜在地震影响，提出了概率性地震滑坡危险性区域划新方法，绘制研究区地震滑坡高危险性边坡分布地图。尽管通过GIS 进行区域边坡危险性分析已经有了较为丰富的研究成果，但是不同工况的危险性发展演化规律还少有文献提到。

无限边坡模型作为相对简单的极限平衡分析方法，在区域边坡危险性分析中应用较为成熟[4-9]。运用该模型计算稳定性系数时，土层厚度通常根据研究区域的地质条件及历史滑坡数据进行粗略估计。如JIBSON 等[10]选取1994 年加州北岭地震引发滑坡典型厚度，绘制Oat Mountain 地区的地震滑坡概率地图；CHEN 等[11]选取2014 年鲁甸地震浅层滑坡的常见厚度，对震中研究区域进行危险性区划分析；RODRÍGUEZ-PECES 等[12]基于2011年洛尔卡地震数据确定土层厚度参数的4 种可能取值，进行地震滑坡危险性概率研究；WANG 等[13]参考甘肃天水地震历史滑坡数据，选用典型土层厚度研究不同Arias 强度下的Newmark 位移。在计算区域内不同位置处，选用相同的土层厚度对边坡稳定性系数进行计算，会导致计算结果出现偏差。

本文以金沙江白鹤滩水电站库区巧家县城北门及附近库岸为研究对象，基于研究区域的钻孔勘察数据，运用普通克里金插值方法建立GIS地层模型，相较于以往方法，能够更准确地反映不同位置处边坡地层厚度的实际分布。随后运用无限边坡分析方法，计算该区域在水库蓄水前、蓄水后、蓄水后强降雨以及蓄水后地震4 种工况下岸坡的稳定性系数。通过岸坡稳定性系数的前后比较，绘制出不同工况下的危险性变化地图，分析土质岸坡的危险性演化规律，可为白鹤滩水电站蓄水后巧家县城附近库岸边坡防护工程提供参考依据。

1 无限边坡模型稳定性分析方法

1.1 蓄水工况稳定性分析方法

对于浅部土体由多层土组成的边坡，各层土的抗剪强度与层厚各不相同，无法基于经验确定其潜在滑移破坏的位置。因此采用无限边坡模型（图1）计算该边坡的稳定性系数时，应首先假定滑移破坏面位于浅部各层土底部。取滑移面位于不同土层底部稳定性系数FSi的最小值作为边坡的天然状态稳定性系数FSstatic。如图1 所示，在水库蓄水前，浅部各土层均处于地下水位以上。

图1 无限边坡计算模型Fig.1 Infinite slope calculation model

式中：ci、ϕi、γi分别为第i层土的黏聚力、内摩擦角、重度；Zi为第i层土顶部到底部的深度；θ为坡度；Zicosθ为第i层土层厚度。

水库蓄水后，对于完全处在库水位以下的边坡，其稳定性受到自身重力与水压的双重影响；对于处在蓄水位以上的边坡，不考虑地下水的影响，仍按天然状态进行计算。分析蓄水工况下的边坡稳定性，首先对浅部各个土层进行蓄水后稳定性系数FSwi计算，最后取FSwi的最小值作为边坡的蓄水工况稳定性系数FSwater：

1.2 降雨工况稳定性分析方法

（1）降雨影响浅部土体稳定性分析方法

瞬态降雨入渗区域边坡稳定模型 TRIGRS（Transient Rainfall Infiltration and Grid Based Regional Slope-Stability Model）是基于GIS 栅格的降雨诱发型边坡稳定性计算程序，其包括入渗、水文、边坡稳定3 个模块。入渗模块采用有透水边界模型，根据边坡土体的初始入渗速率、水力扩散系数、饱和渗透系数等条件，计算不同降雨持续时间、不同深度位置处的地下水压力水头ψ；水文模块则采用质量守恒模型计算降雨所产生的地表径流，该模块运行过程中土体达到饱和时，无法入渗的降雨将全部转化为地表径流，直接流向运算区域下游的其他栅格；边坡稳定模块耦合地下水压力水头的无限边坡稳定性分析公式。综上降雨工况下浅部土体稳定性系数计算公式为：

式中：Z1为浅部土体的深度；c1ʹ、ϕ1ʹ、γ1分别为浅部土体的有效黏聚力、有效内摩擦角、重度；t为降雨持续时间；γw为水的重度。

（2）降雨影响深部土体稳定性分析方法

场地局部存在软弱夹层土时，强降雨过程中，雨水会通过优势通道渗流，使软弱夹层土达到饱和状态。土体强度将会下降，造成边坡安全性降低。当软弱土层位于地下水位以上时，该工况下边坡稳定性系数FSr计算公式为：

中粮集团党组书记、董事长赵双连表示：“按党的十九大提出的构建现代农业产业体系、生产体系、经营体系的要求，中粮正在改变过去坐在家里收粮当‘坐商’的模式，加快从传统农贸企业转变为新型生产服务经营主体，激发农村发展新活力。”

1.3 地震工况稳定性分析方法

地震工况边坡稳定性分析中，拟静力法在工程中应用较多，可用于绘制滑坡危险性地图。该方法将地震动瞬时荷载施加在潜在滑体重心上，基于静态极限平衡理论，确定边坡拟静力稳定性系数[14-15]。对于处在库水位以上的边坡，分析地震工况的边坡稳定性，首先对浅部各土层进行地震工况稳定性系数计算，最后取最小值作为库水位以上边坡的地震工况稳定性系数FSeq：

式中：ci、ϕi、γi分别为第i层土的黏聚力、内摩擦角、天然重度；Zi为第i层土顶部到底部的深度；α为顺坡向地震动加速度峰值。处于库水位及地下水位以下的边坡进行地震工况下稳定性分析时，各土层物理力学参数应采用有效黏聚力、有效内摩擦角、浮重度。

2 白鹤滩库区典型土质岸坡稳定分析

2.1 研究区概况

研究区域为库区云南省巧家县城的北门场地，坐落在金沙江右岸、巧家县城所在的不规则冲洪积缓坡台地，地势东高西低，东面为整体呈南北走向的山体，西面为金沙江。该区域地理范围为东经102°54′～102°56′、北纬26°54′～26°56′，总面积约为8 km2，最低海拔645 m，最高海拔942 m。

研究区域工程地质平面图的底图为比例尺达到1∶2 000 的地形图，在ArcGIS 软件中对其进行地理处理，提取等高线信息及高程点信息，可转化为栅格分辨率为5 m 的高精度数字高程模型（DEM）。如图2 所示，利用ArcToolbox 表面分析工具对DEM 进行处理，得到山体阴影地图，北门场地的地形特征得到良好可视化，同时将典型研究区位置标注其中。北门场地以石灰窑沟为界，划分为南区、北区两片区域。靠近东部巧家县城区域已设置防洪堤，地形以及土体分布有所改变，该区域前期工程建设已完成精细化数值分析，本次研究重点关注天然边坡区域。钻孔在研究区域有大量分布，主要集中在场地东侧的巧家县城附近，钻孔总数为146 个，其深度多分布在20～50 m，部分钻孔深度达80 m 以上。白鹤滩水电站水库死水位765 m，到正常蓄水位825 m 时，将淹没上游北门场地的大部分计算区域，即图2 中蓝色区域。

图2 山体阴影地图Fig.2 Mountain shadow

采用ArcGIS 表面分析工具对DEM 进行坡度分析，得到研究区域的坡度栅格图，如图3 所示。北门场地绝大部分区域坡度比较平缓，坡角小于10° 的区块约占75%，大于20° 的区块约占11%，大于30° 的区块约占4.5%，大于40° 的区块约占3%，其中北区局部存在坡角大于50° 的风化岩层。

图3 坡度图Fig.3 Slope map

2.2 土层分布及力学参数

整理钻孔勘探结果，各土样分布情况有所不同，以下4 层土在场地浅层分布相对广泛，其埋深按照土层编号自上到下：第①层含砾黏土在场地大部分钻孔中均有揭露，在北区局部厚度较大；第②层碎石混合土主要分布在场地泥石流沟底部及沟谷两侧；第③层有机土主要分布在北区的缓坡凹地、垭口等区域；第④层粉土质砂在场地南区分布较普遍，在场地北区小范围分布。除此之外，场地深部为强度高、厚度大的碎石土，因此本文所采用的无限边坡模型假定的滑动破坏位于上方强度较低土层。

场地部分土体为局部分布。编号⑤粉质黏土为区域典型软弱夹层土，在北门南区中部有分布。如表1 所示，该土体在含水量增加后其抗剪强度会急剧下降，黏聚力值由30 kPa 下降到15 kPa。编号⑥洪积碎石土主要分布于大桥沟、窝塘头沟、石灰窑沟等泥石流沟底及沟谷两侧附近；编号⑦堆积碎石土主要分布在窝塘头沟附近的岩质边坡浅层；编号⑧冲积黏土主要分布于南区最南侧坡度较大的山坡。

表1 研究区土层物理力学参数及分布情况Table 1 Physical and mechanical parameters and distribution of soil samples in the study area

本次研究各土样的物理力学参数，来自于华东勘测设计研究院前期进行的室内土工试验以及原位测试，主要包括三轴试验、直剪试验、十字板剪切试验、标准贯入试验等。综合上述试验结果以及地区经验值，最终提出北门场地土体的物理力学参数地质建议值，如表1 所示。

2.3 地层模型构建

ArcGIS地理统计分析（Geostatistic Analysis）模块普通克里金插值可基于钻孔中土层厚度数据的空间相关性构建场地全域的地层模型[16]，进而确定各栅格单元边坡计算时所选取的滑体厚度，同时该方法被验证可以对场地地下水位分布进行空间插值[17]。

普通克里金插值常用的半变异函数有球形模型、指数模型、高斯模型等，交叉验证模块可计算不同半变异函数的插值精度。即假设某一个钻孔的要素值未知，通过其余位置处的数据得出该钻孔的预测值，然后轮流改变钻孔，最后计算所有钻孔的实际观测值与估计值的各项误差。通过比较误差值，选用适合普通克里金插值的半变异函数模型。首先在图2 所示的ArcGIS 钻孔点要素中关联各钻孔点位置处的土层厚度及地下水埋深数据，然后分别采用上述3 种不同的半变异函数进行普通克里金插值，通过交叉验证方法来确定精度最高的插值方法。以第①层含砾黏土为例，其交叉验证结果如表2 所示。计算场地各土层深度及地下水埋深插值的交叉验证数值，根据误差判别标准，最终确定指数模型作为北门计算区域克里金插值所选用的半变异函数。在ArcGIS 地理统计分析模块中，采用指数模型普通克里金插值，可得出浅部4 层土样的深度以及地下水深度插值结果，图4、图5 分别为第①层、第②层土深度（层顶到层底竖向距离）插值地图。

表2 第①层含砾黏土克里金插值交叉验证表Table 2 Cross validation table of Kriging interpolation for pebbly clay layer ①

图4 第①层土深度插值地图Fig.4 Interpolation of the thickness of layer ①soil

图5 第②层土深度插值地图Fig.5 Interpolation of the thickness of layer ② soil

将插值地图与2.2 节中土层实际分布对比，在钻孔分布密集区域或者没有明显的山体、山沟区域，即使钻孔数据量不足，插值结果也与实际相符；在大桥沟、石灰窑沟、窝塘头沟及其附近山坡、北门南区山坡，根据克里金插值出来的结果并未考虑到岩土体的侵蚀搬运作用，与实际分布不符。对表1 中所列局部分布4 类土，结合具体地质情况确定其分布位置及深度。地下水插值结果表明86%计算栅格的地下水位深度大于15 m，绝大部分浅层土体位于地下水位以上。本文假设蓄水前土体为天然状态，边坡稳定性分析采用天然状态指标。

2.4 降雨、地震条件

巧家县地属亚热带季风气候，年内降雨量差异显著。查阅当地气象台站1981—2020 年降雨资料，巧家县多年平均降水量为844 mm，历年最大单日降水量为109.9 mm。基于上述数据，TRIGRS模型所选取降雨持续时间为24 h，每小时降雨量为4.6 mm。已有边坡降雨入渗研究表明，降雨主要对浅层土体造成影响，均质岩土体内雨水入渗深度有限[18-19]，一般不超过5 m。北门浅部土层为第①层含砾黏土，设定为TRIGRS 模型计算土层，如图4 所示局部区域深度大于5 m，因入渗深度有限，取5 m 进行计算，该模型其他计算参数详见表3。当计算区域存在优势渗流通道时，强降雨会造成深部软弱黏土强度降低。本文首先根据钻孔勘察数据确定软弱夹层即表1 中编号⑤粉质黏土分布范围，随后通过普通克里金方法插值出钻孔半径100 m、埋深25 m 以内的软弱夹层土的深度栅格，进行后续稳定性系数计算。

表3 TRIGRS 模型输入参数Table 3 TRIGRS model input parameters

巧家县城具有比较高的区域地震活动背景，在中国地震动峰值加速度区划图[20]中，巧家北门研究区域位于0.2 g 区内，白鹤滩水电站坝址靠近0.2 g 区与0.3 g 区的分界线。本文地震工况分析较为保守，采用白鹤滩水电站场址地震动数据，即超越概率为50 年10%的基岩水平向峰值加速度0.213 g 用于地震工况分析。同时考虑地震动放大效应，采集场地土样进行土质边坡离心机振动台试验以确定研究区地震动放大系数。当振动台面加速度为0.2 g 时，距基岩不同高度的传感器所监测到的峰值加速度放大系数如表4 所示，基于线性拟合的方式，可以得到地震动加速度放大函数。区域分析中，取深部碎石土（图1）顶部加速度为基岩峰值加速度，基于上述加速度放大函数，最终可确定滑体顺坡向地震动加速度峰值。

表4 峰值加速度放大系数Table 4 Amplification coefficient of peakacceleration

2.5 计算结果分析

在ArcGIS 中，将每一个栅格视为独立的坡体，对坡度、土体物理力学参数进行数字化处理，通过栅格计算即可获得稳定性系数，并绘制地图。图6 为天然状态边坡稳定性系数图，灰色栅格稳定性系数大于2，区域占比达到99%以上，为稳定区域；绿色栅格稳定性系数为1.5～2.0，属于较为稳定区域；黄色栅格稳定性系数为1.0～1.5，属于较高风险区域，表明在施加其他荷载之后边坡可能发生失稳；红色栅格主要分布在北门北区的山沟地区，稳定性系数低于1.0，属于高风险区域。

图6 天然状态边坡稳定性系数图Fig.6 Slope stability factor in natural state

图1 所示无限边坡模型假设中，忽略滑体上坡向土体对其施加的下滑力与下坡向土体施加的抗滑力，一般情况下抗滑力大于下滑力，导致稳定性系数计算结果偏低。稳定性系数随着计算厚度的增加而不断减小，当计算厚度为浅层时，下滑力与抗滑力差别不大，对计算结果影响较小；当计算厚度为深层时，下滑力远小于抗滑力，导致稳定性系数结果过于保守。本文重点评估不同工况下稳定性系数的变化规律，用统一的物理模型计算边坡稳定性系数，其变化值相较绝对值更具参考价值。因此图6 中红色栅格不代表边坡处于失稳状态，仅表明在其他因素干扰下，相较其他位置处更易发生失稳。本文重点关注区域为较高风险区域与高风险区域，即图6 中天然状态稳定性系数小于1.5 区域，研究其在不同工况下边坡危险性演化规律。

对重点关注区域绘制蓄水、降雨、地震工况下相较天然状态的稳定性系数变化图。图7 为蓄水工况稳定性系数变化图，蓝色部分为825 m 蓄水位以下区域，黑色部分为蓄水位以上区域，图中未折减栅格仅针对重点关注区域，非重点关注区域稳定性系数变化情况未进行分析。蓄水位以下，绝大部分重点关注区稳定性系数未降低，该区域采用浮重度计算下滑力，而土体抗剪强度下降不大，边坡稳定性上升。绿色栅格和黄色栅格表明蓄水后重点关注区域稳定性系数为蓄水前的75%～100%和 50%～75%，均有不同程度的下降。稳定性系数变化计算结果表明，随着蓄水位由765 m 逐渐上升至825 m，边坡稳定性系数下降区域增加，图7 中黄色栅格与绿色栅格范围相较765 m 蓄水位时增多，大桥沟与窝塘头沟相邻区域整体稳定性下降较为明显，图中红色栅格表明稳定性系数下降比例达到50%及以上，需针对该区域做蓄水工况下的进一步数值模型分析。

图7 蓄水工况稳定性系数变化图Fig.7 Variation of stability factor in water storage condition

降雨工况区域边坡的稳定性分析仅针对未蓄水区域。蓄水位825 m 时，计算面积较小，因此选用水库死水位765 m 进行降雨工况稳定性分析。综合浅部土体TRIGRS 模型以及深部软弱土计算结果，绘制图8 降雨工况稳定性系数变化图。降雨工况下，重点关注区域绝大部分稳定性系数未下降，北门南区局部出现稳定性系数降低25%～50%（黄色栅格），此区域存在深部软弱黏土夹层，降雨后边坡稳定性下降较大，需针对降雨工况进行数值模型分析。

图8 降雨工况稳定性系数变化图Fig.8 Variation of stability factor in rainfall condition

图9 为地震工况稳定性系数变化图，蓝色部分为825 m 蓄水位以下区域，黑色部分为蓄水位以上区域。在地震动作用下，北门重点关注区域稳定性系数都出现了降低，北区部分区域稳定性系数降低较为明显，其中局部稳定性系数值下降达到了50%以上（红色栅格），表明本区域在地震来临时，与其他区域相比存在较大的边坡失稳风险，需针对地震工况进行数值模型分析。

图9 地震工况稳定性系数变化图Fig.9 Variation of stability factor in earthquake condition

3 结 论

本文基于普通克里金插值构建区域地层模型，采用无限边坡模型计算岸坡在不同工况下的稳定性系数及其变化值并绘制相应地图，完成岸坡稳定性的区域分析，得到了以下结论：

（1）研究区域钻孔数据丰富时，本文采用指数函数模型的普通克里金方法所构建的地层模型精度较高，为岸坡稳定性分析提供了相对准确和比较可靠的数据。

（2）蓄水及降雨工况对于巧家县北门岸坡稳定性整体影响较小，含深部软弱黏土夹层岸坡在降雨后稳定性下降较大；地震工况下岸坡稳定性系数整体下降较大，局部下降比例达到了50%以上。同时，稳定性系数分布及变化地图可以甄别出蓄水、降雨及地震各工况下岸坡危险性相对较高区域，如巧家县北门北区局部位置，后续可针对相关区域做进一步数值分析。

（3）在岸坡稳定性的区域分析中，利用无限边坡模型获得的蓄水、降雨和地震各工况下的稳定性系数由于计算假设等原因可能偏低。然而在精细化数值模拟分析无法应用的大尺度区域分析中，采用统一的物理模型定量评估多工况下的边坡稳定性对于区域规划和风险研判有科学性、高效性和实用性等优点。