降雨条件下路基土体湿度场迁移预测模型

Prediction Model of Moisture Field Migration in Subgrade Soil Under Rainfall Conditions

摘要

受降雨的影响，非饱和路基土体的湿度场会发生迁移，改变土体的基质吸力，进而影响土体强度，对边坡稳定性产生威胁。为深入研究降雨条件下非饱和土体的湿度场迁移规律，对土体湿度场迁移的影响因素进行理论分析与公式推导，提出了降雨条件下非饱和路基土体的动态平衡饱和度预测公式、湿润锋迁移预测公式和地下水位波动预测公式。通过自主研发的1维土柱试验模型装置和3维有限元数值模拟，初步揭示了湿度场变化规律，并验证了上述预测公式的可行性。基于预测模型，分析了土体压实度和降雨强度对湿度场迁移的影响规律。研究结果表明：1）土体压实度对湿度场迁移的影响存在一个临界值，当土体压实度低于0.80时，湿度场迁移速度随着压实度的增加而加快，而当土体压实度超过0.80时，湿度场迁移速度反而随着压实度的增加而减慢；2）湿度场迁移速度随着降雨强度的增加而加快，然而，当降雨强度达到15mm/h之后，随着降雨强度的进一步提高，降雨强度对湿度场迁移速度的提升作用逐渐减弱。本文提出了一种湿度场迁移预测公式，可较好地预测非饱和土体水分累积、湿润锋迁移和地下水位上升过程，为进一步完善边坡预警方法、提升边坡稳定性评估与灾害预防能力提供了理论基础。

Abstract

Objective The increasing incidence of rainfall-induced embankment landslides corresponds with the recent rise in extreme rainfall events. During precipitation, moisture field migration occurs within unsaturated embankment soils, which modifies matric suction distribution and soil strength characteristics, threatening embankment stability. Accurately characterizing the spatiotemporal distribution pattern of the water content field and predicting its moisture field migration mechanism is a key prerequisite for advancing slope stability evaluation methods and optimizing geological disaster prevention strategies. Although current research extensively documents macroscopic principles governing moisture field movement, the prediction equation of the moisture field in soil remains scarce and inadequately validated. Prevailing analytical models predominantly assume complete soil saturation behind wetting fronts, but actual field conditions contradict this premise, which influences model accuracy. Therefore, the moisture field migration prediction model is proposed to reveal the migration pattern of the soil moisture field and provide a theoretical basis for improving slope instability warning methods. Methods Firstly, using a self-developed one-dimensional soil column test system, typical embankment filling sand was selected as the research material, and an unsaturated soil column structure was constructed through layered compaction. The test adopted the rainfall intensity control method, with 35 mm/h as the benchmark rainfall intensity, and simulated rainfall conditions using an axial uniform spraying system. Moisture sensors were installed along the depth of the soil column to monitor the dynamic response characteristics of the volumetric moisture content of different soil layers in real time to accurately and continuously capture the moisture migration process. Secondly, using COMSOL finite element software, the modified Cambridge model, unsaturated empirical formulas, and the VG model were introduced to establish one-dimensional soil column models, and the numerical model was validated using the test results. The variation law of the moisture field was preliminarily revealed through both the model test and numerical simulation. Thirdly, based on the indoor test and simulation results, the main influencing factors of the soil dynamic equilibrium saturation degree were analyzed, and a prediction formula for the soil dynamic equilibrium saturation degree was proposed and validated through the numerical simulation model. Then, based on the assumption of a homogeneous soil model, the quantitative relationship between the rainfall volume and rainwater retention in soil was analyzed, the soil dynamic equilibrium saturation prediction formula was introduced, and the prediction formula for the dynamic equilibrium saturation of unsaturated soil was established and verified through indoor tests and numerical simulations. In addition, the quantitative relationship between rainfall volume and the amount of rainwater required to achieve full soil saturation, as well as the prediction formula for groundwater level rise, was established, and the feasibility of these prediction formulas was verified through experiments and numerical models. Finally, based on the above formulas, the effects of wet front migration and groundwater level rise under different compaction degrees (0.75, 0.80, 0.85, and 0.90) and rainfall intensities (5, 10, 15, 20, 25, 30, and 35 mm/s) were discussed, and the migration pattern of the unsaturated soil moisture field under different rainfall intensities was revealed, providing an important theoretical basis for embankment landslide warning. Results and Discussions 1) The results of the model tests and numerical simulations showed that the changes in soil volumetric water content under rainfall exhibit five distinct stages, namely the initial stage, the first growth stage, the dynamic equilibrium stage, the second growth stage, and the complete saturation stage. Among them, the first growth stage occurs due to the influence of rainfall on the soil, during which the wetting front migrates from the soil surface to the bottom and causes a gradual increase in soil volumetric moisture content from top to bottom. After the wetting front passes through the soil, the volumetric water content enters the dynamic equilibrium stage. As the groundwater level rises from below, the volumetric water content enters the second growth stage, and eventually, all soil reaches the fully saturated stage. 2) The feasibility of the dynamic equilibrium saturation prediction formula, the wetting front migration prediction formula, and the groundwater level fluctuation prediction formula for unsaturated soil under rainfall conditions was double verified through indoor tests and numerical simulations. The results showed that the maximum error between the dynamic equilibrium saturation prediction formula and the simulation was 4%, the maximum error between the wetting front migration prediction formula and both the simulation and test was 8.5%, and the maximum error between the groundwater level fluctuation prediction formula and both the simulation and test was 8.5%. The trend of the humidity field change curve calculated by these formulas was consistent with the experimental and numerical simulation results. Therefore, these formulas better predicted the migration behavior of the soil moisture field and were suitable for studying unsaturated soil moisture field migration. 3) The unit pore volume, soil permeability, and soil compaction degree of unsaturated soil showed a clear inverse relationship. The reduction in unit pore volume accelerates the migration of the unsaturated soil moisture field, whereas the decrease in soil permeability slows the migration speed of the moisture field. When the soil compaction degree is less than or equal to 0.80, the migration of the humidity field is mainly controlled by the unit pore volume, and the migration speed increases with increasing compaction degree. When the soil compaction degree exceeds 0.80, the migration of the moisture field is mainly controlled by soil permeability, and the migration speed decreases with a higher compaction degree. 4) The migration speed of the unsaturated soil moisture field showed a positive correlation with rainfall intensity, indicating that the migration speed increases as rainfall intensity increases. However, as rainfall intensity continues to increase, its rising effect on the migration speed of the humidity field gradually weakens. 5) At any given moment, the influence of soil compaction and rainfall intensity on the migration of the humidity field remains consistent and does not change with longer rainfall duration. Conclusions The study demonstrates that the dynamic equilibrium saturation prediction formula, wetting front migration prediction formula, and groundwater level fluctuation prediction formula effectively predict soil moisture field migration. In addition, the effects of soil compaction degree and rainfall intensity on moisture field migration are clarified. The findings are expected to provide a theoretical basis for improving slope instability warning methods and enhancing the capacity for slope stability assessment and disaster prevention.

Graphical abstract

关键词

湿润锋 / 地下水位 / 动态平衡饱和度预测 / 湿度场迁移预测 / 压实度 / 降雨强度

Key words

moist front / groundwater level / dynamic equilibrium saturation prediction / humidity field migration prediction / compaction degree / rainfall intensity

本刊网刊

近年来，因降雨引发的山体与路堤边坡滑坡、坍塌问题已演变成严重且广泛存在的自然灾害。随着极端气候现象的不断加剧，此类灾害的发生愈发频繁，强度也持续攀升。在极端降雨条件下，非饱和路基边坡出现的滑坡与坍塌情况和土体中孔隙水压力场的动态变化紧密相关^[1‒3]。一方面，基质吸力（即负孔隙水压力）在增强土体抗剪强度方面发挥着关键作用，主要通过提高土体的有效内摩擦角及黏聚力提升土质边坡的稳定性^[4]；另一方面，地下水位的上升会导致正孔隙水压力增大，软化土体并加速边坡的蠕变变形，最终导致边坡稳定性大幅降低^[5]。

Liao等^[6]研究表明，当湿润锋穿透土层，湿润锋后的土体将完全饱和，土体基质吸力完全消失，进而产生正孔隙水压力，对边坡稳定性产生不利影响。实际上，在不产生积水的情况下，湿润锋穿过土层后，湿润锋后土体的饱和度将上升，却不会完全达到饱和状态，而是趋于某一动态平衡值。在这种状况下，土体基质吸力虽会下降但不会完全消失，也不会产生正孔隙水压力^[7]。当湿润锋到达地下水位后，随着降雨持续，地下水位将从边坡底部向顶部逐步抬升，使得边坡土体自下而上逐步达到完全饱和状态。此时，土体内的基质吸力完全消散，正孔隙水压力随之产生，进而使边坡土体软化，并加速其蠕变变形，导致边坡稳定性降低^[8‒9]。综上所述，边坡发生滑坡和坍塌是基质吸力减小与地下水位上升耦合作用的结果。因此，深入探究降雨条件下边坡土体中正/负孔隙水压力场的变化规律，分析边坡土体湿度场的演变规律，对于揭示降雨诱发型边坡滑坡的内在机理有至关重要的作用。

目前，众多学者从不同角度对土体湿度场的演化与分布规律开展了较为丰富的研究。在解析方法方面，汪丁建等^[10]以Green‒Ampt模型为基础，推导出一种滑坡降雨入渗函数，旨在理想化地描述土体内湿润锋的迁移进程，Liao等^[6]也建立了湿润锋迁移计算模型，模型均基于降雨强度大于土体饱和渗透系数的假设，且随着入渗过程的持续，湿润锋锋后的土体将完全饱和。实际上，降雨强度大于土体饱和渗透系数的情况并不常见，湿润锋穿过土层后，仅使土体湿润，却难以达到完全饱和的状态，导致解析结果与试验及数值模拟结果相比出现较大误差。大量的边坡滑坡实例也表明，滑坡发生区域的浅层土体通常并未完全饱和^[11‒12]。因此，现有模型难以精确地描述土体湿度场的发展及分布规律。在数值模拟方面，Gallage^[13]和Wang^[14]等分别采用流固耦合MPM材料点法、SEEP/W与SLOPE/W软件对边坡土体湿度场的发展与分布进行了模拟，尽管这些方法在一定程度上能模拟降雨结束后边坡土体的湿度场分布规律，但在模拟降雨过程中土体内湿润锋动态迁移过程方面仍存在局限。在室内试验方面，史振宁^[15]和覃小华^[16]等设计了一种可测量土体体积含水率的降雨入渗试验模型，试验结果表明，随着降雨的不断持续，土体饱和度自上而下逐层提高，直至地下水位上升，导致土体自下而上完全饱和，这一试验结果充分揭示了土体湿度场的迁移规律，但目前尚未提出用于预测土体湿度场迁移的预测公式。

针对现有研究中存在的不足，建立了可模拟降雨条件下土体湿润锋迁移过程的有限元模型，并利用室内试验对模型进行验证。同时，建立了降雨条件下非饱和土体的动态平衡饱和度预测公式、湿润锋迁移预测公式和地下水位波动预测公式，以揭示土体湿度场迁移规律，为进一步完善边坡失稳预警方法提供理论基础，助力提升边坡稳定性评估与灾害预防能力。

1 降雨入渗模型试验及数值模拟

1.1 降雨入渗试验

为揭示降雨入渗条件下土体湿度场的变化规律，本文研制了1维土柱试验模型装置，并开展了降雨入渗模型试验。砂土力学性能稳定、透水性强、分布广泛且施工便捷，其颗粒级配连续性与结构均质特征便于精准控制试验变量，研究结果可外推至粗粒土体系^[17]；同时，砂土的低黏聚力、非膨胀性及弱压缩特性可规避复杂土‒水耦合效应干扰，利于聚焦排水特性主导的湿度场迁移机制。因此，本文选取砂土为研究对象^[18]。图1为1维土柱试验模型装置示意图，该装置由亚克力透明塑料材质制成，高100 cm，直径40 cm。如图1所示，装置底部和顶部分别铺设10 cm厚砾石层，并在底部砾石层上逐层填筑砂土层。

为精准获取砂土的关键物理参数，开展了击实试验、固结试验和常水头试验。经测定，本文所用砂土最优含水率为15.5%，最大干密度为1.60 g/cm³，压实度为0.85，孔隙比为0.694，饱和渗透系数为0.147 4 mm/s。同时，通过筛分法试验，测得本文所用砂土的颗粒级配曲线，如图2所示。为实时观测1维土柱试验模型装置内土体体积含水率的动态变化规律，在装置内布设了由浙江吉欧科技有限公司生产的水分传感器（测量精度为±3%），在布设时，将最下层的水分传感器设置在距桶底10 cm的位置，随后向上以每隔10 cm为间距设置测点（图1）。

近年来，由于极端降雨增加，传统降雨等级划分标准难以指导地质灾害预警，中国气象局公共气象服务中心起草发布了《短时气象服务降雨量等级标准》，如表1所示。但是，降雨对地区的影响受各地地形、土体特征等因素影响，各地对短时强降雨的定义存在较大差异：江苏省界定的极端强降雨阈值为45 mm/h，且不同县市存在±10 mm/h地区差异；青海省将降雨强度大于30 mm/h界定为短时大暴雨；浙江省气象局将可能造成滑坡灾害的降雨强度阈值定为50 mm/h。因此，综合考虑《短时气象服务降雨量等级标准》及各地对于极端降雨的定义，本文将短期极端强降雨下限值设置为35 mm/h，通过对土柱顶部缓慢喷淋，进行土体湿度场迁移试验。试验以土柱内含水率保持不变或土体完全饱和作为试验结束的标志。

1.2 数值模拟

在降雨条件下，非饱和土体的体积含水率在湿润锋附近存在急剧变化现象^[19]。本文通过对比非饱和土体的实时体积含水率与初始体积含水率，判断并描述湿润锋的时空演变规律。对于某一测点，若其某一时间节点i处体积含水率θ(i)大于初始体积含水率θ₀，且在i前的所有时间节点的体积含水率θ(i')均小于或等于初始体积含水率θ₀，则认为在i时湿润锋经过该点的坐标即为i时湿润锋的空间位置。通过对不同测点湿润锋通过时间进行统计并绘制湿润锋时间‒空间曲线图，即可描述湿润锋时空演变规律。

利用商用软件COMSOL Multiphysics建立基于1维土柱降雨入渗试验的土柱3维有限元模型，如图3所示。模型顶部（AB平面）设置通量边界条件，以模拟降雨入渗；根据试验条件，模型侧面（AC和BD）和底部均设置为不排水边界。在室内试验中，采用初始体积含水率为15%的土体开展试验，因此，将数值模型中土体残余体积含水率设置为15%。在室内试验中，装置内初始地下水位位于砾石层表面以下1 cm处，砂土与砾石层交界处受毛细作用所产生的微弱吸力应为-0.1 kPa。在数值模型中，CD平面初始孔隙水压力设置为-0.1 kPa，以表征毛细作用下的初始湿度梯度。根据试验条件，模型侧面应设置为辊支撑，对侧面节点x和y方向进行约束，仅允许其垂直方向上位移，对模型底面设置为固定端约束，不允许土体产生水平和垂直方向位移。模型网格类型采用自由四面体，网格大小根据软件预定义设置为常规网格（最大单元大小为0.080 0，最小单元大小为0.014 4，曲率因子为0.6）。

采用VG（Van Genuchten）模型^[20]对非饱和土体的土水特征进行模拟：

S_{\mathrm{e}}=S_{\mathrm{r}\mathrm{R}}-(\mathrm{1}-S_{\mathrm{r}\mathrm{R}})(\mathrm{1}+{(\alpha s)}^r{)}^{-m}（1）

式中：S_{\mathrm{e}}为饱和度；S_{\mathrm{r}\mathrm{R}}为残余饱和度；s为基质吸力；α、r和m为拟合待定系数，与模拟土体材料有关。

同时，引入非饱和土体渗透系数折减公式^[20]和土体饱和渗透系数与孔隙比关系的经验公式^[21]，以模拟不同土层深度及饱和度下非饱和土渗透特性：

{\kappa }_{\mathrm{r}}=S_{\mathrm{e}}^l(\mathrm{1}-(\mathrm{1}-S_{\mathrm{e}}^{-m}{)}^m{)}^{\mathrm{2}}（2）

k_{\mathrm{s}}=\mathrm{532.4}\frac{e^{\mathrm{3}}}{\mathrm{1}+e}{d}_{\mathrm{9}}^{\mathrm{2}}（3）

式（2）～（3）中：κ_r为渗透系数折减系数；l为拟合待定系数，与模拟土体材料有关，现有研究中l通常取值为0.5^[22]；k_s为土体饱和渗透系数；e为土体孔隙比；d₉为土体有效粒径，其值由柯森公式推导所得^[21,23]。

本文采用修正剑桥模型并引入初始孔隙比计算公式及塑性体积模量λ和对数体积模量κ随基质吸力变化的经验公式^[24]，模拟非饱和土物理力学特性：

c_{\mathrm{u}}=\frac{M}{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{(\mathrm{1}+e_{\mathrm{1}})-(\lambda -\kappa )\mathrm{l}\mathrm{n}\mathrm{2}}{\lambda }-\frac{\mathrm{1}+e_{\mathrm{0}}}{\lambda }\right)（4）

e_{\mathrm{0}}=e_{\mathrm{1}}-(\lambda -\kappa )\mathrm{l}\mathrm{n}{p}_{\mathrm{c}}^{\text{'}}+\kappa \mathrm{l}\mathrm{n}{p}_{\mathrm{0}}^{\text{'}}（5）

\lambda (s)=\lambda (\mathrm{0})-\frac{{\lambda }_{s}s}{p_{\mathrm{1}}+s}（6）

\kappa (s)=\kappa (\mathrm{0})+{\kappa }_{s}s（7）

式（4）～（7）中，c_u为不排水抗剪强度，e₀为各深度初始孔隙比，e₁为INCL线在e-ln p'空间上的截距，λ为塑性体积模量，κ为对数体积模量，M为临界状态线斜率，p'₀为初始平均应力，p'_c为各点对应屈服面大小，λ(s)为不同基质吸力下λ值，λ(0)为基质吸力为0时λ值，λ_s 为表达非饱和土压缩指数随基质吸力大小变化的参数，p₁为大气压，κ(s)为不同基质吸力下κ值，κ(0)为基质吸力为0时κ值，κ_s 为表达非饱和土回弹指数随基质吸力大小变化的参数。

本文通过击实试验、标准固结试验、直接剪切试验^[25]和滤纸法试验^[26]，分别对土体力学模型和水力模型基本参数进行测量，试验所得的土体土水特征曲线（SWCC）如图4所示，土体水力模型和力学模型基本参数如表2所示。

1.3 数值模拟的验证分析

图5为不同埋深下数值模拟与试验体积含水率验证结果。图5中详图为埋深0.2 m处土体在各个阶段的体积含水率曲线，直观呈现了该深度处土体体积含水率在不同阶段的变化情况。由图5可知，在降雨的全过程中，试验结果与数值模拟结果之间的最大误差为14%，而当降雨结束后，二者的最大误差缩小至4%。在试验和数值模拟中，降雨作用下土体体积含水率的变化均呈现出5个明显阶段，分别为初始阶段、首次增长阶段、动态平衡阶段、二次增长阶段及完全饱和阶段。

在首次增长阶段，由于土体受降雨影响，湿润锋自土体表面向底部迁移，造成土体体积含水率呈现出自上而下逐渐增加的态势。待湿润锋经过土体后，土体体积含水率将进入动态平衡阶段。随着地下水位自下而上升，土体体积含水率将进入二次增长阶段，最终，所有土体均进入完全饱和阶段。

综上所述，数值模拟结果与试验结果呈现出高度的一致性，验证了数值模型的可靠性，表明该土柱3维有限元模型在模拟降雨作用下非饱和土体湿度场的演化过程与分布规律方面具备良好的适用性与准确性。

2 非饱和土体湿度场迁移预测

在降雨入渗的作用下，当湿润锋穿过土体后，湿润锋后方土体的体积含水率与饱和度会逐步提升，与此同时，土体的基质吸力与抗剪强度则相应降低。土体饱和度的升高，一方面源于降雨的直接影响，另一方面，在重力作用下，土体内的水分会向深层土体渗透，使得原土层的土体饱和度有所下降。在降雨和重力双重因素影响下的土体饱和度变化对土体的稳定性及相关工程性质有着重要影响。

图6为湿度场迁移示意图。如图6（a）所示，降雨过程中，随着湿润锋的迁移，土体逐渐被湿润。由于土体内水分的下渗作用，湿润锋过后的土体无法达到完全饱和状态，而是进入一种动态平衡状态，此时土体渗透速率k与降雨强度Q相等。同时，在湿润锋前方的传导区域，土体的渗透速率会随着土体饱和度的上升而逐渐增大，最终稳定在某一平衡数值。其中，土体渗透速率k与土体渗透能力和土体饱和度有关，其与饱和渗透速率数量关系系数定义为β。在部分存在水分累积的情况下，土体的渗透系数与土体饱和渗透速率保持一致。

不同情况下土体渗透速率与降雨强度关系如下：

（8）

式中，k_sat为土体饱和渗透速率。

滞留于土体内的雨水体积V_wet由下式计算：

V_{\mathrm{w}\mathrm{e}\mathrm{t}}=V_{\mathrm{f}\mathrm{i}}-V_{\mathrm{i}\mathrm{n}}（9）

式中，V_fi和V_in分别为湿润锋经过后土体体积含水率和土体初始体积含水率。

如图6（b）所示，随着降雨的持续，雨水将逐渐渗透至地下水位处，地下水位由土体底部向顶部逐渐上升，土体由湿润状态达到完全饱和状态的所需雨水体积V_increase可由下式计算：

V_{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}}=V_{\mathrm{s}}-V_{\mathrm{f}\mathrm{i}}（10）

式中，V_s为土体完全饱和后体积含水率。

2.1 非饱和土体动态平衡饱和度预测

为实现对非饱和土体动态平衡饱和度的预测，本文假设：在整个降雨过程中，土体内不会出现局部水分累积的情况。经上述分析可知，土体饱和度仅与降雨强度Q及土体渗透能力相关。本文采用饱和渗透系数k_s来表示土体的渗透能力，并且将降雨强度与土体渗透系数的比值作为动态平衡饱和度预测公式的自变量，将其定义为渗排通量比（D），进而建立起动态平衡饱和度预测公式：

S_{\mathrm{r}\mathrm{b}\mathrm{0}}=a{D}^{\mathrm{0.109}\mathrm{1}}（11）

式中，S_{\mathrm{r}\mathrm{b}\mathrm{0}}为表层土体动态平衡饱和度，a为拟合待定系数，且有：

D=\frac{Q}{k_s}（12）

同时，引入土体饱和渗透系数与孔隙比关系的经验公式（式（3）），得：

S_{\mathrm{r}\mathrm{b}\mathrm{0}}=a{\left(\frac{(\mathrm{1}+e)Q}{\mathrm{532.4}{e}^{\mathrm{3}}{d}_{\mathrm{9}}^{\mathrm{2}}}\right)}^{\mathrm{0.109}\mathrm{1}}（13）

在降雨入渗过程中，随着水分渗透路径的不断延长，土体受降雨作用的影响会逐渐衰减。因此，为使土体动态平衡饱和度预测更为精准，引入深度影响因素，建立修正后的土体动态平衡饱和度（S_rb）预测公式：

S_{\mathrm{r}\mathrm{b}}=b{x}^{\mathrm{2}}+cx+S_{\mathrm{r}\mathrm{b}\mathrm{0}}（14）

式中：x为土体埋深；b、c分别为曲线特征系数，与土体材料性质相关，通过试验测量拟合获得。

将式（13）代入式（14）得：

S_{\mathrm{r}\mathrm{b}}=b{x}^{\mathrm{2}}+cx+a{\left(\frac{(\mathrm{1}+e)Q}{\mathrm{532.4}{e}^{\mathrm{3}}{d}_{\mathrm{9}}^{\mathrm{2}}}\right)}^{\mathrm{0.109}\mathrm{1}}（15）

利用式（15）对不同降雨强度下（5、10、15、20、25、30 mm/h）不同深度的土体动态平衡饱和度进行预测，并通过有限元模型对解析结果进行验证，结果如图7所示。

图7中，解析预测值与数值模拟所得结果最大误差仅为4%。要特别指出的是，除了最接近地下水位的土层（埋深0.7 m）以外，其余各土层的解析预测值与数值模拟结果的误差均低于1%。埋深为0.7 m处误差相对较大，其原因在于：该层距地下水位较近，除了受降雨入渗影响外，还受到地下水位的影响，尤其当降雨强度较小时，动态平衡饱和度受地下水位的影响更为显著，从而导致预测误差较大。然而，随着降雨强度逐渐增加，降雨入渗的影响逐渐占据主导地位，所以预测误差也随之逐渐减小。对于大多数边坡工况，地下水位一般位于坡底以下数米处，可忽视其对土体动态平衡饱和度预测的影响。

综上所述，本文建立的非饱和土体动态平衡饱和度预测公式的解析结果与数值模拟结果可较好吻合，意味着该公式在土体湿度场相关研究以及边坡预警方面具备较高的精准度与可靠性，能够为研究土体在不同降雨条件下的动态平衡饱和度提供可靠依据。

2.2 湿润锋迁移预测

受降雨的影响，湿润锋逐渐从土体表面向深层土体迁移，导致土体基质吸力降低，造成边坡稳定性下降，引发边坡滑坡灾害。为预测湿润锋的迁移过程，本文假设在湿润锋迁移过程中，湿润锋后不存在局部土体水分累积现象。为简化计算模型，假设土柱为均质材料，即不同深度处土体孔隙率相等，且忽略埋深对土体动态平衡饱和度的影响。基于上述假设，降雨入渗的雨水体积与滞留于土体中的雨水体积相等，即：

\frac{V_{\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{a}\mathrm{l}\mathrm{l}}}{V_{\mathrm{w}\mathrm{e}\mathrm{t}}}=\mathrm{1}（16）

式中，V_rainfall和V_wet分别为雨水体积和滞留于土体内雨水体积，并有：

V_{\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{a}\mathrm{l}\mathrm{l}}=tQ{S}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}（17）

V_{\mathrm{w}\mathrm{e}\mathrm{t}}=V_{\mathrm{f}\mathrm{i}}-V_{\mathrm{i}\mathrm{n}}=n{z}_{\mathrm{w}}(S_{\mathrm{r}\mathrm{f}\mathrm{i}}-S_{\mathrm{r}\mathrm{i}\mathrm{n}})S_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}（18）

式（17）、（18）中，S_soil为土体底面积，t为降雨时长，n为土体孔隙率，z_w为湿润区深度，S_{\mathrm{r}\mathrm{f}\mathrm{i}}为湿润区土体饱和度，S_{\mathrm{r}\mathrm{i}\mathrm{n}}为土体初始饱和度。

将式（17）和（18）代入式（16）得：

\frac{Qt}{n(S_{\mathrm{r}\mathrm{f}\mathrm{i}}-S_{\mathrm{r}\mathrm{i}\mathrm{n}})z_{\mathrm{w}}}=\mathrm{1}（19）

即：

z_{\mathrm{w}}=\frac{Qt}{n(S_{\mathrm{r}\mathrm{f}\mathrm{i}}-S_{\mathrm{r}\mathrm{i}\mathrm{n}})}（20）

引入非饱和土体动态平衡饱和度预测公式（式（13））得：

z_{\mathrm{w}}=\frac{Qt}{n\left(a{\left(\frac{(\mathrm{1}+e)Q}{\mathrm{532.4}{e}^{\mathrm{3}}{d}_{\mathrm{9}}^{\mathrm{2}}}\right)}^{\mathrm{0.109}\mathrm{1}}-S_{\mathrm{r}\mathrm{i}\mathrm{n}}\right)}（21）

土体湿润锋迁移的试验结果、模拟结果和解析预测结果如图8所示。通过式（21）计算所得的土体湿润锋迁移曲线与试验结果之间的最大误差为6.0%，而与模拟结果之间的最大误差为8.5%。试验、模拟和解析结果呈现出较好的一致性，共同验证了式（21）的可行性，证明本研究所提出的湿润锋迁移预测公式可预测特定条件下湿润锋迁移过程，为非饱和土体关于湿润锋迁移研究提供了可靠的理论工具。

2.3 地下水位上升预测

为简化计算公式，本文忽略实际工程中可能存在的土体内局部水分累积现象，并假设土体为均质模型，即不同深度处土体孔隙率和动态平衡饱和度相同。图9为地下水位上升预测公式验证结果。如图9所示，随着降雨的持续进行，雨水将逐渐渗透至地下水位所在位置，地下水位将由土体底部向顶部逐渐上升，原处于湿润区的土体将逐渐达到完全饱和状态。湿润区土体从湿润状态转变为完全饱和状态所需的雨水体积可由下式计算：

V_{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}}=Qt{S}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}=V_{\mathrm{s}}-V_{\mathrm{f}\mathrm{i}}=n{z}_{\mathrm{g}\mathrm{w}}\left(\mathrm{1}-S_{\mathrm{r}\mathrm{f}\mathrm{i}}\right)S_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}（22）

即：

z_{g\mathrm{w}}=\frac{Qt}{n\left(\mathrm{1}-S_{\mathrm{r}\mathrm{f}\mathrm{i}}\right)}（23）

式中，z_gw为地下水位上升高度。

通过室内模型试验与数值模拟，对地下水位上升预测公式进行验证，如图9所示，解析结果与数值模拟结果间的最大误差为8.5%，出现在地下水位开始上升后的1 500 s；而解析结果与室内试验结果间的最大误差为7.5%，出现在地下水位开始上升后的2 500 s。这些误差可能源于解析计算过程中未能充分考虑不同深度土体孔隙率的差异性。

综上所述，本文所建立的解析公式与试验及数值模拟结果的吻合程度较高，表明该公式能有效预测特定条件下土体内地下水位的迁移过程，可为非饱和土体湿润锋迁移研究提供可靠的理论工具，有助于深化对土体水分迁移机制的理解。

3 土体湿度场迁移影响因素

3.1 土体压实度对非饱和土体湿度场迁移影响

土体孔隙率随压实度的增加而逐渐减小，导致单位体积内土体孔隙体积减小，进而影响降雨对土体湿度场的影响能力和土体湿度场迁移。

图10为不同压实度下土体湿润锋迁移规律。图10中，差异比值为相邻湿润锋迁移深度的差异比值。如图10所示，当土体压实度较低（0.80及以下）时，随着压实度的提升，土体湿润锋的迁移速度增快。在降雨时长为100 s时，压实度为0.80的土体湿度场迁移深度比压实度为0.75的土体深0.000 66 m，二者湿润锋迁移深度差异比值为7.76%。在降雨时长为7 500 s时，两者湿润锋迁移深度差异仅为0.029 44 m，湿润锋迁移深度比值为7.68%。这表明，随着降雨时长的不断增加，不同压实度土体的湿润锋迁移深度差异会逐渐增大，但其差异的比值却相对稳定，仅呈现出轻微下降趋势，在降雨全过程中差异比值最大变化量仅0.08%。这种现象可以归因于土体孔隙率随压实度的增加而逐渐减小，导致土体达到湿润状态所需的雨水量相应减少，因此在相同的降雨条件下，湿润锋的迁移速度会加快。随着降雨时间的持续增加，不同压实度的土体在湿润锋迁移深度上的差异愈发明显，但降雨时间对于相邻土体压实度湿润锋迁移变化比值影响则相对较小。

当土体压实度处于较高水平（0.8及以上）时，会出现与低压实度情况相反的现象：土体湿润锋的迁移速度随压实度的提高而逐渐减缓。究其原因，土体压实度对土体特性有多方面的影响，不仅影响土体的孔隙率和孔隙体积，还会显著改变雨水在土体内的渗透率。在压实度较高的情况下，土体的孔隙结构变得更为致密，使得雨水在土体内渗透的难度增大，渗透率降低。由于湿润锋的迁移依赖于雨水的渗透，所以渗透率的降低直接导致了湿润锋迁移速度的减缓。此时，土体湿润锋的迁移受渗透率变化的影响更为突出，导致随着压实度进一步增加，湿润锋的迁移速度不升反降。此外，根据现有试验研究与对Green‒Ampt模型^[27]进行的数学推导可知，土体单位孔隙体积（有效孔隙度Δθ）与湿润锋距离L呈负相关。Δθ减小时，土体的储水能力降低，在相同入渗水量条件下，湿润锋需要覆盖更大的体积以容纳水分，导致其迁移距离增加，迁移速度加快^[28]；与此同时，当土体压实度增加时，土体颗粒会堵塞孔隙，降低土体渗透性，从而减缓流体运移^[29]。

最后，随着降雨时长的不断增加，不同压实度土体湿润锋迁移变化比值保持相对稳定，在降雨全过程中仅产生0.26%的波动。这证明了压实度对湿润锋迁移所产生的影响作用，并不会因降雨时长的改变而改变。

对于不同材质的路基土体，均存在使土体湿度场迁移速度产生变化的临界压实度值。但是，不同材质土体的压实度临界值与其物理特性、孔隙结构密切相关。黏性土的临界压实度与其持水特性和孔隙结构密切相关，黏粒含量高导致其临界压实度较低^[30]。砂性土的临界值则受颗粒级配与动荷载响应双重调控，其临界压实度高于黏性土，需通过振动压实增强骨架密实性^[31]。在工程实践中，土体压实度临界值直接影响路基动荷载下的湿度稳定性和路基区域适应性^[32]。若压实度低于临界值，在动荷载作用下土体内孔隙水压力会产生骤增骤减，引发累积塑性变形。若土体压实度达到临界值，湿润锋迁移速度与荷载频率相匹配，则差异沉降风险降低。且土体压实度临界值还能有效平衡水分迁移与持水能力，既可抑制冻区的冻胀与融沉灾害，又可避免干旱区路基因过度干燥引发的收缩裂缝。

图11为不同压实度下地下水位上升规律。如图11所示，土体地下水位的上升速度随压实度的增加而逐渐加快。随着时间的持续推进，不同压实度下土体的地下水位上升高度的差异愈发明显。然而，相邻压实度土体间地下水位上升的差异比值在任何给定时间点都保持相对稳定，波动幅度极小，模拟全过程中最大波动仅为0.01%。这表明，尽管压实度对地下水位上升速度有显著影响，但在不同时间节点上，这种影响的相对变化情况保持一致。

3.2 降雨强度对非饱和土体湿度场迁移的影响

图12为不同降雨强度下土体湿润锋迁移规律。如图12所示，湿润锋的迁移速度与降雨强度呈现正相关，即随着降雨强度的增大，湿润锋迁移速度逐渐加快。然而，从相邻降雨强度条件下湿润锋迁移差异比值曲线分析可知，当降雨强度进一步提升时，降雨强度对湿润锋迁移速度的促进作用逐渐减弱。此外，随着降雨时长的增加，不同降雨强度条件下湿润锋迁移的差异值会不断增大，但相邻降雨强度条件下土体湿润锋迁移的差异比值变化并不显著，在降雨全过程中差异比值的最大波动值仅为0.14%，表明尽管降雨强度对湿润锋迁移速度有着明显影响，但降雨时长的影响相对有限。

图13为不同降雨强度下地下水位上升规律。如图13所示，在降雨时长为100 s时，不同降雨强度条件下土体湿润锋的迁移深度呈现出显著差异。具体来看，当降雨强度分别为5、10、15、20、25、30、35 mm/h时，土体湿润锋的迁移深度分别为0.006 4、0.015 2、0.026 2、0.039 0、0.053 7、0.070 4和0.089 4 m。相邻降雨强度条件下湿润锋迁移的差异比值分别为137.50%、72.37%、48.85%、37.69%、31.10%和26.99%。当降雨时长延长至7 500 s时，各降雨强度下土体湿润锋迁移深度均有显著增加，其中，降雨强度为5、10、15、20、25、30、35 mm/h时，土体湿润锋迁移深度分别达到0.056、0.135、0.230、0.343、0.472、0.620和0.787 m。相应地，相邻降雨强度条件下湿润锋迁移的差异比值分别为139.75%、71.23%、48.94%、37.70%、31.19%和26.97%。这些数据表明，降雨强度与地下水位上升速度紧密相关，随着降雨强度的增加，地下水位上升的速度也逐步加快。同时，这一现象与降雨强度对湿润锋迁移的影响相似，即随着降雨强度的不断增大，其对地下水位上升速度的促进作用逐渐减弱。此外，随着降雨时长的增加，不同降雨强度条件下地下水位上升的差异值不断增大，但相邻降雨强度条件下土体湿润锋迁移的差异比值变化并不显著。这进一步证实了降雨时长对地下水位上升速度的影响相对有限。

4 结 论

本文提出了适用于降雨条件下非饱和路基土体动态平衡饱和度预测公式、湿润锋迁移预测公式和地下水位波动预测公式。基于所提预测公式，揭示了土体压实度和降雨强度对土体湿度场迁移的影响规律。主要结论如下：

1）经室内试验和数值模拟的双重验证，本文所提出的降雨条件下非饱和土体动态平衡饱和度预测公式、湿润锋迁移预测公式和地下水位波动预测公式可较好预测土体湿度场迁移规律，可用于非饱和土体湿度场迁移研究。

2）非饱和土体的单位孔隙体积、土体渗透能力与土体压实度之间存在明确的反比关系。单位孔隙体积的减小会加速非饱和土体湿度场的迁移，而土体渗透能力的降低则会减缓湿度场的迁移速度。当土体压实度小于或等于0.80时，湿度场迁移主要受单位孔隙体积的控制，迁移速度随着压实度的增加而加快；而当土体压实度大于0.80时，湿度场迁移则主要受土体渗透能力的控制，迁移速度随着压实度的增加而减慢。

3）非饱和土体湿度场的迁移速度与降雨强度呈正相关，即随着降雨强度的增大而加快。然而，当降雨强度持续增大时，其对湿度场迁移速度的促进作用会逐渐呈现出减弱的趋势。

4）在任一给定时刻，土体压实度和降雨强度对湿度场迁移的影响作用保持一致，不会随着降雨时长的增加而发生变化。

5）土体压实度临界值通过调控非饱和渗流路径与基质吸力分布，直接影响湿度场迁移速率及边坡稳定性。其差异性受矿物组成、孔隙分形特征及动载下骨架‒孔隙水耦合作用制约。未来需对不同土体的压实度临界值进行研究，揭示其对土体湿度场迁移的影响规律。

参考文献

[1] Baum R L, Godt J W.Early warning of rainfall-induced shallow landslides and debris flows in the USA[J].Landslides,2010,7(3):259‒272. doi:10.1007/s10346-009-0177-0

[2] Chen Xiangyu, Zhang Lulu, Zhang Limin,et al.Modelling rainfall-induced landslides from initiation of instability to post-failure[J].Computers and Geotechnics,2021,129:103877. doi:10.1016/j.compgeo.2020.103877

[3] Yeh H F, Tsai Y J.Effect of variations in long-duration rainfall intensity on unsaturated slope stability[J].Water,2018,10(4):479. doi:10.3390/w10040479

[4] Guo Longxiao, Chen Guangqi, Gong Shilin,et al.Analysis of rainfall-induced landslide using the extended DDA by incorporating matric suction[J].Computers and Geotechnics,2021,135:104145. doi:10.1016/j.compgeo.2021.104145

[5] Li Xiang, Chen Yanni, Handwerger A L,et al.Dynamics of creeping landslides controlled by inelastic hydro-mechanical couplings[J].Engineering Geology,2023,317:107078. doi:10.1016/j.enggeo.2023.107078

[6] Liao Wenwang, Ji Jian.Time-dependent reliability analysis of rainfall-induced shallow landslides considering spatial variability of soil permeability[J].Computers and Geotechnics,2021,129:103903. doi:10.1016/j.compgeo.2020.103903

[7] Gerscovich D M S, Vargas E A, de Campos T M P.On the evaluation of unsaturated flow in a natural slope in Rio de Janeiro,Brazil[J].Engineering Geology,2006,88(1/2):23‒40. doi:10.1016/j.enggeo.2006.07.008

[8] Wang Chen, Hawlader B, Perret D,et al.Modeling of initial stresses and seepage for large deformation finite-element simulation of sensitive clay landslides[J].Journal of Geotechnical and Geoenvironmental Engineering,2021,147(11):04021111. doi:10.1061/(asce)gt.1943-5606.0002626

[9] Duan Guanghui, Song Fei, Wang Huaning,et al.Stability analysis of unsaturated loess slopes subjected to extreme rainfall incorporating creep effects[J].Computers and Geotechnics,2024,169:106231. doi:10.1016/j.compgeo.2024.106231

[10] Wang Dingjian, Tang Huiming, Li Changdong,et al.Stability analysis of colluvial landslide due to heavy rainfall[J].Rock and Soil Mechanics,2016,37(2):439‒445.

[11] 汪丁建,唐辉明,李长冬,等.强降雨作用下堆积层滑坡稳定性分析[J].岩土力学,2016,37(2):439‒445.

[12] Ran Qihua, Hong Yanyan, Li Wei,et al.A modelling study of rainfall-induced shallow landslide mechanisms under different rainfall characteristics[J].Journal of Hydrology,2018,563:790‒801. doi:10.1016/j.jhydrol.2018.06.040

[13] Vahedifard F, Leshchinsky D, Mortezaei K,et al.Effective stress-based limit-equilibrium analysis for homogeneous unsaturated slopes[J].International Journal of Geomechanics,2016,16(6):D4016003. doi:10.1061/(asce)gm.1943-5622.0000554

[14] Gallage C, Abeykoon T, Uchimura T.Instrumented model slopes to investigate the effects of slope inclination on rainfall-induced landslides[J].Soils and Foundations,2021,61(1):160‒174. doi:10.1016/j.sandf.2020.11.006

[15] Wang Di, Wang Bin, Yuan Weihai,et al.Investigation of rainfall intensity on the slope failure process using GPU-accelerated coupled MPM[J].Computers and Geotechnics,2023,163:105718. doi:10.1016/j.compgeo.2023.105718

[16] Shi Zhenning, Qi Shuangxing, Fu Hongyuan,et al.A study of water content distribution and shallow stability of earth slopes subject to rainfall infiltration[J].Rock and Soil Mechanics,2020,41(3):980‒988.

[17] 史振宁,戚双星,付宏渊,等.降雨入渗条件下土质边坡含水率分布与浅层稳定性研究[J].岩土力学,2020,41(3):980‒988.

[18] Qin Xiaohua, Liu Dongsheng, Song Qianghui,et al.Experimental study on one-dimensional vertical infiltration in soil column under rainfall and the derivation of permeability coefficient[J].Chinese Journal of Rock Mechanics and Engineering,2017,36(2):475‒484.

[19] 覃小华,刘东升,宋强辉,等.降雨条件下一维土柱垂直入渗模型试验研究及其渗透系数求解[J].岩石力学与工程学报,2017,36(2):475‒484.

[20] Cui Xianze, Zhou Yuchang, Fan Yong, et al.Experimental study on migration and deposition characteristics of fine particles in sandwich sand under seepage condition[J].Journal of Hydraulic Engineering,2017,36(2):475‒484.

[21] 崔先泽,周钰錩,范勇,等.渗流条件下夹层砂土中细颗粒迁移沉积特性试验研究[J].水利学报,2024,55(8):977‒989.

[22] Bai Qinqin, Li Xinsheng, Yan Mingkang,et al.Experimental study on the effect of water content on the cohesive force of unsaturated sand soil[J].Science Technology and Engineering,2019,19(7):195‒199.

[23] 白琴琴,李新生,严明康,等.含水率对非饱和砂土似黏聚力影响试验研究[J].科学技术与工程,2019,19(7):195‒199.

[24] Zhao Wenhe, Yang Xiujuan, Wang Baozhong,et al.Laws of water migration and settlement at interface in loess filled areas under rainfalls[J].Chinese Journal of Geotechnical Engineering,2022,44(9):1710‒1720.

[25] 赵文赫,杨秀娟,王宝仲,等.降雨作用下黄土填方区交界面水分迁移及沉降规律研究[J].岩土工程学报,2022,44(9):1710‒1720.

[26] van Genuchten M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Science Society of America Journal,1980,44(5):892‒898. doi:10.2136/sssaj1980.03615995004400050002x

[27] Zhu Liyi, Yang Peng, Wensheng Lü,et al.Experimental on influence of pore ratio on tailings sand permeability under different chemical conditions[J].The Chinese Journal of Nonferrous Metals,2020,30(9):2190‒2200.

[28] 诸利一,杨鹏,吕文生,等.不同化学条件下孔隙比对尾矿砂渗透性影响试验[J].中国有色金属学报,2020,30(9):2190‒2200.

[29] Yu Daijin, Huang Qiangbing, Kang Xiaosen,et al.The unsaturated seepage process and mechanism of internal interfaces in loess-filled slopes during intermittent rainfall[J].Journal of Hydrology,2023,619:129317. doi:10.1016/j.jhydrol.2023.129317

[30] Hu Qizhi, Zou Qiang, Li Ming,et al.Study on permeability coefficient of cohesive soil based on three-dimensional ideal model of soil particles[J].Science Technology and Engineering,2020,20(11):4512‒4519.

[31] 胡其志,邹强,李鸣,等.基于土颗粒三维理想模型的黏性土渗透系数研究[J].科学技术与工程,2020,20(11):4512‒4519.

[32] Hu Xiaorong, Lu Xiang, Cai Xiaofeng.The triple-shear unified constitutive model for unsaturated clays based on the generalized effective stress concept[J].Engineering Mechanics,2024,41(5):107‒119.

[33] 胡小荣,路祥,蔡晓锋.非饱和黏性土的广义有效应力三剪弹塑性本构模型[J].工程力学,2024,41(5):107‒119.

[34] 中华人民共和国交通运输部. JTG 3430—2020 公路土工试验规程 [S].北京:人民交通出版社,2020.

[35] Bai Fuqing, Liu Sihong, Yuan Jiao.Measurement of SWCC of Nanyang expansive soil using the filter paper method[J].Chinese Journal of Geotechnical Engineering,2011,33(6):928‒933.

[36] 白福青,刘斯宏,袁骄.滤纸法测定南阳中膨胀土土水特征曲线试验研究[J].岩土工程学报,2011,33(6):928‒933.

[37] Jiang Shuihua, Yuan Zhirong, Liu Xian,et al.Rainfall infiltration model considering spatial variability of multiple layers in transition layer and its application in slope stability analysis[J].Chinese Journal of Geotechnical Engineering,2025,47(2):255‒264.

[38] 蒋水华,元志镕,刘贤,等.考虑过渡层多层空间变异性的边坡降雨入渗模型及稳定性分析应用[J].岩土工程学报,2025,47(2):255‒264.

[39] Cui Xinzhuang, Zhang Xiaoning, Wang Yilin,et al.State-of-the-art review on measurement methods,evolution laws,and control technologies of subgrade moisture field[J].China Journal of Highway and Transport,2024,37(6):1‒33.

[40] 崔新壮,张小宁,王艺霖,等.路基湿度测量方法、演化规律及调控技术研究进展[J].中国公路学报,2024,37(6):1‒33.

[41] Zhu Qin, Su Lijun, Liu Zhenyu,et al.Study of seepage in wide-grading soils with particles migration[J].Rock and Soil Mechanics,2021,42(1):125‒134.

[42] 朱秦,苏立君,刘振宇,等.颗粒迁移作用下宽级配土渗透性研究[J].岩土力学,2021,42(1):125‒134.

[43] Dong Cheng, Yang Xianzhang, Liu Wenjie,et al.Study on compaction degree criterion of red clay highway subgrade in southern Hunan[J].Journal of Highway and Transportation Research and Development,2017,34(2):42‒49.

[44] 董城,杨献章,刘文劼,等.湘南红黏土公路路基压实度标准研究[J].公路交通科技,2017,34(2):42‒49.

[45] Wu Yang, Huang Jinsheng, Cui Jie,et al.Influences of particle shape and degree of compaction on shear response of clinker ash[J].Chinese Journal of Geotechnical Engineering,2021,43(12):2220‒2229.

[46] 吴杨,黄锦盛,崔杰,等.颗粒形状和压实度对炉渣颗粒土力学特性的影响[J].岩土工程学报,2021,43(12):2220‒2229.

[47] Thuy D V, Tiep P D, Hieu N V,et al.Influence of loading frequency and relative compaction on liquefaction behavior of reconstituted sand in cyclic triaxial tests[J].Rock and Soil Mechanics,2024,45(6):1813‒1823. doi:10.26599/rsm.2024.9435673

基金资助

浙江省自然科学基金资助项目(LMS25E080006)

国家自然科学基金项目(52378467)

浙江省大学生科技创新活动计划（新苗人才计划）(2024R403C097)