周苏华,周帅康,张运强,聂志红,雷 瑜
(1.湖南大学土木工程学院,湖南 长沙 410082;2.建筑安全与节能教育部重点实验室,湖南 长沙 410082;3.中南大学土木工程学院,湖南 长沙 410075;4.中国路桥工程有限责任公司,北京 100010)
肯尼亚是我国“一带一路”倡议在非洲的重要支点。近年来,两国在交通、基础设施建设等方面展开了重要的合作。膨胀土(当地俗称“黑棉土”)是肯尼亚境内常见的土壤类型,具有遇水膨胀软化、失水收缩开裂等特性,对当地公路、铁路等项目的建设影响很大[1]。胀缩性是膨胀土最明显的特征,预测膨胀土的胀缩等级是治理工作中的首要问题[2]。工程实践表明,胀缩等级的误判已成为膨胀土工程灾害的主要根源[3−4],低估或高估膨胀土的胀缩等级都可能引起工程安全事故并造成巨大的经济损失。因此,如何准确的预测膨胀土的胀缩等级是目前亟需解决的问题。
现有的膨胀土胀缩等级分类方法研究成果较多,主要有直觉模糊集法[3]、模糊均值法[5]、距离判别法[6]、集对分析法[7]等。不同分类方法会产生不同的分类结果[8],至今还没有一种方法能够被全面的推广。《膨胀土地区建筑技术规范》(GB 50112-2013)[9]按照自由膨胀率大小和地基变形量对膨胀土的胀缩等级进行了划分,以指导膨胀土地区相关工程建设工作。
目前,有关膨胀土的胀缩等级预测研究取得了较为丰硕的成果。黄卫等[10]较早地采用模糊评判方法提出了膨胀土的胀缩等级模糊评判理论模型。之后,一些学者分别采用模糊神经网络[11]、可拓学以及期望理论等[12−13]方法建立了膨胀土胀缩等级评价模型。但上述模型均或多或少有等级划分跳跃、主观性过大以及易遗漏重要约束条件等缺点。近年来,随着计算机学科的发展和相关数学理论的不断完善,为解决膨胀土的胀缩等级预测提供了新思路。VAPNIK[14]在1995年首次提出了支持向量机SVM(Support Vector Machine)模型并最先用于模式识别领域,YANG 等[15]基于SVM 分析原理,利用5 个土体指标建立了膨胀土的分类模型。之后,VAPNIK 等[14]将ε 不敏感损失函数加入到模型中,提出了用于非线性回归估计的支持向量回归机SVR(Support Vector Regression)模型[16]。随着研究的深入,SVR 以良好的性能、可靠的预测结果及核函数的多样性等特点在各个领域被广泛应用。诸如,张骏等[17]将SVR 模型应用于驾驶简单反应时间的预测,结果显示准确率达80%以上;王芬等[18]建立基于MABCSVR 的边坡安全系数预测模型,预测结果最大相对误差为7.62%;林荣安等[19]在SVR 模型的基础上结合粗糙集理论进行地铁隧道工程的地表沉降预测,结果显示相对误差均小于15%等等。
基于此,为了准确预测膨胀土的胀缩等级,并为处理肯尼亚等相关地区膨胀土工程问题提供理论依据,提出了一种基于支持向量回归机SVR 模型的膨胀土胀缩等级预测方法。通过分析肯尼亚“蒙内铁路”沿线膨胀土的土工试验数据,选取土体自由膨胀率作为预测指标,选择相关性显著地参数作为训练参数,以“预测准确率”作为模型预测效果的评价指标,并运用多个核函数构建了基于SVR 模型的膨胀土胀缩等级预测模型。研究结果为肯尼亚及相关地区工程建设提供依据。
本文数据来源于“肯尼亚蒙巴萨—内罗毕标准轨铁路”(简称“蒙内铁路”)(图1)。蒙内铁路正线全长471.65 km,工程建设场地主要含蒙脱石黏土和含高岭土黏土。通过现场实际调研,膨胀土主要集中分布在DK263+300~DK448+000 标段,区域内土体表观较细腻,具有失水易收缩、易干裂、易剥落以及遇水膨胀等特性,一般情况下裂缝宽度在3~5 cm,最大裂缝可达10 cm。图2为现场采集的3 种主要土样,具体特征:(1)土样1 表观为黑色,天然状态下粒径大小不一,大直径土颗粒为小直径土颗粒形成的团聚体,节理裂隙发育,风化碎屑较多,粗颗粒表面粗糙,坚硬,手不易捏碎,干强度高。(2)土样2 表观为灰褐色,天然状态下粒径稍大,呈块状,手摸有滑腻感,风化碎屑较少,粗颗粒表面较光滑,棱角分明,较坚硬,手捏不碎,干强度较高。(3)土样3 表观为灰黑色,天然状态下粒径较小,风化碎屑多,土样中网状裂隙发育,颗粒表面粗糙,较坚硬,手可捏碎(比较困难),干强度较高,含有少许植物根系夹杂物。上述3 种土样均表现为硬塑状态,遇水时强度急剧下降。本文根据现场土工试验,选取63 组数据进行分析,其中包括19 组已知土体自由膨胀率的数据以及44 组未知土体自由膨胀率的数据,土体分类均按照现有标准《铁路工程岩土分类标准》(TB 10077-2019)[20](表1、表2)。由于篇幅限制,表2仅列出部分数据。
图1 肯尼亚地质分布及铁路线路图Fig.1 Geological distribution and railway line map of Kenya
图2 蒙内铁路现场膨胀土土样Fig.2 Expansive soil samples of the Mombasa-Nairobi Railway
表1 19 组已知土体自由膨胀率的土工数据Table 1 Geotechnical data of 19 groups of known soil free expansion rates
表2 44 组未知土体自由膨胀率的土工数据(部分)Table 2 Geotechnical data of 44 groups of unknown soil free expansion rates (Partial data)
支持向量回归机SVR(Support Vector Regression)是支持向量在函数回归领域的应用,是利用支持向量机SVM(Support Vector Machine)解决回归拟合问题的进一步拓展。首先,通过对多个土体参数进行分析,选出与预测量(即土体的自由膨胀率)关系较为密切的指标;然后,将选取的多个指标做为SVR 模型的输入,将预测量作为SVR 模型的输出,建立了一个多变量输入、单变量输出的SVR 预测模型。具体建立过程如下:假设训练样本的数量为l,那么训练数据集合为{(xi,yi),i=1,2,···,l},其中xi是第i个训练样本的输入向量。表示第i个训练样本的N个训练参数,yi是第i个训练样本的土体自由膨胀率。根据从低维空间映射到高维特征空间的关系,建立在高维特征空间中的线性回归函数,如下:
式中:w—高维特征空间权重向量;
Φ(x)—非线性映射函数;
b—阈值。
为了表示回归效果的好坏,引入 ε线性不敏感损失函数
式中:f(x)—回归函数返回的土体自由膨胀率预测值;
y—对应的土体自由膨胀率真实值。
若f(x) 与y二者差值小于等于ε,则损失值等于0。
如何确定权重向量w、阈值b是非常重要的一步,这里引入松弛变量和,建立极小化目标函数,即
式中:C—惩罚因子,C越大表示对训练误差大于 ε的样本惩罚越大。
ε规定了回归函数的误差要求,ε越小表示回归函数的误差越小。
为了求解式,引入Largrange 函数,并将其转化为对偶形式:
式中:K(xi,xj)—K(xi,xj)=Φ(xi)·Φ(xj)为核函数,常用的核函数类型有线性(Linear)函数、多项式(Polynomial)函数、高斯径向基函数RBF(Radial Basis Function)和S形(Sigmoid)函数等,不同核函数对模型性能会有不同的影响。
通过式,可以求得 α 和 α∗。假设 α 和 α∗的最优解分别为便可进一步求得权重向量w、阈值b。于是,回归函数为
式中:f(x)—预测的土体自由膨胀率;
x—输入的相关土体参数。
《膨胀土地区建筑技术规范》(GB 50112-2013)[9]给出了膨胀土胀缩等级的分级规定(表3)。
根据表3,可以建立土的自由膨胀率与胀缩等级之间的关系函数F
表3 膨胀土的胀缩等级分类Table 3 Classification of expansion and contraction grades of expansive soil
使用SVR 模型解决回归拟合问题时,输入参数与输出参数的相关性将会影响SVR 模型的预测结果。为了选取与土体自由膨胀率 δef相关性较为密切的指标参数,将表1中的土体自由膨胀率 δef和其余参数分别进行皮尔逊相关性检验,使用皮尔逊相关系数R(R∈[−1,1]) 来描述各参数与土体自由膨胀率 δef之间的相关性强弱。R的绝对值越大,说明该参数与自由膨胀率的相关性越大;反之,说明该参数与自由膨胀率的相关性越小。相关性检验结果如表4所示。
从表4中可以看出,各参数与土体自由膨胀率均具有一定的相关性,但相关性大小各异。其中,仅有液限、塑限、塑性指数和颗粒组成≤0.075、(20,60]的R值大于零。并且在颗粒组成中,颗粒组成≤0.075、(0.075,0.25]、(0.25,0.5]的R的绝对值较大。特别的,以上三种颗粒组成的R的绝对值虽较为接近,但是颗粒组成(0.075,0.25]、(0.25,0.5]的R的绝对值却大于颗粒组成≤0.075 的R的绝对值。颗粒组成≤0.075 是影响自由膨胀率的主要因素,但是其在相关性分析中的表现却并不如颗粒组成(0.075,0.25]、(0.25,0.5]。为此,再对颗粒组成≤0.075、(0.075,0.25]、(0.25,0.5]进行相关性分析,得到颗粒组成≤0.075 与颗粒组成(0.075,0.25]、(0.25,0.5]之间的相关性分别为−0.950、−0.929,R的绝对值比较接近于−1,在0.01 的显著性水平上相关性极显著。基于上述分析并结合试验数据的情况,认为出现颗粒组成(0.075,0.25]、(0.25,0.5]的R的绝对值大于颗粒组成≤0.075 的R的绝对值的原因为:在本文的试验数据中,三个参数的相关性显著,即颗粒组成≤0.075 的增加和减小伴随着颗粒组成(0.075,0.25]、(0.25~0.5]的减小或增加,故在相关性分析中,颗粒组成≤0.075、(0.075,0.25]、(0.25,0.5]的R的绝对值较接近,且产生了颗粒组成≤0.075 的R的绝对值小于颗粒组成(0.075,0.25]、(0.25,0.5]的R的绝对值的特殊情况。针对该特殊情况,本文在方案设计中将进行特殊的处理。
表4 土的各参数与自由膨胀率之间的相关性Table 4 Correlation between various parameters of soil and free expansion rate
基于上述分析,并考虑到在实际工程中土的类别对自由膨胀率的影响也不容忽视,因此,本文最终确定了两种包含不同输入参数的SVR 模型分别为:①液限、塑限、塑性指数、3 种不同粒径的颗粒含量≤ 0.075、(0.075,0.25]、(0.25,0.5]、土的类型,共计7 个参数;②液限、塑限、塑性指数、粒径< 0.075 的颗粒含量、土的类型,共计5 个参数。
训练样本的选择对模型预测的可靠性影响极大。合理的训练样本数量及样本训练次数可以提高模型预测准确率和计算效率。
根据土体各参数描述中存在定量和定性的区别,可将其分为两类:(1)液限、塑限、塑性指数、3 种颗粒组成≤0.075、(0.075,0.25]、(0.25,0.5];(2)土的类型。第(1)类是定量描述,考虑到其值域的区别较大且量纲不同,故对其进行归一化处理,常用的归一化处理方法有两种,即min-max 标准化和Z-score 标准化方法。本文使用Z-score 标准化方法,采用数据的均值和标准差对数据进行标准化处理,处理后的数据符合标准正态分布,转化函数为
式中:µ—样本数据的均值;
σ—样本数据的标准差。
由于土的类型是定性描述的,无法直接运用计算机语言进行运算处理,本文考虑采用One-Hot 编码进行处理。One-Hot 编码主要是采用N 位状态寄存器来对N 个状态进行编码,每个状态都由他独立的寄存器位,并且在任意时候只有一位有效,即将分类参数变量用一个二进制向量表示。对于本文,土的分类为粉砂、黏土、粉质黏土,对应的二进制向量分别为 [0,0,1]、[0,1,0]和 [1,0,0]。采用这种方式不但可以将土的类型转化为计算机能够识别的信息,而且较好的保留了不同土的类型之间的关系和区别。
式中:Pi—表示第i次随抽取l个样本数据作为训练样本时的预测准确率,i=1,2,···,n;
n—表示随机抽取的次数,一般情况下,n值越大则的可信度越高。
从已知土体自由膨胀率的19 组数据中随机抽取10~18 组样本作为训练样本,剩余组作为验证样本。抽取不同样本数量时可能出现的样本组合数量如表5所示。
表5反应了不同抽取样本数量下的样本组合数,也是下文进行预测的重要基础。针对预测次数的概念,例如19 个数据中抽取10 个训练样本(每个土样的参数属于一个整体,训练样本是某土样的所有参数的一个整体)的情况,由于每次预测均采用电脑进行随机抽取,所以虽然都是属于19 个数据中抽取10 个训练样本的情况,但其实对于每一次的预测,训练样本及预测样本几乎均是不同的。当然,这只是对不同土样之间的随机抽取,并不是对每个土样下的各参数的随机组合(即土样A 的所有参数都是一直属于土样A 的,并不会与其他土样的参数进行样本组合)。
表5 不同抽取样本数量下的样本组合数Table 5 Number of sample combinations under different sample sizes
为了消除随机不同样本组合可能带来的影响,下面分析预测次数对预测结果产生的影响。分别采用Linear 函数、Polynomial 函数、Sigmoid 函数、RBF 函数作为模型的核函数,利用MATLAB 编程计算预测准确率随预测次数n的变化趋势。如图3所示,为随机抽取10 组数据时,随预测次数n的变化趋势图,由图可知,当预测次数n较小时,波动较大,随着预测次数的增加,逐渐趋于稳定。图4给出了随机抽取16 组数据时随 预测次数n的变化趋势图,由图4可知,其变化规律与图3相同,即随着n的增加,逐渐趋于稳定。此外,观察图3、图4并结合其余的随机抽取若干组数据时随预测次数n的变化趋势图可知,当预测次数达到1 000 次时,预测准确率已基本趋于稳定,故本文在分析过程中均采用n=2 000,可满足精度要求。
图3 随机抽取10 组训练样本时的预测准确率Fig.3 Prediction accuracy when randomly selecting 10 sets of training samples
图4 随机抽取16 组训练样本的预测准确率Fig.4 Prediction accuracy when randomly selecting 16 sets of training samples
图5 预测准确率随训练样本数量增加时的变化规律Fig.5 Variation of prediction accuracy as the number of training samples increases
综上,从上述分析可知,方案一(7 个输入参数)的预测效果好于方案二(5 个输入参数)的预测效果,并根据准确率变化趋势可以认为当采用19 组试验数据为训练样本时,RBF 函数和Linear 函数为核函数时模型的预测效果较好,其次是Sigmoid 核函数模型,效果最差的则是Polynomial 核函数模型。需要指出的是,由于本文训练样本数量有限,这可能会对预测准确率结果产生一定影响,但从目前的训练结果来看,当采用RBF、Linear 和Sigmoid 核函数时,预测准确率基本处于70%(方案一)和67%(方案二)以上,证实了本文所提方法的可行性。此外,该模型还具有以下2 个优点:(1)随着样本数据的进一步积累,预测准确率也会相应地提高;(2)本文模型并不局限于肯尼亚蒙内铁路工程,对于其他地区膨胀土的胀缩等级预测同样适用。
采用上述方法,将表1中19 组土工数据作为训练样本,对表2中自由膨胀率未知的44 组土样进行胀缩等级预测。考虑到上述4 种函数的预测效果,本节采用linear 函数、RBF 函数和Sigmoid 函数3 个函数作为模型的核函数进行预测(图6)。
图6 三种核函数模型预测的土体自由膨胀率Fig.6 Free expansion rate of soils predicted by three kernel function models
为了更直观的显示预测得到的44 组土样的胀缩等级,将44 组土样的预测等级用柱状图展现出来(图7)。对于方案一,从预测结果来看,采用Linear 函数时,胀缩等级为“强”、“中”、“弱”和“非”占比分别为29.5%、40.9%、27.3%和2.3%;当采用RBF 函数时,占比依次为29.5%、22.7%、43.2%和4.6%;采用Sigmoid 函数时,预测结果与采用Linear 函数时保持一致。通过统计可以发现,三种核函数模型预测得到的土样胀缩等级均相同的数量为32 组,占比为73%;其余12 组土样胀缩预测等级或相同或相邻,不存在“越级”现象。相比于RBF 函数,采用Linear 函数和Sigmoid 函数时,土样胀缩等级预测结果更为保守(即胀缩等级预测值较大)。对于方案二,采用Linear 函数时,占比分别为31.8%、22.7%、40.9%和4.5%;当采用RBF 函数时,占比分别为34.1%、27.2%、25%和13.6%;采用Sigmoid 函数时,占比分别为29.5%、25%、38.6%和6.8%。通过统计发现,三种核函数模型预测得到的土样胀缩等级均相同的数量为30 组,占比为68.2%;其余14 组土样胀缩预测等级或相同或相邻,不存在“越级”现象。
图7 三种核函数模型预测的土体胀缩等级Fig.7 Grades of soil swelling predicted by three kernel function models
对于胀缩预测等级不同的土样,从安全性方面考虑,可以取保守值,即取3 种函数中胀缩等级较大的预测值。表6给出了44 组土体胀缩等级预测结果,其中胀缩等级为“强”、“中”、“弱”和“非”的占比分别为29.5%、47.7%、20.5%和 2.3%(方案一)、36.4%、27.3%、31.8%和4.5%(方案二)。对于两种预测方案,可见对于胀缩等级为“中”的土样的预测结果差异较大,从而导致与胀缩等级为“中”临近的“强”、“弱”土样数量也有一定的差异。
为了验证本文方法的合理性,使用模糊层次分析法对44 组土样进行胀缩等级的评价。选择了液限、塑限、塑性指数、自由膨胀率、颗粒组成≤0.075 共计五个指标作为评价依据,其中自由膨胀率为三种核函数模型预测结果的平均值。选用熵权法和层次分析法进行结合,并选取了基于模糊分布原理下的x次抛物线形或半抛物线形分布函数,使用MATLAB 软件编程计算,得到相应的评价结果如表6所示。模糊评价结果与方案一预测结果的吻合度为50%,在结果不同的情况中基本是模糊评价的结果要低一级,不存在“越级”现象;模糊评价结果与方案二预测结果的吻合度为77%,在结果不同的情况中基本是模糊评价的结果要低一级,不存在“越级”现象。
表6 44 组土样的胀缩等级预测及评价结果Table 6 Prediction results and evaluation results of swelling and shrinking grades of 44 groups of soil samples
综上,本文虽然是通过训练多指标参数得到单一指标而确定膨胀土的胀缩等级,但是通过与多指标综合评价方法对比可以发现本文方法得到膨胀土胀缩等级相对于工程施工而言更加“安全”。
(1)基于肯尼亚“蒙内铁路”沿线膨胀土的土工试验数据,选择土体自由膨胀率作为预测指标,通过皮尔逊相关性分析及工程经验根据训练参数的不同设计了两种方案:液限、塑限、塑性指数、3 种不同颗粒组成≤0.075、(0.075,0.25]、(0.25,0.5]和土的类型;液限、塑限、塑性指数、颗粒组成≤0.075 和土的类型。采用Zscore 标准化方法对定量指标进行归一化处理以消除不同量纲的影响,采用One-Hot 编码将定性指标转化为亚变量使其能够被计算机识别和储存,从而建立了基于SVR 的膨胀土胀缩等级预测模型。
(2)提出以“预测准确率”作为模型预测效果的评价指标,采用4 种不同核函数的SVR 模型所得到的预测准确率虽存在明显差异,但根据该评价指标随着随机抽取训练样本次数的变化情况,当预测次数达到1 000 次时,不同核函数的SVR 模型得到的预测准确率均已基本趋于稳定,本文选择预测次数为2 000 次,可满足精度要求。
(3)不同核函数SVR 模型之间的预测准确率差异为核函数的选择提供了重要的参考依据,其中采用RBF 函数和Linear 函数时的预测效果较好,其次是Sigmoid 函数,最差的是Polynomial 函数。采用前三者函数时,方案一的预测准确率基本在70%以上,方案二的预测准确率基本在68%以上,证实了本文所提方法的可行性,且随着样本数据的进一步积累,预测准确率会相应地提高。
(4)通过已知自由膨胀率的19 组土工数据作为训练样本对自由膨胀率未知的44 组土样进行胀缩等级预测,结果表明采用Linear、Sigmoid 和RBF 函数作为核函数的SVR 模型时,方案一中土样预测的胀缩等级均相同的数量占比为73%,其余12 组土样胀缩预测等级相同或相邻,不存在“越级”现象,方案二中土样预测的胀缩等级均相同的数量占比为68.2%,也不存在“越级现象”。而两种预测方案的结果中,对胀缩等级为“中”的土样的预测结论差别较大。
(5)将本文预测结果与模糊层次分析法评价结果对比,进一步反映了本文预测方法结果的安全性。
(6)本文所提膨胀土的胀缩等级预测模型及方法,并不局限于肯尼亚蒙内铁路工程,对于其他地区的膨胀土胀缩等级预测同样适用。