0.   引言

路基回弹模量是路面结构设计中不可或缺的设计参数之一,直接影响了路面设计的厚度等^[1-3]。它表征了路基土在弹性变形阶段内,在垂直荷载作用下,抵抗竖向变形的能力。通常,回弹模量可以通过室内试验或现场试验直接获得。然而,由于原状软土土样的获取较为费时费力,因此,通过大量试验获得软土路基的回弹模量并非易事。工程实践中,工程师通常根据路基回弹模量与其他容易获得的土性参数（如静力触探试验CPT数据、含水率、干密度、液限、塑限等）之间的关系,间接估计土体的回弹模量,进而进行设计与分析^[4]。

Mohammad等^[5]研究发现,软土路基回弹模量M_r与CPT锥尖阻力q_c、侧摩阻力f_s、土体含水率w以及干重度${\gamma }_{\text{d}}$有着密切联系。因此,本文重点在于构建回弹模量M_r与上述四个指标之间相关关系的数学模型。与其他土体参数相比,CPT数据几乎是完全连续的^[6-10],且土体含水率w以及干重度${\gamma }_{\text{d}}$数据的获取也较为简单。最为常见的经验模型之一即假定M_r与其他土体参数的方程形式,通过最小二乘法获得对应的参数。虽然该方法直接、简单、可操作性强、,但由于方程形式假定的局限性,在不同工程场地的预测效果往往难以尽如人意^[5,11]。近年来,国内外研究人员提出基于多维正态概率分布的模型以期解决此类问题。目前,该方法在岩石、砂土、黏土甚至边坡风险预测等领域均有所应用,且取得良好效果^[12-15]。然而,该方法对数据量的要求较高,而工程实践中数据量往往相对较少,难以满足要求。此外,该方法要求不同的岩土参数试验数据能够转换成近似服从正态分布的数据,且转换后的试验数据满足线性相关或近似线性相关。该假设通常与实际工程数据有所出入。

为解决这一问题,本文提出基于高斯过程回归（gaussian process regression,简称GPR）的方法对江苏省9个城市的路基回弹模量与上述四个土体参数进行建模,从而实现在给定静力触探锥尖阻力q_c、侧摩阻力f_s、土体含水率w以及干重度${\gamma }_{\text{d}}$的条件下,对路基回弹模量进行合理估计,并计算估计的标准差,量化估计的可靠性。该方法一方面对数据量的需求相对较小,在小样本数据下即可获得较高的训练精度以及较强的泛化能力；另一方面,该模型具有较强的解译性,能够准确地反映出回弹模量与不同土体参数间依赖关系的强弱,体现了数据驱动的特点。此外,高斯过程回归不需要将输入参数进行非正态空间到正态空间的转换以及保证两两参数之间的线性关系。本文通过已有的江苏部分地区的软土路基参数数据库,系统进行了大量的计算与检验,验证了基于高斯过程回归的模型在预测软土路基回弹模量上的准确性、高效性以及鲁棒性。

1.   软土路基参数特征分析

为全面反映江苏省整个区域的软土路基特性,本节收集了文献[4]中江苏省9个城市、15个工程场地的M_r,q_c,f_s,w以及${\gamma }_{\text{d}}$数据,共计124组。其中,M_r根据美国AASHTO标准测量,q_c和f_s由现场圆锥静力触探试验获得,${\gamma }_{\text{d}}$由室内试验常规方法获得。表1总结了五个软土路基参数的均值、标准差、最大值和最小值。最大值、最小值反映了依据此数据集构建方法的适用范围,即12.5 MPa≤M_r≤95.8 MPa,0.22 MPa≤ q_c≤3.93 MPa,0.03 MPa≤f_s≤0.14 MPa,6.90%≤ w≤78.1%,10.50 kN/m³≤${\gamma }_{\text{d}}$≤19.90 kN/m³。此外,图1绘制了每个参数的直方图,以及任意两参数之间关系的散点图。如图1所示,大多数参数并不严格服从正态分布,且两两之间亦非全部线性相关（如回弹模量M_r与含水率w的散点图）。这也解释了文献[4]采用Johnson’s system distribution或者Box-Cox变换将非正态分布转换为正态分布的原因。但是,上述的转换方法无法保证转换后变量之间的线性相关性,如前文所述。为了解决参数间的非线性问题,q_c,f_s,w,${\gamma }_{\text{d}}$的非正态性以及对数据量的依赖性,本文采用高斯过程回归的方法对已有少量数据进行建模,详见第2节。

表  1  软土回弹模量（M_r）、锥尖阻力（q_c）、侧摩阻力（f_s）、含水率（w）和干重度（${\gamma }_{\text{d}}$）的统计特征

Table  1.  Statistics of resilient modulus (M_r), tip resistance (q_c), sleeve friction (f_s), water content (w) and dry unit weight (${\gamma }_{\text{d}}$) of subgrade soft soils

软土参数	均值	标准差	最小值	最大值
Mr/MPa	46.09	17.28	12.50	95.80
qc/MPa	1.76	0.85	0.22	3.93
fs/MPa	0.09	0.03	0.03	0.14
w/%	32.00	14.22	6.90	78.10
${\gamma }_{\text{d}}$/(kN·m-3)	15.86	2.07	10.50	19.90

下载: 导出CSV 

| 显示表格

图  1  软土回弹模量（M_r）、锥尖阻力（q_c）、侧摩阻力（f_s）、含水率（w）、干重度（${\gamma }_{\text{d}}$）的直方图与散点图

Figure  1.  Histograms of resilient modulus (M_r), tip resistance (q_c), sleeve friction (f_s), water content (w) and dry unit weight (${\gamma }_{\text{d}}$) and scatter plot among them

下载: 全尺寸图片 幻灯片

2.   高斯过程回归建模

2.1   基本原理

高斯过程回归方法是一种非线性的、基于贝叶斯思想的无参推断方法,对于高维数据、小样本数据等具有良好的适应性^[16-17]。近年来,该方法在许多领域得到了成功应用。高斯过程回归方法的基本原理是假定因变量的先验分布,如本文中的软土路基回弹模量M_r = f($q_{\text{c}}$,f_s,w,${\gamma }_{\text{d}}$),服从多维正态分布,即f(X)～GP(μ(X),k(X,X'))。其中,f ($\cdot$)表示未知函数。本文中X = (q_c,f_s,w,${\gamma }_{\text{d}}$)；μ(X)表示f(X)的均值向量,k(X,X')表示f(X)与f(X')之间的协方差矩阵。通常,可通过预处理使μ(X) = 0。k(X,X')可以表示为^[17]

式中 X₁ =$q_{\text{c}}$,X₂ =$f_{\text{s}}$,X₃ = w,X₄ =${\gamma }_{\text{d}}$；${\lambda }_i^{}$表示X_i (i = 1,2,3,4) 坐标方向的相关长度,反映了f(X)在X_i坐标方向上${\lambda }_i^{}$范围内具有较高的相关性。

基于高斯过程的基本特性,f(X)的子集的先验分布仍然服从多维正态分布。假设列向量y表示f(X)的m个观测值,对应m组q_c,f_s,w,${\gamma }_{\text{d}}$的已有数据。那么,根据上述讨论,y的先验分布同样服从GP(0,$K(X_{\text{y}},{X^{\prime }}_{\text{y}})+{\sigma }_{\epsilon }^2{\rm I}$),即均值为0,协方差为$K(X_{\text{y}},{X^{\prime }}_{\text{y}})+{\sigma }_{\epsilon }^2{\rm I}$的高斯分布。${\sigma }_{\epsilon }^2$表示y的试验误差。令X^*表示待预测点,那么,y和f(X^*)的联合分布为

根据高斯过程的解析特性和贝叶斯框架,在给定y的条件下,f(X^*)的后验分布仍然服从高斯分布,其均值和协方差分别为

式（3）表示f(X^*)的预测均值${\mu }_{f({\text{X}}^{*})}$可看作已有试验数据y的权重叠加。权重系数α的分配取决于待预测数据X^*与已知测量数据X_y之间的距离。由式（1）可知,距离越近,相似度越高,对应的权重系数越大。若待预测X^*只是一组(q_c,f_s,w,${\gamma }_{\text{d}}$)数据,那么$Va{r}_{f({\text{X}}^{*})}$则表示预测回弹模量的方差。方差越大,表示预测结果越不准确,反之,则说明预测结果较为合理。若X^*表示多组数据,那么${\mu }_{f({\text{X}}^{*})}$为均值向量,$Va{r}_{f({\text{X}}^{*})}$为协方差矩阵,其对角线表示预测结果的方差。毫无疑问,预测结果取决于式（1）中的超参数${\lambda }_i$ (i = 1,2,3,4)、${\sigma }_f^2$以及${\sigma }_{\epsilon }^2$。超参数不同,预测结果有所不同。因此,如何确定最为合理的超参数值至关重要,详细讨论见2.2节。

2.2   高斯过程回归超参数的估计

为实现超参数最优值的估计,本文采用最大似然估计的方法解决这一问题。最大似然估计的核心是寻找最优的超参数使获得测量数据y的可能性最大。获得测量数据y的可能性可以由其概率密度函数反映,可表示为：

因此,可以通过最大化式（4）或者等价的最小化-lnP(y|X_y)估算出超参数的最优值,即^[17]

为了寻找式（5）的最小值,可以采用多种优化算法,如牛顿法、共轭梯度法等。此外,为了降低初始值对最终优化结果的影响,可以对多组初始值进行试运算,最终选取一组对训练结果估算最为准确的变量作为优化的初始值。

3.   工程案例应用

3.1   高斯过程回归模型构建

为了说明高斯过程回归在江苏软土路基回弹模量预测上的适用性,从文献[4]收集的124组数据中随机抽取30组数据以训练高斯过程回归模型。为使该模型预测更加准确、适用范围更广,在构建模型前,首先对30组X_i (i = 1, 2, 3, 4) 数据进行归一化处理,即

式中,${\widehat{\mu }}_{X_i}$和${\widehat{\sigma }}_{X_i}$分别表示m=30组X_i数据的均值与标准差,$X_i{}^{\prime }$对应X_i转换后的数据。式（6）是一种线性转换,不会改变与回弹模量M_r的相关性。因此,回弹模量M_r仍然与归一化锥尖阻力$X_1{}^{\prime }$ = q_c和软土含水率$X_3{}^{\prime }$ = w相关性更大。然后,利用m=30组M_r和归一化后的(q_c,f_s,w,${\gamma }_{\text{d}}$)数据优化式（5）,得到高斯过程回归模型${\lambda }_i$的最优值分别为${\lambda }_1$= 5.43,${\lambda }_2$= 31.07,${\lambda }_3$= 5.58,${\lambda }_4$ = 17.65。由式（1）可知,k(X,X)主要取决对应于$\lambda$较小的参数。因此,$\lambda$越小,理论上对M_r的预测贡献愈明显。由训练结果可知,${\lambda }_1$（对应$q_{\text{c}}$）和${\lambda }_3$（对应w）相对较小,即q_c与w对M_r的预测贡献更大。这与图1所示的信息一致,间接反映了高斯过程回归模型的数据驱动特性。

3.2   高斯过程回归模型验证

为了检验参数估计是否合理,将上述参数代入式（2）重新预测30组训练数据的M_r,所得结果与M_r试验值的对比如图2（a）所示,对应的决定系数R² = 0.97。为了检验模型的泛化能力,即给定新的X_i时M_r预测的准确性,这里基于剩余的94组X_i数据（需要归一化处理）,可得到回弹模量M_r的预测均值,如图2（b）所示。M_r的预测值与试验值的决定系数R² = 0.96。无论在训练集或是在验证集中,这组试验中决定系数R²均大于0.95,这表示基于高斯过程回归模型的预测具有较高的准确度。为了更好地对比预测结果,图3绘制了M_r预测值与试验值的对比图,实线表示从小到大M_r的试验值,长虚线表示由高斯过程回归模型得到的M_r的预测值,短虚线表示预测均值±1.96预测标准差,对应了M_r预测值的95%置信区间。由图3可知：①M_r预测值与试验值在绝大多数情况下非常匹配；②M_r试验值绝大多数（88/94= 93.6%）落在由高斯过程回归模型所得的95%置信区间内。这些结果均表示基于高斯过程回归模型回弹模量的预测比较合理,且量化的不确定性符合预期。

图  2  软土回弹模量预测值与真实值的散点图

Figure  2.  Scatter plot between predicted and measured resilient moduli

下载: 全尺寸图片 幻灯片

图  3  软土回弹模量预测值与测量值的对比

Figure  3.  Comparison between predicted and measured resilient moduli in an ascending order

下载: 全尺寸图片 幻灯片

然而,上述结果（图2,图3）只是基于一次随机试验。为了更好地验证本文中所用高斯过程回归模型的鲁棒性,重复上述随机试验100次：即随机将124组数据分成m=30的训练集以及m=94的验证集,然后每次训练后可以计算验证集对应的R²。因此,可得对应于100次验证集的100个R²。图4绘制了100个R²的箱型图。由图4可以看出,基于30组数据训练的高斯过程回归模型在验证集上的预测值对应的R²基本全部高于0.92。这说明高斯过程回归模型具有较高的鲁棒性,且可以准确地基于（q_c,f_s,w,${\gamma }_{\text{d}}$）数据预测软土地基回弹模量。

图  4  100次随机试验中验证数据集中高斯过程回归方法所得决定系数R²的箱型图

Figure  4.  Box-plot of R² for validation dataset under four different scenarios with 100 random experiments

下载: 全尺寸图片 幻灯片

此外,为了讨论训练数据量m对预测结果的影响,分别取m = 45,60,75对高斯过程回归模型进行训练。具体过程与前文m = 30一致,均随机进行100次试验。对应的结果仍然以箱型图的形式呈现,如图4所示。由图4可知,随着训练集数据m的增加,高斯过程模型泛化能力逐渐提高,预测结果愈加准确,这与预期完全一致。同时,训练集数量增加后,高斯过程回归模型获得了更为准确的超参数,因此,预测结果更为合理、鲁棒性更好。由图4可以看出,当m ≥ 60时,R²的最小值大于0.95；且相比m=30和m=45,100次随机试验结果的变异性更小。

4.   结论

本文在已有江苏软土地基数据集的基础之上,提出了高斯过程回归的建模方法。通过对不同样本数目的训练集的建模分析,验证了高斯过程回归方法在江苏软土路基回弹模量预测上的准确性、高效性以及鲁棒性。主要得出以下结论：

（1）高斯过程回归方法可以准确地反映江苏软土回弹模量与不同土参数间的相互关系,体现了高斯过程回归的数据驱动特性。

（2）小样本条件下,高斯过程回归模型针对江苏软土路基回弹模量的预测值与试验值之间的拟合程度（R²）最高可达到0.97,最低亦可达0.93。

（3）高斯过程回归模型的准确度、鲁棒性均随着训练样本数目的增大而不断提高。当训练数据由30增长到75时,100次随机试验中,验证集中预测值R²的均值由0.957提升至0.975,标准差由0.011降低到0.005。