基于粒径分布曲线的非饱和砂土土水特征曲线概率预测模型*

张 洁 阳 帅 谭泽颖 李 旭

(①同济大学地下建筑与工程系，上海 200092，中国)(②北京交通大学土木建筑工程学院，北京 100044，中国)

0 引 言

土水特征曲线是研究非饱和土学的重要手段，定义了非饱和土的基质吸力和含水量之间的关系(李志清等，2007)，反映了土的持水能力特性。土水特征曲线广泛应用于岩土工程，环境工程以及农业工程。土水特征曲线可用来推断非饱和土的渗透系数、抗剪强度(Vanapalli et al.，1996; 高游等，2017)等性质，对涉及到非饱和土的工程具有重要意义(包承纲，2004)，目前对于土水特征曲线的影响因素开展了许多研究(石振明等，2018；李同录等，2019)。为方便工程应用，很多学者提出不同的拟合方程来描述土水特征曲线，如VG模型(van Genuchten，1980)、FX模型(Fredlund et al.，1994)、Gardner模型(Gardner et al.，1958)等。其中：VG模型由于参数较少、形式简单，是目前使用最广泛的土水特征曲线模型之一。

非饱和土土水特征曲线可通过室内试验获得(Li et al.，2014)，但相关试验一般比较费时且费用高昂，对非饱和土力学的工程应用带来了挑战。另一方面，由于土壤粒径分布曲线和土水特征曲线形状相似，一些学者发现可基于土体粒径分布曲线来预测土体的土水特征曲线。其中：Arya和Paris在1982年(Arya et al.，1982)首次尝试通过粒径分布曲线参数预测土水特征曲线(A & P模型)，将颗粒累积分数转化为含水率，将颗粒半径通过经验系数α转化为孔隙半径，并通过毛细方程计算吸力。一些研究者还在A & P模型的基础上提出一些物理或概念模型对土水特征曲线进行预测(如Mohammadi et al.，2011)。这些研究结果表明，可以基于一些易于获得的物理指标，如粒径分布，通过一定的转换方法来估算土的土水特征曲线，从而实现用较低的代价来获得非饱和土的主要工程性质。上述方法虽然为快速预测土体土水特征曲线提供了一个很好的思路，但本质上是一种经验方法，不可避免地存在预测误差。目前，现有研究中极少对上述误差进行估算，导致不少工程人员直接将由粒径分布曲线估算获得的土水特征曲线作为真实的土水特征曲线，无法考虑土水特征曲线估算误差的影响。

本文将以VG模型为例，提出一种基于土体粒径分布曲线的土水特征曲线概率预测方法。与传统方法相比，本文提出的方法可以全面考虑土水特征曲线的估算误差以及由此带来的不确定性。本文主要分为以下几个部分：首先建立非饱和砂土的数据库，利用极大似然估计法拟合出VG模型的参数；然后，利用粒径分布曲线参数和孔隙率等参数建立VG模型参数的线性方程；进而，通过回归分析建立土水特征曲线的经验公式；最后，对建立的概率分布模型进行验证。本文提出的方法对非饱和土土水特征曲线的预测和非饱和土可靠度分析均具有重要的意义。

1 土水特征曲线模型

土水特征曲线定义了非饱和土的基质吸力和含水量之间的关系，反映了土的持水能力特性。在VG模型(van Genuchten，1980)中，其表达式为：

(1)

式中：S为土的体积饱和度；ψ为吸力；a、m、n均为拟合参数。一般而言，a与土的进气状态有关，n与土的孔径分布有关，参数m与土体特征曲线的整体对称性有关。为简便计算，m可通过参数n基于如下公式进行估算(van Genuchten，1980)：

(2)

基于式(2)，VG模型的模型参数可进一步减少为2个，即a和n。土水特征曲线模型参数可通过最小二乘法等方法来拟合。本文将利用最小二乘法对参数进行拟合。

图1给出了某土样实测土水特征曲线值。基于上述最小二乘法，其最优模型参数值为a=0.032、n=2.15、m=0.46，模型的拟合优度为0.99。图中也同时给出了基于上述参数获得的预测土水特征曲线。由图可知，实测曲线与由VG模型获得的拟合曲线十分接近，由此根据拟合参数a、n、m确定土水特征曲线的参数相比于传统图解法能更精确(潘登丽等，2020)。

图1 土水特征曲线拟合图

2 基于粒径分布曲线的土水特征曲线参数概率预测方法

2.1 土水特征曲线及粒径分布曲线参数的相关性分析

已有研究表明，土体土水特征曲线参数与粒径分布曲线参数具有高度的相关性(Cosby et al.，1984；Saxton et al.，1986)。Li et al.(2014)指出土体土水特征曲线与土体粒径分布曲线的d10、d30、d60、不均匀系数Cu、曲率系数Cc等特征参数具有密切关系。

根据《岩土工程勘察规范》(中华人民共和国国家标准编写组，2009)定义，砂土指粒径大于2 mm的颗粒质量不超过总质量的50%且粒径大于0.075 mm的颗粒质量超过总质量的50%的土。为研究砂土土水特征曲线与土体粒径分布曲线之间的关系，本文首先从美国农业安全部UNSODA数据库(Leij，1996)中选取了100组砂土数据的脱湿土水特征曲线及土壤物理性质参数。图2和图3给出了这100组砂土数据的土水特征曲线及土粒径分布曲线。

图2 土水特征曲线

图3 土样粒径分布曲线

为建立砂土的土水特征曲线经验预测模型，论文首先采用最小二乘法对上述100组土水特征曲线进行了拟合。表1给出了100组砂土的土体粒径分布曲线拟合参数与VG模型参数的相关性分析。

表1 VG模型参数相关性分析

2.2 基于最大似然法的回归分析

首先介绍lna的回归方程标定方法。为方便表述，令X={lnd10，lnd30，lnd60}，令Y=lna。假设Y和X1，X2，X3的关系可用如下方程表示：

Y=b0+b1X1+b2X2+b3X3+ε

(3)

式中：ε为一均值为0、标准差为σ的正态分布随机变量。上述方程中，未知参数包括B={b0，b1，b2，b3，σ}。当B已知的条件下，基于式(3)，Y的均值和标准差分别为：

μY=b0+b1X1+b2X2+b3X3

(4)

σY=σ

(5)

令X1i，X2i，X3i代表第i组土样的lnd10，lnd30，lnd60值，令δi代表第i个土样的Y值。给定B条件下，观测到δi的概率为：

(6)

式中：φ为一标准正态分布的概率密度函数。假设存在s组土样，令δ={δ1，δ2，…，δs}。当B已知的条件下，观测到所有数据的概率为：

(7)

上述表达式即为回归方程参数B的似然函数。根据最大似然性原理(赵军圣等，2010)，使得式(5)取值最大的B即为模型参数的最优值。实际优化过程中，式(5)中由于连乘的存在会导致似然性函数的数值绝对值很小，在优化中容易受到计算误差的影响而导致算法难以收敛。此时，可对似然函数取对数，通过对似然函数对数值取最大值而获得模型参数的最优值。

采用相同的方法，也可对lnn的预测方程进行标定。

2.3 VG模型参数的概率分布预测

下面将采用上述最大似然估计法对lna、lnn进行预测。为方便表述，令θ={lna，lnn}。令μ1、μ2分别代表lna、lnn的均值，令σ1、σ2分别代表lna、lnn的标准差。首先采用最大似然法基于式(7)对lna的回归方程进行标定，可得B={0.32、-0.27、1.5、-0.19、1.49}。将上述参数带入方程(4)，可得μ1的计算公式为：

(8)

基于式(5)，可得σ1=1.49。μ1代表了lna的最可能取值，或最优预测值，σ1衡量了预测lna过程中的变异性。图4a给出了采用式(8)对数据库中100组土样预测lna与实际lna的对比图。通过计算，两者的相关系数为0.79，表示计算值与预测值具有较好的相关性。

图4 预测值与实际值对比

采用最大似然法对lnn进行标定，可得其均值预测公式为：

(9)

lnn的标准差为σ2=0.43。图4b给出了基于式(9)获得的lnn的预测值与实测值的对比图。由图中可知，两者之间的相关系数为0.58。相比而言，lnn预测值与实测值之间的相关性要低一些，表明lnn更难预测。

通过上述分析获得了lna和lnn的边缘概率分布。为获得lna和lnn之间的联合概率分布，还需获得lna和lnn之间的相关性。为获取上述相关性，可对100个土样lna、lnn预测值与实测值之间的残差进行分析。令δ1i、δ2i分别代表第i个土样的lna、lnn的实测值，令e1i、e2i分别代表第i个土样关于lna和lnn的预测误差。第i个土样关于lna和lnn的残差可按下式计算：

(10)

图5a、图5b分别给出了100组土样lna、lnn的预测残差直方图。图中同时给出了根据直方图拟合获得的正态分布曲线。Kolmogorov-Smirnov检验(侯澍旻等，2007)表明，在0.05显著性水平下，lna、lnn的预测残差均服从正态分布。

图5 模型参数残差直方图

图6给出了lna、lnn残差的散点图。由图可知，lna、lnn残差之间呈正相关关系，两者之间的相关系数为0.41，即ρ12=0.41。

图6 ln a-ln n残差散点图

令μθ代表θ的均值，即令Cθ代表θ的协方差。Cθ可按下式计算：

(11)

由于lna、lnn分别服从正态分布，θ={lna，lnn}服从多元正态分布。基于多元正态分布特性，θ的联合概率密度函数为：

(12)

当θ={lna，lnn}服从多元正态分布时，{a，n}服从多元对数正态分布。在实际应用中，可先根据式(12)生成θ={lna，lnn}的样本，再通过指数变换获得{a，n}的样本，由此获得对应的土水特征曲线。

在式(12)中，μθ={μlna，μlnn}代表了θ的最可能取值，由此生成的土水特征曲线为最可能的土水特征曲线；Cθ给出了θ的估算误差，界定了θ可能的取值范围。因此，本文提出的方法不但可提供土水特征曲线的最可能位置，还能提供土水特征曲线可能的取值范围。

3 模型验证

为进一步验证方法的可靠性，本文搜集了文献外30组砂土的脱湿土水特征曲线及其粒径分布曲线，采用本文方法对其土水特征曲线进行了预测。图7给出了其中4组土样的粒径分布曲线，图8给出了对应土样的实测土水特征曲线。

图7 土样粒径分布曲线

图8 土样实测土水特征曲线

以土样1为例，根据该土样粒径分布曲线，其d10、d30、d60，Cu、Cc参数分别为0.09 mm、0.18 mm、0.24 mm、2.72、0.70。将上述参数带入式(8)、式(9)，可得μθ={-1.36，0.47}。该分布的协方差为式(11)，Cθ={2.231，0.266；0.266，0.188}。根据土水特征曲线参数的均值，图9a中给出了土样1的最可能土水特征曲线。由图可知，最可能土水特征曲线与实际土水特征曲线具有相似性，但不完全一致，给定吸力水平下，最可能土水特征曲线与实测特征曲线饱和度最大差别为0.091。

图9 验证土样土水特征曲线及预测曲线

将μθ、Cθ代入式(12)，即可获得土样土水特征曲线参数 {lna，lnn}的分布。根据该分布，获得了 {lna，lnn}的50个样本，由此生成了50条土水特征曲线。作为比较，图中也给出了上述50条土水特征曲线。由图可知，由于lna、lnn预测误差的存在，土水特征曲线也存在多种可能性，根据粒径分布曲线无法唯一确定土水特征曲线。

另一方面，与图2中100条砂土曲线的土水特征曲线相比，图9a中土水特征曲线的离散性有显著降低，且实测土水特征曲线位于预测土水特征曲线范围之内，说明粒径分布曲线显著降低了土水特征曲线的不确定性。

对于图9a中的每一个吸力，50条预测曲线分别对应50个饱和度。根据上述50个饱和度数值，可分别获得该吸力水平下90%的置信区间。改变不同的吸力值，由此可获得不同吸力水平下饱和度的90%置信区间。图9a中比较了实测的土水特征曲线以及不同吸力水平下饱和度的90%置信区间。由图可知，实测土水特征曲线位于预测土水特征曲线的90%置信区间内。

采用类似的方法，图9b、图9c和图9d分别对其他3种土样的实测土水特征曲线与预测土水特征曲线进行了比较。由图可知，另外3种土样的实测土水特征曲线虽与预测最可能土水特征曲线不完全一致，但都位于预测土水特征曲线90%置信区间范围内。

为对本文提出的方法进行进一步的验证，采用本文方法对验证数据库中30个土样的土水特征曲线进行了预测。为方便比较，对每条实测土水特征曲线按其90%的上下界进行了归一化，具体方法介绍如下。令φi代表某一具体的吸力水平，令Si，SiU和SiL分别代表其实测饱和度、实测饱和度90%置信区间上界、实测饱和度90%置信区间下界。与φi对应的归一化饱和度为：

(13)

采用类似的方法，也可对该吸力水平下饱和度的90%置信区间上下界进行标准化。由式(13)可以看出，归一化后饱和度90%置信区间上界的值为1、下界的值为0。按上述方法，图10给出了上述30组土样的归一化的土水特征曲线。由图可知，所有土水特征曲线的归一化吸力值均位于0～1之间，表明所有实测土水特征曲线均位于预测土水特征曲线90%置信区间内。

图10 土水特征曲线置信区间

总体而言，由于土体的孔隙结构对土水特征曲线的形态具有绝对性的影响，土的持水能力是其微观孔隙分布的宏观表现(胡冉等，2013)，因此土的粒径分布为土水特征曲线的预测提供有用的信息，可用来减少土水特征曲线的不确定性。另一方面，即使级配相同的土体，其在不同的应力、干湿循环历史作用下，会形同不同的孔隙结构和储水状态，对土水特征曲线也存在一定的影响(张雪东等，2010)，导致土水特征曲线无法由粒径分布曲线唯一确定，采用粒径分布曲线对土水特征曲线进行预测时，不可避免地存在模型预测误差。

与通过非饱和土试验获得土水特征曲线相比，基于土体粒径分布曲线估算土水特征曲线效率较高，但其伴随的预测误差不可忽视。如前所述，文献中已提出了不少基于粒径分布曲线的土水特征曲线预测方程，但极少对预测方程的误差进行研究，导致基于粒径分布曲线的土水特征曲线预测方程存在被滥用的危险。与已有研究相比，本文提出的方法不但能提供土水特征曲线参数的最优取值，还能给出土水特征曲线的变异范围，为分析土水特征曲线预测误差对非饱和岩土体性能评价的影响、评定预测误差是否位于可接受范围内、以及是否应采纳其他更为精确的土水特征曲线测试和预测方法提供了基础。例如，在边坡稳定性分析中，如土水特征曲线变异范围内边坡始终处于稳定或不稳定状态，表明土水特征曲线的预测误差对边坡是否需要进行加固的决策影响不大，该预测误差可视为可接受；如土水特征曲线在其预测误差范围内变化时，边坡是否需要加固的结论发生显著变化，说明土水特征曲线预测误差对工程决策具有重要影响，此时可考虑采用更精确的方法对土水特征曲线进行预测或测量。

土水特征曲线变异性对非饱和土工程系统的影响也可采用可靠度理论来进行定量分析。如可靠度分析表明土水特征曲线的不确定性不是岩土体性能预测的主导不确定性因素，则预测误差可视为可接受；当土水特征曲线的不确定性对非饱和土系统的性能预测具有重大影响时，可考虑采用更精确的土水特征曲线测试或预测方法。

4 结 论

本文基于100组砂土的土水特征曲线和粒径分布曲线试验数据，分析了土水特征曲线VG模型参数与粒径分布曲线参数的相关性，采用线性回归方法获得了基于粒径分布曲线对VG模型参数的预测方程，通过残差分析建立了土水特征曲线参数的概率预测模型，并基于实测数据对提出的模型进行了验证。研究表明，基于粒径分布曲线无法唯一确定土体的土水特征曲线，采用粒径分布曲线对土水特征曲线进行预测时不可避免地存在模型误差。与已有研究相比，本文提出的方法不但能提供土水特征曲线参数的最优取值，还能给出土水特征曲线的变异范围，为明晰土水特征曲线预测误差对非饱和土工程力学性能分析与评价的影响奠定了基础。