隧道掘进机（tunnel boring machine, TBM）是一种大型掘进装备，相比于传统的钻爆法，具有开挖速度快、施工质量高和对周围环境扰动小等优点，因而广泛应用于山岭隧洞、引水隧洞、城市地铁等的开挖中^[1-2]. 然而，TBM对地质条件的适应性较差^[3]，如果遭遇不良的地质条件或者设置的操作参数与地质条件不匹配，将导致掘进速度慢、利用率低的情况.

TBM有4个性能参数：净掘进速度(penetration rate, PR)，表示机器处于掘进状态时的推进速度；利用率(utilization factor, U)，表示TBM 掘进时间与总当班时间的百分比；施工进度(advance rate, AR)，为净掘进速度与利用率的乘积；滚刀寿命，表示单把刀开挖的岩土体积或开挖距离^[4]. 它们表征着TBM的运行成本及施工效率，是TBM隧道建设的重要经济性指标. 当前，对TBM性能参数的研究主要集中在净掘进速度预测^[5-6]和刀具磨损评估^[2]上，对利用率的关注相对较少. 主要有2方面原因：1）利用率影响因素众多，精确评估面临困难；2）各因素导致的停机数据较难完整地获取. 实际上，如果净掘进速度较大而利用率低，比如软岩/土情况下，TBM的施工进度也不会高；反之，如果净掘进速度较小而利用率高，比如硬岩情况下，TBM的施工进度也可能不低. 因此，准确预测TBM的利用率具有重要的意义.

在利用率的定义中，总当班时间包括掘进、换步、支护/拼装、刀具检查与更换、设备维修保养等时间. 据统计，TBM施工项目中掘进机利用率一般为5%~65%，且不同工程掘进机利用率差别较大^[4]. 段晓晨等^[7]基于西康铁路秦岭隧道统计数据，指出加强管理、提高工人水平、加强机械维保能提高利用率. 龚秋明等^[8]基于引汉济渭岭北隧洞统计数据，指出岩体条件或者操作参数的细小变化会对利用率产生较大影响. Alber等^[9]指出地质条件对利用率影响很大，因为岩石性质会影响隧道支护时间和刀盘刀具更换时间. Rostami^[4]指出利用率对地质条件非常敏感，断层破碎带会影响开敞式TBM的利用率，而不良地质导致的掌子面坍塌或者“卡盾”则会降低盾构的利用率. 综上所述，TBM利用率的影响因素众多，其中地质条件和操作参数是2个重要的影响因素.

在国外，CSM与NTNU模型是2个最经典的模型，它们通过将TBM 的停机时间分配成不同的活动来评估利用率^[10]，但是这2个模型提出的时间较早，在当下应用中精度不高^[11]. Simoes 等^[12]以机器直径、岩体质量等级、地下水流入量等为输入，利用模糊逻辑方法建立了 TBM 利用率的预测模型. Frough等^[13]应用岩石工程系统来计算与地质和岩体相关的停机时间指标，结果表明其相比传统的统计方法具有更高的精度. Noori等^[14]分别利用多元线性回归、人工神经网络和人工神经网络-粒子群优化算法建立利用率预测模型，模型的输入考虑地质条件、推力与扭矩等，结果表明人工神经网络-粒子群优化算法表现最佳. 在国内，利用率预测模型的开发落后于国外，主要停留在定性分析的层面.

以上研究，要么是在规划阶段对利用率进行预估，要么是对现场施工管理提出一些建设性意见，并没有研究在现有的地质状况、设备规格、施工管理水平等条件下应如何实时预测利用率并优化司机操作以提高TBM的施工进度. 因此，本研究以新加坡某双层隧道项目为背景，着重分析在现有的设备规格、施工管理水平下地质类型和操作参数对利用率的影响，建立动态的利用率预测模型，并在此基础上对操作参数进行优化以最大化施工进度.

SVR适用于解决小样本、非线性问题，在TBM隧道施工中有许多应用，如刀盘扭矩预测^[15]. 给定训练样本D={(x₁, y₁), (x₂, y₂) $,\cdots , $ (x_m, y_m)}，其中x_i∈Rⁿ, y_i∈R. 对于样本(x, y)，SVR能容忍模型输出f(x)与真实值y之间最多有ε的偏差，当f(x)与y之间差别的绝对值大于ε时才计算损失，如图1所示. 图中，ε为模型能容忍的最大偏差，ω为模型的权重，b为模型的偏置.

SVR可以表示为如下数学形式：

(1) $ \mathop {\min }\limits_{{\boldsymbol{\omega}} ,b}\; {\text{ }}\frac{1}{2}{\left\| {\boldsymbol{\omega}} \right\|^2}+C\sum\limits_{i = 1}^m {{l_\varepsilon }\left( {f\left( {{{\boldsymbol{x}}_i}} \right) - {y_i}} \right)} {\text{ }}{\text{.}} $

式中：C为正则化常数， $ {l_\varepsilon } $为ε-不敏感损失函数. $ {l_\varepsilon } $表达式如下：

(2) $ {l_\varepsilon }\left( z \right) = \left\{ \begin{array}{*{20}{l}} {0,}& {\left| z \right| \leqslant \varepsilon ;} \\ { \left| z \right| - \varepsilon ,} & {\left| z \right| > \varepsilon .} \\ \end{array} \right. $

引入松弛变量 ${\xi _i}$和 ${\hat \xi _i}$，可以将式（1）转化为如下形式：

(3) $ \left. \begin{array}{l} \mathop {\min }\limits_{{\boldsymbol{\omega}} ,b,{\xi _i},{{\hat \xi }_i}} {\text{ }}\dfrac{1}{2}{\left\| {\boldsymbol{\omega }} \right\|^2}+C\displaystyle \sum\limits_{i = 1}^m {\left( {{\xi _i}+{{\hat \xi }_i}} \right)}; \\ \quad\quad{\text{ }}{\rm{s.t.}}{\text{ }}f\left( {{{\boldsymbol{x}}_i}} \right) - {y_i} \leqslant \varepsilon +{\xi _i}, \\ \quad\quad{\text{ }}{y_i} - f\left( {{{\boldsymbol{x}}_i}} \right) \leqslant \varepsilon +{{\hat \xi }_i}, \\ \quad\quad{\text{ }}{\xi _i} \geqslant 0,\;{{\hat \xi }_i} \geqslant 0,\;i = 1,2,\cdots ,m{\text{ }}{\text{.}} \end{array} \right\} $

通过引入拉格朗日乘子 ${\mu _i} \geqslant 0$， ${\hat \mu _i} \geqslant 0$， ${\alpha _i} \geqslant 0$， ${\hat \alpha _i} \geqslant 0$，由拉格朗日乘子法可以得到式（3）的拉格朗日函数：

(4) $ \begin{split} & L\left( {{\boldsymbol{\omega}} ,b,{\boldsymbol{\alpha}} ,\hat {\boldsymbol{\alpha}} ,\xi ,\hat \xi ,{\boldsymbol{\mu}} ,\hat {\boldsymbol{\mu}} } \right) = \frac{1}{2}{\left\| {\boldsymbol{\omega}} \right\|^2}+C\sum\limits_{i = 1}^m {\left( {{\xi _i}+{{\hat \xi }_i}} \right)} - \\ &\quad\quad \sum\limits_{i = 1}^m {{\mu _i}{\xi _i}} - \sum\limits_{i = 1}^m {{{\hat \mu }_i}{{\hat \xi }_i}} +\sum\limits_{i = 1}^m {{\alpha _i}\left( {f\left( {{{\boldsymbol{x}}_i}} \right) - {y_i} - \varepsilon - {\xi _i}} \right)} + \\ &\quad\quad\sum\limits_{i = 1}^m {{{\hat \alpha }_i}} \left( {{y_i} - f\left( {{{\boldsymbol{x}}_i}} \right) - \varepsilon - {{\hat \xi }_i}} \right){\text{ }}{\text{.}} \end{split} $

再令 $L\left( {{\boldsymbol{\omega}} ,b,{\boldsymbol{\alpha}} ,\hat {\boldsymbol{\alpha}} ,\xi ,\hat \xi ,{\boldsymbol{\mu}} ,\hat {\boldsymbol{\mu}} } \right)$对ω、b、 ${\xi _i}$和 ${\hat \xi _i}$的偏导为0，并将求出的方程代入式（4），即可求得SVR的对偶问题：

(5) $ \begin{split} & \mathop {\max }\limits_{{\boldsymbol{\alpha}} ,\hat {\boldsymbol{\alpha}} } {\text{ }}\sum\limits_{i = 1}^m {\left( {{y_i}\left( {{{\hat \alpha }_i} - {\alpha _i}} \right) - \varepsilon \left( {{{\hat \alpha }_i}+{\alpha _i}} \right)} \right)} - \\ &\quad\quad {\text{ }}\frac{1}{2}\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^m {\left( {{{\hat \alpha }_i} - {\alpha _i}} \right)\left( {{{\hat \alpha }_j} - {\alpha _j}} \right)} } {\boldsymbol{x}}_i^{\text{T}}{{\boldsymbol{x}}_j}; \\ &\quad\quad {\rm{s.t.}}{\text{ }}\sum\limits_{i = 1}^m {\left( {{{\hat {{\alpha}} }_i} - {{{\alpha}} _i}} \right)} = 0,{\text{ 0}} \leqslant {\alpha _i},{{\hat \alpha }_i} \leqslant C{\text{ }}{\text{.}} \end{split}$

上述过程须满足KKT条件，最终得到SVR的解或决策函数如下：

(6) $ f\left( {\boldsymbol{x}} \right) = \sum\limits_{i = 1}^m {\left( {{{\hat \alpha }_i} - {\alpha _i}} \right){\boldsymbol{x}}_i^{\text{T}}{\boldsymbol{x}}} +b, $

(7) $ b = {y_i}+\varepsilon - \sum\limits_{j = 1}^m {\left( {{{\hat \alpha }_j} - {\alpha _j}} \right)} {\boldsymbol{x}}_j^{\text{T}}{{\boldsymbol{x}}_i}{\text{ }}{\text{.}} $

若考虑特征的非线性映射（ ${\boldsymbol{x}}{\text{ }} \to \phi \left( {\boldsymbol{x}} \right)$），则SVR的解或决策函数如下：

(8) $ f\left( {\boldsymbol{x}} \right) = \sum\limits_{i = 1}^m {\left( {{{\hat \alpha }_i} - \alpha {}_i} \right)} {\text{ }}k\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_i}} \right)+b{\text{ }}{\text{.}} $

式中： $k\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_i}} \right) = \phi {\left( {{{\boldsymbol{x}}_i}} \right)^{\text{T}}}\phi \left( {{{\boldsymbol{x}}_i}} \right)$为核函数，常用的核函数有径向基函数、线性函数和多项式函数.

PSO是一种全局优化算法^[16]，它源于对鸟群捕食行为的研究，具有结构简单、收敛速度快、适应性强等优点，在参数优化^[17]领域有较好的应用.

在PSO中, 每个粒子都代表着寻优空间的1个可能解，粒子有其位置和速度，粒子位置坐标对应的目标函数f(x)值即为该粒子的适应度. 在每1轮迭代进化过程中, 粒子通过自身历史最优位置p_best和群体的历史最优位置g_best来更新当前的速度和位置，直至收敛到最优解. 假设X=[X₁, X₂$,\cdots , $ X_n]表示大小为n的种群，每个粒子的搜索空间为d维，X_i=[x_i1, x_i2$,\cdots , $ x_id]和V_i=[v_i1, v_i2$,\cdots , $ v_id]分别表示第i个粒子的位置和速度，则粒子的速度和位置更新公式如下：

(9) $ {\boldsymbol{V}}_i^{t+1} = \omega {\boldsymbol{V}}_i^t+{c_1}r_{i1}^t\left( {{\boldsymbol{p}}_{{\text{best}},i}^t - {\boldsymbol{X}}_i^t} \right)+{c_2}r_{i2}^t\left( {{\boldsymbol{g}}_{{\text{best}}}^t - {\boldsymbol{X}}_i^t} \right), $

(10) $ {\boldsymbol{X}}_i^{t+1} = {\boldsymbol{X}}_i^t+{\boldsymbol{V}}_i^{t+1};{\text{ }} {i = 1,2,\cdots ,n} {\text{ }}{\text{.}} $

式中： ${\boldsymbol{V}}_i^{t+1}$表示第 $i$个粒子在第t+1次迭代时的速度，初始化速度为0；ω为惯性权重，用于控制全局搜索和局部寻优能力；c₁为认知常数，用于实现自我认知，通常取值为2；c₂为社会常数，用于模拟与社会种群的交互，通常取值也为2； $r_{i1}^t$和 $r_{i2}^t$为正态分布的随机数，范围为[0，1.0]，用于增加粒子的随机性，避免陷入局部最优； ${\boldsymbol{p}}_{{\text{best}},i}^t$为截止到第t次迭代，第i个粒子的历史最优位置； ${\boldsymbol{g}}_{{\text{best}}}^t$为截止到第t次迭代，种群的历史最优位置； ${\boldsymbol{X}}_i^t$为第t次迭代时第i个粒子的位置.

PSO求解带约束的优化问题，须在原始的目标函数f(x)中加入与约束条件有关的惩罚项^[18]. 改进的目标函数F(x)计算如下：

(11) $ F\left( {\boldsymbol{x}} \right) = f\left( {\boldsymbol{x}} \right)+h\left( t \right)H\left( {\boldsymbol{x}} \right){\text{ }}, $

(12) $ h\left( t \right) = \sqrt t {\text{ }}或{\text{ }}t\sqrt t {\text{ }},$

(13) $ H\left( {\boldsymbol{x}} \right) = \sum\limits_{i = 1}^m {\theta \left( {{q_i}\left( {\boldsymbol{x}} \right)} \right)} {q_i}{\left( {\boldsymbol{x}} \right)^{\gamma \left( {{q_i}\left( {\boldsymbol{x}} \right)} \right)}}{\text{ }}, $

(14) $ {q_i}\left( {\boldsymbol{x}} \right) = \max \;\left( {0,{\text{ }}{g_i}\left( {\boldsymbol{x}} \right)} \right){\text{ }}{\text{,}} $

(15) $ \theta \left( {{q_i}\left( {\boldsymbol{x}} \right)} \right) = \left\{ \begin{array}{*{20}{l}} {10,} &{{q_i}\left( {\boldsymbol{x}} \right) < 0.001;} \\ {20,}&{0.001 \leqslant {q_i}\left( {\boldsymbol{x}} \right) \leqslant 0.1;} \\ {100,}&{0.1 < {q_i}\left( {\boldsymbol{x}} \right) \leqslant 1.0;} \\ {300,}&{其他.} \\ \end{array} \right. $

(16) $ \gamma \left( {{q_i}\left( {\boldsymbol{x}} \right)} \right) = \left\{ \begin{array}{*{20}{l}} {1,}&{{q_i}\left( {\boldsymbol{x}} \right) < 1.0 ;} \\ {2,} &{其他.} \\ \end{array} \right. $

式中：h(t)为动态更新的惩罚系数，与迭代次数t相关，通常可使用 $t\sqrt t $，如果问题规模小，可以使用 $\sqrt t $；H(x)为约束惩罚项；q_i(x)为相对约束惩罚函数；g_i(x)为约束条件；θ(q_i(x))为分段赋值函数；γ(q_i(x))为惩罚指数.

本研究提出基于PSO-SVR的TBM施工进度优化方法，流程如图2所示.

首先，对原始数据进行预处理，包括：1）筛选出含地质标签的环的推力、扭矩、刀盘转速和推进速度；2）进一步筛选出掘进状态对应的数据，将一环数据划分为上升段与稳定段，并选出稳定段数据；3）利用3σ准则检测数据中异常值并进行插值处理，对处理后的数据按环求均值. 注意：由于刀盘转速会有正负值（对应正反转），对于刀盘转速的预处理是取绝对值后再求平均.

其次，构建数据集，利用上隧道数据建立利用率预测模型U=f₀(T, F, n, v, geo)，利用下隧道数据验证模型的有效性.

再次，建立施工进度优化方程：

(17) $ \left. \begin{array}{l} \max \; f\left( {n,v} \right) = {\text{AR}} = v {f_0}\left( {T,{\text{ }}F,{\text{ }}n,{\text{ }}v,{\text{ }}{\bf{geo}}} \right){\text{ }}{\text{;}} \\ {\rm{s.t.}}\; F = {g_1}\left( {n,{\text{ }}v,{\text{ }}{\bf{geo}}} \right){\text{ }}{\text{,}} \\ \quad\quad T = {g_2}\left( {n,{\text{ }}v,{\text{ }}{\bf{geo}}} \right){\text{ }}, \\ \quad\quad T n \leqslant P{\text{ }}{\text{,}} \\ \quad\quad {n_{\text{l}}} \leqslant n \leqslant {n_{\text{u}}}{\text{ }}{\text{,}} \\ \quad\quad {v_{\text{l}}} \leqslant v \leqslant {v_{\text{u}}}{\text{ }}{\text{.}} \end{array} \right\} $

式中：AR为施工进度，n为刀盘转速，v为推进速度，f₀(T, F, n, v, geo)为利用率预测模型，T为扭矩，F为推力，geo为地质类型，P为刀盘最大功率，n_l和n_u为刀盘转速的上、下界，v_l和v_u为推进速度的上、下界. 最后，利用PSO优化式（17），得到特定地质类型下最优的操作参数（n^*,v^*）.

项目是新加坡地铁系统的一部分，待研究标段长约710 m，从珊顿大道站到麦克斯韦路站，为堆叠式双层隧道，如图3所示. 上、下隧道分别采用相同型号土压平衡盾构(earth pressure balance shield, EPB)开挖，其主要技术规格如表1所示. 表中，D为开挖直径，T_max为最大扭矩，F_max为最大推力，v_max为最大推进速度，R为环号. 这2条隧道均采用预制钢筋混凝土管片衬砌，掘进环编号为1~500.

隧道沿线的地质信息主要来源于新加坡政府文件和现场勘探. 隧道沿线会遭遇多种地质构造，包括卡兰组、古冲积层、坎宁堡卵石层和句容组. 整体上，有4~15 m厚的填充层覆盖在卡兰组上方，卡兰组下方依次是古冲积层和句容组. 隧道沿线布设有15个钻孔，其中14个穿过上下隧道. 如表2所示为钻孔的地质信息和相应的环号. 可以看出，在这些钻孔中，句容组最常见，其次是坎宁堡卵石层. 隧道断面除了单一的地质构造类型外，还存在复合类型，总体地质条件较复杂. 值得注意的是，距离较近的钻孔更有可能具有相似或相同的地质类型，如10号和11号钻孔、12号和13号钻孔.

EPB的数据采集系统每5 s采集1446个现场参数，包括推力、扭矩、刀盘转速、推进速度等. EPB是一环一环地掘进的，一个完整的掘进循环包括掘进状态、等待状态和拼装状态. 其中，等待状态主要由土箱车运载能力有限导致，拼装状态主要用于管片拼装以及设备检修. 如图4所示为上隧道第90环推力随着系统状态的变化. 图中，S为系统状态（0、1、2分别表示拼装、等待和掘进状态），s为采样序号. 可以看出，推力在等待状态和拼装状态为0. 如图5 (a)~(d)所示分别为上隧道EPB掘进状态下推力、扭矩、刀盘转速和推进速度的分布. 图中，c为频数.

当前的研究主要是从一定掘进距离或者工作日的尺度来计算利用率^{[12, 14]}，而本研究则将利用率的计算细化到掘进环尺度. 将一环掘进状态持续时间除以所有状态持续时间便得到该环的利用率，而1~500环利用率的平均值即为整个标段的平均利用率. 如图6所示为上、下隧道各环利用率的变化情况. 可以看出，各环利用率在0~70%波动，上隧道平均利用率为41.5%，下隧道平均利用率为37.1%. 由引言中前人的统计结果^[4]可知，TBM利用率一般在5%~65%，而本项目少量环的利用率超过了该范围. 利用率过低主要由非掘进状态耗时严重引起，具体原因为等待土箱车时间过长或者设备故障导致的停机检修时间过长；利用率过高主要是由掘进时间过长引起的，在一环掘进长度固定的情况下，这意味着净掘进速度过低. 另外，还可以发现上、下隧道沿着掘进方向利用率的变化相似，这可以从侧面证明利用率与地质条件有关，因为上、下隧道的地质条件较相似. 由于本研究主要关注地质条件和司机操作参数对利用率的影响，因此对利用率的“异常值”进行预处理，用前后2环的利用率的均值对“异常值”进行插值，并将处理后的利用率也绘制在图6中. 如图7所示为上、下隧道不同地质类型下利用率的统计结果. 图中，JR- V、JR-IV、FCBB、FC和MC分别表示句容组(V)、句容组(IV)、坎宁堡卵石层、河流黏土和海洋黏土. 可以发现：1）整体上，EPB在坎宁堡卵石层下具有最高的利用率，其次是句容组(V&IV)和河流黏土，EPB在海洋黏土下利用率最低；2）即使在相同的地质条件下，EPB也会具有不同的利用率.

如图8 (a)、(b)所示为操作参数刀盘转速和推进速度对利用率的影响. 可以看出：1）刀盘转速对利用率的影响不如推进速度那么明显；2）整体上，随着推进速度的加快，利用率有逐渐下降的趋势. 由于施工进度为利用率与推进速度的乘积，而利用率又随着推进速度的加快而减小，必然存在一个最优的推进速度使得施工进度最大.

操作参数除了对利用率有影响外，还对载荷有影响. 如图9 (a) 、(b)所示分别为不同的地质类型和贯入度下扭矩和推力的变化，同时还在图中列出了相应的拟合曲线及其R²值. 图中，PR为贯入度，它是推进速度与刀盘转速的比值. 可以发现：1）在句容组(IV)和河流黏土下，扭矩随着贯入度的增大而增大，其他地质类型则是呈现相反的趋势；2）推力基本上随着贯入度的增大而减小，句容组(V)是例外. 从地质具有连续性^[19]的认识出发，本研究假设钻探点所在的环的前后2环与该环具有相同的地质类型，而表2也能在一定程度上证实该假设. 该假设有助于丰富数据集，进而开展更为深入的分析.

综上所述，地质类型和操作参数均对EPB的利用率有很大的影响. 在实际施工中，司机一般凭借经验判断地质类型并据此设置相应的操作参数，但这不能保证操作参数是最优的. 如何设置合适的操作参数，使得盾构掘进得又快又稳，这是一个值得研究的问题，也是本研究力图解决的问题.

由表2可知，上隧道钻探点对应的地质类型有5种，分别是句容组(V)、句容组(IV)、坎宁堡卵石层、河流黏土和海洋黏土，其对应的地质标签分别设置为0、1、2、3、4. 进一步，对数字标签进行向量化处理，这有助于建模成功，如标签0转换为[1, 0, 0, 0, 0]，标签1转换为[0, 1, 0, 0, 0]，其余依次类推. 下隧道的地质条件较复杂，除了有上隧道中出现的句容组(V)、句容组(IV)和坎宁堡卵石层外，还有句容组、古冲积层和河流黏土的复合地层. 由于复合地层复杂多变且数据十分有限，本研究暂不考虑复合地层这种情况. 接着，按照1.3节介绍的数据预处理流程处理数据，利用上隧道的数据构造训练集，利用下隧道数据构造测试集. 为了调节并获得较优的模型参数，对训练集按照7∶3进一步细分为训练集和验证集. 模型的输入为扭矩(T)、推力(F)、刀盘转速(n)、推进速度(v)和地质类型(geo，向量形式)，输入的维度为9；模型的输出为利用率(U)，输出的维度为1. 如表3所示为训练集和测试集的信息.

采用2个指标来评估模型的性能，分别为均方误差(mean square error, MSE)和决定系数(coefficient of determination, R²)，其计算公式如下：

(18) $ {\text{MSE}} = \frac{1}{N}{\sum\limits_{i = 1}^N {\left( {{y_i} - {{\hat y}_i}} \right)} ^2}{\text{ }}{\text{,}} $

(19) $ {{{R}}^{\text{2}}} = 1 - {{\displaystyle \sum\limits_{i = 1}^N {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}\Bigg/{{\displaystyle \sum\limits_{i = 1}^N {{{\left( {{y_i} - {y_{\text{m}}}} \right)}^2}} }}{\text{ }}{\text{.}} $

式中：N为样本量，y_i为回归变量的真实值， ${\hat y_i}$为y_i对应的预测值，y_m为回归变量真实值的均值.

采用网格搜索法来获取最优的模型参数. 对于SVR，首先选择核函数类型，常见的核函数有线性核、多项式核和高斯核；然后设置惩罚系数C和核系数γ的寻优范围. 此处，核函数选为线性核，ln C 与ln γ的寻优范围均设置为[−7.0, 0]，步长为0.2. 如图10所示为验证集下，不同C与γ组合的MSE的变化. 可以看出：1）在线性核条件下，γ的取值对模型的性能无影响；2）当C=0.002时，模型在验证集上的MSE最小，为0.0038，此时对应的模型参数是最优的. 如图11 (a) 、 (b)所示为SVR模型在验证集与测试集上的表现. 图中，U_t、U_p分别表示利用率的真实值与预测值. 可以看出，模型在验证集与测试集上的R²分别为0.729和0.625. 因此，预测结果与真实值较吻合，模型在另一条隧道上具有良好的泛化能力.

为了进一步证明SVR模型的优越性，采用多元线性回归(multiple linear regression, MLR)、岭回归(ridge regression, RR)及传统的机器学习算法，如决策树(decision tree, DT)、k最近邻(k-nearest neighbors, KNN)、随机森林(random forest, RF)、AdaBoost、XGBoost等，建立每环利用率预测模型. 由于训练样本有限，本研究不考虑采用深度学习算法建模. 如图12所示为各算法模型在测试集上的表现. 可以看出，SVR模型表现最佳，其次是RR模型和MLR模型，DT模型表现最差. 如表4所示为各算法模型在验证集与测试集上的性能参数，表中还列出了各模型最优的模型参数. 可以看出：1）SVR模型在验证集与测试集上均表现最佳，其在验证集和测试集上的MSE分别为0.003 8、0.004 5；2）DT模型与RF模型虽然在验证集上表现良好，但是在测试集上却表现较差，表明模型过拟合了. 另外，值得注意的是，RF、AdaBoost、XGBoost等在许多其他应用领域中表现良好的算法在测试集上的表现还不如线性模型MLR和RR. 训练样本少，算法模型无法有效地学习到内在的非线性映射关系，是这些算法模型精度不高的一个重要原因. 这也说明，SVR在小样本建模上有其特有优势.

根据图5，选定刀盘转速的优化范围为1.0~3.5 r/min，推进速度的优化范围为5.0~40.0 mm/min. 设置PSO的种群大小为300，认知常数c₁和社会常数c₂均为2.0，最大进化代数为100. 如图13所示为句容组(V)下，PSO寻优时目标函数值的变化. 图中，N_g为进化代数. 可以看出，PSO进化到20代就已经找到最优解了，最大的施工进度为8.89 mm/min. 如图14所示为不同地质类型下，通过枚举法绘制的AR=f(n,v)三维曲面，同时将PSO寻优的结果（绿色球表示）也展示在图中. 可见：1）相比于推进速度，刀盘转速对施工进度的影响较小；2）当推进速度较慢时，施工进度会随着推进速度的加快而增大，而当推进速度快到一定程度后，施工进度会有下降的趋势；3）PSO能准确地寻找到最优的操作参数，使得施工进度最大. 如表5所示为不同地质类型下PSO寻优的结果. 表中，AR^*为最优施工进度. 句容组(V)、句容组(IV)、坎宁堡卵石层、河流黏土和海洋黏土下最大施工进度分别为8.89、9.34、8.78、7.49 、6.45 mm/min.

本研究提出基于PSO-SVR的EPB施工进度优化方法，结论如下：

（1）施工进度优化研究的关键是首先建立利用率预测模型，然后以此建立施工进度优化方程；

（2）将利用率的计算细化到掘进环尺度，能丰富数据集，并为掘进环内的操作优化奠定基础；

（3）地质类型、操作参数与载荷对利用率有重要影响，且随着推进速度的加快，利用率有逐渐下降的趋势；

（4）SVR适合小样本建模，建立的利用率预测模型在测试集上的R²为0.625，优于多元线性回归、决策树、k最近邻、随机森林、AdaBoost和XGBoost模型；

（5）PSO具有全局寻优能力，能准确地找出不同地质类型下最优的操作参数，使得施工进度最大.

综上，所提方法能准确地预测EPB的利用率并优化其操作，使其安全、快速地掘进，具有较好的工程应用价值. 然而，所建立的基于SVR的利用率预测模型的泛化能力只是在同一项目的另一条相近的隧道中展开了验证，未来考虑将模型推广应用到其他具有相似地质条件的隧道中.