Positioning method for roadheaders based on fusion of LiDAR and inertial navigation
-
摘要:
煤矿掘进机精准定位是智能掘进的基础,但井下低光照、高粉尘等恶劣作业环境导致单一定位方法精度低、稳定性差。为提高掘进机在恶劣环境中的定位精度,提出了一种基于误差状态卡尔曼滤波(ESKF)的激光雷达与惯导融合的掘进机定位方法。首先,以悬挂在巷道顶部的球靶中心为巷道坐标系原点,设计基于密度的噪声鲁棒空间聚类(DBSCAN)算法和基于形状特征的球靶点云提取算法,解决传统依靠反射强度区分球靶的方法在粉尘堆积时易失效的问题,结合坐标变换方法构建雷达位置测量系统以获得融合定位基准。其次,利用惯导积分得到掘进机的位置和姿态信息。然后,基于一阶高斯马尔可夫过程进行误差状态建模,采用误差状态卡尔曼滤波算法融合雷达和惯导的输出,得到掘进机在巷道中的融合定位结果,并将融合定位结果反馈给惯导,以校正其累计误差,从而获得精准的定位结果。定位试验结果表明:在掘进机静止状态下,不同位置和姿态角下雷达定位系统的位置误差小于10 cm,惯导定位系统的位置误差小于70 cm;在掘进机运动状态下,融合系统的位置误差为5.8 cm,相比雷达系统的位置误差降低了12.1%。基于激光雷达与惯导融合的掘进机定位方法可以在复杂掘进工况中满足煤矿掘进机自动截割时的定位需求。
Abstract:Accurate positioning of roadheaders in coal mines is fundamental to intelligent tunneling. However, harsh working conditions, such as low illumination and high dust levels in underground mines, often degrade the accuracy and stability of single-source positioning methods. To improve the positioning accuracy of the roadheaders in these harsh conditions, a new positioning method based on the fusion of LiDAR and inertial navigation using error state kalman filter (ESKF) was developed. First, the center of the spherical target suspended from the tunnel roof was defined as the origin of the tunnel coordinate system. A density-based spatial clustering of applications with noise (DBSCAN) and a shape-feature-based spherical target point cloud extraction algorithm were designed to address the problem that conventional methods relying on reflection intensity for distinguishing spherical targets fail in environments with dust accumulation. The coordinate transformation method is then used to build a radar position measurement system to obtain a reference for the fusion positioning. Next, position and attitude information of the roadheader were obtained through inertial navigation integration. Subsequently, an error-state model was formulated based on a first-order Gaussian-Markov process, and the ESKF algorithm was applied to fuse the outputs of LiDAR and the inertial navigation, providing the fusion positioning results of the roadheader within the tunnel. The fusion positioning results were then fed back into the inertial navigation to correct accumulated errors, achieving precise positioning. Experimental results demonstrated that, under static conditions, the position error of the LiDAR-based positioning system remained below 10 cm across different positions and attitude angles, and the inertial navigation system exhibited a position error of less than 70 cm. In dynamic conditions, the fusion positioning system achieved a position error of 5.8 cm, reducing the LiDAR system's error by 12.1%. The proposed LiDAR and inertial navigation fusion-based roadheader positioning method meets the positioning requirements for automated cutting operations of roadheaders in complex tunneling conditions.
-
0. 引言
煤矿巷道掘进是煤炭开采的重要环节,然而受复杂的地质条件和人工主导的掘进工艺等因素的限制,巷道掘进效率低,精度不足[1]。因此无人化、远程化和机器人化的智能掘进技术是技术发展的必然趋势。掘进机的精准定位是实现智能掘进的基础。
根据传感信息源的不同,掘进机定位技术主要包括惯导[2-3]、全站仪[4-5]、超宽带(Ultra−Wideband, UWB)无线通信[6-7]、视觉[8-9]和激光定位[10-11]等。对于低光照、强干扰等复杂的巷道工况,基于单传感信息的定位技术精度难以达到要求,且定位结果鲁棒性低。为提高定位准确性,融合多传感信息的掘进机定位方法成为当前的研究热点。马宏伟等[4]基于卡尔曼滤波,融合光纤惯导和全站仪信息,在模拟巷道中实现了高精度定位,全站仪分支定位精度高,但全站仪成本高昂、依赖人工操作,难以实现自动定位。史苏阳[6]借助井下巷道特征提取掘进机运动惯性规律,提出了基于卡尔曼滤波的惯性测量单元(Inertial Measurement Unit, IMU)和UWB融合定位算法。然而,巷道壁表面会使UWB信号发生多次反射,造成多径效应,影响定位系统的准确性和稳定性。Lü Hongbo等[7]建立了基于牛顿迭代的UWB位姿解算算法,提出了一种基于机器学习的融合定位方法,利用IMU积分获得的高精度姿态角补偿UWB位姿解算获得的姿态角,提高了掘进机姿态角测量的精度,UWB分支利用超短脉冲信号的时间差获取掘进机位置,对信号的多路径干扰具有较强的抑制能力,但在煤矿井下,电磁干扰及UWB信号衰减和穿透力的限制导致该类方法不适用于长距离定位。毛清华等[8]利用相机采集固定在巷道顶部的激光标靶图像,通过PnP(Perspective−n−Point)算法解算出掘进机位置,结合惯导获得更准确的位姿。然而,井下的煤尘、烟雾和液体等会覆盖并污染相机镜头,造成图像模糊,影响定位精度。Wan Jicheng等[9]建立了基于三激光标靶的掘进机位姿视觉测量模型,同时通过视觉信息校正了惯导的累计误差,视觉分支依赖相机捕捉的环境特征进行位置估计,可以提供丰富的巷道信息且无定位距离限制,但容易受到巷道中的低亮度和光照变化等环境因素的影响。孙凌飞等[10]通过激光雷达扫描掘进机后的反光贴获得行进距离,结合光纤惯导和激光感知系统得到横向偏移信息,但该方法无法估计斜坡上的掘进机高度变化。刘亚超[11]基于同步定位与建图(Simultaneous Localization and Mapping, SLAM)技术配准巷道点云,利用激光雷达输出的高精度点云作为观测信息与惯导融合进行掘进机定位,然而,巷道环境单一、纹理相似,缺乏结构元素和显著的点云特征,导致难以区分不同区域,位姿估计容易偏移。
结合惯导与其他传感器的方法不仅能校正惯导的累计误差并提高定位精度,也能在其他传感器被遮挡或干扰时实现稳定定位。其中,基于激光雷达与惯导的融合定位方法不受定位距离和光照条件的限制,但不显著的巷道特征仍制约着该技术的进步。因此,笔者借助特征明显的球靶,提出了一种基于激光雷达与惯导融合的掘进机定位方法。在巷道悬挂球靶并利用雷达扫描巷道,提取球靶点云以实现基于雷达的掘进机定位,作为融合定位基准,通过对惯导输出的加速度和角速度积分获得基于惯导的机身位姿,使用误差状态卡尔曼滤波(Error State Kalman Filter, ESKF)算法融合雷达和惯导的输出得到掘进机在巷道中的精准定位结果,为实现截割过程中自动移机等操作奠定基础。
1. 基于激光雷达与惯导融合的掘进机定位方法
1.1 掘进机融合定位系统组成与坐标系设置
掘进机融合定位系统包括悬臂式掘进机、光纤惯导、球靶、3个激光雷达、嵌入式处理器,如图1所示。光纤惯导安装在机身的惯导箱中,与机体平行;球靶位于机体的后方,通过钢丝绳悬挂于巷道顶部中线处,在球靶上喷涂材料为铝酸锶铕的反光漆,以使其点云具有更高的反射强度;3个激光雷达分别通过螺栓固定在掘进机顶板的左、中、右侧,用于扫描球靶,联合3个雷达可以提供更密集的扫描线;嵌入式处理器用于处理输入的激光雷达和惯导数据。
为了准确检测掘进机的位姿,建立融合定位系统坐标系,如图2所示。
1) 雷达坐标系OlXlYlZl。以雷达光心为坐标原点Ol,Yl轴与雷达尾部线槽方向平行并指向前方,Zl轴与雷达中轴线平行并指向上方,Xl轴与Yl,Zl轴构成右手正交系并指向右方。
2) 掘进机坐标系 O2X2Y2Z2。以机身主体架的最小外接长方体的中心为坐标原点O2,X2轴垂直于机身侧壁并指向右方,Y2轴与侧壁平行并指向掘进机前进方向,Z2轴与X2,Y2轴构成右手正交系并指向上方。
3) 惯导坐标系O3X3Y3Z3。以惯导底面中心为坐标原点O3,X3,Y3,Z3方向与掘进机坐标系3个对应坐标轴方向一致。
4) 导航坐标系(东北天坐标系)O4X4Y4Z4[12]。以掘进机重心为坐标原点O4,X4轴指向地理北方,Z4轴垂直于地面向上,Y4轴与X4,Z4轴符合右手正交系规则,指向地理东方,形成东北天坐标系。
5) 地心地固坐标系O5X5Y5Z5[13]。以地球质心为坐标原点O5,X5轴在赤道平面内,位于格林尼治子午圈中,Z5轴与地球旋转方向一致,Y5轴与X5,Z5轴符合右手正交系规则。
6) 巷道坐标系O6X6Y6Z6。以球靶中心为坐标原点O6,X6轴垂直于巷道侧壁并指向右方,Y6轴与巷道走向平行并指向截割面,Z6轴与X6,Y6轴构成右手正交系并指向上方。
1.2 基于激光雷达的掘进机定位方法
在球靶点云提取的基础上,通过求取球靶中心坐标并进行坐标变换,以确定掘进机的位置。颜色、反射强度和形状特征通常用来提取目标点云[14-15]。然而,激光雷达点云缺乏颜色特征,且久置在巷道中的球靶上粉尘堆积严重,导致球靶反射强度下降,无法根据反射强度区分球靶和其他巷道点云。因此,为了提高分割效率,提出基于形状特征的球靶点云提取方法,首先将雷达坐标系中的点云变换到以雷达原点为坐标原点、巷道坐标系的方向为方向的坐标系下,设定距离阈值,滤除占比最高的巷道围岩,然后在各扫描线上结合基于密度的噪声鲁棒空间聚类(Density-Based Spatial Clustering of Applications with Noise,DBSCAN)算法和形状特征滤除曲率远小于球靶的点云,最后保留与球靶曲率最接近的点云作为球靶并获得球靶中心坐标,进而得到掘进机在巷道坐标系下的位置。
雷达坐标系下的巷道点云坐标范围随着掘进机姿态的变化而显著变化,将其转换到巷道坐标系下可有效减小坐标范围的变化程度,有助于滤除巷道围岩的点云。使用3个激光雷达扫描巷道点云,借助雷达标定[16]、惯导输出和测绘提供的雷达坐标系、掘进机坐标系、惯导坐标系、东北天坐标系和巷道坐标系间的旋转姿态角,将巷道点云中第$ k' $个点($ k' = 1,2, \cdots ,n'' $,$ n'' $为巷道点云中点的总数)在雷达坐标系中的坐标$\left( {x_{{\mathrm{pOl}}}^{k^\prime },y_{{\mathrm{pOl}}}^{k^\prime },{\textit{z}}_{{\mathrm{pOl}}}^{k^\prime }} \right)$转换到以雷达原点为坐标原点、巷道坐标系的方向为方向的坐标系,得到坐标$\left( {x_{{\mathrm{pOlt}}}^{k^\prime },y_{{\mathrm{pOlt}}}^{k^\prime },{\textit{z}}_{{\mathrm{pOlt}}}^{k^\prime }} \right)$:
$$ \left[\begin{array}{l}x_{\mathrm{pOlt}}^{k'} \\ y_{\mathrm{pOlt}}^{k'} \\ \text{z}_{\mathrm{pOlt}}^{k'}\end{array} \right]=\boldsymbol{R}_4^6\boldsymbol{R}_4^{3\mathrm{T}}\boldsymbol{R}_2^3\boldsymbol{R}_2^{1\mathrm{T}}\left[ \begin{array}{l}x_{\mathrm{pOl}}^{k'} \\ y_{\mathrm{pOl}}^{k'} \\ \text{z}_{\mathrm{pOl}}^{k'}\end{array}\right] $$ (1) 式中$ {\boldsymbol{R}}_4^6,{\boldsymbol{R}}_4^3,{\boldsymbol{R}}_2^3,{\boldsymbol{R}}_2^1 $分别为东北天坐标系旋转到巷道坐标系的旋转矩阵、东北天坐标系旋转到惯导坐标系的旋转矩阵、掘进机坐标系旋转到惯导坐标系的旋转矩阵、掘进机坐标系旋转到雷达坐标系的旋转矩阵。
坐标系间的旋转矩阵[17]为
$$ \begin{split} \boldsymbol{R}_V^U= & \left[ \begin{array}{*{20}{c}}\text{cos}\text{ }\beta_V^U\text{ } & 0 & -\text{ }\mathrm{sin}\text{ }\beta_V^U\text{ } \\ 0 & 1 & 0 \\ \text{sin}\text{ }\beta_V^U\text{ } & 0 & \text{cos}\text{ }\beta_V^U\text{ }\end{array} \right]\left[ \begin{array}{*{20}{c}}1 & 0 & 0 \\ 0 & \text{cos}\text{ }\alpha_V^U\text{ } & \text{sin}\text{ }\alpha_V^U\text{ } \\ 0 & -\text{ }\mathrm{sin}\text{ }\alpha_V^U\text{ } & \text{cos}\text{ }\alpha_V^U\text{ }\end{array} \right]\times \\ &\left[ \begin{array}{*{20}{c}}\cos\gamma_V^U & \text{sin}\text{ }\gamma_V^U & 0 \\ -\text{ }\mathrm{sin}\text{ }\gamma_V^U\text{ } & \text{cos}\text{ }\gamma_V^U\text{ } & 0 \\ 0 & 0 & 1\end{array} \right] \end{split} $$ (2) 式中${\boldsymbol{R}}_V^U,\alpha _V^U,\beta _V^U,\gamma _V^U$分别为坐标系V旋转到坐标系U的旋转矩阵、俯仰角、横滚角和偏航角。
在X轴和Z轴方向上设置过滤距离阈值(0.3 m),以滤除巷道侧壁、巷道顶部和巷道底部的点云,将其变换回雷达坐标系,分别得到3个雷达的过滤点云。
在滤除巷道围岩的基础上,进一步滤除曲率远小于球靶各截面曲率的点云。通过基于密度的DBSCAN算法[18]对各雷达每条扫描线的点云进行聚类以滤除孤立点,同时利用奇异值分解对各类进行平面拟合,投影点云到拟合的平面以消除雷达扫描误差引入的点云波动对后续曲率计算的影响。通过将各类点云拟合成圆弧来计算点云曲率,进而滤除小曲率点云类。
拟合圆弧时,点云、圆弧中心坐标(xpcOl, ypcOl, zpcOl)和半径r满足以下关系:
$$ {\left( {x_{{\mathrm{pOl}}}^i - {x_{{\mathrm{pcOl}}}}} \right)^2} + {\left( {y_{{\mathrm{pOl}}}^i - {y_{{\mathrm{pcOl}}}}} \right)^2} + {\left( {x_{{\mathrm{pOl}}}^i - {{\textit{z}}_{{\mathrm{pcOl}}}}} \right)^2} = {r^2} $$ (3) 式中$\left( {x_{{\mathrm{pOl}}}^i,y_{{\mathrm{pOl}}}^i,{\textit{z}}_{{\mathrm{pOl}}}^i} \right)$为点云中第i个点的坐标,$ i = 1,2, \cdots ,n' $,$ n' $为点云中点的总数。
将式(3)展开并代入点云中各点坐标,可得矩阵形式的方程组:
$$\begin{split} &\left[ {\begin{array}{*{20}{c}} {x_{\rm{pOl}}^1}&{y_{\rm{pOl}}^1}&{{\textit{z}}_{\rm{pOl}}^1}&1 \\ {x_{\rm{pOl}}^2}&{y_{\rm{pOl}}^2}&{{\textit{z}}_{\rm{pOl}}^2}&1 \\ \vdots & \vdots &\vdots & \vdots \\ {x_{\rm{pOl}}^{n'}}&{y_{\rm{pOl}}^{n'}}&{{\textit{z}}_{\rm{pOl}}^{n'}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} A \\ B \\ C \\ D \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - (x{{_{\rm{pOl}}^1})^2} - (y{{_{\rm{pOl}}^1})^2} - ({\textit{z}}{{_{\rm{pOl}}^1})^2}} \\ { - (x{{_{\rm{pOl}}^2})^2} - (y{{_{\rm{pOl}}^2})^2} - ({\textit{z}}{{_{\rm{pOl}}^2})^2}}\\ \vdots\quad\quad\quad \vdots \quad\quad\quad \vdots\\ { - (x{{_{\rm{pOl}}^{n'}})^2} - (y{{_{\rm{pOl}}^{n'}})^2} - ({\textit{z}}{{_{\rm{pOl}}^{n'}})^2}} \end{array}} \right] \end{split} $$ (4) $$ \left\{ {\begin{array}{*{20}{l}} {A = - 2{x_{\rm{pcOl}}}} \\ {B = - 2{y_{\rm{pcOl}}}} \\ {C = - 2{{\textit{z}}_{\rm{pcOl}}}} \\ {D = x_{\rm{pcOl}}^2 + y_{\rm{pcOl}}^2 + {\textit{z}}_{\rm{pcOl}}^2 - {r^2}} \end{array}} \right. $$ (5) 式中A,B,C,D分别为与圆弧中心坐标xpcOl, ypcOl, zpcOl和半径r相关的中间变量。
利用最小二乘法求解式(4),可得圆弧中心坐标和半径:
$$ \left\{ {\begin{array}{*{20}{l}} {{x_{\rm{pcOl}}} = - {A}/{2}} \\ {{y_{\rm{pcOl}}} = - {B}/{2}} \\ {{{\textit{z}}_{\rm{pcOl}}} = - {C}/{2}} \\ {r = \sqrt {x_{\rm{pcOl}}^2 + y_{\rm{pcOl}}^2 + {\textit{z}}_{\rm{pcOl}}^2 - D} } \end{array}} \right. $$ (6) 利用标定获得雷达安装角度和安装位置[16],合并各雷达的剩余点云到掘进机坐标系。基于DBSCAN算法对合并点云进行聚类,设置聚类半径为球靶直径尺寸,将各类点云拟合成球面,球面和圆弧拟合公式相同,区别在于前者使用空间分布的点,而后者需使所有点在同一平面上,然后滤除拟合效果较差的点云。
$$ \left|\frac{1}{n}\sum\limits_{k=1}^n \left\| o_k-\text{ }c \right\| _2-r'\right| > J $$ (7) 式中:$ k $为每类点云中的点数;n为每类点云中点的总数;ok为点坐标;c为拟合的球心坐标;$ r' $为拟合的球径;J为球体阈值,值越小表示点云越接近球形。
选择拟合的球径最接近球靶半径的类作为球形标靶的点云,进而获得球靶中心在掘进机坐标系下的坐标。球靶点云提取流程如图3所示,虚线框内为各阶段球靶点云。
球靶中心是巷道坐标系原点,对该点在掘进机坐标系下的坐标(xtOm, ytOm, ztOm)进行变换,得到其在以掘进机原点为坐标原点、巷道坐标系的方向为方向的坐标系中的坐标(xtOmt, ytOmt, ztOmt):
$$ \left[\begin{array}{l}x_{\rm{tOmt}} \\ y_{\rm{tOmt}} \\ \text{z}_{\rm{tOmt}}\end{array}\right]=\boldsymbol{R}_{\mathrm{4}}^6\boldsymbol{R}_{\mathrm{4}}^{3\mathrm{T}}\boldsymbol{R}_{\mathrm{2}}^3\left[\begin{array}{l}x_{\rm{tOm}} \\ y_{\rm{tOm}} \\ \text{z}_{\rm{tOm}}\end{array}\right] $$ (8) 将式(8)中坐标变换为掘进机原点在巷道坐标系中的坐标$\left( {x_{\rm{mOt}}^{{\rm{Lidar}}{\text{ }}},y_{\rm{mOt}}^{{\rm{Lidar}}{\text{ }}},z_{\rm{mOt}}^{{\rm{Lidar}}{\text{ }}}} \right)$,即为机身在巷道中的位置:
$$\left[ \begin{array}{l} x_{\rm{mOt}}^{\rm{Lidar}}\\ y_{\rm{mOt}}^{\rm{Lidar}}\\ {\textit{z}}_{\rm{mOt}}^{\rm{Lidar}} \end{array} \right]=-\left[ \begin{array}{l} x_{\rm{tOmt}}\\ y_{\rm{tOmt}}\\ {\textit{z}}_{\rm{tOmt}} \end{array} \right] $$ (9) 1.3 基于惯导的掘进机定位方法
惯导实时输出其原点在地心地固坐标系下的经纬高坐标,结合巷道原点的经纬高坐标及惯导、东北天和掘进机坐标系间的变换关系获得机身在巷道坐标系中的位置。地心地固坐标系下的经度、纬度、高度坐标(lon, lat, alt)可以转换为笛卡尔坐标(x, y, z):
$$ \left\{\begin{array}{*{20}{l}}x=(W+a_{\mathrm{lt}})\cos\; l_{\mathrm{at}}\cos\; l_{\mathrm{on}} \\ y=(W+a_{\mathrm{lt}})\cos\; l_{\mathrm{at}}\sin\; l_{\mathrm{on}} \\ \text{z}=\left[W\left(1-e^2\right)+a_{\mathrm{lt}}\right]\sin\; l_{\mathrm{at}}\end{array}\right. $$ (10) 式中:W为地球的曲率半径;e为地球偏心率。
已知巷道原点Ot (lontOe, lattOe, alttOe)和当前位置的惯导原点 Oi (loniOe, latiOe, altiOe),根据式(10)计算点Ot (xtOe, ytOe, ztOe)和Oi (xiOe, yiOe, ziOe),则巷道原点在以惯导原点为坐标原点、东北天坐标系的方向为方向的坐标系下的坐标(xtOin, ytOin, ztOin)为
$$ \left[\begin{array}{l}x_{\rm{tOin}} \\ y_{\rm{tOin}} \\ \text{z}_{\rm{tOin}}\end{array}\right]=\boldsymbol{S}\left(\left[\begin{array}{l}x_{\rm{tOe}} \\ y_{\rm{tOe}} \\ \text{z}_{\rm{tOe}}\end{array}\right]-\left[\begin{array}{l}x_{\rm{iOe}} \\ y_{\rm{iOe}} \\ \text{z}_{\rm{iOe}}\end{array}\right]\right) $$ (11) 式中S为从地心地固笛卡尔坐标系到东北天坐标系的旋转矩阵[19]。
通过惯导原点在掘进机坐标系中的坐标 (xiOm, yiOm, ziOm),得到掘进机原点在以惯导原点为坐标原点、东北天坐标系的方向为方向的坐标系下的坐标 (xmOin, ymOin, zmOin):
$$ \left[\begin{array}{l}x_{\rm{mOin}} \\ y_{\rm{mOin}} \\ \text{z}_{\rm{mOin}}\end{array}\right]=-\boldsymbol{R}_{\mathrm{4}}^{3\mathrm{T}}\boldsymbol{R}_{\mathrm{2}}^3\left[\begin{array}{l}x_{\rm{iOm}} \\ y_{\rm{iOm}} \\ \text{z}_{\rm{iOm}}\end{array}\right] $$ (12) 由式(11)和式(12)可得掘进机原点在以巷道原点为坐标原点、东北天坐标系的方向为方向的坐标系下的坐标 (xmOtn, ymOtn, zmOin):
$$\left[ \begin{array}{l} x_{\rm{mOtn}}\\ y_{\rm{mOtn}}\\ {\textit{z}}_{\rm{mOtn}} \end{array} \right]=\left[ \begin{array}{l} x_{\rm{mOin}}\\ y_{\rm{mOin}}\\ {\textit{z}}_{\rm{mOin}} \end{array} \right]-\left[ \begin{array}{l} x_{\rm{tOin}}\\ y_{\rm{tOin}}\\ {\textit{z}}_{\rm{tOin}} \end{array} \right]$$ (13) 以巷道原点为坐标原点、东北天坐标系的方向为方向的坐标系与巷道坐标系原点重合,结合测绘部门给出的东北天坐标系与巷道坐标系间的姿态角,可得掘进机原点在巷道坐标系下的坐标$\left( {x_{\rm{mOt}}^{{\text{Ins}}},y_{\rm{mOt}}^{{\text{Ins}}},{\textit{z}}_{\rm{mOt}}^{{\text{Ins}}}} \right)$:
$$ \left[\begin{array}{*{20}{l}}x_{\rm{mOt}}^{\rm{Ins}} \\ y_{\rm{mOt}}^{\rm{Ins}} \\ \text{z}_{\rm{mOt}}^{\rm{Ins}}\end{array}\right]=\boldsymbol{R}_{\mathrm{4}}^3\left[\begin{array}{*{20}{l}}x_{\rm{mOtn}} \\ y_{\rm{mOtn}} \\ \text{z}_{\rm{mOtn}}\end{array}\right] $$ (14) 1.4 基于激光雷达与惯导融合的掘进机定位方法
惯导和激光雷达提供2种互补的位置信息,惯导不受外部环境影响,能够实时输出定位数据,但误差累计影响定位精度;雷达虽然处理速度较慢,但可以为惯导提供精确的校正信息。将掘进机状态误差建模为一阶高斯马尔可夫过程,基于ESKF算法建立融合激光雷达和惯导信息的掘进机位置测量模型,通过状态方程和观测方程描述非线性系统的动态特性和测量信息,利用雷达数据对惯导误差连续校正,得到更准确和鲁棒的掘进机状态估计。
基于惯导的加速度计误差进行掘进机误差建模,位置误差δp随时间的变化率${\dot {\boldsymbol{\delta}} _{\mathrm{p}}}$等于速度误差δv,速度误差随时间的变化率$ {\dot {\boldsymbol{\delta }}_v} $等价于加速度误差δa和其常值零位误差δa0:
$$ {\dot {\boldsymbol{\delta}} _{\mathrm{p}}} = {{\boldsymbol{\delta}} _{\mathrm{v}}} $$ (15) $$ {\dot {\boldsymbol{\delta }}_{\mathrm{v}}} = {{\boldsymbol{\delta }}_{\mathrm{a}}} + {{\boldsymbol{\delta }}_{{\mathrm{a}}0}} $$ (16) 常值零位误差随时间的变化率$ {\dot \delta _{{\mathrm{a}}0}} = 0 $。在一阶高斯马尔科夫过程中,系统下一时刻的状态只与当前时刻有关,而且该随机过程具有负反馈机制,不会发散。基于一阶高斯马尔可夫过程近似得到掘进机的加速度误差:
$$ {\dot {\boldsymbol{\delta }}_{\mathrm{a}}} = \tau {{\boldsymbol{\delta }}_{\mathrm{a}}} + \sigma {\boldsymbol{w}}(t) $$ (17) 式中:τ为时间相关系数;σ为随机系数;w(t)~N(0,1)为标准正态分布;t为时间变量。
根据掘进机的位置p、速度v、加速度a在XYZ 3个维度上的误差,以及加速度常值零位误差建立误差状态矢量X = [δp δv δa δa0],利用一阶高斯马尔可夫过程构建系统方程:
$$ \dot {\boldsymbol{X}}={\boldsymbol{\varPhi X}}+{\boldsymbol{G}} $$ (18) 式中:Φ为动态矩阵;G为驱动项。
以雷达输出位置坐标的时间间隔为数据处理的最小时间单元M,得
$$ {{\boldsymbol{X}}_t} = {{\boldsymbol{X}}_{t - 1}} + M{\dot {\boldsymbol{X}}_{t - 1}} $$ (19) 式中Xt为t时刻的状态向量。
将式(18)代入式(19)构建状态方程:
$$ {{\boldsymbol{X}}_t} = {\boldsymbol{F}}{{\boldsymbol{X}}_{t - 1}} + {\boldsymbol{\varepsilon}} $$ (20) 式中:F为状态转移矩阵,F = MΦ + I,I为单位矩阵;ε为过程噪声,ε = MG。
过程噪声ε的协方差Q为
$$ {\boldsymbol{Q}} = {\mathrm{E}}\left( {\left( {{\boldsymbol{\varepsilon}} - {{0}}} \right){{\left( {{\boldsymbol{\varepsilon}} - {{0}}} \right)}^{\mathrm{T}}}} \right) = {M^2}{\mathrm{E}}({{{\boldsymbol{G}}}}{{{{\boldsymbol{G}}}}^{\mathrm{T}}}) $$ (21) 式中E为求矩阵的期望运算。
根据雷达和惯导输出的位置结果,得到观测值Z为
$$ {\boldsymbol{Z}} = \left[ {\begin{array}{*{20}{c}} {x_{\rm{mOt}}^{\rm{Lidar}} - x_{\rm{mOt}}^{\rm{Ins}}} \\ {y_{\rm{mOt}}^{\rm{Lidar}} - y_{\rm{mOt}}^{\rm{Ins}}} \\ {{\textit{z}}_{\rm{mOt}}^{\rm{Lidar}} - {\textit{z}}_{\rm{mOt}}^{\rm{Ins}}} \end{array}} \right] + {\boldsymbol{\eta}} $$ (22) 式中η为雷达的测量噪声。
用误差状态矢量X表示观测值Z,得到观测方程:
$$ {{\boldsymbol{Z}}_t} = {\boldsymbol{HX}}_t +{\boldsymbol{\eta}} $$ (23) 式中H为观测转移矩阵。
在对掘进机状态误差建模的基础上,基于ESKF算法融合激光雷达和惯导测量的掘进机位置。ESKF算法[20]是一种递归算法,通过状态反馈来估计系统的状态。在状态方程的基础上,利用上一时刻的误差状态修正值$ {\boldsymbol{X}}_{t - 1}^ + $和状态协方差修正值$ {\boldsymbol{P}}_{t - 1}^ + $分别预测当前时刻的误差状态$ {\boldsymbol{X}}_t^ - $和状态的协方差$ {\boldsymbol{P}}_t^ - $:
$$ {\boldsymbol{X}}_t^ - = {\boldsymbol{FX}}_{t - 1}^ + + {\boldsymbol{\varepsilon}} $$ (24) $$ {\boldsymbol{P}}_t^ - = {\boldsymbol{FP}}_{t - 1}^ + {{\boldsymbol{F}}^{\mathrm{T}}} + {{{\boldsymbol{Q}}}} $$ (25) 利用以上预测值和观测方程,可得卡尔曼增益Kt、误差状态修正值$ {\boldsymbol{X}}_t^ + $和状态协方差修正值$ {\boldsymbol{P}}_t^ + $:
$$ \boldsymbol{K}_t=\boldsymbol{P}_t^-\boldsymbol{H}^{\mathrm{T}}\left(\boldsymbol{HP}_t^-\boldsymbol{H}^{\mathrm{T}}+\boldsymbol{R}\right)^{-1} $$ (26) $$ {\boldsymbol{X}}_t^ + = {{\boldsymbol{K}}_t}{{\boldsymbol{Z}}_t} - {{\boldsymbol{K}}_t}{\boldsymbol{HX}}_t^ - + {\boldsymbol{X}}_t^ - $$ (27) $$ {\boldsymbol{P}}_t^ + = {\boldsymbol{P}}_t^ - - {{\boldsymbol{K}}_t}{\boldsymbol{HP}}_t^ - $$ (28) 将掘进机误差状态修正值与惯导输出的掘进机位置数据相加,得到融合的掘进机位置结果。将该结果反馈给惯导可以遏制误差累计,Fuse为输出的融合定位结果, Ins为惯导输出的定位结果,Correct为修正后的定位误差。
$$ {\left[ {\begin{array}{*{20}{c}} {\boldsymbol{p}} \\ {\boldsymbol{v}} \\ {\boldsymbol{a}} \end{array}} \right]_{{\mathrm{Fuse}}}} = {\left[ {\begin{array}{*{20}{c}} {\boldsymbol{p}} \\ {\boldsymbol{v}} \\ {\boldsymbol{a}} \end{array}} \right]_{{\mathrm{Ins}}}} + {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\delta }}_p}} \\ {{{\boldsymbol{\delta }}_v}} \\ {{{\boldsymbol{\delta }}_a}} \end{array}} \right]_{{\mathrm{Correct}}}} + {\left[ {\begin{array}{*{20}{c}} 0 \\ 0 \\ {{{\boldsymbol{\delta }}_{a0}}} \end{array}} \right]_{{\mathrm{Correct}}}} $$ (29) 1.5 基于激光雷达与惯导融合的掘进机定位流程
基于激光雷达和惯导融合的掘进机定位流程包括雷达定位、惯导定位、融合定位,如图4所示。首先,在巷道中悬挂3D球靶,使用激光雷达采集巷道点云,利用DBSCAN算法和形状特征分割巷道点云,得到掘进机坐标系下3个雷达在球靶上的联合点云,进而确定球靶中心在掘进机坐标系下的坐标;选择球靶中心为巷道坐标系原点并对其进行坐标变换,得到基于雷达的机体原点在巷道坐标系下的坐标,作为融合定位系统的位置基准。然后,在初始经纬高坐标的基础上,惯导通过对陀螺仪和加速度计分别采集的角速度和加速度积分,输出自身在当前位置的经纬高坐标,综合当前位置和球靶处的经纬高变换坐标,获得基于惯导的掘进机原点在巷道坐标系的坐标。最后,基于ESKF算法融合激光雷达和惯导数据,得到机身原点在巷道坐标系下的精准坐标。机身在巷道坐标系下的姿态通过分析惯导输出的姿态角和煤矿测绘提供的巷道方位角确定,此外,惯导在运行过程中不定时进行自对准,以保持对姿态角的高精度测量。
2. 基于激光雷达与惯导融合的掘进机定位试验
2.1 掘进机定位试验设计
为了验证本文提出的掘进机融合定位系统的有效性,搭建悬臂式掘进机位姿测量试验平台,包括模拟巷道、EBZ200S悬臂式掘进机、3个16线的GUL50矿用本安型激光雷达、SIRIUS−F002光纤捷联惯导和1个球靶。设计球靶直径为30 cm,使球靶上至少有2条激光线,确保球靶中心坐标计算的准确性,使用在煤矿井下久置后粉尘堆积的球靶,放置在掘进机后方的巷道中线处,下方悬挂铅锤,用于确定其中心在地面的投影点。模拟巷道在东北天坐标系的偏航角、俯仰角和翻滚角分别为270.13°, 0.02°和0.25°,内部为低光环境,底部为平坦的厚钢板,巷道内划分边长为1 m的正方形网格,网格与巷道平行。为了便于进行真实位置坐标测量和误差计算,在机体尾部选择第二掘进机原点,在该点悬挂铅锤以确定其在地面的投影点,借助该点和球靶中心的投影点及网格确定第二掘进机原点在巷道坐标系下的X轴和Y轴真实坐标,利用球靶和第二掘进机原点下的铅垂线长度获得其Z轴真实坐标,通过比较真实坐标和定位系统输出的第二掘进机原点坐标得到位置误差。
2.2 基于激光雷达的掘进机定位试验
保持掘进机静止并在模拟巷道中移动球靶,以球靶和掘进机间的相对运动模拟机体在巷道中移动。将球靶移动到巷道不同位置进行定点测试,球靶相对第二掘进机原点的距离分布在3.37~7.96 m。考虑掘进机在工作过程中通常会出现旋转和俯仰运动,保持球靶静止,通过水平偏转机身方向、伸缩铲板和后支撑等依次改变掘进机偏航角和俯仰角,同时使其余姿态角与巷道方位角近似相等,偏航角和俯仰角分别分布在−14.6~15.42°和−5.89~6.82°。使用雷达扫描巷道并提取球靶点云,对球靶中心坐标进行变换,获得第二掘进机原点在巷道坐标系中的位置,对比第二掘进机原点在巷道坐标系中的真实位置坐标,获得定位误差。掘进机在不同距离、偏航角和俯仰角下的定位误差在10 cm以下,表明雷达定位具有较高的精度,可以作为融合定位系统的位置基准。
2.3 基于惯导的掘进机定位试验
将掘进机移动到巷道中不同位置进行定点测试,第二掘进机原点相对球靶的距离分布在2.69~8.62 m。首先进行惯导对准,然后分别改变掘进机偏航角和俯仰角,保持其余姿态角与巷道方位角近似相等,偏航角和俯仰角分别分布在−10.31~12.26°和−4.03~5.83°。对比各个位置下惯导输出的第二掘进机原点在巷道坐标系中的位置坐标与真实位置坐标,发现掘进机在不同距离、偏航角和俯仰角下的定位误差在70 cm以下,惯导输出的掘进机位置精度低,无法直接用于掘进机定位。
2.4 基于激光雷达与惯导融合的掘进机定位试验
在融合试验开始前对惯导进行对准,然后控制掘进机在巷道内以S形前进,分别通过雷达和惯导同步测量第二掘进机原点的位置坐标,基于ESKF算法融合雷达和惯导数据获得位置坐标。掘进机沿Y轴方向前进约8.5 m。将一个匀速滴水的装置固定在第二掘进机原点处,通过测量地面上水滴的坐标形成机体真实运动轨迹。掘进机运动轨迹及雷达和融合方法的定位轨迹如图5所示。可看出基于雷达的定位结果和融合定位结果与真实坐标变化趋势一致,而融合定位方法的曲线与真实位置曲线更接近,位置误差更小。
在融合试验中,利用平均Hausdorff距离来度量估计轨迹和真实轨迹间的距离。若$L $和${{L'}} $分别为估计轨迹和真实轨迹的有限点集,则其平均Hausdorff距离$ {d_{\mathrm{H}}}({{L}},{{L'}}) $为
$$ {d_{\mathrm{H}}}({{L}},{{L'}}) = \max \left( {\mathop {{{\mathrm{mean}}} }\limits_{f \in {\boldsymbol{L}}} \mathop {\min }\limits_{b \in {{L'}}} d(f,b),\mathop {{{\mathrm{mean}}} }\limits_{b \in {{L'}}} \mathop {\min }\limits_{f \in {\boldsymbol{L}}} d(f,b)} \right) $$ (30) 式中d(f,b)为点f和点b的欧氏距离。
将三维空间的估计轨迹与真实轨迹向Y−X平面、Y−Z平面投影,并根据式(30)得到Y−X平面、Y−Z平面估计轨迹与真实轨迹间的平均Hausdorff距离,见表1。可看出与雷达定位相比,融合定位方法在Y−X平面、Y−Z平面和三维空间的估计轨迹与真实轨迹间的平均Hausdorff距离分别降低了17.4%,7.5%和12.1%,融合定位方法抑制了惯导定位的误差累计,提升了系统定位精度,可以满足掘进机定位精度要求。
表 1 基于激光雷达与惯导融合的掘进机定位试验的平均Hausdorff距离Table 1. Mean Hausdorff distance in roadheader positioning experiments based on LiDAR and inertial navigation fusion轨迹空间 平均 Hausdorff 距离/cm 雷达定位 融合定位 Y−X平面 4.6 3.8 Y−Z平面 4.0 3.7 三维空间 6.6 5.8 3. 结论
1) 为提高掘进机定位精度,提出了基于误差状态卡尔曼滤波的激光雷达和惯导融合的定位方法。结合球靶点云提取算法和坐标变换获得基于激光雷达的掘进机位置测量结果,作为融合定位系统的定位基准。对惯导输出的加速度和角速度积分,得到基于惯导的掘进机定位结果。采用一阶高斯马尔可夫过程进行误差建模,在ESKF迭代过程中将融合定位结果反馈给惯导,修正了惯导自身的累计误差,提升了融合定位系统的精度。
2) 在EBZ200S掘进机上部署融合定位系统,基于雷达测量球靶与掘进机在不同距离及掘进机在不同姿态角下机身的位置,位置误差小于10 cm,惯导定位误差小于70 cm。与雷达定位相比,融合定位方法在Y−X平面、Y−Z平面和三维空间的估计轨迹与真实轨迹间的平均Hausdorff距离分别降低了17.4%,7.5%和12.1%,融合定位方法满足低光照、高粉尘等复杂工况下巷道内精准定位的需要。
-
表 1 基于激光雷达与惯导融合的掘进机定位试验的平均Hausdorff距离
Table 1 Mean Hausdorff distance in roadheader positioning experiments based on LiDAR and inertial navigation fusion
轨迹空间 平均 Hausdorff 距离/cm 雷达定位 融合定位 Y−X平面 4.6 3.8 Y−Z平面 4.0 3.7 三维空间 6.6 5.8 -
[1] 谢和平,王金华,鞠杨,等. 煤炭革命的战略与方向[M]. 北京:科学出版社,2018. XIE Heping,WANG Jinhua,JU Yang,et al. Coal industry reform:strategies and directions[M]. Beijing:Science Press,2018.
[2] SHEN Yang,WANG Pengjiang,ZHENG Weixiong,et al. Error compensation of strapdown inertial navigation system for the boom-type roadheader under complex vibration[J]. Axioms,2021,10(3) . DOI: 10.3390/axioms10030224.
[3] 马宏伟,王鹏,张旭辉,等. 煤矿巷道智能掘进机器人系统关键技术研究[J]. 西安科技大学学报,2020,40(5):751-759. MA Hongwei,WANG Peng,ZHANG Xuhui,et al. Research on key technology of intelligent tunneling robotic system in coal mine[J]. Journal of Xi'an University of Science and Technology,2020,40(5):751-759.
[4] 马宏伟,毛金根,毛清华,等. 基于惯导/全站仪组合的掘进机自主定位定向方法[J]. 煤炭科学技术,2022,50(8):189-195. MA Hongwei,MAO Jingen,MAO Qinghua,et al. Automatic positioning and orientation method of roadheader based on combination of ins and digital total station[J]. Coal Science and Technology,2022,50(8):189-195.
[5] 张旭辉,刘博兴,张超,等. 掘进机全站仪与捷联惯导组合定位方法[J]. 工矿自动化,2020,46(9):1-7. ZHANG Xuhui,LIU Boxing,ZHANG Chao,et al. Roadheader positioning method combining total station and strapdown inertial navigation system[J]. Industry and Mine Automation,2020,46(9):1-7.
[6] 史苏阳. 煤矿复杂环境移动目标精确定位技术研究[D]. 徐州:中国矿业大学,2023. SHI Suyang. Research on accurate positioning technology of moving target in complex environment of coal mine[D]. Xuzhou:China University of Mining and Technology,2023.
[7] LYU Hongbo,ZHENG Xinyue,QI Yuhao,et al. UWB-IMU pose estimation for roadheader based on machine learning[C]. International Conference on Intelligent Control,Measurement and Signal Processing,Chengdu,2023:1153-1156.
[8] 毛清华,周庆,安炎基,等. 惯导与视觉信息融合的掘进机精确定位方法[J]. 煤炭科学技术,2024,52(5):236-248. DOI: 10.12438/cst.2023-1003 MAO Qinghua,ZHOU Qing,AN Yanji,et al. Precise positioning method of tunneling machine for inertial navigation and visual information fusion[J]. Coal Science and Technology,2024,52(5):236-248. DOI: 10.12438/cst.2023-1003
[9] WAN Jicheng,ZHANG Xuhui,ZHANG Chao,et al. Vision and inertial navigation combined-based pose measurement method of cantilever roadheader[J]. Sustainability,2023,15(5). DOI: 10.3390/SU15054018.
[10] 孙凌飞,刘亚,彭继国,等. 基于惯性技术的掘进机组合定位方法[J]. 煤炭科学技术,2024,52(12):300-310. DOI: 10.12438/cst.2023-1648 SUN Lingfei,LIU Ya,PENG Jiguo,et al. Integrated positioning method of roadheader based on inertial technology[J]. Coal Science and Technology,2024,52(12):300-310. DOI: 10.12438/cst.2023-1648
[11] 刘亚超. 面向掘进机的三维激光同步定位与建图[D]. 北京:北方工业大学,2024. LIU Yachao. Three-dimensional laser synchronous positioning and mapping for tunneling machines[D]. Beijing:North China University of Technology,2024.
[12] 李朕阳,邹鹏,刘振海,等. 同平台双偏振仪器地理定位及校正方法[J]. 光子学报,2021,50(2):150-162. LI Zhenyang,ZOU Peng,LIU Zhenhai,et al. Geolocation and correction method for dual polarization instrument on same platform[J]. Acta Photonica Sinica,2021,50(2):150-162.
[13] 司垒,王忠宾,谭超,等. 基于差分式惯性传感组件的采煤机位姿解算法[J]. 振动. 测试与诊断,2021,41(2):220-227,406. SI Lei,WANG Zhongbin,TAN Chao,et al. Position and attitude calculation algorithm of shearer based on differential inertial sensors[J]. Journal of Vibration,Measurement & Diagnosis,2021,41(2):220-227,406.
[14] 石勇. 基于三维激光雷达的掘进机实时位姿纠偏系统[J]. 煤矿机械,2023,44(5):64-66. SHI Yong. Real-time pose correction system of roadheader based on three-dimensional lidar[J]. Coal Mine Machinery,2023,44(5):64-66.
[15] QU Yuanyuan,YANG Teng,LI Tao,et al. Path tracking of underground mining boom roadheader combining BP neural network and state estimation[J]. Applied Sciences,2022,12(10). DOI: 10.3390/APP12105165.
[16] 蔡春蒙,刘京,周江涛,等. 激光雷达的标定方法、系统和可读存储介质:CN202411239096.7[P]. 2024-11-29. CAI Chunmeng,LIU Jing,ZHOU Jiangtao,et al. The calibration method,system,and readable storage medium of lidar:CN202411239096.7[P]. 2024-11-29.
[17] 郭伦锋,郭一楠,蒋康庆,等. 掘进机姿态参数测量及解算方法[J]. 工矿自动化,2021,47(12):46-54. GUO Lunfeng,GUO Yinan,JIANG Kangqing,et al. Measurement and calculation method of attitude parameters of roadheader[J]. Industry and Mine Automation,2021,47(12):46-54.
[18] HAHSLER M,PIEKENBROCK M,DORAN D. Dbscan:fast density-based clustering with R[J]. Journal of Statistical Software,2019,91. DOI: 10.18637/jss.v091.i01.
[19] SZABOVA M,DUCHON F. Survey of GNSS coordinates systems[J]. European Scientific Journal,2016,12(24):33-42.
[20] BRIGADNOV I,LUTONIN A,BOGDANOVA K. Error state extended Kalman filter localization for underground mining environments[J]. Symmetry,2023,15(2). DOI: 10.3390/SYM15020344.