坑道近水平回转钻进是保障煤炭资源安全高效开采的重要手段之一。在钻进过程中，司钻调整操作参数时往往依赖自身经验，难以实现最优钻进效率同时节省钻进成本。为平衡钻进效率、钻进成本、钻进安全等多个指标，研究钻进过程多目标优化，给出操作参数的最优推荐值，是保证钻进性能的关键。

由于钻进过程充满复杂性和耦合性，单目标操作参数优化在实际应用中往往存在局限。为平衡钻进过程的经济成本，学者们利用机理方法建立了评价钻速、钻头寿命及钻进效率的多种钻进多参数多目标优化模型^[1-2]；还有学者运用混合蝙蝠算法、支持向量机等数据驱动方法建立了钻速预测模型^[3-4]。基于这些目标评价模型，学者们运用粒子群算法^[5]、改进蚁群算法^[6]、非支配排序遗传(NSGA-II)算法^[7]等智能优化算法得到了钻进过程操作参数的最优推荐结果。利用随钻测量技术获知孔底钻进数据为进行操作参数优化提供了便利^[8]。但由于随钻测量设备成本较高且多以轨迹参数测量为主，其在坑道回转钻进中应用较少。为分析钻头的运动状态和钻柱系统的稳定性，有学者根据钻柱系统结构建立了钻柱动力学模型^[9-10]。还有学者将钻柱动力学模型、目标函数和操作变量相结合，提出了一种基于模型的操作参数优化方法^[11]。因此，相比于较高成本的孔底信息测量，建立钻柱动力学模型进行孔底信息估计，是分析孔底钻柱运动状态的一个重要方式。

笔者以坑道近水平回转钻进为研究对象，建立了轴向和扭转维度的钻柱动力学模型，设计全维状态观测器，建立孔口−孔底钻柱运动状态映射关系；提出了机械钻速和钻头磨损的优化目标评价模型，基于孔口数据和孔底状态估计，利用NSGA-II算法优化转速和给进压力这两个参数；最后利用淮南某煤矿实钻数据进行验证与分析，以验证所提方法的有效性。

1   煤矿井下坑道回转钻进过程与特性

在煤矿井下钻进过程中，坑道钻机的液压系统在扭转维度上通过液压泵将液压油输送至液压马达驱动动力头回转，进而带动钻杆旋转运动；在轴向维度上通过输送液压油至液压缸，使液压缸推动钻杆给进或起拔。钻头不断旋转和推进，从而实现碎岩钻进。在实际钻机操作中，司钻往往通过遥控器下发控制指令调整动力头转速和给进系统给进压力，以适应不同地质条件和钻孔要求^[12]。

随着煤矿坑道钻进不断深入，地层分布愈加复杂，可能存在褶皱、断层等地质构造，容易造成卡钻埋钻事故，安全钻进难度大。为了确保煤矿坑道钻进的安全和效率，需要对钻机操作参数进行严格的监控和调整，以保持煤矿钻进过程稳步进行。煤矿井下坑道钻进过程的复杂性、非线性使得钻进参数受到多种因素的影响，研究钻进参数多目标优化方法能够提高钻进效率和成孔质量。

2   钻进参数多目标优化

煤矿井下钻进参数多目标优化通过考虑多个目标，综合评估不同参数组合的优劣，从而选择最优的钻进参数，这样可以提高钻进效率、降低成本、减少风险，从而提高煤矿的生产效益和安全性。针对本文的研究问题，设计的技术路线如图1所示。

图  1  钻进参数多目标优化技术路线

Figure  1.  Roadmap of multi-objective optimization technology for drilling parameters

下载: 全尺寸图片 幻灯片

首先将孔口信息输入到钻柱运动模型获取状态信息，进而通过非线性状态观测器得到孔内状态估计信息；然后利用获取的孔内信息，通过非线性多元回归方法辨识机械钻速模型的煤岩层系数、门限钻压和转速指数3个系数。最后利用NSGA-II算法进行多目标优化求解，得到转速和钻压等操作参数，并对优化的结果进行分析。

2.1   基于孔口数据驱动的钻柱状态估计

针对煤矿坑道近水平回转钻进中孔内数据难以获得、操作参数优化缺少依据等问题，进行基于轴向和扭转钻柱模型的非线性状态观测器设计。

所使用的动力学模型采用集中质量法对钻柱系统各部分动力学特性进行建模。其中，轴向钻柱动力学模型采用下式表示：

式中：j_al为各集中质量单元的质量，l取1，2，3，下标a为轴向钻柱模型；z_l、$\dot {\textit{z}_l}$、$\ddot {\textit{z}_l}$分别为各单元的轴向位移、轴向速度和轴向加速度，l取1，2，3；h_al为各单元轴向阻尼，l取1，2，3；k_al为各单元轴向刚度，l取1，2；W₀为给进力；W_b为轴向钻头−煤岩相互作用力。

选择${[{{\textit{z}}_1},{{\textit{z}}_2},{{\textit{z}}_3}]^{\text{T}}}$作为状态变量可构成模型系数矩阵J_a、H_a、K_a，表示为：

扭转钻柱动力学模型采用下式表示：

式中：j_tl、h_tl、k_tl分别为钻柱系统各单元的转动惯量、扭转阻尼、扭转刚度，j_tl和h_tl中l取1，2，3，k_tl中l取1，2，下标t为扭转钻柱模型；$ {\theta _l} $、$ {\dot \theta _l} $、$ {\ddot \theta _l} $分别为钻柱系统各单元的转动位移、转速和扭转加速度；T₀为动力头输入扭矩；T_b为钻头−煤岩相互作用带来的负载扭矩。选择$ {[{\theta _1},{\theta _2},{\theta _3}]^{\text{T}}} $作为状态变量可构成模型系数矩阵J_t、H_t、K_t，表示为：

为后续进行状态变量观测，需要重构状态变量，选择$ {[{\dot {\textit{z}}_1},{\dot {\textit{z}}_2},{\dot {\textit{z}}_3},{{\textit{z}}_2} - {{\textit{z}}_1},{{\textit{z}}_3} - {{\textit{z}}_2}]^{\text{T}}} $作为轴向钻柱模型的状态变量x_a，选择$ {[{\dot \theta _1},{\dot \theta _2},{\dot \theta _3},{\theta _2} - {\theta _1},{\theta _3} - {\theta _2}]^{\text{T}}} $为扭转钻柱模型的状态变量x_t。

由于轴向钻柱模型和扭转钻柱模型的状态空间方程结构相同，结合所重构状态向量，通过矩阵代换，采用以下通式对轴向和扭转钻柱动力学模型的状态空间方程进行描述：

式中：i=a为轴向钻柱模型，i=t为扭转钻柱模型；x_i为钻柱模型的状态变量；y_i为输出变量；A_i为状态转移矩阵；B_i为输入矩阵；R_i为钻头−煤岩作用输入矩阵；C_i为输出矩阵；n为模型的自由度；J_i、H_i、K_i为动力学系数矩阵，由已建立的轴向、扭转动力学模型系数确定，具体形式可参考文献[13]；u_i为钻柱模型输入(给进力或扭矩)；w_i(x_i)为钻柱模型钻头−煤岩作用，||w_i(x_i)||<M_i；M_i为钻头−煤岩相互作用产生的最大值，并且${{\boldsymbol{S}}_1} = {[1,{0_{1 \times (n - 1)}}]^{\text{T}}}$，$ {{\boldsymbol{S}}_2} = {[{0_{1 \times (n - 1)}}, - 1]^{\text{T}}} $。

后续对状态观测器进行李雅普诺夫稳定性分析时，会用到以下几个引理。

引理1：当V(t)是连续函数，w(t)为指数有界函数时，如果存在$ \alpha $>0和b>0，那么对于所有t>$ {t_0} $有以下不等式成立：

式中：t为时间；t₀为初始时间；b为标量。

进一步可知：

引理2：当系统具有指数衰减率$ \alpha $>0并存在指数最终有界时，若有常数$ \beta $≥1，则：

其中：

式中：$ \kappa $为$ ||x(t)|{|^2} $的最终上界；$ {\eta _2} $为给定标量，代表$ \theta $所存在的范围大小。

状态空间方程的优点是能够更直观地描述系统的行为，也更容易进行分析和控制。它可以用来描述线性和非线性系统，也可以用来描述时变和时不变系统。状态空间方程在控制工程、信号处理、机器人控制、电力系统等领域都有广泛的应用。

状态观测器是一种用于估计系统状态的工具。它的作用是通过测量系统的输出和输入，以及系统的动态模型，来估计系统的状态变量，包括未被直接测量的状态变量。基于孔口已知信息，通过状态观测器估计孔内信息，状态观测器原理如图2所示，图中s表示积分。

图  2  孔内状态观测器原理

Figure  2.  Principle of the in-hole state observer

下载: 全尺寸图片 幻灯片

式(5)中引入反馈矩阵L，并加入观测器估计值与系统的被估计量的误差$ {y_i} - {\hat y_i} $，构造如下状态观测器：

式中：$ {\hat x_i} $为输入估计状态变量；$ {\hat y_i} $为输出估计状态变量；${\hat w}_i(x_i)$为估计的钻头−煤岩作用。

基于建立的模型状态空间方程，输入给进力或扭矩，利用可测量参数(动力头转速和动力头钻速)，通过观测器估计煤矿井下孔底所需参数(钻头转速和钻头钻速)。

引入状态估计误差量$ {r_i} = {x_i} - {\hat x_i} $来分析观测器的收敛性及收敛速度，式(5)与式(11)相减可得：

式(12)中引入状态估计误差量r_i可得：

已知||w_i(x_i)||<M_i，||$ {\hat w_i}({x_i}) $||<M_i，M_i为钻头−煤岩作用最大值，设定${\varphi _i}({x_i}) = {w_i}({x_i}) - {\hat w_i}({x_i})$，则$||{\varphi _i}({x_i})||\leqslant2{M_i}$。矩阵A_i、B_i和C_i由已建立模型参数确定，所以系统状态估计误差能否衰减为零，主要取决于反馈矩阵L，若能选取合适的反馈矩阵L，就能使系统状态估计误差具有良好的衰减率。

针对的状态估计误差r_i，构造李雅普诺夫函数V：

式中：P为正定对称矩阵。对式(14)求导，并在等式两边加上$ 2\alpha V $可得：

其中，$ \alpha > 0 $且为实数。

由上述条件可知，存在$ {\varphi _i}{({x_i})^{\text{T}}}\varepsilon {\varphi _i}({x_i}) < 4{M_i}^2\varepsilon $，其中$ \varepsilon $为设置大于0的中间变量，后续可根据线性矩阵不等式计算得出，代入式(15)可得：

设$ \sigma (t) = {[{r_i},{\varphi _i}({x_i})]^{\text{T}}} $，基于线性矩阵不等式性质设置反馈增益矩阵$ {\boldsymbol{L}} = \mu {{\boldsymbol{P}}^{ - 1}}{{\boldsymbol{C}}_i}^{\text{T}} $，$ \mu $为代求常数，代入式(16)可得：

其中，$\Sigma = {\boldsymbol{P}}{{\boldsymbol{A}}_i} + {{\boldsymbol{A}}_i}^{\text{T}}{\boldsymbol{P}} - 2\mu {{\boldsymbol{C}}_i}^{\text{T}}{{\boldsymbol{C}}_i} + 2\alpha {\boldsymbol{P}}$$\left[ {\begin{array}{*{20}{c}} \Sigma &{{\boldsymbol{P}}{{\boldsymbol{R}}_i}} \\ {{{\boldsymbol{R}}_i}^{\text{T}}{\boldsymbol{P}}}&{ - \varepsilon } \end{array}} \right] < 0$成立时，意味$ \dot V + 2\alpha V < 4{M_i}^2\varepsilon $，由引理1可知存在：

式中：${\lambda _{\min }}({\boldsymbol{P}})$为P的最小特征值；$ \bar c $为一个不小于$ {\lambda _{\max }}({\boldsymbol{P}}) $的常数，由式(18)可得：

根据引理2可知，系统实际上是指数稳定，并且$ ||{r_i}(t)|{|^2} $的最终上界为$\dfrac{{2\varepsilon {M_i}^2}}{{{\lambda _{\min }}({\boldsymbol{P}})\alpha }}$。

利用Matlab中的LMI工具箱求解线性矩阵不等式$ \left[ {\begin{array}{*{20}{c}} \Sigma &{{\boldsymbol{P}}{{\boldsymbol{R}}_i}} \\ {{{\boldsymbol{R}}_i}^{\text{T}}{\boldsymbol{P}}}&{ - \varepsilon } \end{array}} \right] < 0 $可获取正定矩阵P、实数$ \mu $、标量$ \varepsilon $，获得观测器增益$ {\boldsymbol{L}} = \mu {{\boldsymbol{P}}^{ - 1}}{{\boldsymbol{C}}_i}^{\text{T}} $，代入式(11)完成非线性观测器设计。

将给进力和扭转分别作为轴向钻柱和扭转系统模型的输入，通过观测器式(11)获取状态变量$ {\hat x_i} $的值，即可获取钻头钻速估计，为下一步多目标优化问题提供孔内信息。

2.2   钻进参数优化方法

由于煤矿环境的复杂性和多样性，在实际煤矿钻进工程中，需要把如机械钻速、钻头磨损等目标综合考虑。传统的多目标优化问题是转化成单目标优化问题来解决，这种方法受人为影响，优化结果具有较大误差。随着智能优化算法的发展，国内外学者将这些算法应用到钻进参数优化中，本文尝试用NSGA-II算法来解决多目标优化问题。

在煤矿井下实际钻进中，钻进速度和钻头磨损对煤矿生产效率、安全和环保等方面都有着重要的影响，希望机械钻速和钻头磨损都能达到最优的指标。以最大钻速和最小钻头磨损为优化目标，建立能定量评估钻进效率同时能反映钻进参数对钻进过程影响的多目标优化模型。修正的杨格机械钻速模型^[14]是一种典型的物理模型，综合考虑了多种因素对钻速的影响：

式中：v为机械钻速，m/h；k_d为地层可钻性系数；C_p为压差影响系数；C_H为水力净化系数；p₁为钻压，MPa；p₂为门限钻压，MPa；N为转速，r/min；λ为转速指数；C₂为切削齿磨损系数；h为钻头切削齿相对磨损量。其中C_P和C_H的值在本文中设定为1，k_d、p₂、λ的值与实际钻进现场工况有关，通过回归分析得到这些值。为便于回归分析求取模型参数，将钻速模型进行简化：

式中：A₁=k_d为煤岩层系数；A₂=p₁为门限钻压；A₃=λ为转速指数。通过回归分析可得到A₁、A₂和A₃。其中，钻速模型中的1/(1+C₂h)是切削齿磨损对钻速的影响，该影响在下文的钻头磨损模型中已经体现，故在钻速模型中将这部分忽略。该模型是在杨格钻速模型基础上，综合考虑钻压、转速、钻井液性能以及水力因素对钻速的影响，并结合鲍格因等人的研究成果归纳建立^[15]。

本文中的钻头磨损主要考虑切削齿磨损，切削齿磨损的速度模型如下：

式中：D₁、D₂为钻压影响系数；Q₁、Q₂为转速影响系数；A_f为地层研磨性系数；C₁为切削齿磨损减慢系数。相关系数的具体值通过查阅相关资料在下文的实验部分给出。

钻头转速和钻压是影响钻进效率和质量的两个重要因素，同时也是钻进过程中易于控制的两个变量。以钻头转速N和钻压W为决策变量，综合考虑机械钻速和钻头磨损，建立多目标优化模型，优化问题可用下式的数学模型表示：

回归分析是用于研究自变量和因变量之间关系的一种统计学方法，利用孔口数据结合状态观测器估计的孔底信息，可以通过回归分析求得机械钻速模型简化形式(式(21))中的系数A₁、A₂和A₃。具体的步骤如下。

(1) 收集数据：收集的数据包括状态观测器观测的孔底钻速，孔口转速和钻压。

(2) 数据处理：使用均值滤波对初始数据进行预处理。

(3) 估计参数：利用非线性多元回归估计钻速模型中的参数。

(4) 检验模型：对模型的拟合程度进行判断，检验回归模型是否合适，包括检验模型的显著性和拟合程度等。

上述优化问题(式(23))提出了两个优化目标，需要利用有效的算法寻得最优解。NSGA-II算法在解决双目标优化问题时表现出良好的效果，已广泛应用于实际工程问题中，该方法适用于本文的研究。

NSGA-II主要克服了NSGA算法的一些缺点。在NSGA-II中，使用了快速非支配排序算法，针对多个优化目标无法同时达到最优，将种群根据相应的支配关系划分为不同级别的Pareto解，算法的目标是找到一组尽可能接近Pareto最优域的解，即对算法的收敛性有要求，同时要保证这组解的多样性。NSGA-II同时采用精英保留策略，父代优秀的个体被保留下来的概率大大增加。NSGA-II算法采用拥挤度比较算子作为种群中个体间的比较标准，克服了NSGA算法需要人为指定共享参数的缺陷。算法流程如下。

(1)初始化种群：随机生成一定数量的个体作为初始种群，这里的个体是一组机械钻速和钻头磨损的值。

(2)计算适应度：对于每个个体，计算其在多个目标函数下的适应度值。

(3)非支配排序：根据个体之间的支配关系，将种群中的个体划分为不同的等级，即非支配解、一级支配解、二级支配解等。

(4)计算拥挤度：对于每个等级中的个体，计算其在目标空间中的拥挤度，即其周围个体的密度。

(5)选择操作：根据非支配排序和拥挤度计算，选择一定数量的个体作为下一代种群。

(6)进化操作：通过交叉和变异等进化操作，生成新的个体，加入下一代种群中。

(7)终止条件判断：判断是否满足终止条件，如达到最大迭代次数或目标函数值达到一定精度等。

(8)输出结果：输出最终的帕累托前沿解集。

其中种群大小、迭代次数、交叉和变异的概率，根据实际问题进行设置。从得到的Pareto解中，根据现场的实际情况和钻进需求，挑选出满意解作为钻进参数的优化值，为决策者在实际工程问题中提供有效参考。

3   仿真分析与实验

为了验证上述方法的有效性和可行性，基于状态观测器所获得的数据进行仿真实验，对相关钻进参数进行优化，上述算法基于Python 3.8实现。

3.1   回归分析与优化结果

通过非线性多元回归对机械钻速模型式(18)进行拟合，所得的3个参数值为A₁=12.9392，A₂=8.4479，A₃=−0.3277。通过回归分析常用评价指标对结果进行分析见表1。

表  1  回归分析评价参数

Table  1.  Regression analysis evaluation parameters

模型	均方误差	平均绝对误差	R2
钻速模型	4.2×10−4	0.01193	0.9881

下载: 导出CSV 

| 显示表格

均方误差是观测值与真值偏差的平方和与观测次数的比值，均方误差的值越小，说明预测模型描述实验数据具有更好的精确度。平均绝对误差能更好地反映预测值误差的实际情况。R²介于0到1之间，可以理解为拟合程度的好坏，R²为0.9881，可以认为模型能解释98%的输出变化。

通过非线性多元回归确定钻速模型的系数后，利用机械钻速模型得出的钻速和状态观测器估计的孔底钻速值对比如图3所示。机械钻速模型计算的钻速和状态观测器估计的钻速相差较小，可认为非线性多元回归拟合的钻速方程效果较好，具有实用价值。

图  3  计算钻速与孔底估计钻速对比

Figure  3.  Comparison between ROP calculated by ROP model and downhole ROP estimated by state observer

下载: 全尺寸图片 幻灯片

钻头磨损方程中的相关系数可根据实际工况和所使用的钻头情况来确定^[16]。D₁=0.014 8，D₂=6.38，Q₁=1.5，Q₂ = 6.53×10⁻⁵，C₁=5，A_f = 0.002 28。优化算法中的相关参数设置根据经验来确定，种群大小为100，进化代数为120，图4—图6表示不同遗传代数解的分布情况。

图  4  未优化时初始解的分布情况

Figure  4.  Distribution of initial solution without optimization

下载: 全尺寸图片 幻灯片

图  5  循环次数为60时解的分布情况

Figure  5.  Distribution of solution under 60 cycles

下载: 全尺寸图片 幻灯片

图  6  优化完成时解的分布情况

Figure  6.  Distribution of solutions under completed optimization

下载: 全尺寸图片 幻灯片

从图4—图6可以看出，初始时，种群散乱地分布在解的空间中，当种群进化到60代时，种群个体已经分布得较为均匀，当完成最终的进化时，种群个体均匀地分布在解空间中，获得的Pareto曲线较为平滑，最终解的分布情况证实了NSGA-II算法的有效性，选取3组Pareto最优解分析对比，见表2。

表  2  3组Pareto最优解对比

Table  2.  Comparison of three sets of Pareto optimal solutions

分组	钻速/(mm·s−1)	钻头磨损速度/(%·s−1)
第1组	2.06	0.013 10
第2组	6.05	0.013 14
第3组	10.17	0.013 19

下载: 导出CSV 

| 显示表格

由表2可知，当钻速增大时，钻头磨损速度加快，可见两者是相互制约的关系，一个目标性能的改善将引起另一个目标性能的降低。在实际煤矿钻进的现场中，对钻速有一定的要求，首先考虑满足钻速在一定的范围内，然后去考虑钻头磨损速度，满意解的选取还是需要根据现场的实际工况来确定。

3.2   实测数据实验分析

为了进一步说明钻柱状态估计对优化结果的提升，利用淮南某煤矿实测数据进行了实验分析。所使用的数据为淮南顾北煤矿采集到的井下钻进数据，此时钻进深度为18.75 m，所使用的钻杆长度为0.75 m。在下文的实验中，钻头磨损模型相关系数与上文中进行仿真数据使用一致，机械钻速模型中的参数通过实测数据进行回归拟合得到。文献[17]仅利用孔口的实测数据，未使用状态观测器获取孔底信息，并利用NSGA-II算法进行仿真优化。本文方法利用状态观测器估计的孔底信息，并使用NSGA-II算法进行优化。将两种方法得到的结果同时与现场实际操作参数进行对比，并计算预计钻速，实验结果见表2和表3。由于机械钻速的增大，另一个优化目标钻头磨损速度也会增大，表2中优化后的目标值选择方案是保持钻头磨损速度变化较小，让钻速较大的增大，对比情况见表3。

表  3  实际煤矿钻进数据与优化后的数据对比

Table  3.  Comparison between actual coal mine drilling data and optimized data

实际数据与方法	钻压/MPa	转速/(r·min−1)	钻速/(mm·s−1)
实际采集数据	2.03	146	19.68
文献[17]	2.12	100	22.64
本文方法	2.25	100	26.07

下载: 导出CSV 

| 显示表格

分析表3可得出如下结论：使用文献[17]中的方法，基于实测孔口数据进行优化，优化后的钻速为22.64 mm/s，相比于实际采集数据中的钻速19.68 mm/s，提升了15.04%；基于钻柱状态估计孔底信息进行优化后，钻速为26.07 mm/s，相比于实际采集的数据钻速提升了32.47%。表明构造状态观测器估计孔底信息，能让优化结果提升更高，证实了本文所提方法的有效性。

4   结 论

a. 分析了煤矿井下钻进过程存在的优化问题，综合考虑机械钻速模型和钻头磨损模型这两个优化目标，通过NSGA-II算法求解优化问题，实现了最优钻压和转速的匹配，用实测数据验证了该方法的有效性。

b. 通过建立钻柱动力学模型和设计全维状态观测器，建立了孔口−孔底钻柱运动状态映射关系，实现了孔底不可测量数据的估计。实验结果表明，基于钻柱状态估计孔底信息进行优化后的钻速优于仅利用孔口的实测数据优化的钻速。

c. 提出的钻速模型计算钻速与状态观测器估计孔底钻速拟合程度好，钻柱状态估计对优化结果的提升显著，但还需要通过实践验证优化结果的准确性。