1.1　钻进数据预处理

在工程实际的钻进工作中，受钻进设备、采集设备、钻进环境、人员操作等客观因素的影响，所采集的钻进数据不可避免地混杂一些漂移、残缺、失真及处于特殊状态的异常数据。为提高钻进数据的质量以更准确地反映围岩破碎情况，首先需要对获取的钻进数据进行预处理，具体步骤包括异常数据处理、数据降噪、数据插值和段落划分4项。

在钻进初始阶段，由于受到钻机空转、开孔及掌子面喷浆等因素的影响，所采集的钻进数据并不具备解译价值，因此，应将钻进初始阶段的异常钻进数据直接剔除。在正常钻进过程中，由于受换杆操作的影响，钻压、转速、钻速、扭矩和钻进深度在一定时间内均会出现异常变化，同样将这些数据剔除。在剔除上述异常钻进数据后，所对应的钻进时间因缺失而出现不连续的现象(图1)，为便于后续的分析，应对此类时间数据进行衔接处理。此外，受钻机设备稳定性及钻进参数采集设备精度的影响，钻进数据会产生不规则性和粗糙性。小波分析可以从多尺度、分辨率方面对信号进行分析从而发掘其中的隐含信息，有效抑制随机噪声。通过选取一个合适的阈值，可以使信号趋向于真实值^[20]。因此，本文选择小波阈值降噪方法对钻进数据进行降噪处理。小波阈值降噪的基本原理是信号经过小波变换得到的小波系数含有重要时频信息，真实信号的小波系数较大，噪声的小波系数较小，选取一个合适的阈值即可保留信号并滤除噪声。

图1全屏查看

异常值剔除后的时间不连续段落示例

钻进参数采集设备的采样频率并不固定，为便于后续对各钻进参数的二级指标波动信号特征进行时域、频域分析及特征提取，采用插值方法将原始钻进信号的采样频率进行固定处理。本文采用简单直观的线性插值方法将采样频率统一预处理为1 Hz。在隧道及地下工程中，围岩完整性均是按段落进行评价，因此，在对钻进数据完成降噪处理后，需根据实际需求进一步对钻进数据进行等距分割，以便分段研究钻进数据特征及其对应的围岩特征。本文中，考虑到所依托工程隧道单循环开挖长度为1 m左右，因此，将各钻进数据划分间距定为1 m。

1.2　钻进数据特征指标解译与优选

1.2.1　钻进数据特征定量解译

钻进参数信号由多个具有不同调幅-调频特征的信号分量组成。如何从预处理完成后的钻进参数信号中解译出与围岩完整性密切相关的动态变化特征，是实现围岩完整性随钻评价的关键。预处理后的钻进信号可作为一种时域信号，通过时频分析方法又可将时域信号转换到频域进行观察，本文拟从时域和频域这2个角度对钻进参数波动特征进行全面解译和提取。

1) 钻进信号时域特征提取。

时域特征反映了信号在时间域上的变化规律，是动态信号最基本、最直观的表达。时域信号中主要提取的时域特征可分为2种：一种是有量纲参数，常用于描述信号的强度特征，本文选择了最大值、最小值、峰峰值、平均值、均方根值、标准差、信号总值7个指标；另一种是量纲一参数，常用于表征信号的波形和变化特征，本文选择峭度、偏度、峰值因子、脉冲因子、波形指标、裕度因子6个指标。以下是相关指标的计算方法，其中，xi为信号值(i=1，2，3，…，n)，N为信号样本总数。

最大值x_max表示同一段钻进信号的最大瞬时幅值，

(1)

最小值x_min表示同一段钻进信号的最大瞬时幅值，

(2)

峰峰值x_peak表征信号整体振幅，

(3)

平均值表示段内各采样点的平均幅值，

(4)

均方根值x_rms反映信号能量平均水平，

(5)

标准差x_std用于刻画信号波动强度，

(6)

信号总值x_E用于衡量总能量大小，

(7)

峭度K用于描述波形尖锐程度、异常尖峰敏感性，

(8)

偏度S用于描述波形对称性，

(9)

峰值因子P用于衡量最大峰值相对能量大小，

(10)

脉冲因子I用于衡量单点尖峰程度，

(11)

波形指标W用于衡量波形平坦度，

(12)

裕度因子M用于综合反映信号裕量：

(13)

2) 钻进信号频域特征提取。

频域分析可从频域角度了解信号幅值在频率范围内的分布情况，从而进一步反映信号的复杂变化。本文拟通过经典的时频分析方法——快速傅里叶变换(fast fourier transform, FFT)，将时域钻进信号转换为频域，获取相应的频谱图，再对信号进行频域特征提取。本文采用频域幅值平均值、重心频率、均方频率、频率方差对信号进行频域特征描述^[21-22]，各频域特征的计算方法如下。

频域幅值平均值A_m表示平均谱幅，

(14)

重心频率F_c反映能量中心，

(15)

均方频率M_sf用于刻画频谱能量分布，

(16)

频率方差V_f用于度量频域离散性；

(17)

式中：fi为频率；pi为相应的幅值(钻速或扭矩等)。

1.2.2　钻进数据特征指标优选

合理选择特征指标作为围岩完整性评价模型的输入参数是机器学习的关键步骤。为了在不降低围岩完整性评价质量的前提下筛选出反映围岩完整性的关键因素，本文采用RReliefF算法并结合Pearson相关系数来完成指标特征优选。RReliefF算法是一种经典的过滤式特征选择算法，专门用于解决回归模型的特征选择问题^[23]，该算法通过引入2个相异实例的概率建立模型：

(18) (19) (20)

式中：为实例在特征A上的差异值；为实例的预测差异值；为相近实例；P_diffA为在相近实例中在特征A上的差异概率；P_diffC为在相近实例中的预测值差异概率；P_diffC|_diffA为在特征A上有差异的相近实例中的预测值差异概率。根据贝叶斯公式，特征A的权值W(A)表示为：

(21)

结合上述公式可以估算出各输入特征的权值，算法的实施流程和伪代码如图2所示。

图2全屏查看

RReliefF算法伪代码

为减少各特征指标之间因较大相关性而对数据集造成的重复干扰，确定最优特征指标组合，本文进一步引入Pearson相关系数对选出的特征指标进行二次筛选。2个特征指标之间的Pearson相关系数计算式如下。

(22)

式中：r为Pearson相关系数，取值范围为[-1，1]；为x的均值；为y的均值；n为样本总数。

当r的绝对值较大时，2个变量间的相关性较强。|r|大于0.8时表明参数之间存在强相关性，小于0.2则表明参数之间存在弱相关性^[24]。在实际检验过程中，选择将|r|>0.8的特征指标中的一组进行剔除。

RReliefF-Pearson特征指标优选方法的流程见图3。具体步骤如下：首先基于所提取的钻进曲线时域特征和频域特征构建初始特征指标数据集，并采用Min-Max标准化对数据集进行预处理；使用RReliefF算法寻找各特征指标与围岩完整性系数之间的相似度和差异度，以权重表示特征指标对围岩完整性系数的重要程度并进行排序；通过对各特征指标的权重进行对比分析确定权重阈值，进而剔除小于权重阈值的特征指标。为去除RReliefF算法下的冗余特征指标，引入Pearson相关系数对选出的特征指标进行二次筛选求解特征指标之间的相关性，并剔除相关性较高的特征指标，最后得到最优特征指标，完成特征指标的降维与选择。

图3全屏查看

RReliefF-Pearson特征指标优选流程

2.1　主要算法介绍

随机森林算法(random forest，RF)是一种以决策树为基学习器并结合Bagging集成思想的非线性统计集成算法，广泛用于分类和回归分析^[25]。随机森林通过构建多棵决策树对样本进行训练，然后，通过集成决策树的结果来提高模型的准确性和稳定性。算法原理如图4所示。随机森林有2个重要的随机性来源：一是数据的随机性，二是特征的随机性。对于数据的随机性，随机森林算法使用自助采样法从原始数据集中随机选择n个样本，组成一个新的训练数据子集，这个新的数据子集被用来训练一个新的决策树。对于特征的随机性，随机森林算法在每个决策树的节点上随机选择m个特征，并从这m个特征中选择最佳的特征用于分裂。

图4全屏查看

随机森林算法原理示意图

2.2　模型训练与构建

在机器学习模型中，超参数的设定对模型的运算效率和泛化性能有重大影响。超参数优化就是寻找这些超参数的最优设定，将模型在验证集上的泛化性能最大化。贝叶斯优化算法(bayesian optimization，BO)是一种基于高斯过程和贝叶斯定理的全局优化算法，该算法利用已有的样本点构造代理函数近似替代黑箱目标函数，降低了计算目标函数的时间成本，并引入了采集函数，对损失函数最小的评估点进行有效探索和利用，能在较少的评估次数下取得复杂目标函数的最优解^[26-27]。为寻求预测性能最优的模型，采用BO算法对RF的超参数进行优化，通过设置参数空间对指定的超参数进行搜索，以获取最优超参数组合，RF拟优化的超参数组合及寻优范围见表1。超参数寻优迭代次数定为200次，并选取均方根误差E_RMSE作为模型每次迭代后的评价分数。在寻优过程中，通过迭代尝试不同的超参数组合，记录并保存超参数组合和模型评价分数。在五折交叉验证下经过200次迭代，预测模型总共被建立和评价5×200=1 000次，最终通过选取模型最佳分数来确定最优的超参数组合。

本文采用Jupyter Notebook作为Python语言开发环境，依托Scikit-learn机器学习库作为预测模型建立的基础工具，实现数据加载切分、数据预处理、模型建立及训练结果保存等步骤。模型训练与构建过程如下。

1) 数据加载及切分。指定数据库文件所在路径，加载数据，对加载后的数据进行切分处理，完成模型输入特征和输出标签的分离与确定。

2) 数据预处理。对数据集中的异常值进行处理，检查是否含空值或无穷大值。为消除数据集中各输入特征之间量纲的限制、提高模型迭代收敛的速度等，对数据集进行Z-Score标准化处理，将样本数据集均值转换为0、标准差为1，使得数据分布符合标准正态分布。

3) 数据集划分。为确保训练集和测试集的数据之间具有相互独立性，对数据集进行随机划分。常见的训练集和测试集划分比例为7꞉3和8꞉2^[28-30]，经实际验证对比后，本文所选择的训练集和测试集的划分比例为8꞉2。

4) 模型构建及预测。通过调用Scikit-learn库中的Random Forest Regressor函数，设置模型初始超参数，构建初始模型，采用5折交叉验证对随机森林模型进行训练。基于表1执行贝叶斯算法的超参数优化，并使用最优超参数组合进行后期的训练与预测。在模型训练完成后，再分别对训练集和测试集进行预测。

2.3　模型评估

为全面评估预测模型的预测性能，选择均方根误差E_RMSE、决定系数R²和平均绝对百分比误差E_MAPE这3种指标，相应计算公式见表2。其中，R²表征回归模型对观测数据拟合的准确程度，其值越接近1表明模型的拟合效果越好；E_RMSE表征回归模型预测值与真实值的偏离程度；E_MAPE为预测值和真实值之间绝对误差的平均值。

表2全屏查看

预测模型的评价指标

评价指标	计算公式	评判标准
R2		数值越大，模型的预测性能越好
ERMSE		数值越小，模型的预测性能越好
EMAPE		数值越小，模型的预测性能越好

展开更多

3.1　工程概况及现场试验

为验证所提围岩完整性评价方法的有效性，依托我国正在建设的某铁路隧道的辅助坑道2号斜井工程，现场采集了钻进数据及其所对应的围岩完整性数据，开展了围岩完整性评价案例分析与验证。该工程位于我国高原地区，海拔高程在 4 100 m以上，属构造剥蚀高山-极高山地貌，最大埋深超过800 m，为高原深埋隧道。隧址区主要发育有燕山期正长花岗岩、二长花岗岩，多呈弱风化到强风化，节理裂隙较发育，工程地质条件从一般到较差。受构造运动及后期风化作用的影响，该段围岩岩体完整性差，围岩类别主要为IV类和Ⅴ类。钻探数据收集于2号斜井中的5次超前地质钻探工作，每次超前探孔施作长度为20~30 m不等，钻进参数采集系统的采样频率约为90条/m，共收集原始钻进数据10 000余条，累计覆盖隧道长度约120 m。在超前钻探作业完成后，对探孔进行了孔内声波测试，以便定量化获取隧道围岩的完整情况。受时间和成本的约束，现场每次超前钻探作业仅施作1孔。因此，使用非金属声波检测仪的一发双收探头对超前探孔进行声波测试(采样间隔为0.2 m)，进而获取隧道围岩的纵波波速。基于后续隧道开挖过程中的岩块采样及钻孔取芯，用于开展岩石声波测试和饱和单轴抗压强度测试，隧道围岩的完整性系数K_v计算式^[31]如式(23)所示。

(23)

式中：V_pm为岩体纵波波速，km/s；V_pr为岩石纵波波速，km/s。

最终120 m超前钻孔采集岩体纵波波速约600组，经过现场取样及后续的钻孔取芯，获取完整岩石的纵波波速112组，故而计算出围岩K_v共112组，其统计结果见图5。从图5可以看出：现场隧道围岩的完整性系数在0.10到0.45之间，参照表3中K_v与围岩完整程度的对应关系^[31]，所研究隧道段落的围岩完整性主要分为极破碎、破碎和较破碎共3种，占比分别为5.88%、69.75%及24.37%，以较破碎围岩为主。

图5全屏查看

围岩完整性系数K_v统计图

3.2　钻进数据定量解译结果分析

基于前文所提出的钻进数据预处理及指标提取方法，以2号斜井X2DK2+261.1~X2DK2+235.2里程段(长度为25.8 m)的钻压和钻速数据为例，将原始钻进数据与异常值剔除、降噪和插值预处理后的钻进数据对比如图6和图7所示。

图6全屏查看

钻进数据曲线随深度变化对比

图7全屏查看

钻进数据曲线随时间变化对比

从图6可以看出进行预处理后的数据在一定程度上滤除了因钻具振动等因素所产生的噪声，因而钻进参数能够更加准确地反映出围岩的真实情况。从图7可以看出，进行预处理后的钻进数据有效剔除了因换杆停机等因素导致的非正常时间记录，真实地反映了各钻进参数在时间上的变化特征，为后续的钻进数据时域分析奠定基础。此外，图6和图7显示钻压与转速受操作和地质变化影响显著。因此，在钻进数据进行段落划分处理后，对各段落的钻压、转速、钻速和扭矩的波动特征进行时频域分析和特征指标提取，将其作为围岩完整性评价模型构建的初选特征。以2号斜井X2DK2+261.1~X2DK2+235.2里程的26段钻速为例，其时频特征指标见表4和表5。

表4全屏查看

2号斜井X2DK2+261.1~+235.2里程处钻速的时域特征指标

段落	xmax	xmin	xpeck		xrms	xstd	K	S	xE	W	P	I	M
1	84.30	16.79	41.57	67.51	157.40	12.55	6.49	1.82	43.40	1.04	1.56	1.62	1.66
2	84.86	34.84	45.63	50.03	102.73	10.14	4.79	1.32	46.73	1.02	1.07	1.10	1.11
3	75.74	13.59	35.98	62.15	84.47	9.19	6.41	1.40	37.12	1.03	1.67	1.73	1.75
4	80.44	21.37	32.91	59.07	123.96	11.13	8.79	2.35	34.73	1.06	1.70	1.79	1.83
5	65.98	10.38	29.61	55.59	41.59	6.45	14.89	2.61	30.31	1.02	1.83	1.88	1.90
6	85.75	22.36	34.89	63.39	166.03	12.89	6.23	1.84	37.17	1.07	1.71	1.82	1.87
7	73.64	18.67	36.12	54.97	67.54	8.22	6.58	1.37	37.03	1.03	1.48	1.52	1.54
8	123.66	21.33	54.03	102.34	609.67	24.69	3.97	1.40	59.33	1.10	1.72	1.89	1.98
9	112.26	26.03	51.93	86.24	354.32	18.82	3.13	0.77	55.19	1.06	1.56	1.66	1.71
10	99.13	30.88	49.94	68.26	254.79	15.96	3.79	1.23	52.39	1.05	1.30	1.37	1.40
11	128.78	28.60	60.35	100.18	527.79	22.97	3.55	0.97	64.51	1.07	1.55	1.66	1.72
12	152.90	1.65	23.11	151.26	132.49	11.51	100.59	9.19	25.80	1.12	5.86	6.54	6.73
13	69.51	19.29	28.28	50.22	112.16	10.59	6.25	1.98	30.18	1.07	1.66	1.78	1.83
14	89.97	17.69	48.80	72.27	408.89	20.22	1.58	0.32	52.81	1.08	1.37	1.48	1.55
15	180.37	20.59	61.28	159.78	1 337.70	36.57	4.06	1.15	71.22	1.16	2.24	2.61	2.83
16	56.45	7.218	20.47	49.23	71.22	8.44	12.86	3.22	22.13	1.08	2.22	2.40	2.47
17	85.19	19.82	37.02	65.36	206.44	14.37	3.95	1.13	39.69	1.07	1.65	1.77	1.83
18	118.48	10.15	26.94	108.32	286.89	16.94	12.58	2.85	31.79	1.18	3.41	4.02	4.29
19	63.32	10.30	25.89	53.02	158.03	12.57	3.53	1.24	28.77	1.11	1.84	2.05	2.15
20	135.55	19.79	46.39	115.76	589.82	24.29	5.56	1.69	52.29	1.13	2.21	2.50	2.64
21	138.42	13.00	26.03	125.42	581.88	24.12	9.87	2.58	35.43	1.36	3.54	4.82	5.51
22	63.09	9.27	30.99	53.83	150.54	12.27	2.39	0.44	33.31	1.07	1.62	1.74	1.81
23	87.38	10.05	19.79	77.34	90.44	9.51	25.47	4.21	21.95	1.11	3.52	3.91	4.05
24	27.93	9.86	13.26	18.07	5.62	2.37	15.50	2.87	13.47	1.02	1.34	1.36	1.37
25	58.13	8.81	12.39	49.32	24.68	4.97	36.75	5.28	13.35	1.08	3.70	3.98	4.07
26	96.22	7.57	13.49	88.65	38.41	6.20	130.52	10.06	14.84	1.10	5.97	6.57	6.73

展开更多

最终提取特征指标组合包括4种钻进参数，共计68项，分别标记为X1~X68。其中，X1~X17为钻速的时频特征指标，X18~X34为扭矩的时频特征指标，X35~X51为钻压的时频特征指标，X52~X68为转速的时频特征指标。对2号斜井的5次超前钻孔的钻进参数进行指标提取后，获取钻进特征指标组合119组。由于现场获取的K_v仅有112组，最终整理出包含112组数据的数据集，用于构建K_v预测模型。采用RReliefF算法对4种钻进参数的特征指标组合进行优选，并根据权重从大到小进行排序，特征权重值越大，表明该特征对于回归模型越重要，为保证所训练模型的准确度，适当放宽特征筛选的数量，本文以0.005作为特征筛选的阈值，最终筛选出23个特征，如图8所示。

图8全屏查看

围岩完整性评价特征指标权重图

特征筛选后，利用Pearson相关系数求解初步筛选后各特征之间的相关性。综合考虑特征权重，剔除具有强相关性(r>0.8)的特征指标，保留权重较大且与其他特征相关性较弱的特征，获取最优的围岩完整性评价指标组合，如图9所示。具体包括钻速最大值、钻速最小值、钻速标准差、钻速均方根、钻速平均频率、扭矩平均值、扭矩均方根、扭矩峰值因子、扭矩脉冲因子、扭矩平均频率、钻压峭度、钻压平均频率、转速最大值、转速最小值、转速平均值、转速峭度、转速均方根、转速平均频率共18个钻进特征指标。

图9全屏查看

用于围岩完整性性评价的最优指标组合相关性热力图

3.3　围岩完整性预测模型构建及结果分析

为对比验证本节提出的BO-RF模型的优越性以及围岩完整性评价指标体系的优越性，从输入特征划分和超参数优化2个方面设计了模型训练方案。方案1为本文的所提方法，其输入指标为优选的18个钻进特征指标(F₁)，预测模型为BO-RF。方案2的输入参数为原始4个钻进参数的各段均值(F₂)，预测模型为BO-RF；对比两者验证围岩完整性评价指标体系的优越性。方案3的输入参数为优选的18个钻进特征指标(F₁)，预测模型为RF，与方案1对比用于验证所提BO-RF模型的优越性。方案1的寻优结果见表6，不同方案条件下模型的评估结果见表7。从表7可以看出：方案1在训练集和测试集上的R²分别为0.959 7和0.925 4，E_RMSE分别为0.02和0.03，E_MAPE分别为6.35%和8.03%。该方案在训练集与测试集上的各项评价指标均表现最好，表明该模型不仅在训练集上有较强的学习能力，在测试集上也表现出了较好的泛化性能。方案2在训练集和测试集上的R²分别为0.872 9和0.842 5，E_RMSE分别为0.09和0.11，E_MAPE分别为11.04%和17.26%。经对比可发现无论是在训练集还是测试集，以围岩完整性评价指标体系作为输入特征的PSO-RF模型在围岩完整性预测性能上有了显著提升。方案3在训练集和测试集上的R²分别为0.928 5和0.908 6，E_RMSE分别为0.06和0.08，E_MAPE分别为9.72%和12.94%。方案3在训练集和测试集上的性能均比方案1的低，说明采用BO算法对RF的部分超参数进行寻优可有效提升模型的预测性能。以上结果表明，本文提出的以围岩完整性评价指标体系作为输入特征的PSO-RF模型在围岩完整性预测上具有显著优越性。

通过预测值-实测值散点对角线图系统比较了训练方案1在训练集与测试集上的回归表现及泛化差异，如图10所示，在理想情形下散点应集中分布于45°对角线上。从图10可以看出：无论是训练集还是测试集，基于钻进特征指标所构建的K_v预测模型中的数据点多集中在等值线附近，表明该模型中的预测值和实际值的一致性较好，具有较高预测精度，能够实现现场隧道围岩完整性系数的有效预测。

图10全屏查看

隧道围岩K_v预测值与真实值对比

3.4　工程验证

根据建立的基于钻进数据定量解译的岩体完整性系数预测方法，结合表3中K_v与围岩完整程度的对应关系，即可实现对隧道围岩完整性的评价。鉴于在后续隧道开挖过程中，仅对预测里程的前18 m进行了连续跟踪和地质分析，因此，仅展示18段K_v预测结果，最终预测出隧道围岩的完整性系数K_v如图11所示。从图11可以看出：预测里程范围内，围岩完整性系数呈现一定的波动性，完整性系数最大为0.41，最小为0.17，平均值为0.30，岩体完整程度为较破碎到破碎。其中，在0~3 m、9~12 m的深度内，围岩完整性系数不低于0.35，围岩较破碎；在3~9 m、12~18 m深度范围内，围岩完整性系数均在0.15~0.35之间，围岩破碎。

图11全屏查看

隧道围岩各段落K_v预测结果

为验证隧道围岩完整性随钻判识结果，在后续隧道开挖过程的掌子面进行连续跟踪，图12为钻深18 m范围内不同开挖面的照片。从图12(a)可以看出：钻深3 m和钻深10 m处掌子面的花岗岩整体呈微风化到弱风化，节理裂隙较发育，节理间为强风化，完整程度属于较破碎，现场围岩分级定为Ⅳ级。本文预测该段围岩K_v分别为0.35和0.39，同样为较破碎，预测结果与实际情况同样较相符。从图12(c)可以看出：钻深15 m处掌子面花岗岩呈弱风化，岩石强度较低，另发育密集节理，完整性较差，掌子面呈点滴状出水与线流状出水，伴随掉块及坍塌现象，现场定为Ⅳ级围岩。本文预测该段围岩K_v为0.17左右，为破碎围岩。预测结果与实际结果基本符合。综上可知，本文所提出的隧道围岩完整性随钻识别方法具有有效性，可为现场施工与支护设计提供一定的指导。

图12全屏查看

隧道开挖后的掌子面照片

本文提出的隧道围岩完整性评价方法在裂隙发育的花岗岩山区铁路隧道中验证效果良好，钻进参数预处理、时频特征解译、RReliefF-Pearson优选、BO-RF预测的流程和思路对其他隧道同样具有可行性。然而，该方法应用于其他工况时主要受两方面约束：首先，模型依赖于钻进参数与岩体完整性之间的统计相关性，当钻机类型、传感器精度或采样频率发生显著变化时，特征分布将偏离训练域，需进行小样本微调或重新训练；其次，在不同工况下，岩石对钻进参数敏感度不同，优选出的完整性预测指标体系可能发生变化。为确保预测精度和工程可靠性，应基于新工程背景重新执行特征优选与后续流程。