考虑渗透系数各向异性的降雨入渗岩质边坡稳定性分析

沙 鹏 何群霄 刘冬艳

(绍兴文理学院 土木工程学院，浙江 绍兴312000)

0 引言

山体滑坡是一种破坏力强、影响范围广的地质灾害．造成山体滑坡的原因多种多样，其中水的作用不可忽视．在《中国典型滑坡》一书中列举出的众多滑坡实例，有90%以上的滑坡与降雨或地下水渗流有着紧密的联系［1］．不同于土质边坡，岩体边坡内部中的渗流主要由裂隙控制渗流大小与方向．水在裂隙岩体中的渗流情况比在孔隙介质中复杂得多，岩体结构的各向异性和非均质性以及裂隙分布的不连续性，导致了渗流的各向异性、不连续性和非均质性．因此研究分析岩体边坡中降雨入渗转化为地下水的过程并合理确定地下水位，对判断其稳定性有着重要意义．

近年来，有不少国内外学者采用数值方法对岩体边坡的地下水位和稳定性进行研究．李晓等［2］基于有限差分法研究库水位涨落与降雨联合作用下滑坡的地下水位变化;卢刚，尹小涛等［3－4］应用非饱和渗流理论，对降雨作用下的地下水入渗机制进行了探讨．一种简单有效的解析方法对实际工程的帮助也很大;国外Kacimov、Obnosov 和Warric 等［5－8］分别建立了两种不可压缩达西流体之间的界面突变的水力模型、瞬态过程的数学模型、三维拉普拉斯或平流色散方程场等压线模型(等压线的重构形状为贝壳形)，利用多种求解方式对浸润线的解析进行了系统研究;CHEN R H 等［9］研究了岩质边坡中裂隙水压力与稳定性的关系，提出降雨强度与时间历程是边坡失稳的最大影响因素．吴琼等［10］推导了库水位升降联合降雨作用下的浸润线公式，并通过工程实例进行了验证;王飞跃等［11］根据浸润线的影响因素，建立了浸润线矩阵，进行了拟合和求解;张利博等［12］依据同一条浸润线不同渗流截面的渗流量相等以及不同浸润线同一截面的渗流量差值为库水位的增量建立平衡方程，求得浸润线方程的解;刘燕等［13］引入了渗透力的计算方法，以黄蜡石边坡为例定量分析了地下水对边坡稳定性的作用．

然而目前的研究还存在诸多问题．首先，为方便浸润线计算往往假设岩体各向同性，这与实际情况存在较大差异．其次，在边坡稳定性验算中，未分清边坡岩体的透水特性进行水压力作用分析易造成计算混乱［14］．

针对上述问题，根据统计岩体力学理论对传统的浸润线方程进行修正，考虑了裂隙岩体中渗透系数的各向异性特征，探讨了降雨作用下的地下水渗流过程，以期为岩体边坡浸润线的研究提供一种新思路．在上述基础上，采用C#语言编写了基于不平衡推力法的边坡稳定性计算程序，并运用于实际工程进行验证．

1 裂隙岩体边坡的渗透性

岩体在生成过程中形成了各种结构面，水在岩体中渗流的实质是在裂隙网络中的渗流．因此边坡岩体的渗流特性主要取决于岩体结构特性，特别是结构面的产状、张开度和水力学连通性．

岩质边坡的渗透性可用裂隙岩体的渗透系数来表示．岩体裂隙水渗流服从达西定律，基于统计岩体力学方法［15］，采用等效连续介质模型推导了渗透张量计算公式，表示裂隙岩体各向异性的渗透性．

渗透张量一般可表示为:

式中Kf为沿裂面的渗透系数，v 为运动黏度，t 为隙宽，n 为结构面法线的方向余弦，δij为单位张量．

r 为结构面的视平均间距最小值，由渗透性可叠加原理，渗透张量可表示为，

式中可见岩体的渗透张量与裂隙张开宽度的三次方成正比，说明岩体的渗透性能对隙宽更为敏感．另外n 为结构面法线的方向余弦，因此岩体的渗透性呈现显著的各向异性．

由于边坡卸荷的定向性，一般会导致边坡上部张应力区中陡倾坡内或坡外的结构面张开，这将大大提升岩体该方向的渗透性能，显著加剧岩体渗透性能的各向异性．此外，由于边坡应力场的特殊性，裂隙张开度随卸荷变形程度向坡内递减，导致岩体渗透性能也向坡内减弱．为简化计算过程，本文对裂隙张开度变化不做详细讨论．

2 降雨入渗对地下水位的影响

随降雨强度的不同，边坡区域会发生不同程度的降雨入渗．入渗的水量将转化为地下水流，改变渗流场和地下水位，从而影响岩体的变形与稳定性．在边坡岩体中，降雨入渗和转化的地下水渗流是两个不同的阶段，分别具有不同的特征规律．下面分别讨论降雨入渗和地下水渗流．

2.1 降雨入渗

对于裸露的基岩，降雨沿岩体裂隙网络入渗，毛细作用和基质吸力等因素的影响可以不考虑．入渗量主要受降雨强度、降雨历时、岩体入渗能力等因素的综合影响．

现设边坡岩体的入渗强度为W，降雨强度为F，降雨强度按日降雨量确定，例如，暴雨的降雨强度为50～100 mm/d．当F＞W 时，降雨不能完全入渗，将在坡面形成径流;当F＜W 时，降雨将全部入渗．我们近似认为边坡岩体中的入渗沿铅直方向进行，铅直方向水力梯度J3=1，而其他方向为0．于是根据Darcy 定律和式(3)，饱和入渗强度为:

因此，边坡中的降雨入渗强度可写为:

2.2 地下水渗流

上述入渗水量将转化为地下水，水力坡度的存在使地下水在边坡岩体间也存在着渗流．由于岩体渗透张量的各向异性，降雨入渗方向和地下水流动方向的渗透性能是不同的．根据《地下水动力学》［16］中的浸润线方程

可知，当K 越小，降雨入渗作用下的浸润线变化越大．这个K 的实质是水力坡度方向上的渗透系数KJ．上述现象产生的原因为KJ越小，地下水在边坡岩体间渗流越困难，地下水位越易抬高．

综上，由于岩体渗透张量的各向异性，我们有必要对浸润线方程进行修正．

2.3 浸润线方程修正

现以图1 岩体边坡为例，做如下基本假设:

(1)对称的存在两条完整切割潜水含水层的平行河流，令h2处为分水岭．

(2)潜水含水层均质，各向异性，隔水底板水平方向不透水．

(3)为方便计算为稳定渗流，我们假设潜水流为一维流．

(4)河水位保持不变．

取水力坡度方向上的渗透系数KJ．分析地下水的浸润曲线．数学模型为:

图1 边坡降雨入渗示意图

对式(7)做定积分，得

由于h2实为分水岭位置，根据分水岭定义，浸润线方程可表示为:

考虑到降雨入渗强度W 式(5)可将上式分别写为:

上列各式中，水力坡度方向上的渗透系数可近似取横向渗透系数，即可参照式(4)的方法取KJ=K11．

以大渡河瀑布沟水电站经变形松动的变质玄武岩边坡的渗透系数为例进行计算．现取东西向的渗透系数K11为水力坡度方向的渗透系数KJ，K11为4.55 m/d，K33为14.68 m/d，K为9.95 m/d，取h1为20 m，h2为50 m，l 为500 m，F 为50 mm/d，F 明显小于K33．具体计算结果如图2 所示．可见两条浸润线存在明显的差异，因此根据岩体的各向异性渗透系数对原浸润线公式进行修正是有必要的．

3 裂隙岩体边坡稳定性分析

目前用于定量分析边坡稳定性的方法主要为极限平衡法和数值分析法．对于实际工程而言，刚体极限平衡法由于其技术成熟且精度满足要求，是目前规范中评价边坡稳定的主要方法．

图2 修正浸润线和原浸润线计算结果对比

这里主要研究裂隙化程度较高的岩坡，故将选用不平衡推力法来分析岩坡的稳定性［17］．不平衡推力法是我国独创的边坡稳定分析方法，适用于任意形状的滑坡．现根据上文计算的浸润线方程，在不平衡推力法中，采用水土合算的形式加入水的影响．条块受力分析见图3．

图3 条块受力分析图

不平衡推力法的计算公式为:

式中Pi为条块的剩余推力;Ti为下滑力;Ri为抗滑力;Fs为稳定性系数;ψi为力的传递系数;Fi－1为剩余推力．其中相关参数计算公式如表1 所示．式中αi为条块的滑面倾角;βi为水位线坡度角;φi为岩体摩擦角;ci为岩体黏聚力;li为各条块底边长;PGA 为地震系数．

对于(11)的计算公式通常有两种解法，即强度储备法和超载法［18］．强度储备法需要进行迭代计算，为优化运行速度，采用超载法进行程序编写．超载法的核心计算公式［19］为:

在超载法中，式(11)右边需要乘以Fs，其中，

基于上述理论和公式，本文采用C#语言进行了降雨入渗作用下岩坡稳定性计算程序的编写．程序界面如图4 所示．图中三线从上到下依次为坡面线、水位线和滑面线，色块从蓝色渐变为红色代表边坡稳定性降低或条块下滑力增大．

图4 降雨入渗岩坡稳定性计算程序

该程序的主要功能包括:(1)通过读取Excel的数据实现任意条分下的稳定性计算，同时也可修改窗口中的数据进行计算绘图;(2)直观的体现了降雨作用下水位线的变化情况;(3)通过色块图能直观了解各条块的稳定状态，点击条块能获得该条块处的剩余推力和稳定系数．

4 工程实例

本文以江西省萍乡市山口岩大坝附近一处岩质边坡作为研究对象．利用实测浸润线对浸润线的近似解析解进行验证，并使用自编程序进行了边坡的稳定性分析．此外进一步分析了降雨强度、河水位、渗透系数、内聚力和内摩擦角等因素对地下水位及边坡稳定性的影响．

4.1 滑坡概况

该滑坡位于山口岩大坝下坝址轴线上游400 m 左右的右岸，为一处古滑坡体，该古滑坡体在坡面上的出露面积约为28 350 m2，其平均厚度约为19 m，方量约为5.3×105m3．该滑坡体分布高程为190～340 m，地形边坡坡角在316 m高程以下为30°～35°，在316 m 以上地形坡度变缓，坡角为15°～25°．滑坡体地层岩性分布主要分为三层，表层为夹碎石、块石粉质黏土;中间岩层为千枚岩;基岩为千枚状变质砂岩．其剖面形状如图5 所示．

图5 古滑坡体剖面图

该滑坡体所处位置距离大坝较近，在降雨条件下，其各岩层的物理力学性质会发生变化，可能会导致二次滑坡的发生．因此须对该滑坡体进行降雨条件下的稳定性分析，验证其是否稳定．

4.2 实测浸润线对近似解析解的验证

本文根据当地水文站提供的降水资料，采用历年实测出的最大三天的降水量情况和地下水位线，对浸润线的近似解析解进行验证分析．最大三天的降水量为289.6 mm，因此取降雨强度为96 mm/d，根据式(3)，计算的岩体渗透系数K11为0.72 m/d，K33为0.86 m/d，降雨强度明显小于竖向渗透系数K33，因此取W 为96 mm/d．

图6 为180 m 水位实测浸润线与近似解析解预测结果的对比图．结果表明两者基本吻合，现将对比结果分为三段进行分析．在靠近岸坡的部分，预测结果略低于实测值，这是因为在预测浸润线计算过程中假定河水位高度不变，但实际降雨过程中河水位存在抬高;在滑坡的中间部分，预测结果略高于实测值，这是因为预测浸润线计算过程中假定坡体为裸露的岩体边坡，但实际该坡体上部覆盖了一层厚度约2 ～4 m 的土体，存在部分上层滞水;而在远离岸坡的部分，预测结果略低于实测值，这是因为边坡应力场的特殊性，岩体渗透性能随卸荷变形程度而向坡内减弱．

图6 实测浸润线与预测浸润线的对比

4.3 边坡的稳定性计算

获得合理浸润线的目的是为了能更准确地判断边坡的稳定性．根据上述工况信息以及该滑坡体的剖面几何形状和物理力学参数(见表2)，采用上文提及的基于不平衡推力法的稳定性计算程序进行该边坡的稳定性分析．

现根据滑坡体的几何参数按20 m 一段进行条分，条分数为14，输入上述计算参数，程序运行结果如图7 所示．

在上述条件下，根据程序运行结果可知在条块11、12 处的剩余推力最大，这与文献中的数值模拟结果是一致的．同时，程序运行计算得到的边坡稳定系数为1.220，小于安全系数1.250，这与文献［20］中利用修正的极限平衡法求得的稳定系数1.212 和考虑各向异性的数值模拟求得的稳定系数1.231 极为接近，该边坡属于欠稳定状态．上述结果说明本文提出的稳定性计算程序在便捷实现的同时，满足工程计算要求．

表2 计算参数赋值列表

图7 边坡稳定性计算程序运行结果

4.4 边坡稳定性的影响因素分析

从上文中可知程序计算涉及大量参数，这些参数不同程度上影响浸润线的变化和稳定系数的大小．现以上节中的计算参数为基准，以降雨强度、河水位、横向渗透系数、内聚力、内摩擦角为例，通过单独改变某一参数数值，分析其对浸润线和边坡稳定性的影响．对比结果见表3．

根据上述数据进行各参数对边坡稳定性影响的敏感性分析．采用参数变化量/原参数值表示参数变化程度，采用稳定系数变化量表示稳定系数变化程度，分析见图8．

根据表3 和图8，可知降雨强度、河水位、横向渗透系数主要通过改变边坡地下水位来影响坡体的稳定性．其中，降雨强度、河水位变化与稳定系数变化呈负相关，随着降雨强度增大、河水位抬高，地下水位越高，坡体稳定系数越小，边坡越不稳定．而横向渗透系数变化与稳定系数变化呈正相关，横向渗透系数增大，意味着水力坡度方向上的水力通道较为通畅，降雨入渗的水体在边坡内易于流动，致使地下水不易抬高，即横向渗透系数越大，降雨作用下的地下水位抬高越小，坡体稳定系数变化越小，边坡越稳定．这三个参数对边坡稳定性影响的敏感性为渗透系数＞河水位＞降雨强度．

表3 5#条块参数影响对比

图8 各参数对边坡稳定性影响的敏感性对比

而内聚力和内摩擦角代表的是边坡岩体的岩性，且都与稳定系数呈正相关，即内聚力、内摩擦角越大，表示岩体岩性越好，坡体稳定系数越大，边坡越稳定．这两个参数对边坡稳定性影响的敏感性为内聚力＞内摩擦角．值得一提的是，根据图3 中的条块受力分析图，可知在计算中内聚力和摩擦角的作用受到滑动面角度的影响，当滑动面角度变小至一定值时，对边坡稳定性影响的敏感性为内摩擦角＞内聚力．

5 结论

(1)从渗流基本原理出发，根据统计岩体力学理论推导了各向异性渗透系数的计算公式，探明了降雨入渗及地下水径流两个阶段中渗透系数起到的作用，表明采用各向异性的渗透系数更符合实际情况．

(2)根据各向异性的渗透系数对浸润线方程进行修正，并基于不平衡推力法采用C#语言进行了边坡稳定性计算模块的编写．以山口岩大坝附近一处岩质边坡为例，证明了浸润线的近似解析解的合理性以及稳定性计算程序的可应用性．

(3)在工程实例的基础上进一步分析了降雨强度、河水位、渗透系数、内聚力和内摩擦角等因素对地下水位及边坡稳定性的影响，并进行了敏感性分析，结果为内聚力＞内摩擦角＞渗透系数＞河水位＞降雨强度．

(4)本文引入了渗透系数的各向异性，在一定程度上使计算结果接近于实际岩体状态．然而，在推导过程中的部分假设，如不考虑坡体风化程度、潜水渗透为一维流等仍偏于理想．由于稳定流时的边坡稳定性更趋于坡体极限状态，现实降雨条件下边坡基本处于非稳定流状态．因此，为克服工程应用中的局限性，该状态下的浸润线计算有待进一步探讨研究．