在热带和亚热带地区，大部分滑坡是由于降雨渗入风化的残积土斜坡所引起的^[1]。降雨入渗将增大边坡土体含水率，导致负孔隙水压力 (即土体基质吸力) 锐减，甚至产生正的孔隙水压力，降低土体抗剪强度^[2-3]，从而增大滑坡风险。

由于地质作用的差异以及气候条件不断变化，自然土体的孔隙结构具有空间时间变异性，是土体形成和发展过程中固有的特性，因此，采用可靠度分析方法进行边坡稳定性分析已成为岩土工程界的共识。土体是由不同固体颗粒组成的多孔介质，其孔隙结构与颗粒材料决定了诸如饱和渗透系数、黏聚力、内摩擦角等土体性质，而微观尺度上土体颗粒结构与孔隙分布的空间变异性也决定了宏观特性如饱和渗透系数等的空间变异性。通常，短期降雨条件下滑坡风险评估均假设土体饱和渗透系数不随时间变化，忽略其时间变异性，并采用随机场理论模拟分析土体饱和渗透系数的空间变异性。Zhu等^[4]采用快速傅里叶变换方法模拟渗透系数函数空间变异性，通过蒙特卡罗法 (MCS) 研究了中雨强度下饱和渗透系数的空间变异性对边坡孔隙水压力、地下水位及安全系数的影响规律。Cho^[5]采用Karhunen-Loeve (K-L) 展开方法离散土体饱和渗透系数随机场，通过MCS法研究了恒定降雨入渗条件下渗透系数的空间变异性对无限边坡体内孔隙水压力和稳定可靠性的影响。蒋水华等^[6]采用K-L展开方法离散土体饱和渗透系数、有效黏聚力和内摩擦角三参数相关非高斯随机场，通过稳定渗流分析和边坡稳定性分析，研究了坡面降雨强度、土体特性参数变异性和参数间互相关性对边坡体内孔隙水压力以及稳定可靠度的影响。上述文献均假设土体抗剪强度不随时间变化，侧重研究降雨条件下饱和渗透系数空间变异性对边坡体内水力响应及稳定可靠性的影响。事实上，土体抗剪强度随降雨过程中土体含水率的变化而变化，忽略抗剪强度的时间变异性对边坡稳定可靠性分析结论会产生较大影响^[7]。而现有文献对考虑空间变异性与时变性耦合的土体抗剪强度研究有所欠缺。

为此，文中将对降雨雨型、雨强和持时对边坡体内孔隙水压力分布的影响展开研究，并通过土体饱和渗透系数随机场的模拟，结合随机渗流分析，研究降雨入渗条件下饱和渗透系数的空间变异性对孔隙水压力、含水率、抗剪强度参数等空间分布的影响，从而建立局部安全系数的空间分布随机场，为进一步研究降雨条件下边坡稳定可靠性提供依据。

1 饱和渗透系数随机场模拟

一般地，饱和渗透系数y为服从对数正态分布的随机场，则饱和渗透系数对数样本空间服从正态分布N[μ_Y，σ_Y；λ_Y]，其中Y=lny，μ_Y和σ_Y为均值和标准差，可根据对数变换得到，λ_Y为自相关系数。根据Miller and Miller相似介质假设，认为所研究的多孔介质 (如土壤等) 为相似介质，具有相似的孔隙形状和分布或相似的水力函数，不同位置处的土壤性质可通过其内部几何特征尺寸的变化来实现，即不同位置处的土壤性质可通过对土壤性质的参考值进行缩放来简化地描述其空间变异性。因此饱和渗透系数对数Y总样本可表述为如下形式：

$ {{Y}_{ij}}={{\mu }_{Y}}+{{\sigma }_{Y}}\cdot {{\varepsilon }_{ij}}={{\mu }_{Y}}\cdot \left( 1+{{\varepsilon }_{ij}}\cdot {{\sigma }_{Y}}/{{\mu }_{Y}} \right)={{\mu }_{Y}}\cdot {{A}_{ij}}, $

(1)

式中：A_ij称为缩放因子，表示某点的饱和渗透系数对数值与饱和渗透系数对数均值的相对大小；ε_ij为服从标准正态分布N[0，1；λ_Y]的随机场。饱和渗透系数对数随机场的模拟就转化为对标准高斯随机场ε_ij的模拟，可采用如下基于相关系数与相关函数的谱密度函数方法：

$ \varepsilon \left( u \right)={{\left( \frac{2}{N} \right)}^{1/2}}\sum\limits_{i=1}^{N}{\cos \left( {{u}^{\text{T}}}{{Y}_{i}}{{W}_{i}}+\mathit{{{\Theta }}_{i}} \right)}, $

(2)

式中，N为谐波项数；u为欧氏空间坐标向量，对于二维空间u=(x，z)，(x_i, z_j) 为空间坐标；Y为单位球体表面上等分布的随机变量，对于二维空间为单位圆周上等分布的随机变量，故Y={cos γ，sin γ}^T，γ在[0，2π]均匀分布；W为与相关函数有关的随机变量，Θ为在[0，2π]均匀分布的随机变量；各随机变量相互独立。对于二维随机场，具体为

$ {{\varepsilon }_{ij}}=\varepsilon \left( {{x}_{i}}, {{z}_{j}} \right)={{\left( 2/N \right)}^{1/2}}\sum\limits_{m=1}^{N}{\cos \left[{{\omega }_{m}}\left( {{x}_{i}}\cos \ {{\gamma }_{m}}+{{z}_{j}}\sin \ {{\gamma }_{m}} \right)+{{\varphi }_{m}} \right]}。 $

(3)

由于土体水力参数的相关函数对计算结果影响不大^[9]，文中采用简单的指数函数形式ρ(τ)=e^-λ_Yτ，自相关系数λ_Y与波动尺度δ互为倒数关系。该相关函数对应的径向谱分布函数为

$ G\left( \omega \right)=1-1/{{\left[1+\left( {{\omega }^{2}}/\lambda _{Y}^{2} \right) \right]}^{1/2}}, 0\le \omega \le \infty 。 $

(4)

由于方程u=F(y)(u为[0,1]内均匀分布的随机变量序列) 的解{y}可组成累积分布函数为F(y) 的序列。故：

$ {{\omega }_{m}}={{\lambda }_{Y}}{{\left\{ {{\left\{ 1/\left[1-G\left( {{\omega }_{m}} \right) \right] \right\}}^{2}}-1 \right\}}^{1/2}}, $

(5)

式中，G(ω_m) 为在0~1之间均匀分布随机数。

当水平波动尺度与竖向波动尺度不同时，可以根据水平波动尺度与竖向波动尺度的比值对空间坐标进行缩放来考虑波动尺度的各向异性。生成Y_ij后经过指数变换即可得到满足对数正态分布的饱和渗透系数随机场。

2 降雨入渗和稳定性数值分析 2.1 渗流分析

降雨入渗情况下，根据达西定律和质量守恒定律推导多孔介质渗流控制微分方程：

$ \frac{\partial }{\partial x}\left[{{K}_{x}}\left( \theta \right)\frac{\partial H}{\partial x} \right]+\frac{\partial }{\partial z}\left[{{K}_{z}}\left( \theta \right)\frac{\partial H}{\partial z} \right]={{m}_{\text{w}}}{{\gamma }_{\text{w}}}\frac{\partial H}{\partial t}, $

(6)

上式即为Richards方程^[8]。其中，K_x(θ), K_z(θ) 为x，z方向的渗透系数，为依赖于空间位置的量，文中假设x，z方向的渗透系数相等，即K_x(θ)=K_z(θ)；H为总水头，即位置水头h_z与孔隙水压力水头h_w之和，其中h_w=u_w/γ_w，u_w为孔隙水压力，γ_w为水的重度，可设为常数；m_w为比水容量，即体积含水率与基质吸力水头曲线的斜率，可由土水特征曲线得到，如下所示van Genuchten-Mualem模型^[9]：

$ \theta \left( h \right)=\left\{ \begin{align} &{{\theta }_{\text{r}}}+\frac{{{\theta }_{\text{s}}}-{{\theta }_{\text{r}}}}{{{\left[1+{{\left| \alpha h \right|}^{n}} \right]}^{m}}}, h < 0 \\ &{{\theta }_{s}}, h>0 \\ \end{align} \right\}, $

(7)

其中：m=1-1/n，n＞1；h为基质吸力水头，空气气压力为零时，基质吸力水头等于负孔隙水压力水头-h_w；θ为土壤体积含水率；θ_s为土壤饱和含水率；θ_r为土壤残余含水率；α、n为形状参数，对于某一类型土体，可设α、n为常数。研究边坡土体为壤土，van Genuchten-Mualem模型参数采用文献[10]的土水特征曲线，如表 1所示。

渗透系数函数为

$ {{K}_{i}}\left( \theta \right)={{K}_{\text{s}}}{{\left( \frac{\theta-{{\theta }_{\text{r}}}}{{{\theta }_{\text{s}}}-{{\theta }_{\text{r}}}} \right)}^{\frac{1}{2}}}{{\left\{ 1-{{\left[1-{{\left( \frac{\theta-{{\theta }_{\text{r}}}}{{{\theta }_{\text{s}}}-{{\theta }_{\text{r}}}} \right)}^{\frac{1}{m}}} \right]}^{m}} \right\}}^{2}}, i=x, z, $

(8)

式中，K_s为饱和渗透系数，为服从对数正态分布的随机场。

当降雨强度大于土壤渗透系数时采用定水头边界条件，多余的雨水将形成地表径流从边坡表面排走；当降雨量小于土壤渗透系数时采用定流量边界条件，雨水将全部入渗边坡体内。坡面及坡顶处采用大气边界条件，模型两侧按零流量边界处理；底面为渗出面排水边界。初始条件为地下水位以上土体的基质吸力成线性分布，其中底面为0，顶面为-100 kPa，如图 1所示。忽略蒸发。

2.2 含水率对抗剪强度的影响

含水率的变化是导致土体抗剪强度参数粘聚力c′和内摩擦角φ′发生变化的最直接因素。国内学者通过试验研究了含水率ω与c′、φ′之间的关系，指出c′、φ′随含水率ω的变化表现出良好的指数关系^[11-13]。采用文献[14]中的试验数据，取样点位于重庆市主城区排水管道A管线187#黑石子段岸坡体内。c′、φ′随含水率ω的变化如图 2所示。经数值拟合，得到c′、φ′与含水率ω的关系式如下：

$ \left\{ \begin{array}{*{35}{l}} {c}'=2\ 629.3{{\text{e}}^{-0.217\ 9\omega }}, {{R}_{{{c}'}}}=0.984\ 9; \\ {\varphi }'=343.58{{\text{e}}^{-0.223\ 6\omega }}, {{R}_{{{\varphi }'}}}=0.938\ 0。 \\ \end{array} \right. $

(9)

其中，R_c′，R_φ′分别为2条曲线的拟合优度。

2.3 稳定性分析

以往边坡稳定性分析是在假设潜在完整滑移面的基础上计算边坡稳定安全系数，整个边坡仅有1个安全系数，无法描述滑坡的始发点。然而由于降雨条件和孔隙结构的空间变异性，坡体内每点的抗剪强度不相同，每点的相对安全程度也各有差异。文中采用文献[15]中提出了“局部安全系数”的概念来描述边坡体内每点的安全性，如图 3和式 (10) 所示。

$ \text{LFS=}\frac{{{\tau }^{*}}}{\tau }=\frac{2\cos \ \varphi '}{\sigma {{'}_{1}}-\sigma {{'}_{3}}}\left[c'+\left( \frac{\sigma {{'}_{1}}+\sigma {{'}_{3}}}{2}-{{\sigma }^{s}} \right)\tan \ \varphi ' \right], $

(10)

式中:τ^*为潜在的库仑应力; τ为当前状态的库仑应力；σ^s为吸应力，其表达式如式 (11) 所示；α、n为van Genuchten-Mualem模型参数；u_a、u_w分别为孔隙气压力和孔隙水压力。

$ {{\sigma }^{s}}=\left\{ \begin{align} & -\left( {{u}_{\text{a}}}-{{u}_{\text{w}}} \right),{{u}_{\text{a}}}-{{u}_{\text{w}}}\le 0; \\ & -\frac{\left( {{u}_{\text{a}}}-{{u}_{\text{w}}} \right)}{{{\left( 1+{{\left[ \alpha \left( {{u}_{a}}-{{u}_{w}} \right) \right]}^{n}} \right)}^{\left( n-1 \right)/n}}},{{u}_{a}}-{{u}_{w}}\ge 0。 \\ \end{align} \right. $

(11)

因此，基于随机场理论的降雨致滑坡稳定性分析步骤如下：

1) 输入饱和渗透系数的均值、方差以及波动尺度，采用上述基于谱密度函数的随机场生成方法，模拟饱和渗透系数随机场，得到饱和渗透系数随机场的一次实现；

2) 结合多孔介质渗流控制微分方程和土壤的土水特征曲线模型，求解降雨条件下边坡体内的孔隙水压力、含水率随机场空间分布样本；

3) 根据含水率抗剪强度参数经验关系曲线和步骤2) 所得含水率的空间分布样本，得到降雨条件下抗剪强度参数的空间分布样本；

4) 通过有限元方法求解在自重情况下边坡体内每点的应力状态，结合上文得到的抗剪强度参数空间分布，即可获得局部安全系数空间分布样本。

5) 重复步骤1)~4)，获得降雨条件下的边坡局部安全系数随机场。

3 降雨特性对孔压的影响 3.1 降雨雨型对孔压的影响

伴随着降雨的发生、发展和停止，不同降雨类型具有不同的雨量分配。假设48 h内总降雨量相同的情况下，即小雨、中雨、大雨、暴雨和大暴雨的总降雨量分别为10、35、75、150、350 mm，研究降雨雨型对孔隙水压力的影响。按降雨分配特点划分了3种降雨类型：均匀型、单峰型和阶梯型。三类降雨强度与时间的关系示意如图 4所示。

不同降雨强度三类雨型下边坡特征点的孔压如图 5所示。从图 5(a)可以看出，小雨时，均匀型降雨12 h导致边坡表面的孔压降低最大，单峰型和阶梯型降雨的孔压降低量相同且小于均匀型降雨；降雨24 h、48 h及之后，三类雨型下边坡表面的孔压降低相同；均匀型降雨36 h导致的孔压降低最小，单峰型和阶梯型降雨的孔压降低量相同且大于均匀型降雨，与12 h时恰好相反。这是因为小雨条件下，土体的饱和渗透系数大于降雨强度，雨水全部入渗边坡体内，边坡表面的孔压降低与累积降雨量直接相关，雨型对孔压的影响较大。随着降雨强度增大，三类雨型下边坡表面的孔压降低量差异逐渐减小，其原因在于当降雨强度较大时，土体的饱和渗透系数接近或小于降雨强度，此时边坡表面土体逐渐饱和，孔压接近于0，多余的雨水形成地表径流，从坡面排走，因此，三类雨型下边坡表面的孔压变化基本一致。对于边坡体内的孔压，当降雨强度较小时，降雨入渗速度较慢，三类雨型下边坡体内特征点处的孔压不变；随着降雨强度增大，降雨入渗速度逐渐增大，特征点处孔压发生变化的时间越来越快；除大暴雨0~24 h外，降雨雨型对边坡体内孔压影响较小，主要原因在于当小雨、中雨、大雨、暴雨时，降雨持时较小，降雨湿润锋还没有达到特征点位置，各雨型对特征点的孔压影响很小；随着降雨入渗进行，特征点处的孔压逐渐减小，但三类雨型下的孔压变化相差较小，这是因为三类雨型下实际入渗量均较小且相差不大的缘故；当大暴雨时，均匀型降雨24 h前导致边坡体内特征点处孔压降低最大，单峰型次之，阶梯型最小，24 h后三类雨型对孔压的影响基本一致。究其原因在于边坡体内的孔压变化主要取决于实际累积入渗量以及土体的饱和状态。当降雨24 h前，土体处于非饱和状态，均匀型降雨的实际累积入渗量最大，导致的孔压降低量也最大，单峰型次之，阶梯型最小。当降雨24 h后，土体接近饱和状态 (孔压为0)，此时不管何种降雨类型下，只要不断有降雨入渗，孔压将一直保持不变。

3.2 降雨强度对孔压的影响

主要考虑阶梯型降雨下5种降雨强度对边坡孔压的影响。根据我国气象部门^[16]所采用的降雨等级划分标准，分别以各区间中点值代表各种降雨强度，即0.020 8 cm/h (小雨)、0.072 9 cm/h (中雨)、0.156 3 cm/h (大雨)、0.312 5 cm/h (暴雨)、0.729 2cm/h (大暴雨)。5种降雨强度下边坡特征点的孔压如图 6所示。可见随着降雨强度增加，相同降雨持时下特征点处孔压降低幅度逐渐增大；当土体达到饱和状态后，降雨强度对特征点处孔压的影响较小。

3.3 降雨持时对孔压的影响

以均匀型降雨总降雨量为75 mm为例研究降雨持时对边坡孔压的影响 (见图 7)，从小雨到大暴雨所对应的降雨持时分别为：361、103、48、24、10 h。

从图 7可以看出，在总降雨量相同的情况下，降雨持时越大，边坡表面的孔压降低越少，体内的孔压降低越多。原因在于当总降雨量相同时，降雨持时较小，即降雨强度较大，边坡表面的雨水不能及时入渗，土体将处于饱和状态，多余的雨水将以径流的方式从表面排走，此时表面孔压接近于0。降雨持时越大，意味着降雨强度越小，边坡表面土体饱和度越小，故表面的孔压降低越少；当降雨强度小于边坡的入渗能力时，雨水全部入渗到边坡内部，导致边坡体内的饱和度增加，孔压降低越多。

4 饱和渗透系数的空间变异性对孔压和稳定性的影响 4.1 饱和渗透系数缩放因子随机场

由于自然界中土体的矿物构成、沉积条件、应力历史、风化以及其他地质作用的差异，土性参数具有明显的各向异性，同一土层内不同方向的波动尺度是不同的。通常，水平波动尺度比竖向大一个数量级^[17]。文中研究不同波动尺度对饱和渗透系数的影响。饱和渗透系数缩放因子变异系数假设为100%，波动尺度的选取如表 2所示。

采用上述随机场模拟方法生成饱和渗透系数缩放因子随机场如图 8所示，δ_h和δ_v分别为水平和竖向波动尺度，由红色过渡到浅蓝色表示土体渗透系数由大到小。图 8(a)~(c)为相同的δ_h下，不同δ_v对饱和渗透系数缩放因子的影响，随着水平波动尺度的增大，同一水平位置处的饱和渗透系数的波动性逐渐减少，饱和渗透系数的水平层状分布愈见明显，表明土体的分层分类情况。图 8(c)~(e)为相同的δ_v下，不同δ_h对饱和渗透系数缩放因子的影响，随着竖向波动尺度的减小，饱和渗透系数沿竖向的波动性逐渐增大，相同高度的边坡出现较小 (或较大) 饱和渗透系数的次数增多。

4.2 孔隙水压力及含水率分布

设边坡土体的饱和渗透系数均值为36 cm/h，研究大雨条件下 (持续240 h、10 d) 饱和渗透系数的波动尺度对孔隙水压力及含水率分布的影响，图 9(a)~(c)与图 10(a)~(c)分别为相同的δ_h下，不同δ_v对孔隙水压力及含水率分布的影响，图 9(c)~(e)与图 10(c)~(e)分别为相同的δ_v下，不同δ_h对孔隙水压力及含水率分布的影响如图 9和图 10所示。随着水平波动尺度的增大，边坡体内的负孔隙水压力和低含水率范围逐渐减小，还可能在坡脚附近形成一水平饱和区，这是由于竖向波动尺度不变时，随着水平波动尺度的增大，同一水平处相近两点的渗透系数表现出很强的相关性，雨水在边坡体内更容易流动；反之，当水平波动尺度较少时，同一水平处相近两点的渗透系数波动性较大，将阻碍水流的流动，雨水入渗的深度随之减小，边坡体内的负孔隙水压力的范围减小较小，在渗透系数较小的区域雨水积聚，导致含水率较大。从图 9和图 10可以看出，水平波动尺度对孔隙水压力及含水率的变异性影响小于竖向波动尺度，与阻水土坝相反；对于阻水土坝，其渗流过程主要是由高水位向低水位进行，水流以水平方向为主，水平波动尺度越大，同一水平处的饱和渗透系数越接近，有利于水流水平发展。而对于边坡降雨入渗来说，其渗流过程主要受重力的影响，水流以竖向方向为主，竖向波动尺度越大，同一竖向截面处的饱和渗透系数波动性较小，有利于雨水的竖向渗透。

4.3 抗剪强度参数的空间分布

根据含水率分布以及c′、φ′与含水率ω的关系式即可获得粘聚力c′和内摩擦角φ′的空间分布，以δ_h=0.5 m, δ_v=0.5 m的情况为例，如图 11所示。从图 11可以看出，抗剪强度参数的空间分布特点与含水率分布基本一致；由于降雨条件和孔隙结构的空间变异性的存在，土体抗剪强度将在边坡表面产生局部软弱区，对边坡稳定性极为不利。

4.4 边坡稳定性分析

分析结果如图 12所示，图 12(a)~(c)为相同的δ_h下，不同δ_v对局部安全系数分布的影响，图 12(c)~(e)为相同的δ_v下，不同δ_h对局部安全系数分布的影响。临界滑面 (即局部安全系数为1的等值线) 基本与坡面平行且深度较浅，与实际观测到的降雨致滑坡情况相符^[18-19]。这是由于降雨入渗导致与边坡表面等距离的土体的抗剪强度几乎同步下降，抗剪强度等值线几乎平行于坡面，而由自重产生的剪应力等值线也平行于坡面，当抗剪强度下降到不能抵抗当前剪应力时，将发生滑坡。随着水平波动尺度的增加，临界滑面深度随之增大，主要原因是当水平波动尺度较小时，渗透系数表现出很强的变异性，出现低渗透性的概率增多，雨水不能继续渗透，雨水积聚导致边坡土体抗剪强度局部降低，产生局部小滑坡。随着水平波动尺度增大，渗透系数表现出很强的相关性，雨水在边坡体内更容易流动，抗剪强度降低区域范围更广，形成整体大滑坡。随着竖向波动尺度的增加，渗透系数在一定范围内的相关性逐渐增大，有利于雨水持续向底部下渗，坡体上部积水较小，土体抗剪强度降低较小，总体而言，坡体上部局部安全系数随竖向波动尺度的增加而减小。

5 结论

文中基于土体饱和渗透系数随机场，通过随机渗流分析，研究了降雨入渗条件下饱和渗透系数的空间变异性对孔隙水压力、含水率、抗剪强度参数等空间分布的影响，以及降雨雨型、雨强和持时等降雨特性对孔隙水压力的影响，并结合应力场与抗剪强度参数的空间分布，建立了局部安全系数随机场，为进一步研究降雨条件下边坡稳定可靠性奠定了基础。

1) 当降雨强度较小时，边坡表面的孔压降低与累积降雨量直接相关，雨型对孔压的影响较大。随着降雨强度增大，三类雨型下边坡表面的孔压降低量差异逐渐减小。边坡体内的孔压变化主要取决于实际累积入渗量以及土体的饱和状态，当降雨强度较小时，三类雨型下孔压差异较小，随着降雨强度增大，均匀型降雨的实际累积入渗量最大，导致的孔压降低量也最大，单峰型次之，阶梯型最小。随着降雨强度增加，相同降雨持时下特征点处孔压降低幅度逐渐增大；当土体达到饱和状态后，降雨强度对特征点处孔压的影响较小。在总降雨量相同的情况下，降雨持时越大，边坡表面的孔压降低越少，体内的孔压降低越多。

2) 随着水平波动尺度的增大，边坡体内的负孔隙水压力和低含水率范围逐渐减小，还可能在坡脚附近形成一水平饱和区。受渗流主方向影响，降雨条件下水平波动尺度对孔隙水压力及含水率的变异性影响小于竖向波动尺度，与阻水土坝相反。

3) 临界滑面 (即局部安全系数为1的等值线) 基本与坡面平行且深度较浅。随着水平波动尺度的增加，临界滑面宽度随之增大；总体而言坡体上部局部安全系数随竖向波动尺度的增加而减小。