福建-广东地区花岗岩残积土地下水模型参数的敏感性分析 谢永宁、 李坤、 陈志博

《Eng》：Sensitivity Analysis of UH Model Parameters for Granite Residual Soils in the Fujian–Guangdong Region Yongning Xie, Kun Li and Zhibo Chen

　　摘要 本研究收集了来自福建-广东地区的155组花岗岩残余土的测试数据，并应用卡方检验（chi-square test）分析了八个常见物理和力学参数的分布情况。随后使用统一硬化（Unified Hardening, UH）模型对排水三轴试验进行了模拟，并基于土壤性质的分布进行了Sobol全局敏感性分析。结果表明，天然密度和粘聚力大致遵循韦伯分布（Weibull distribution）；孔隙比、液限和塑限遵循对数正态分布（lognormal distribution）；含水量和内摩擦角遵循正态分布（normal distribution）；而塑性指数遵循Gumbel分布。Sobol分析表明，临界状态偏应力主要取决于临界状态应力比（critical state stress ratio, M），临界状态体积应变则同时受到M和正压曲线斜率（λ）的影响。偏应力的整体变化主要取决于M，而体积应变的整体变化主要取决于λ，泊松比（Poisson’s ratio, ν）对土壤应力-应变响应的影响较小。这些发现为福建-广东地区花岗岩残余土的参数选择和数值模拟提供了参考。

1. 引言 花岗岩残余土在中国南部广泛分布。在福建、广东等省份，由于风化条件的差异，导致物理和力学参数存在显著的空间变化性，从而增加了岩土工程设计和分析的不确定性。已有许多研究探讨了岩土参数的统计特性。de Oliveira等人[1]使用多元分析和地质统计方法研究了土壤物理和化学性质的空间分布特征。Paterson等人[2]提出了一种多分数框架来描述岩土参数的空间变异性和随机性。Ching和Phoon[3]构建了基于粘土参数的十维多元概率分布。L?fman和Korkiala-Tanttu[4]整理了芬兰粘土土的多元数据库，并开发了土壤性质转换模型。Kirts等人[5]利用佛罗里达州各地土壤的619组固结试验结果，采用数据驱动方法估算了压缩指数和再压缩指数。尽管取得了这些进展，但专门针对中国南部花岗岩残余土的研究仍然相对较少。Dai等人[6]通过现场调查数据比较了不同地区的物理和力学参数差异。Liang等人[7]对广州东部花岗岩残余土常用的工程参数进行了统计和相关性分析，并建立了参数之间的线性回归模型。Zhang等人[8]基于大量的实验室测试数据，研究了物理和力学参数之间的相关性，并使用两种不同的方法建立了力学参数的多元对数正态分布模型。Yu[9]对广东省花岗岩残余土的物理和力学参数进行了统计分析和分布测试。这些研究主要基于广东省的数据。因此，需要结合多个地区的数据进一步统计分析物理和力学参数，以更全面地了解中国南部花岗岩残余土的工程特性。随着数值分析在岩土工程中的广泛应用，它已成为解决复杂工程问题最重要的基本方法之一[10]。数值模拟的可靠性很大程度上取决于合适的本构模型的选择。例如，新加坡委员会得出结论，Nicoll公路倒塌的一个原因是设计中使用了不合适的土壤本构模型，导致倒塌位置的土壤强度被过高估计[11]。统一硬化（UH）模型非常适合中国的多种粘土类型；它可以有效地模拟土壤行为，包括剪切膨胀、临界状态特性和应力路径依赖性行为。先前的研究表明，在不同加载条件下，使用UH模型计算的结果与实验观测结果吻合良好[12,13,14,15]。然而，UH模型很少用于模拟花岗岩残余土。由于空间变异性，岩土参数难以准确确定，这导致了数值模拟的不确定性。敏感性分析提供了一种定量方法，用于评估参数变化对模型输出的影响，并评估各个参数的相对重要性[16]。Shi等人[17]研究了厦门地区的典型土层，并采用了HS-small模型。通过对某挖掘项目施工过程的有限元模拟，他们对小应变参数对地表沉降和围护墙弯矩的影响进行了敏感性分析。Xie等人[18]研究了HS-small模型性能对关键参数（包括初始剪切模量和临界剪切应变）的敏感性。Guo[19]研究了HS-small模型中高级参数对土壤韧性的敏感性，但仅考虑了单参数变化。Fang[20]以隧道支护设计为例，构建了一个全局敏感性分析的潜在模型。Xu[21]提出了一种基于Sobol序列和最小角回归（Least Angle Regression, LARS）算法的坡度稳定性全局敏感性分析方法。Guo[22]提出了一种结合克里金替代模型（Kriging surrogate model）、蒙特卡洛模拟（Monte Carlo simulation）、全局敏感性分析和一阶可靠性方法（first-order reliability method）的程序，以帮助工程师判断极限平衡方法何时足够准确，以及何时需要进行概率分析时的强度折减方法。Pandit等人[23]对四种具有不同潜在破坏模式的岩石边坡案例进行了敏感性分析，发现全局敏感性方法通常比局部方法更准确。Kumar和Tiwari[24]提出了一种不精确的、与矩无关的全局敏感性分析方法，并通过岩土工程案例证明了其适用性。总体而言，大多数关于本构模型参数敏感性的现有研究都集中在单参数分析上，而考虑参数相互作用的多参数敏感性分析仍相对有限。

3. 敏感性分析方法 3.1. UH模型 Yao等人[30]开发了统一硬化（Unified Hardening, UH）模型，该模型是在改进的Cam-Clay（MCC）模型基础上发展而来的。UH模型引入了一种变换应力方法和一个不依赖于应力路径的统一硬化参数[31]。基于从SMP准则[32]推导出的变换应力概念，UH模型可以表征一般三维应力状态下的土壤强度和变形。统一硬化参数使模型能够再现不同加载条件下的多种土壤力学行为，包括剪切膨胀、超固结行为、应变硬化和软化以及应力路径依赖性行为。UH模型参数列于表3中。

4. 物理和力学参数的分布 4.1. 花岗岩残余土的物理和力学参数统计 由于成因相似，中国东南沿海地区的花岗岩残余土的参数变化相对较小[9]。因此，本研究重点关注来自福建-广东地区的花岗岩残余土。本研究共收集了155组来自福建-广东地区的花岗岩残余土数据集。数据集包括常用的工程参数，如天然含水量??、天然密度ρ、孔隙比??0、液限????、塑限????、塑性指数IP、粘聚力??、内摩擦角??和压缩指数Cc。相应的统计结果见表1。表1显示，天然密度的变异系数（COV）为0.07，表明其变异性较低。相比之下，其余八个参数的COV均超过0.15，表明变异性较大。与Zhang等人[25]关于马来西亚花岗岩残余土的已发布结果相比，本研究的物理和力学参数在均值和COV上存在差异。这些差异可能与区域风化程度和地质条件的差异有关，突显了特定区域参数收集和统计描述的必要性。因此，本研究的统计结果有助于表征福建-广东地区花岗岩残余土的基本工程特性和区域特征，并为后续的UH模型分析提供了岩土工程背景。 4.2. 分布拟合的卡方检验 根据Cheng等人的方法[26]，我们选择了正态分布（Normal distribution）、Gumbel分布、韦伯分布（Weibull distribution）和对数正态分布（Lognormal distribution）来拟合收集的数据。Karl Pearson[27]提出的卡方（??2）拟合优度检验用于评估每个候选分布的适用性。该检验量化了观察到的样本频率与从假设概率分布得出的理论频率之间的差异，从而确定样本数据是否符合假设的分布。为了减少区间选择的主观性，使用Sturges规则[?? =1 +3.3?log???]来确定 class intervals 的数量。用????和????分别表示第??个区间内的观察频率和预期频率，卡方统计量表示为： ??2=∑(?????????)2???? 对于给定的显著性水平，临界值??2??可以从统计表中获得。根据岩土工程可靠性分析的常规做法[28,29]，本研究采用了?? =0.05的显著性水平。当??2 ≤??2??时，接受假设的概率密度函数；否则拒绝该函数。如果有多个候选分布满足接受标准，则选择??2/??2??值最小的分布作为最终分布。 4.3. 分布拟合结果 各候选概率密度函数的??2、??2??和接受结果总结在表2中。相应的直方图、核密度估计（KDE）和采用的分布图显示在图1中。值得注意的是，只有43个样本可用于压缩指数的分析，这限制了统计的可靠性。因此，压缩指数的统计分析结果仅用于参考。 4.4. 卡方检验的适用性和局限性 本研究使用概率分布拟合来表征福建-广东地区花岗岩残余土的物理和力学参数的概率特性，并为后续的全局敏感性分析中的拉丁超立方抽样（Latin Hypercube Sampling, LHS）提供输入分布函数。因此，主要要求是拟合分布能够合理捕捉每个参数的总体统计特性。需要注意的是，不同的拟合优度检验可能会因为基于不同的统计假设和检验统计量而得出略有不同的结果。卡方检验基于分组频率差异，其结果可能会受到样本分组和区间划分的影响。相比之下，Kolmogorov–Smirnov（K-S）和Anderson–Darling（A-D）检验基于累积分布函数。出于本研究的目的，卡方检验被用作识别合适概率分布的实用主要方法，而补充的K-S检验结果和比较信息则提供在附录A中。索博尔灵敏度方法（Sobol Sensitivity Method）是一种广泛使用的全局灵敏度分析技术，它源自方差分析（ANOVA）框架内的方差分解方法。该技术通过评估各个输入变量对模型输出变异性的影响程度来确定参数的灵敏度。该方法中的主要灵敏度度量是一阶灵敏度指数（first-order index）和总效应指数（total-effect index）。

考虑一个表示为 \(Y = f(X_1, X_2, \ldots, X_n)\) 的模型（公式1），其中 \(Y\) 是模型输出，\(X = (X_1, X_2, \ldots, X_n)\) 表示 \(n\) 个输入变量的向量。总输出方差可以分解为：
\[Var(Y) = \sum_{i=1}^n W_i + \sum_{i < j} W_i W_j + \ldots + W_{1,2, \ldots, n}\]
（公式2）
其中 \(W_i\) 表示 \(X_i\) 的单独效应的贡献，\(W_{ij}\) 表示 \(X_i\) 与 \(X_j\) 之间的交互作用的贡献，更高阶项表示多个变量之间的相互作用。

第 \(i\) 个变量 \(X_i\) 的一阶灵敏度指数定义为：
\[S_{i} = \frac{W_i}{Var(Y)}\]
（公式3）
该指数量化了 \(X_i\) 对输出方差的直接贡献，排除了交互作用效应。

第 \(i\) 个变量 \(X_i\) 的总效应指数定义为：
\[S_{iT} = 1 - \frac{Var(Y - X_i)}{Var(Y)}\]
（公式4）
其中 \(Var(Y - X_i)\) 表示在 \(X_i\) 固定的情况下，所有其他变量变化时的模型输出方差。\(S_{iT}\) 测量了 \(X_i\) 的总体贡献，包括其主要效应及其与其他变量的相互作用。

较大的 \(S_i\) 指数表明参数对模型响应有较强的独立影响。相反，相对较小的 \(S_{iT}\) 暗示参数的总体影响有限，可能认为其对模型输出的重要性较低。

3.3 参数分布
在对UH模型进行索博尔灵敏度分析之前，需要定义每个输入参数的概率分布。分布类型的选取基于表1和表2中总结的统计结果。
\(M\) 的计算公式为：
\[M = \frac{6s_i n\mu_0}{\sigma_0^2 - s_i n\mu_0^2}\]
并进行了分布测试。\
\(\mu_0\) 可通过 \(\mu_0 = 0.434\chi_c\) 获得[33]。根据Gao [34] 的研究，典型粘土的 \(\mu_0\) 与 \(\sigma_0\) 的比值介于0.1到0.2之间。出于采样目的，本研究中的 \(\mu_0\) 取值范围大约为 \([0.15\mu_0 m, 0.15\mu_0 a]\)，而不是将 \(\mu_0\) 视为每个采样实现的固定值0.15\(\lambda\)。在此基础上，统计参数 \(\kappa\) 取值为 \(\kappa = 0.15\mu_0\) 和 \(\sigma_2 = 0.0225\sigma_0^2\)。\(\nu\) 的取值参考了Foti等人[35]的研究。假设 \(\lambda\)、\(\kappa\) 和 \(\nu\) 遵循正态分布。UH模型灵敏度分析中采用的参数范围及其相应的分布列在表4中。

3.4 参数采样方法
在进行全局灵敏度分析时，必须根据模型参数的分布生成一组具有代表性的输入样本。为了提高采样效率并更好地覆盖参数空间，采用了拉丁超立方采样（Latin Hypercube Sampling，LHS）[36]方法。
在LHS方法中，每个参数的范围被划分为等概率且不重叠的区间，然后从中随机抽取样本。这种分层采样策略确保了对每个变量整个范围的系统探索，与传统的随机采样相比，在高维参数空间中产生更均匀分布的样本集。

4. 结果与讨论
4.1 基线数值模型
本研究使用Python 3.10对UH模型进行了数值实现。为了验证实现的正确性及其适用性（特别是对于花岗岩残余土），将数值预测结果与花岗岩残余土的三轴试验数据[37]进行了比较。如图2所示，试验数据和数值预测在应力-应变关系和体积应变演变方面有很好的一致性。基于上述验证，可以认为UH模型适用于模拟花岗岩残余土的力学行为。作为灵敏度分析的基准数值模型，采用了在 \(\p_0 = 100\ kPa\) 和 \(OCR = 4\) 条件下的花岗岩残余土的排水三轴压缩试验。输入参数取为收集数据的平均值，表5总结了模型中使用的参数值。基准模型对应的模拟应力-应变关系如图3所示。

4.2 样本量收敛性分析
全局灵敏度分析的结果与样本量 \(N\) 密切相关。由于不同的分析方法和模型问题可能需要不同的样本量，因此有必要对灵敏度结果进行收敛性分析。
在本研究中，选择了均方误差 \(\|\hat{y}_i - y_{ref}_i\|^2\) 作为代表性的模型输出来进行收敛性分析，因为这个输出受到输入参数之间相对强相互作用的影响，在索博尔分析中对应较大的样本量要求。初始样本量设置为500，然后以500为增量递增，以检查灵敏度指数的收敛行为。如图4所示，当 \(N = 1500\) 时，\(S_i\) 和 \(S_{iT}\) 均基本稳定，进一步增加样本量也不会观察到显著变化。这表明灵敏度结果已收敛且足够稳定，适用于当前的分析。因此，本研究选择 \(N = 1500\) 作为索博尔灵敏度分析的样本量。

4.3 灵敏度分析结果
本研究选择不同的物理量作为模型输出，并对UH模型参数进行了多次索博尔灵敏度分析，以研究在不同加载条件下参数灵敏度的变化情况。
4.3.1 以临界状态偏应力 \(q_c^s\) 作为模型输出
在以强度为主要考虑因素的工程场景中（如边坡稳定性和隧道挖掘），\(q_c^s\) 特别重要。UH模型参数的相应索博尔灵敏度指数如图5所示。\(M\) 的 \(S_i\) 和 \(S_{iT}\) 分别为0.9684和0.9843。相比之下，其他参数的 \(S_i\) 和 \(S_{iT}\) 接近于零。

4.3.2 以临界状态体积应变 \(\Delta_v,_c^s\) 作为模型输出
在变形主导性能的工程场景中（例如循环或动态加载条件下），\(\Delta_v,_c^s\) 特别重要。相应的索博尔灵敏度指数如图6所示。\(M\) 的 \(S_i\) 和 \(S_{iT}\) 分别为0.5262和0.7311；\(\lambda\) 的 \(S_i\) 和 \(S_{iT}\) 分别为0.3978和0.5418；\(\epsilon_0\) 和 \(\kappa\) 的 \(S_i\) 和 \(S_{iT}\) 接近于零。

4.3.3 以偏应力均方误差 \(MSE\) 作为模型输出
均方误差（MSE）用于量化模拟输出与参考结果之间的总体差异。其定义如下：
\[MSE = \frac{1}{N} \sum_{i=1}^N (y_{ref}_i - y_{sim}_i)^2\]
（公式5）
其中 \(y_{ref}_i\) 和 \(y_{sim}_i\) 分别表示第 \(i\) 个参考值和模拟值，\(N\) 是数据点的总数。
在此分析中，\(y_{ref}_i\) 是从基准数值模型获得的第 \(i\) 个 \(q\) 值。使用 \(q\) 的MSE作为模型输出可以反映整体应力偏差随参数变化的情况，以及它们与基准模型的偏差程度。得到的索博尔灵敏度指数如图7所示。\(M\) 的 \(S_i\) 和 \(S_{iT}\) 分别为0.8977和0.9245，而 \(\mu_0\) 和 \(\kappa\) 的 \(S_i\) 和 \(S_{iT}\) 分别为0.0812和0.1181。其他参数的相应指数接近于零。

4.3.4 以体积应变 \(\Delta_v\) 的均方误差 \(MSE\) 作为模型输出
与4.3.3节一致，采用MSE作为模型输出；在这种情况下，\(y_{ref}_i\) 是从基准数值模型获得的第 \(i\) 个 \(\Delta_v\) 值。使用 \(\Delta_v\) 的MSE作为模型输出可以反映整体体积应变随参数变化的情况，以及它们与基准模型的偏差程度。得到的索博尔灵敏度指数如图8所示。\(M\) 的 \(S_i\) 和 \(S_{iT}\) 分别为0.7113和0.9244，而 \(\epsilon_0\) 和 \(\kappa\) 的 \(S_i\) 和 \(S_{iT}\) 分别为0.1632和0.2269。\(E_0\) 和 \(\kappa\) 的灵敏度指数较低，\(\nu\) 的 \(S_i\) 和 \(S_{iT}\) 接近于零。

4.4 结果讨论
(1) 当以 \(q_c^s\) 作为模型输出时，\(M\) 显示出最高的灵敏度。根据临界状态土力学，我们知道 \(q_c^s = M\cdot p_c^s\)；在排水三轴条件下，我们有 \(\Delta q = 3\Delta p\)，因此可以确定临界状态下的约束压力 \(p_c^s = 3p_0/(3 - M)\)。因此，一旦 \(p_0\) 和 \(M\) 被确定，\(q_c^s\) 也就确定了。这种基本关系解释了索博尔灵敏度分析中识别出的 \(M\) 的主导贡献。
(2) 当以 \(\Delta_v,_c^s\) 作为模型输出时，根据临界状态土力学，临界状态孔隙比通过临界状态线（CSL）与 \(p_c^s\) 唯一相关。CSL 几乎与 NCL 平行，因此 \(\lambda\) 和 \(M\) 共同影响 \(\Delta_v,_c^s\)。
(3) 当以 \(q\) 的MSE作为模型输出时，\(M\) 仍然是最敏感的参数。MSE衡量了曲线的整体偏差，这反映了从初始状态到临界状态的整个过程中应力-应变曲线的演变差异。\(M\) 决定了屈服面的形状，从而控制了加载过程中的 \(q\) 变化。因此，\(M\) 显著影响这个模型输出。
(4) 当以 \(\Delta_v\) 的MSE作为模型输出时，\(\lambda\) 显示出最高的灵敏度。\(M\) 和 \(\kappa\) 决定了体积模量，从而控制了土壤从初始状态到临界状态的体积变形特性。由于这种演变主要是由塑性变形主导的，因此 \(\lambda\) 对模型输出有更强的影响。
(5) 在四个模型输出中，\(\nu\) 的 \(S_i\) 和 \(S_{iT}\) 都接近于0，表明它对模型输出几乎没有影响。这是因为模型输出主要由塑性变形产生大的变形，而 \(\nu\) 主要控制弹性变形。
(6) 灵敏度结果对UH模型的校准和实际应用有直接的影响。在校准过程中应优先考虑高灵敏度的参数 \(M\) 和 \(\lambda\)，因为这些参数的微小变化可能导致预测输出的显著变化。具体来说，对于偏应力响应，应更加仔细地校准 \(M\)；而在体积变形是主要关注点时，应对 \(\lambda\) 进行更严格的约束。相比之下，对于所有考虑的输出，\(\nu\) 的灵敏度接近于零，表明从测量数据中可靠识别它的信息有限。因此，\(\nu\) 可以固定在一个经验值或文献推荐的值，而不是作为关键的校准参数。这种基于灵敏度的策略有助于减少需要密集校准的参数数量，降低校准不确定性，并提高UH模型在实际工程应用中的预测可靠性。

5. 结论
(1) 对 Fujian–Guangdong地区花岗岩残余土的八个常用物理和力学参数的统计评估表明，在常规显著性水平 \(\alpha = 0.05\) 下，大多数参数偏离正态分布。天然密度和粘聚力大致遵循Weibull分布；孔隙比、液限和塑限遵循对数正态分布；含水量和内摩擦角遵循正态分布；塑性指数遵循Gumbel分布。这些结果突显了花岗岩残余土的非正态性和参数变异性。这种区域统计特征为Fujian–Guangdong地区的参数选择提供了依据。
(2) 应注意，不同的拟合优度测试可能会因为强调分布拟合的不同方面而产生略微不同的分布选择。在本研究中，使用卡方检验（chi-squared test）作为主要方法来描述区域参数的整体概率特征，并为后续基于LHS的分析定义输入分布。其适用性和局限性在第2.4节中进行了讨论，补充的K-S检验结果见附录A。(3)敏感性分析表明，当UH模型应用于福建-广东地区的花岗岩残积土时，参数??和λ在控制模拟结果中起主导作用。??主要影响???????以及偏应力的整体演变，而λ主要控制体积应变的演变。这两个参数共同决定了????,?????。因此，为了提高数值预测的可靠性，应根据主导的工程条件来调整??和λ的校准精度要求。相比之下，ν的敏感性可以忽略不计，在实际工程应用中可以使用传统的经验值进行赋值。