滑坡是自然界分布广泛的地质灾害之一，其中降雨是滑坡的主要触发因素之一^[1-2]。目前，国内外学者对于边坡降雨入渗的研究主要包含理论^[3-4]、试验研究^[5-6]以及数值模拟^[7-8]等方面。相较于理论分析而言，模型试验和数值模拟的计算成本偏高^[9]。现有理论研究中，由于经典Green-Ampt模型参数简明，计算便捷，被广泛应用于非饱和土的降雨入渗分析^[10-11]，但其假定湿润锋面以上均为饱和区域，且未考虑土体渗透系数变化对雨水入渗规律的影响，具有一定的局限性。对此，众多学者对GA模型进行改进，并进一步将其用于分析边坡的稳定性。

文献[12]通过对土体的饱和渗透系数进行随机正态分布，以此对GA模型进行修正，并进一步将其应用于边坡稳定性分析，得到了土体渗透系数再分布过程中边坡破坏概率随降雨时间的变化规律。文献[13]以黄土为例，在分层假定的基础上考虑了入渗过程中黄土渗透系数随深度的变化，以研究黄土地区的雨水入渗规律。文献[14]在BOUWER的研究基础上考虑了土体初始含水率非均匀分布的影响，将土体渗透系数修正为0.5倍的饱和渗透系数，但其在研究中并未考虑土体渗透系数随降雨过程的变化。文献[15]分别采用Richards和GA模型对土体的降雨入渗结果进行对比，发现当土体渗透系数为0.7倍的饱和渗透系数时拟合效果最好。文献[16]利用Van Genuchten模型拟合土水特征曲线，得到不同含水率下的非饱和渗透系数，并采用AHFO技术定位湿润锋的迁移位置。上述文献均在研究土体降雨入渗的过程中对土体渗透系数进行了修正，但在确定土体渗透系数时，大多采用假设或拟合计算的方法，仅适合某些特定土体，未考虑土体颗粒形状、颗粒级配分布等因素对渗透系数的影响。

文献[17]分析了黏土孔隙及其弯曲度与渗透系数的关系，并基于颗粒间微观力的影响，模拟了黏土的固结试验，采用离散元法计算了渗透系数。文献[18]分别对原状土和重塑土进行了固结渗透试验，认为渗透系数与土颗粒大小有关，并根据相关试验结果提出了孔隙比与渗透系数的关系模型。土体孔隙具有较为明显的分形特征，采用分形理论研究非饱和土的力学参数可以加强对非饱和土的认识^[19]。文献[20]用分形维数描述了土体孔隙比，并在此基础上对黏土的渗透系数进行了修正，以减小渗透系数在压缩变形中的误差。

此外，土体的渗透系数还需结合降雨强度、地层特征以及土体非饱和区等因素的影响考虑。许多学者在研究GA模型时假定湿润区完全饱和，这会高估湿润区土体的渗透系数和水力梯度，低估初期湿润锋的形成时间^[21]。文献[22]发现土体含水率在降雨过程中有较为明显的层状分布，从上至下可大致分为饱和带、饱和过渡带、天然土层。对此可将湿润区土体划分为饱和带土层和过渡带土层两部分。文献[23]在分析降雨过程中边坡的入渗规律及稳定性影响时，考虑了非饱和过渡带土层的存在，并假定其占整个湿润区的一半。文献[24]研究了非饱和湿润区的特性，分析模型参数并基于1/4椭圆曲线，提出了一种采用分段法简化GA模型的方法。因此，在研究非饱和土入渗问题中，需确定土体水分剖面模型，考虑土体渗透系数随入渗区过渡带占比的动态变化。

鉴于此，本文从土体孔隙分布特征的角度出发，基于分形理论提出了土体饱和渗透系数的分析模型，考虑入渗区过渡带随降雨的实时变化关系，推导了不同降雨历时下土体的水分剖面模型，从而对入渗区等效渗透系数进行修正，得到了改进的GA入渗模型；进一步考虑土体孔隙分布特征对抗剪强度的影响，建立了修正后的稳定性分析模型，使本文理论具有一定的普遍性。通过COMSOL数值仿真和已有试验结果验证了模型的准确性。并以深圳市龙岗区良白路的路基边坡为依托工程，探究土体孔隙分布特征及初始含水率对边坡降雨入渗规律及稳定性的影响，以帮助理解非饱和土的降雨入渗机制，为实际工程提供参考。

1 改进Green-Ampt模型 1.1 土体饱和渗透系数修正

土体的孔隙分布直接影响了土体的孔隙度和孔隙连通性。较高的孔隙度和良好的孔隙连通性有利于土体水分入渗，对降雨入渗过程中的水分流动速度和路径选择起着重要作用。但对于土体孔隙这种无序的物体，很难通过欧式几何的方法进行描述。由于土体孔隙表面具有分形特征，通过应用分形几何，可以更好地理解和描述土体孔隙的复杂性。

非饱和土的基质吸力是表示土体力学性质，用于了解非饱和土入渗规律的一个重要参数。非饱和土基质吸力与土体孔隙半径r之间的关系可以通过Young-Laplace方程确定(ψ=(2Tcos a)/r，T为表面张力，a为接触角，取为0)。基于土体孔隙表面具有分形特征，可得由非饱和土孔隙孔径增量导致的土体相对含水率变化关系：

$ \left\{\begin{array}{l} \mathit{\Lambda }=\theta-\theta_{\mathrm{r}}=\frac{4 {\rm{ \mathsf{ π} }} C}{V_{\mathrm{t}}(3-D)} r^{3-D} \\ \mathrm{~d} \mathit{\Lambda }=4 {\rm{ \mathsf{ π} }} K(A, r) r^2 \mathrm{~d} r / V_{\mathrm{t}} \end{array}\right. $

(1)

式中：r为假定覆盖土体孔隙表面的圆球半径，C为比例系数，K(A，r)为圆球数目，V_t为孔隙体积，θ为非饱和土的体积含水量，θ_r为残余体积含水率，D为孔隙分形维数。

当土体为饱和状态时，相对体积含水量可表示为

$ \mathit{\Lambda }_{\mathrm{s}}=\theta_{\mathrm{s}}-\theta_{\mathrm{r}}=\frac{4 {\rm{ \mathsf{ π} }} C}{V_{\mathrm{t}}(3-D)} L^{3-D} $

(2)

式中：L为最大孔隙半径，θ_s为土体的饱和体积含水率。

由Young-Laplace方程可得土体粒径分布描述的SWCC模型：

$ S_{\mathrm{e}}=\frac{\mathit{\Lambda }}{\mathit{\Lambda }_{\mathrm{s}}}=\left(\frac{r}{L}\right)^{3-D}=\left(\frac{\psi}{\psi_{\mathrm{c}}}\right)^{D-3} $

(3)

式中：S_e为土体的饱和度，ψ为土体的基质吸力，ψ_c为土体的进气值。

文献[25]根据土体孔隙分布规律得到了土体渗透系数和基质吸力的表达式：

$ \left\{\begin{array}{l} k(\theta)=k_{\mathrm{s}} S_{\mathrm{e}}^{3+2 /(3-D)} \\ \psi(\theta)=\psi_{\mathrm{c}} S_{\mathrm{e}}^{-1 /(3-D)} \end{array}\right. $

(4)

式中k_s为土体的饱和渗透系数。

文献[26]用孔隙-土体分形模型(PSF)计算分形维数时，发现以lg(d_i/d_max)为横坐标，lg[(M(r < d_i)/M₀]为纵坐标，基于最小二乘法进行线性回归拟合，拟合的斜率λ=3-D，即可求出各土样的分形维数D，表达式为

$ (3-D) \lg \frac{d_{\mathrm{i}}}{d_{\max }}=\lg \frac{M\left(r <d_{\mathrm{i}}\right)}{M_0} $

(5)

式中：M(r < d_i)/M₀为粒径小于d_i的土颗粒质量百分数，d_i为土体颗粒直径，d_max为土体中最大颗粒直径。该模型适用于已知土体级配分布求解分形维数。

文献[27]将土-水特征曲线作为反映土体孔隙通道的指标，发现土体的饱和渗透系数是大体积孔隙通道渗透系数的总和，结合基于经典流体力学，提出了饱和渗透系数模型：

$ k_{\mathrm{s}}=k_{\mathrm{c}} \int_{\theta_{\min }}^{\theta_{\max }} \frac{\mathrm{d} \theta}{\psi^2(\theta)} $

(6)

式中k_c=γT_s²cos² α/(2p_iμ)，对于同一土样该值为常数。

当基质吸力大于土体进气值(ψ＞ψ_c)，忽略残余体积含水率影响时，可进一步将式(3)修正为

$ \mathrm{d} \theta=(D-3) \psi_{\mathrm{c}}^{3-D} \psi^{D-4} \frac{e}{1+e} \mathrm{~d} \psi $

(7)

式中: e为土体的孔隙比; θ_s为土体的饱和体积含水率, θ_s=e/(1+e)。

假定残余体积含水率θ_r对应最小的土体孔隙(θ_min=θ_r)，土体进气值ψ_c对应最大的土体孔隙，则联立式(6)、(7)可得

$ k_{\mathrm{s}}=\frac{e k_{\mathrm{c}}(D-3)}{(1+e)(D-5) \psi_{\mathrm{c}}^2}\left[1-\left(\frac{\psi_{\mathrm{a}}}{\psi_{\mathrm{c}}}\right)^{D-5}\right] $

(8)

式中ψ_a为残余体积含水率对应的基质吸力。

在三维空间中，维数D介于2和3之间，且ψ_a远大于ψ_c，则

$ 1-\left(\psi_{\mathrm{a}} / \psi_{\mathrm{c}}\right)^{D-5} \approx 1 $

(9)

可进一步将式(8)简化为

$ k_{\mathrm{s}}=\frac{\lambda e k_{\mathrm{c}}}{(1+e)(2+\lambda) \psi_{\mathrm{c}}^2} $

(10)

式中λ=3-D。

1.2 强降雨作用下改进GA模型的建立

在强降雨条件下，土体的降雨入渗过程可分为降雨强度控制和土体入渗速率控制阶段。这是由于在降雨初期，非饱和土的入渗强度大于降雨强度，雨水全部被土体所吸收。随着降雨时间的推移，到达某一临界时刻t_p，降雨强度开始大于土体雨水入渗速率，雨水不能完全入渗，坡面开始出现积水^[28]。故在强降雨工况下降雨入渗速率可表示为

$ i=\left\{\begin{array}{l} q \cos \alpha, t \leqslant t_{\mathrm{p}} \\ k(\theta) \frac{z_{\mathrm{f}} \cos \alpha+\psi_{\mathrm{f}}+h}{z_{\mathrm{f}}}, t>t_{\mathrm{p}} \end{array}\right. $

(11)

式中：k(θ)为土体入渗区等效渗透系数，z_f为湿润峰深度，α为边坡坡角，ψ_f为湿润峰处土体的基质吸力水头值，可由式(4)推导得出。

降雨入渗初期，土体中形成的湿润区域不断扩大，经典GA模型在分析中假定湿润锋面以上土体饱和，湿润锋面以下土体的体积含水率不变。然而实际情况下，在饱和带和天然土层之间应该存在一个过渡带，即在降雨后土体的含水率分布应该呈现为饱和带、过渡带以及天然土层3个部分。文献[29]研究Richards入渗方程，基于试验和数值拟合，得到了过渡带与湿润区的厚度比例随入渗过程的变化规律。文献[24]以1/4椭圆曲线作为土体入渗过渡带的分布模型分析土体的入渗规律。

本文同样采用椭圆形曲线描述过渡带体积含水率沿边坡土体深度的分布规律，如图 1所示。在降雨入渗初期，土体水分剖面呈椭圆形分布，见图 1(b)，椭圆过渡带的平均含水率θ可表示为

$ \left(\bar{\theta}-\theta_{\mathrm{i}}\right) \times z_{\mathrm{t}}=\frac{{\rm{ \mathsf{ π} }}}{4}\left(\theta_{\mathrm{s}}-\theta_{\mathrm{i}}\right) \times z_{\mathrm{t}} $

(12)

式中：z_t为过渡带深度；θ_i为边坡初始体积含水率。

根据水量守恒定律，可知t_p时间内土体的入渗总量I：

$ I=\frac{{\rm{ \mathsf{ π} }}}{4}\left(\theta_{\mathrm{s}}-\theta_{\mathrm{i}}\right) \times z=q t \cos \alpha $

(13)

临界时刻t_p及对应的临界湿润锋深度z_p的函数关系为

$ \left\{\begin{array}{l} t_{\mathrm{p}}=\frac{I_{\mathrm{p}}}{q \cos \alpha} \\ z_{\mathrm{p}}=\frac{k(\bar{\theta}) \psi_{\mathrm{f}}}{[q-k(\bar{\theta})] \cos \alpha} \end{array}\right. $

(14)

在降雨入渗后期，土体过渡带占总湿润区的比例随降雨历时不断降低^[30]。且当降雨强度不变时，入渗过程中过渡带土层的厚度将趋于稳定。因此本文假设在降雨初期已形成的过渡带土层厚度不再发生改变，如图 1(c)所示，土体水分剖面出现分层现象，过渡带以上土体含水率饱和，其中饱和带深度为z_w，过渡带深度为z_t。降雨后期边坡土体含水率的分布模型可修正为

$ \theta(z)=\left\{\begin{array}{l} \theta_{\mathrm{s}}, 0 \leqslant z \leqslant z_{\mathrm{w}} \\ \theta_{\mathrm{i}}+\left[\theta_{\mathrm{s}}-\theta_{\mathrm{i}}\right] \sqrt{\left[1-\left(z-z_{\mathrm{w}}\right)^2 / z_{\mathrm{t}}^2\right]}, z_{\mathrm{w}} \leqslant z \leqslant z_{\mathrm{f}} \\ \theta_{\mathrm{i}}, z \geqslant z_{\mathrm{f}} \end{array}\right. $

(15)

式中z_f为湿润锋深度，z_f=z_w+z_t。

由于土体的渗透系数与含水率存在明显的正相关关系，整个湿润区内的等效渗透系数可以表示为

$ k(\theta)=k_{\mathrm{s}} \cdot(1-R)+k(\bar{\theta}) \cdot R $

(16)

式中R=z_t/z_f。

对于浅层边坡，在临界时刻t_p之后，结合式(10)、(11)、(16)可将土体的雨水入渗速率i_s修正为

$ \begin{aligned} i_{\mathrm{s}}= & \frac{\lambda k_{\mathrm{c}} \theta_{\mathrm{s}}}{(2+\lambda) \psi_{\mathrm{c}}^2}\left[(1-R)+R S_{\mathrm{e}}(\bar{\theta})^{3+2 / \lambda}\right] \cdot \\ & \frac{z_{\mathrm{f}} \cos \alpha+\psi_{\mathrm{f}}+h}{z_{\mathrm{f}}} \end{aligned} $

(17)

当t>t_p时，坡面径流产生，此时累计入渗量为

$ I=t_{\mathrm{p}} q \cos \alpha+\int_{t_{\mathrm{p}}}^t i_{\mathrm{s}} \mathrm{~d} t $

(18)

基于以上假定与推导，各阶段累积入渗总量可表示为

$ I=\left\{\begin{array}{l} t q \cos \alpha=\frac{{\rm{ \mathsf{ π} }}}{4}\left(\theta_{\mathrm{s}}-\theta_{\mathrm{i}}\right) z_{\mathrm{f}}, t \leqslant t_{\mathrm{p}} \\ \int_0^t i \mathrm{~d} t=\left(\theta_{\mathrm{s}}-\theta_{\mathrm{i}}\right) z_{\mathrm{w}}+\frac{{\rm{ \mathsf{ π} }}}{4}\left(\theta_{\mathrm{s}}-\theta_{\mathrm{i}}\right)\left(z_{\mathrm{f}}-z_{\mathrm{w}}\right), t>t_{\mathrm{p}} \end{array}\right. $

(19)

2 边坡稳定性分析

降雨入渗条件下，对于基岩型边坡，除了基岩处的危险点外，也容易发生平行于边坡表面破坏的浅层滑坡，故常采用无限长边坡模型分析浅层滑坡的稳定性^[31-33]，图 2为边坡稳定性分析模型。

当考虑土体孔隙分布特征λ的影响时，可将土体有效应力公式进行修正，其表达式为

$ \sigma^{\prime}=\left(\sigma_{\mathrm{n}}-u_{\mathrm{a}}\right)+\gamma_{\mathrm{w}} \psi_{\mathrm{c}}^\lambda \psi^{1-\lambda} $

(20)

式中：σ_n-u_a为单位土条底部净法向力，u_a通常取0，γ_w为水的重度，其他参数符号与上文同义。

降雨过程中，对于饱和带的土体，原有的负孔隙水压力消失，可忽略基质吸力影响；但对于边坡过渡带土层以及天然土层，则必须考虑基质吸力影响，由基质吸力引起的剪切强度τ_s可表示为

$ \tau_{\mathrm{s}}=\gamma_{\mathrm{w}} \psi_{\mathrm{c}}^\lambda \psi^{1-\lambda} \tan \varphi^{\prime} $

(21)

非饱和土的剪切强度τ_f可进一步修正为

$ \tau_{\mathrm{f}}=c^{\prime}+\left(\sigma_{\mathrm{n}}-u_{\mathrm{a}}\right) \tan \varphi^{\prime}+\gamma_{\mathrm{w}} \psi_{\mathrm{c}}^\lambda \psi^{1-\lambda} \tan \varphi^{\prime} $

(22)

式中：c′为土体的黏聚力，φ′为土体的内摩擦角。

由降雨入渗过程中坡体含水率变化可得坡体各点位容重γ_m(z)的表达式为

$ \left\{\begin{array}{l} \text { 饱和土层: } \gamma_{\mathrm{m}}(z)=\gamma_{\mathrm{s}}, 0 \leqslant z \leqslant z_{\mathrm{w}} \\ \text { 过渡土层: } \gamma_{\mathrm{m}}(z)=\gamma_{\mathrm{d}}+\gamma_{\mathrm{w}} \theta(z), z_{\mathrm{w}} \leqslant z \leqslant z_{\mathrm{f}} \\ \text { 天然土层: } \gamma_{\mathrm{m}}(z)=\gamma_{\mathrm{d}}+\gamma_{\mathrm{w}} \theta_{\mathrm{i}}, z \geqslant z_{\mathrm{f}} \end{array}\right. $

(23)

式中：γ_d为土体的干重度；z为坡内任意位置距离坡顶的深度；其他参数符号与式(15)同义。

在式(23)的基础上可得边坡滑动面以上土体重量：

$ W=\int_0^{z_{\mathrm{i}}} \gamma_{\mathrm{m}}(z) \mathrm{d} z $

(24)

式中z_i为坡体滑动面所处点位深度。

进一步计算滑动面底部正应力σ_n和土体破坏切应力τ_m：

$ \left\{\begin{array}{l} \sigma_{\mathrm{n}}=\int_0^{z_{\mathrm{i}}} \gamma_{\mathrm{m}} \mathrm{~d} z \cos ^2 \alpha \\ \tau_{\mathrm{m}}=\int_0^{z_{\mathrm{i}}} \gamma_{\mathrm{m}} \mathrm{~d} z \sin \alpha \cos \alpha \end{array}\right. $

(25)

由修正的有效应力公式和Mohr-Coulomb强度准则可得坡内各点位的安全系数表达式：

$ \begin{aligned} F_{\mathrm{s}}=\frac{\tau_{\mathrm{f}}}{\tau_{\mathrm{m}}}=& \frac{c^{\prime}+\int_0^{z_{\mathrm{i}}} \gamma_{\mathrm{m}} \mathrm{~d}z \cos ^2 \alpha \tan \varphi^{\prime}+\gamma_{\mathrm{w}} \psi_{\mathrm{c}}^\lambda \psi^{1-\lambda} \tan \varphi^{\prime}}{\int_0^{z_{\mathrm{i}}} \gamma_{\mathrm{m}} \mathrm{~d} z \sin \alpha \cos \alpha}= \\ & \frac{\tan \varphi^{\prime}}{\tan \alpha}+\frac{c^{\prime}+\gamma_{\mathrm{w}} \psi_{\mathrm{c}}^\lambda \psi^{1-\lambda} \tan \varphi^{\prime}}{\int_0^{z_{\mathrm{i}}} \gamma_{\mathrm{m}} \mathrm{~d} z \sin \alpha \cos \alpha} \end{aligned} $

(26)

综上，基于改进GA模型及分形理论对边坡稳定性分析模型进行修正以研究强降雨作用下边坡的失稳破坏规律，具体分析步骤如下：1)确定土体的物理力学参数，明确边坡初始含水率的分布模型；2)通过式(14)确定临界降雨时刻t_p以及对应的临界湿润锋深度z_p，得出边坡降雨入渗及径流边界条件，根据式(15)和式(19)，得出边坡湿润锋随降雨历时的变化规律，进一步考虑过渡带的影响，得出降雨后坡体的含水率分布；3)根据式(26)计算边坡各点位的稳定性系数，取其中的最小值作为边坡的安全系数，边坡的最危险滑动面则为最小值对应的滑动面。

3 模型验证 3.1 水分剖面分布模型的验证

为验证本文改进GA模型的有效性，以非饱和无限长边坡为例，将本文入渗模型的计算结果与Richards方程数值解和文献理论模型^[34-35]计算结果进行对比分析。边坡上覆厚度约为3 m风化残积土，下部为不透水基岩层，且平行于坡面，边坡角度为30°，风化残积土计算参数为深圳市龙岗区良白路的路基边坡上覆土体实际勘测数据，取值如表 1所示。

采用COMSOL软件中的Richards方程接口对边坡的降雨过程进行数值模拟分析，并利用式(19)计算不同时刻的降雨入渗深度。图 3为本文降雨入渗理论结果及其他文献理论结果^[34-35]与数值模拟解析解的对比。从图 3可以看出，降雨后坡体的含水率分布呈现为饱和带、过渡带以及天然土层3个部分，湿润锋以上并非全为饱和区域，相较于文献理论结果，本文理论结果与Richards方程结果更为一致。此外，由于临界时刻t_p之前为降雨强度控制阶段，此时雨水全部入渗，而文献[35]在分析中未考虑降雨过渡带的影响，从而低估初期湿润锋的形成时间，导致降雨初期本文理论模型得出的湿润峰推移速率要高于潘永亮模型^[35]结果；随降雨继续进行，整个湿润区的等效渗透系数开始变化，而文献[35]仍采用饱和渗透系数计算，所以在降雨后期湿润峰推移速率要大于本文模型。且降雨后期坡体的入渗速率开始缓慢减小，并逐渐趋于稳定，这与Richards方程结果较为一致，可以验证本文模型的有效性。

3.2 剪切强度分析模型的验证

验证分形理论修正下边坡稳定性分析模型的准确性，重点需明确基质吸力引起的剪切强度的准确性。根据文献[36]土-水特征曲线的试验结果可得不同固结压力(0、25、100、200 kPa)下对应的土体进气值分别为27、30、62、86 kPa，将其试验结果通过式(3)进行拟合，得到非饱和土孔隙分布的分形维数约为2.77，见图 4(a)。文献[36]中非饱和土的有效内摩擦角φ′为23°，通过式(21)可得由基质吸力引起的剪切强度变化曲线，图 4(b)所示为不同初始条件下式(21)的理论结果与试验结果的对比分析。从图中可以看出，随着基质吸力的增加，由基质吸力引起的剪切强度τ_s逐渐趋近平稳，本文的理论计算结果与试验结果具有较好的一致性，可以证明本文理论的有效性。

4 参数分析与讨论

深圳市龙岗区良白路的路基边坡上覆土体多为粉质黏土和风化残积土，将其简化为图 5所示的非饱和无限长边坡，对本文改进的入渗模型及稳定性模型进行分析与讨论。坡体上覆厚度约为2 m粉质黏土，下部为不透水基岩层，且平行于坡面，边坡角度为33.7°。研究主要采用控制变量法进行，以探究土体孔隙分布特征及初始含水率对非饱和土边坡降雨入渗及稳定性的影响。根据当地气象资料，选取降雨强度为3 mm/h的强降雨工况进行分析，假设降雨为均匀雨型形式。表 2所示为粉质黏土颗粒尺度与质量分布的实测数据，表 3所示为粉质黏土计算参数的现场勘测数据。

将表 2中粉质黏土的颗粒分布数据代入粒径分布(PSF)分形模型(式5)中，并进行线性回归拟合，处理结果如图 6所示，线性拟合的平均斜率约为0.285，得到粉质黏土的分形维数D=2.715。

4.1 渗透系数随降雨时间及土体孔隙特征的变化

图 7显示了表 3中粉质黏土在3种不同级配分布下的等效渗透系数与饱和渗透系数之比(k(θ)/k_s)随R值的变化，并在此基础上考虑了两种不同的初始含水率(0.2θ_s和0.4θ_s)。从图 7中可以看出，对于相同类型的土体，k(θ)/k_s与R呈单调递减关系；随着R的减小，k(θ)/k_s逐渐接近1。对于相同类型的土体，初始含水量对k(θ)/k_s-R曲线的影响较大，这是由于过渡带中平均饱和度受初始含水率影响较大，进一步导致了过渡带中平均渗透系数的显著变化，从而影响了整个湿润区的等效渗透系数。

在降雨入渗过程中，湿润区内过渡带的占比随着入渗时间的增加而不断减小，因此，R的减小过程反映了降雨时间的增加过程。图中的分析表明，在入渗过程中，土体等效渗透系数小于土体饱和的渗透系数，尤其是当R较小时。因此，在入渗过程中，特别是在入渗的初始阶段，假设整个湿润区的渗透系数与土体饱和时的渗透系数相同，将产生较大的误差。经典Green Ampt模型假设湿润区的渗透系数为土体饱和渗透系数，因此高估了土体中实际入渗速率，从而对边坡稳定性的评估造成不利影响。图 8所示为土体的饱和渗透系数随孔隙分布特征的变化，从图中可以看出两者直接存在较为明显的正相关关系。

4.2 土体孔隙分布特征λ对入渗及稳定性的影响

在所讨论的模型参数中，参数λ与土体颗粒孔隙分布有关，影响了土体中水分的流动路径及土体的剪切强度，并进一步影响坡体的入渗规律及稳定性。对于同一土体由于地质原因以及其所处深度不同, 其级配分布规律也有所差别，造成其渗透性能有较为明显的差异，但现有研究大多没有考虑这一影响因素^{[14-15, 22]}。本文从土体颗粒级配分布的角度出发，假设了几种不同参数λ，以探究土体孔隙分布特征对边坡渗流规律及稳定性的影响，结果如图 9所示。

图 9(a)所示为边坡湿润区等效渗透系数k(θ)及入渗速率i的变化示意图，从图中可以看出土体的渗透能力随着参数λ的增加而增加，这是由于对于同一类型土体，土体的饱和渗透系数与参数λ成正相关关系，而饱和渗透系数的增长会直接影响土体的实际入渗能力，导致相同降雨时间内参数λ越大，湿润峰推移深度越深，如图 9(b)。

图 9(c)所示为边坡水分剖面示意图，土体的颗粒粒径分布会直接影响入渗过程中边坡湿润区内过渡带的占比，相同降雨时间，参数λ越大过渡带与入渗区的厚度比例越大。从图 9(d)可以看出，土体孔隙分布特征对边坡的稳定性影响较为明显，随参数λ的增加，边坡的稳定性进一步降低。这是由于一方面，湿润峰的推移速率受到土体孔隙分布特征的影响，致使相同降雨时间对应的坡底正应力增加；另一方面，土体孔隙分布特征也会影响了坡体水分剖面随降雨时间的变化规律，改变了土体的基质吸力，并进一步影响了土体的剪切强度τ_s，从而对边坡稳定性造成较大的影响。

4.3 初始含水率对入渗及稳定性的影响

边坡初始含水率的大小是影响边坡入渗规律及稳定性的重要因素之一，但很多研究往往假定边坡初始含水率为一定值，未深入分析边坡初始含水率分布的影响。对此本文假设了4种不同初始含水率大小，分析结果如图 10所示。

对于相同湿润区深度，土体初始含水率越大，湿润区整体的等效渗透系数越大。但土体初始含水率越大，也意味着相同降雨条件下边坡越容易达到饱和，边坡湿润峰的推移速率也越快，致使相同降雨时间内初始含水率越大，土体渗透速率i反而下降越快，如图 10(a)、10(b)。对于过渡带土层，由于剪切强度τ_s受土体含水率的影响显著，在过渡带区域变化较大，容易致使坡体在过渡带区域发生失稳破坏，所以在进行边坡稳定性分析时不仅需考虑基岩处的边坡稳定性，边坡过渡带内的稳定性也是需要考虑的重点，如图 10(c)、10(d)。相较与文献[22, 35]本文模型可以更直观的看出坡体内各点位稳定性随降雨时间的变化规律。

5 结论

土体入渗能力的确定是研究土体降雨入渗规律的难点之一，本文从土体孔隙分布特征的角度出发，基于分形理论对边坡的降雨入渗模型及稳定性分析模型进行修正，并通过数值仿真和已有试验对本文理论模型进行验证。以良白路路基边坡为依托工程，进一步探究了土体孔隙分布特征及初始含水率对坡体入渗特征及稳定性的影响规律。主要结论如下：

1) 本文提出的改进模型考虑了土体入渗过渡带以及整个湿润区等效渗透系数随降雨过程的实时变化，故能准确模拟出降雨过程中的土体入渗状态及剪切强度变化，计算精度高于传统模型，更贴近客观实际。且相比于传统模型，本文模型能够准确识别边坡最危险滑动面位置，更有利于工程上的应用。

2) 工程实际中，边坡的降雨入渗过程受到众多环境因素的影响，其中坡体的初始含水率极大影响了降雨后坡体的水分剖面分布，显著改变了土体的基质吸力大小，并进一步影响了土体的剪切强度τ_s，从而对边坡稳定性造成了较大的影响。

3) 土体孔隙分布特征决定了水分在土体中的流动性，影响土体实际入渗能力；土体颗粒粒径越大，孔径分布越粗糙，土体允许入渗能力越大，坡表产生径流时间越晚。此外，过渡带土层的剪切强度受土体孔隙分布特征及初始含水率影响尤为明显，极大影响了边坡在过渡带稳定性，可为揭示降雨诱发的浅层滑坡提供理论依据。