基于钻进参数实时预测土体力学性质的Stacking集成模型

摘要

岩土体物理力学参数对工程勘察、设计、施工等作业不可或缺，但常规取样试验或原位检测均存在明显的精度误差。据此本文提出基于勘察钻探的实时钻进参数，建立基于机器学习的随钻土体物理力学参数模型。通过采集位于珠海市国家高新技术产业开发区内20 m勘探孔的真实数据，将EP-200G型钻机实时随钻采集的钻压、扭矩和三轴振动作为输入数据，将全孔土体粘聚力、内摩擦角、含水量与弹性模量试验数据作为输出。基于建模数据分析，证明使用单算法的3类机器学习模型（支持向量机、神经网络和决策树）的预测精度最高仅为0.78，而基于Stacking理念的集成模型可将预测精度提升至最高0.98。结合该模型，进行了随钻参数与土体参数间的敏感性分析，证实当不同土体参数发生变化时，不同随钻参数会发生明显变化，证明了随钻参数预测土体参数的可靠性与适用性。

关键词

土体参数; 钻进参数; 实时预测模型; 敏感性分析; 机器学习; Stacking理念; 工程勘察

0 引言

获取准确的岩土体物理力学参数是工程勘察的重要目标之一，也是后期设计、施工过程的重要保障。当前获取岩土体物理力学参数的常规方法是通过钻探取样进行室内土工试验以及现场原位试验两种方式，但这两种方式均存在较为明显的不足。对于钻探取样土工试验，取样和送样的过程对样品会产生扰动，以及样品脱离原始地应力环境会对试验值带来误差；对于现场原位检测，不仅需要满足严格的试验条件，并且大多原位检测的实际检测参数并非拟测量土体参数，检测参数与实际土体参数间的转换只能依靠经验公式。在未建立经验公式的新勘察区域，原位检测获得的岩土体物理力学参数仍将存在较大的误差。

由于岩土体物理力学参数的随机性与不确定性，学界对岩土体参数的取值研究更多的侧重于基于统计建立岩土体参数的概率分布模型。如徐卫亚等^［1］引入模糊不确定性，建立了对模糊特征统计量的计算公式。基于随机场理论，张继周等^［2］提出了图形三σ法的统计工具，吴长富等^［3］结合土体相关距离对岩土参数标准值进行了计算。陈炜韬等^［4］利用贝叶斯更新的方法对第四系残坡积粘性土隧道围岩的抗剪强度进行了优化。苏永华等^［5］通过多项式逼近实现了对岩土体随机参数概率密度函数的逼近。李小勇等^［6］结合实验数据的可靠性检验、概率模型有限比较法拟合优度和分布参数的推广贝叶斯法，在统计意义上对岩土参数的概率分布进行了优化分析。

随着机器学习算法在非线性拟合方向的高精度表现，神经网络、支持向量机、随机森林等机器学习算法也越来越多的被引入到对岩土体参数的预测分析中^［7-8］。如Najjar等^［9］通过神经网络建立了粘土渗透性的预测模型；Nath等^［10］实现了基于神经网络的地震波反演；Kaunda^［11］建立了岩石单轴抗压强度和各应力状态间的神经网络模型。阮永芬等^［12-13］构造了以小波分析理论为核函数的支持向量机模型，实现了引入粒子群、灰狼等优化算法的岩土力学参数反演模型。李文斌等^［14］基于随机森林建立了预测结构面轮廓的JRC模型。Zhang等^［15］融合XGBoost和随机森林算法分析软粘土抗剪强度与其他参数间的关系。

基于上述研究现状，可见当前学界对岩土体参数的测定侧重于获取大量数据后的统计分析以及建立机器学习模型分析不同参数间的关系上，但获取岩土体参数仍停留在传统的取样实验室分析以及原位检测上。本文拟通过在工程勘察钻机上实时监测钻进参数，借助集成算法实现对原位岩土体参数的准确预判。同时，基于准确的预测算法，本文还实现了地下参数变化时对随钻参数的敏感性分析。

1 建模数据获取

1.1 土体物理力学参数

本次试验区域位于广东省珠海市国家高新技术产业开发区，场地面积约为70 m²。试验使用图1（a）所示的EP-200G型钻机进行钻进，并在钻机上安装传感器获取包含钻压、扭矩及振动3类钻进参数。试验钻进取心共钻进3个孔，一字排开，每个孔之间相距1~2 m（如图1b所示）。每个孔钻进深度约为20 m，并分别在3个孔的不同深度处进行取心，共取得11组岩心（如图1c所示）。打捞岩心后，将岩心封装保存运往实验室做土工试验，得到土体物理力学参数，包括粘聚力c、内摩擦角$\phi$和弹性模量E。

图1  现场试验示意

对采集样品试验测定的结果如图2所示，在试验深度范围内自上而下可见粘性土层和砂土层两层土体，换层深度13 m左右。在粘性土层中，粘聚力c（图2a）在5~9 m深度范围内一直减小，从20.6 kPa变为15.3 kPa，但从9 m开始短暂增大，在11 m处增大至20 kPa，之后在11~13 m深度范围内骤减，在13 m处减至最小值12.3 kPa，与5 m处粘聚力相比减小了40.3%。内摩擦角$\phi$（图2b）趋势与粘聚力较为一致，整体下降，5 m处为22.1°，在13 m处减小为18.8°，减小幅度14.9%；同样在11 m处存在峰值27.5°。弹性模量（图2c）则整体稳定在70 MPa左右。

图2  土体物理力学性质测试结果

在13~21 m的砂土层中，粘聚力在17 m前呈现先缓慢增长至15.4 kPa，后急剧增加至28 kPa（19 m左右），最后波动下降至18.4 kPa。砂土层内摩擦角呈现典型先增后减的单峰值趋势，峰值出现在17 m处的28.9°。就弹性模量而言，砂土层的弹性模量与内摩擦角类似，呈现明显的单峰值趋势，在17 m达到最高的96.7 MPa。

1.2 随钻工艺参数的采集

试验采集的随钻工艺参数包括钻压、扭矩和振动。其中，基于EP-200G型液压钻机，对油缸和动力头的压力进行监测和计算可直接获得钻压和扭矩信号。由图3可见，钻压和扭矩值受土体强度和取心工艺共同影响，其中取心段的钻压和扭矩明显高于非取心段，当钻进深度进入强度更大的砂土层后，同类钻进工艺（取心或不取心）的钻压和扭矩也较上部粘土层钻进略微升高。

图3  钻压和扭矩实测数据

具体而言，在未取心情况下，钻进过程中粘性土层所产生的钻压与砂土层所产生的钻压差不多，前者与后者钻压均值的比值为108%。而在粘性土层中非取心段与取心段所产生的钻压相差极大，前者与后者的均值比值为60.8%，砂土层中非取心段与取心段所产生的钻压均值比值为18.2%。综上可知，土体强度变大时，钻压不会产生太大的变化；钻进取心时，钻进过程所产生的钻压急剧增大，其中砂土层取心时带来的钻压增加量是粘性土层取心时带来的钻压增加量的5.3倍。类似的，对于钻进扭矩，粘性土层（5.24~13.24 m）中土体强度小于砂土层（13.24~21.24 m）中土体强度，因此在未取心情况下，钻进过程中粘性土层所产生的钻进扭矩小于砂土层所产生的钻进扭矩，前者与后者钻进扭矩均值的比值为36.9%；而在粘性土层中非取心段与取心段所产生的钻进扭矩相差极大，前者与后者的均值比值为17.1%，砂土层中非取心段与取心段所产生的钻进扭矩均值比值为36.1%。综上可知，土体强度变大与取心时，钻进过程所产生的钻进扭矩均会增大，但取心时带来的钻进扭矩增加量约为土体强度增大带来钻进扭矩增加量的10倍，其中砂土层取心时带来的钻进扭矩增加量是粘性土层取心时带来的钻进扭矩增加量的1.6倍。

振动信号可为钻进过程中发生取心、岩土体强度变化（是否分层）提供准确的判据，故本次试验将三轴振动传感器固定于动力头上进行振动信号的采集（图4）。考虑动力头上下移动、回转自身也会产生振动，故试验时先后采集了钻进作业的振动信号与动力头空转并上下移动的振动信号。将钻进作业的振动信号减去空转信号后获得的3个试验孔的采集振动信号如图5所示。

图4  振动传感器安装位置

图5  振动监测数据（用红色边框表示取心段）

结合图5，取心段X、Y、Z轴振动信号的振动幅度远大于非取心段；粘性土层中非取心段与取心段所产生的三轴振动信号相差显著，对于X轴前者与后者的均值比值为25.6%，Y轴为34.5%，Z轴为44.3%；砂土层中非取心段与取心段所产生的三轴振动信号同样相差极大，对于X轴前者与后者均值比值为12.5%，Y轴为27.2%，Z轴为38.6%。相对于取心造成的振动影响，不同深度的土体物理力学性质对振动产生的影响可忽略不计。

2 随钻土体物理力学参数模型

2.1 单算法基本模型

本次试验获得的钻进参数、土体物理力学性质随深度的变化趋势较为复杂，均为非线性变化；且由于钻进过程中取心与非取心共存，钻进监测参数的数量和范围均差异明显。综上，在建模预测算法的选择上，应该更加倾向于适合处理非线性问题、能够处理较大的输入数据的算法，据此选择支持向量机、神经网络和决策树作为建模单算法，拟建立的模型结构如图6所示。

图6  单算法基本模型结构示意

针对拟使用的支持向量机、神经网络以及决策树3种算法，每一种算法将分别对3项土体参数进行建模预测。使用如式1所示的决定系数R²（用于评估模型的解释能力）与均方误差根RMSE（用于评估模型的准确性）作为训练模型性能的评价指标。各算法预测结果如图7所示，精度对比如表1所示。3类机器学习算法对3项土体参数（粘聚力c、内摩擦角φ以及弹性模量E）都有着一定的预测精度。3类算法纵向对比看，决策树对3项土体物理力学参数的评价指标R²最高（均大于0.7），RMSE最低，表征预测精度最高。4项物理力学参数横向对比看，3类算法均对内摩擦角和弹性模量具有较高的预测精度，对含水量的预测精度最低。

表1  单算法预测精度对比

算法	粘聚力c		内摩擦角φ		弹性模量E
	R2	RMSE	R2	RMSE	R2	RMSE
支持向量机	0.41	0.3207	0.44	0.3827	0.42	0.3813
神经网络	0.65	0.591	0.68	0.547	0.76	0.505
决策树	0.73	0.2178	0.76	0.25	0.78	0.2346

图7  基于单算法的参数预测结果

$$\left\{\begin{array}{l}{R}^{\mathrm{2}}=\mathrm{1}-{\displaystyle \sum _{i=\mathrm{1}}^n}{\left(y_i-\stackrel{\wedge }{y_i}\right)}^{\mathrm{2}}/{\displaystyle \sum _{i=\mathrm{1}}^n}{\left(y_i-\overline{y_i}\right)}^{\mathrm{2}}\\ RMSE=\sqrt[]{{\displaystyle \sum _{i=\mathrm{1}}^n}{\left(\stackrel{\wedge }{y_i}-y_i\right)}^{\mathrm{2}}/n}\end{array}\right.$$

(1)

2.2 集成模型架构

由表1可知，单一机器学习模型对3项土体参数的评价指标均不高，最高的决策树模型针对弹性模量的预测精度也仅为0.78，无法满足实际工程应用的需求。为提升预测精度，建立了基于Stacking思想的双层分析集成体系融合模型（如图8所示）。其核心原理在于首先构建一级学习器（基学习器）并训练获得一级预测结果，其次将一级预测结果作为输入参数引入二级学习器（元学习器），训练获得二级预测结果，从而进一步提升初级基学习器的极限精度。结合前文单一算法结果，本次训练中基学习器采用决策树模型，元学习器仍从支持向量机、神经网络和决策树3种经典单模型中进行选择建立。

图8  基于Stacking的集成算法模型架构

依次将支持向量机、神经网络和决策树算法带入如图8所示集成模型中的元学习器，对集成模型进行训练，得到的训练精度如图9和表2所示。由图9和表2可见，对于元学习器的选择而言，仍然是决策树能够带来最高的预测精度，以其为元学习器时，对3项土体物理力学参数的评价指标R²均超过0.95，最高位内摩擦角达0.98。

表2  基于Stacking的集成算法精度

算法	粘聚力c		内摩擦角φ		弹性模量E
	R2	RMSE	R2	RMSE	R2	RMSE
支持向量机	0.91	0.1215	0.9	0.1638	0.92	0.1425
神经网络	0.94	0.0927	0.95	0.104	0.95	0.106
决策树	0.97	0.0725	0.98	0.0762	0.96	0.0997

图9  基于决策树的集成模型预测效果

2.3 不同模型对比

横向对比单算法模型与集成模型如图10所示，可见基于Stacking的集成算法模型相对于传统的单算法模型预测精度显著提高。其中就决定系数R²而言，相对于单预测模型最高不足0.8的预测精度，集成模型3种不同元学习器的效果均超过了0.9，精度提升幅度最小为23%（使用决策树对弹性模量E的单预测模型精度为0.78，使用集成模型后提升至0.96），最大为308%（使用神经网络对含水率的单模型预测精度仅为0.24，使用集成模型后提升至0.98）。

图10  集成学习模型与单算法模型精度对比

将本试验建立模型与同类模型进行对比。张洁等^［16］基于500多组花岗岩残积土的试验数据采用经验法与主成分分析法建立了联合概率密度函数模型，对粘聚力、内摩擦角的预测性能指标R²分别为0.27、0.52；杨鹏等^［17］基于100多组土工试验数据，建立了预测粘聚力、内摩擦角的多元回归预测模型，决定系数R²分别为0.66、0.662；李澄清等^［18］建立了BP人工神经网络反演系统，对粘聚力、内摩擦角的预测精度分别为0.89、0.91；陈庆贺^［19］通过引入人工神经网络对粘聚力、内摩擦角进行了预测，预测指标R²分别为0.92、0.95。而本文建立的Stacking融合模型，对土体的4项参数预测指标R²均在0.96及以上，各类模型的平均预测精度如图11所示，本文的联合模型较之前的预测模型在预测种类以及精度上都有了一定提升，可在实际工程中推广应用。

图11  本文Stacking融合模型对比其他模型精度

3 钻进工艺参数对土体物理力学参数的敏感性

3.1 敏感性计算方案

定义本文建立的Stacking融合模型如式2所示，可通过调整输入参数（钻压、扭矩、振动）观察输出参数（粘聚力、内摩擦角和弹性模量）的变化，进而分析实时勘探作业过程中土体性能变化时对实时钻探参数的潜在影响。

$$(c,\phi ,E)=F_s(WOB,T,V_X,V_Y,V_Z)$$

(2)

对拟引入灵敏度的参数进行分析，令$\mathrm{\Delta }x$为钻进参数的变化量，$\mathrm{\Delta }y$为土体物理力学参数的变化量，则钻进参数对土体物理力学参数变化的灵敏度$N$为：

$$N=\frac{\mathrm{\Delta }x}{\mathrm{\Delta }y}$$

(3)

为消除不同参数间的数量级与单位，在进行敏感性分析前对计算模型的输入和输出参数均进行了归一化处理如式（4）所示：

$$x\mathrm{\text{'}}=\frac{x-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}$$

(4)

式中：$x\mathrm{\text{'}}\mathrm{、}x$——归一化后与归一化前的样本值；$x_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{、}{x}_{\mathrm{m}\mathrm{a}\mathrm{x}}$——归一化前样本集最小值、最大值。

将合理范围内的钻进参数依次取4个不同值，可分别获得指定钻进参数与指定土体物理力学参数之间的关系图，结合式（3）所示的灵敏度计算方程，可如图12所示依次计算每一段钻进参数变化时对应预测土体物理力学参数变化的灵敏度。最终灵敏度取3段灵敏度的加权平均值如式（5）所示：

$$\overline{N}=w_{\mathrm{1}}{N}_{\mathrm{1}}+w_{\mathrm{2}}{N}_{\mathrm{2}}+w_{\mathrm{3}}{N}_{\mathrm{3}}$$

(5)

其中权重系数$w_i$取各段钻进参数与原位参数（横轴）的变化量$x_i$在总长度$x$中的占比：

$$w_i=\frac{x_i}{x}$$

(6)

图12  土体参数与钻进参数灵敏度计算示意

3.2 敏感性分析结论

基于式（2）的模型与式（5）的计算方法，依次计算本次试验中各钻进参数与土体物理力学参数间的敏感性如表3和图13所示，可见不同钻进参数与不同土体物理力学参数间存在明显不同的灵敏度。从灵敏度的大小上看，明显钻进参数对弹性模量的灵敏度更高，即弹性模量的变化会明显诱发本次试验所有钻进参数发生变化。相对而言，粘聚力的变化只会引起钻压和扭矩的变化，而内摩擦角的变化则更多的只会引起振动发生变化。

图13  土体参数与钻进参数敏感性矩阵

表3  土体力学参数敏感性计算

土体参数	钻进参数	组1		组2		组3		平均灵敏度
		Δy	Δx	Δy	Δx	Δy	Δx
粘聚力c	钻压	0.304	0.33	0.696	0.33	0.957	0.33	0.65
	扭矩	1	0.33	0.33	0.33	—	0.33	0.44
	振动	0.86	0.41	1	0.222	—	0.222	0.25
内摩擦角φ	钻压	1	0.33	—	0.33	—	0.33	0.11
	扭矩	1	0.33	—	0.33	—	0.33	0.11
	振动	0.25	0.41	0.75	0.222	—	0.222	0.74
弹性模量E	钻压	0.27	0.33	0.53	0.33	0.2	0.33	1.15
	扭矩	1	0.33	0.25	0.33	0.2	0.33	1.09
	振动	0.51	0.41	0.36	0.222	0.13	0.222	0.85

注：  —表示Δy无变化。

4 结论

通过收集并测定珠海市国家高新技术产业开发区试验区域的土体物理力学性质与钻进取样过程中的随钻参数，本文建立了随钻参数预测土体物理力学参数的集成预测模型。通过本文的研究，主要成果与结论如下：

（1）使用随钻参数建立预测土体物理力学参数的智能模型，分别建立了单算法直接预测模型与基于Stacking融合理念的集成模型，通过对比预测准确性，集成算法准确性显著高于单算法模型，且通过不同学习器的精度对比，可知当集成模型的基学习器和元学习器均使用决策树算法时预测精度最高，针对粘聚力、内摩擦角和弹性模量，可至少达到0.96的预测精度。

（2）基于集成模型，推导获得了钻进参数与土体物理力学参数之间的敏感性计算方案，并对试验区域的敏感性进行了计算对比，结果显示不同钻进参数与不同土体物理力学参数间存在明显不同的灵敏度，进而证实可通过观测随钻参数判断地层参数变化的可能，同时也证实使用随钻参数准确预测土体物理力学参数的可能。