Similar simulation experiment of water loss and settlement in thick loose aquifer
-
摘要: 目前缺乏对厚松散含水层地质采矿条件下覆岩破断及变形规律的深入研究。以淮南矿区潘四东煤矿11111工作面为工程背景，构建相似材料模型，采用数字摄影测量提取位移法记录模型开挖过程中覆岩破断过程及覆岩变形情况。分析了含水层失水沉降原因：覆岩在W型剪切应力拱作用下形成2条纵向的主导水裂隙带，导水裂隙带的进一步发育引起含水层失水固结，在厚松散层重力作用下进一步压实，随着覆岩破断运动的加剧，在弯曲带和覆岩共同挤压下形成О型剪切应力拱，压缩薄层空间，导致地表下沉量增大。分析了失水状态下覆岩损伤情况：工作面开采工作完成且覆岩达到稳态后，前垮落角为57°，后垮落角为62°，导水裂隙带高度为63 m，开切眼及终采线上方覆岩在应力集中作用下断裂，产生纵向裂隙，开切眼及终采线上方垮落带区域内覆岩产生横向离层裂隙，纵向裂隙和横向离层裂隙加剧了覆岩与含水层间的水力联系。给出了失水状态下覆岩动态运动规律：随着开采工作面的推进，各观测线覆岩下沉量逐渐增大，接近开采工作面的观测线覆岩下沉量最大，工作面上方覆岩的观测线下沉量曲线走势基本类似且跳变一致，含水层上方的观测线下沉量曲线走势基本吻合且跳变同步，工作面上方与含水层上方的观测线下沉量跳变异步，表明含水层对覆岩移动变形具有重要作用。Abstract: The in-depth research on the breaking and deformation law of overburden rock under the geological and mining conditions of thick loose aquifer are lacking at present. Taking 11111 working face of Pansidong Coal Mine in Huainan mining area as the engineering background, the similar material model is constructed, and the digital photogrammetry extraction displacement method is used to record the overburden rock breaking process and overburden rock deformation during the model roadway heading. The causes of water loss and settlement of aquifer are analyzed. The overburden rocks form two main longitudinal diversion fissure zones under the action of W-type shear stress arch. The further development of the diversion fissure zone causes water loss and consolidation of the aquifer, and the aquifer is further compacted under the action of gravity of the thick loose layer. With the intensification of the overburden rock breaking movement, О type shear stress arch is formed under the joint extrusion of bending zone and overburden rock, which compresses the thin space and leads to the large amount of surface subsidence. The damage of overburden rock under water loss condition is analyzed. After the roadway heading work of the working face is completed and the overburden rock reaches a steady state, the front caving angle is 57°, the rear caving angle is 62°, and the height of the diversion fissure zone is 63 m. Under the action of stress concentration, the overburden rock above the open-cut hole and the stop-mining line is broken to produce longitudinal fissure, and the overburden rock in the area of the collapse zone above the open-cut hole and the stop-mining line produces lateral separation fissure. The longitudinal fissures and lateral separation fissures intensify the hydraulic connection between overburden rock and the aquifer. The dynamic movement law of overburden rock under water loss state is given. With the advance of mining face, the overburden settlement of each observation line increases gradually, and the overburden settlement of the observation line close to the working face is the largest. The trend of the subsidence curves of the observation lines in the overburden rock above the working face is basically similar, and the jump of the subsidence curves is consistent. The trend of the subsidence curves of the observation lines above the aquifer is basically consistent, and the jump of the subsidence curves is synchronous. The jump of the subsidence curves of observation lines in the overburden rock above the working face and the one of the observation line above the aquifer are asynchronous, indicating that the aquifer plays an important role in the movement and deformation of the overburden rock.
-
0. 引言
尾矿库作为矿山安全生产的重大危险源之一，存在极大的潜在风险 [1-2]。对此，国家出台的《尾矿库安全监督管理规定》[3]和GB 39496−2020《尾矿库安全规程》[4]中要求定期组织尾矿库专项检查，对发现的事故隐患及时进行治理。尾矿库排洪隧洞主要用于泄洪和回收尾矿澄清水，其排洪能力和结构性能直接关系到矿山安全，因此必须对排洪隧洞进行安全巡检。目前排洪隧洞巡检以人工为主，因巡检线路较长，巡检过程存在安全风险，提出以探测机器人代替巡检人员完成巡检任务。探测机器人巡检需要在长度超过1 km的排洪隧洞内部精准定位隐患，所以有必要对探测机器人在巡检过程中的自主定位进行研究。
近年来，学者们对探测机器人自主定位进行了深入研究，提出了多种算法提高探测机器人实时位姿计算精度。文献[5-6]使用测程法进行定位，利用UMBmark算法修正探测机器人轮径和轴径两类系统误差对里程计定位造成的影响，提高了测程法定位精度。测程法易于实现且适用性广，但其鲁棒性和稳定性较差，且存在累计误差。文献[7]提出了一种基于概率估计的里程计模型。文献[8]使用有限记忆卡尔曼方法，以逐级滤波的方法融合多传感器信息，推算矿井机器人的航向和运动轨迹，减小了单一传感器测量误差的影响，但仍未解决累计误差问题，不适合长距离定位。文献[9]使用惯性模块（Inertial Measurement Unit, IMU）、里程计、全球导航卫星系统（Global Navigation Satellite System，GNSS）等传感器进行机器人定位，通过扩展卡尔曼滤波（Extended Kalman Filter, EKF）进行数据滤波，实现多传感器融合，提升了系统鲁棒性和稳定性，但累计误差的消除需依赖卫星提供的绝对位置信息，当卫星信号微弱时该类算法易因校正失败出现累计误差，受使用环境限制。文献[10]利用Fast−SLAM及包含激光雷达在内的多种传感器，利用基于粒子滤波的自适应蒙特卡洛算法构建环境栅格地图，实现了探测机器人在矿洞室内环境的精确定位，但Fast−SLAM算法在弱纹理场景中易出现粒子收敛缓慢情况，影响定位效果。文献[11-13]均利用图优化算法实现了探测机器人多传感器融合定位，并利用子图匹配实现回环检测，有效消除了定位过程的累计误差，但所用算法在弱纹理场景中匹配错误率较高，未考虑弱纹理场景的特殊性。文献[14-15]利用二维码阵列进行室内定位与导航，使得机器人能够在缺少纹理的室内环境实现精确定位，但均需要提前布置密集的二维码，且维护工作量极大，不适合大场景定位。文献[16]利用实验对比分析了目前常用标签的检测精度和效率，结果表明ArUco码的检测速度有明显优势，适用于对标签检测速度要求高的场合。
尾矿库排洪隧洞环境属于狭长室内环境，且隧洞内卫星信号弱，特征点少，现有探测机器人定位算法存在弱纹理室内场景定位失效的问题，难以实现此类环境中的精确定位和累计误差消除。本文提出了一种适用于尾矿库排洪隧洞环境的探测机器人多传感器融合定位算法。该算法将排洪隧洞环境的远距离定位以ArUco码为边界点，简化为多段短距离定位，利用EKF算法融合里程计和IMU信息实现短距离精确定位，利用预设的ArUco码作为路标消除EKF定位带来的累计误差，并实现位姿优化，最终实现狭长空间弱纹理场景下探测机器人的远距离高精度定位。该算法有效利用了多种传感器数据，避免了单一传感器定位误差，提升了定位鲁棒性。为验证探测机器人多传感器融合定位算法效果，设计了尾矿库排洪隧洞探测机器人系统实验平台，并在尾矿库排洪隧洞实际场景中运行和关闭多传感器融合定位算法，分别进行10组行进40 m的重复定位实验，结果证明了该算法具有较高的稳定性和精确性。
1. 多传感器融合定位算法
1.1 融合定位算法流程
探测机器人多传感器融合定位算法如图1所示。首先选取UMBmark算法进行里程计标定，有效消除探测机器人轮径和轴径两类系统误差。然后利用EKF算法融合IMU数据和里程计数据，并利用测程法实现运动过程中探测机器人位置和姿态信息计算。最后选用ArUco码作为路标并固定于隧洞内部，探测机器人携带经张正友标定法[17]标定后的相机对其进行识别与处理，利用图优化方式实现位姿优化，并根据ArUco码携带的信息实现累计误差消除，从而实现狭长空间弱纹理场景下探测机器人的高精度定位。
![]()
图 1 探测机器人多传感器融合定位算法Figure 1. Detection robot multi-sensor fusion positioning algorithm探测机器人车身上安装的里程计、IMU、单目相机数据作为算法输入，其中里程计需利用UMBmark算法进行标定，以消除系统误差。
工控机通过串口获取到传感器数据后，将里程计与IMU数据进行时间戳对齐，并且判断是否检测到ArUco码。未检测到ArUco码时，将里程计和IMU数据通过EKF算法实现融合定位。当检测到ArUco码时，对当前探测机器人位姿状态利用图优化方法进行局部优化。随后将优化结果和参考位姿信息结合，修正融合定位结果，消除累计误差影响，最后将定位结果输出。
1.2 里程计与IMU数据融合
考虑单一传感器容易出现测量错误，且依靠里程计推算得到的偏航角容易出现较大误差，因此将里程计数据和IMU数据利用EKF算法进行融合，避免单一传感器测量误差影响定位结果。
EKF算法的实现包括预测和更新2个步骤。设
$\gamma $ 时刻状态变量为$ {{\boldsymbol{x}}_\gamma } = {\left[ {{X_\gamma }\;{Y_\gamma }\;{\theta _\gamma }} \right]^{\rm{T}}} $ ，$ {X_\gamma },{Y_\gamma } $ 分别为$\gamma $ 时刻探测机器人二维坐标系的横纵坐标，$ {\theta _\gamma } $ 为$\gamma $ 时刻探测机器人偏航角。根据轮式里程计圆弧运动模型进行状态预测：$$ {{\mathop {\boldsymbol{x}}\limits^ \vee} _{{\gamma}} } = {{\mathop {\boldsymbol{x}}\limits^ \wedge } _{{{\gamma}} - {{1}}}} + \left[ {\begin{array}{*{20}{c}} {\Delta {d_{\gamma - 1,\gamma }}\cos ({\theta _{\gamma - 1}} + \Delta {\theta _{\gamma - 1,\gamma }}/2)} \\ {\Delta {d_{\gamma - 1,\gamma }}\sin ({\theta _{\gamma - 1}} + \Delta {\theta _{\gamma - 1,\gamma }}/2)} \\ {{\theta _{\gamma - 1}} + \Delta {\theta _{\gamma - 1,\gamma }}} \end{array}} \right] $$ (1) 式中：
$ {{\mathop {\boldsymbol{x}}\limits^ \vee} _{{\gamma}} } $ 为$\gamma $ 时刻状态变量的先验估计；$ {{\mathop {\boldsymbol{x}}\limits^ \wedge} _{{{\gamma}} -{{ 1}}}} $ 为$\gamma - 1$ 时刻状态变量的后验估计；$ \Delta {d_{\gamma - 1,\gamma }} $ 为$\gamma - 1$ 时刻到$\gamma $ 时刻的位移测量值；$ {\theta _{\gamma - 1}} $ 为$\gamma - 1$ 时刻探测机器人偏航角；$ \Delta {\theta _{\gamma - 1,\gamma }} $ 为$\gamma - 1$ 时刻到$\gamma $ 时刻探测机器人偏航角变化。状态变量的协方差预测方程为
$$ {{\mathop {\boldsymbol{P}}\limits^ \vee } _{{\gamma}} } = {{\boldsymbol{F}}_{{\gamma}} }{{\mathop {\boldsymbol{P}}\limits^ \wedge } _{{{\gamma}} - {{1}}}}{\boldsymbol{F}}_\gamma^{\rm{T}} + {\boldsymbol{R}} $$ (2) 式中：
$ {{\mathop {\boldsymbol{P}}\limits^ \vee} _{{\gamma}} } $ 和$ {{\mathop {\boldsymbol{P}}\limits^ \wedge } _{{{\gamma}} - {{ 1}}}} $ 分别为$\gamma $ 时刻状态变量的先验协方差矩阵和$\gamma - 1$ 时刻后验协方差矩阵；$ {{\boldsymbol{F}}_{{\gamma}} } $ 为$\gamma $ 时刻运动方程对状态变量求导得到的协方差矩阵；$ {\boldsymbol{R}} $ 为系统噪声协方差矩阵。$$ {{\boldsymbol{F}}_{{\gamma}} } = \left[ {\begin{array}{*{20}{c}} 1&0&{ - \Delta {d_{\gamma - 1,\gamma }}\sin ({\theta _{\gamma - 1}} + \Delta {\theta _{\gamma - 1,\gamma }}/2)} \\ 0&1&{\Delta {d_{\gamma - 1,\gamma }}\cos ({\theta _{\gamma - 1}} + \Delta {\theta _{\gamma - 1,\gamma }}/2)} \\ 0&0&1 \end{array}} \right] $$ (3) 将IMU数据积分得到的位置信息
$({X_{\gamma }},{Y_{\gamma }})$ 、偏航角信息${\theta _{\gamma }}$ 作为$\gamma $ 时刻观测值，则有如下观测方程：$$ {{\boldsymbol{z}}_{{\gamma}} } = \left[ {\begin{array}{*{20}{c}} {{X_{\gamma }}} \\ {{Y_{\gamma }}} \\ {{\theta _{\gamma }}} \end{array}} \right] + {{\boldsymbol{V}}_\gamma } $$ (4) 式中
${{\boldsymbol{V}}_\gamma }$ 为观测噪声。根据式(3)和式(4)进行状态更新，卡尔曼增益计算公式为
$$ {{\boldsymbol{\kappa}} _{{\gamma}} } = {{\mathop {\boldsymbol{P}}\limits^ \vee} _{{\gamma}} }{\kern 1pt} {{\boldsymbol{\varTheta}}^{\bf{T}}}{\left( {{\boldsymbol{\varTheta}}{\kern 1pt} {{\mathop {\boldsymbol{P}}\limits^ \vee }_{{\gamma}} }{{\boldsymbol{\varTheta}}^{\bf{T}}} + {{\boldsymbol{Q}}_{{\gamma}} }} \right)^{ - 1}} $$ (5) 式中：
$ {{\boldsymbol{\kappa}} _{\boldsymbol{\gamma}} } $ 为$\gamma $ 时刻计算得到的卡尔曼增益；${\boldsymbol{\varTheta}}$ 为观测方程对观测值求偏导得到的雅克比矩阵；$ {{\boldsymbol{Q}}_{{\gamma}} } $ 为观测噪声的协方差矩阵。根据式（5）可得状态变量的后验估计：
$$ {{\mathop {\boldsymbol{x}}\limits^ \wedge} _{\boldsymbol{\gamma}} } = {{\mathop {\boldsymbol{x}}\limits^ \vee } _{\boldsymbol{\gamma}} } + {{\boldsymbol{\kappa}} _{\boldsymbol{\gamma}} }\left[ {\begin{array}{*{20}{c}} {{X_{\gamma }} - {{\mathop {X}\limits^ \vee }_\gamma }} \\ {{Y_{\gamma }} - {{\mathop {Y}\limits^ \vee }_\gamma }} \\ {{\theta _{\gamma }} - {{\mathop {{\theta}} \limits^ \vee }_\gamma }} \end{array}} \right] $$ (6) 式中
${{\mathop {{X}}\limits^ \vee} _\gamma }$ ，${{\mathop {{Y}}\limits^ \vee } _\gamma }$ ，${{\mathop {{\theta}} \limits^ \vee} _\gamma }$ 分别为$\gamma $ 时刻状态变量的先验估计值。状态变量的后验协方差矩阵为
$$ {{\mathop {\boldsymbol{P}}\limits^ \wedge } _{{\gamma}} } = \left( {{\boldsymbol{I}} - {{\boldsymbol{\kappa}} _{{\gamma}} }{\boldsymbol{\varTheta }}} \right){{\mathop {\boldsymbol{P}}\limits^ \vee } _{{\gamma}} } $$ (7) 式中
$ {\boldsymbol{I}} $ 为单位矩阵。1.3 相对位姿关系
ArUco码到探测机器人坐标系需要经过多次坐标变换，假设以洞口位置探测机器人的车体中心为原点的三维世界坐标系下某点p的齐次坐标为
${{\boldsymbol{p}}_{{{{\rm{wq}}}}}} = $ $ {\left[ {{X_{\rm{w}}}\;\;{Y_{\rm{w}}}\;\;{Z_{\rm{w}}}\;1} \right]^{\rm{T}}}$ ，且其对应的相机坐标系下坐标为$ {{\boldsymbol{p}}_{{{\rm{c}}}}} = $ $ {[{X_{\text{c}}}\;\;\;{Y_{\text{c}}}\;\;\;{Z_{\text{c}}}]^{\rm{T}}} $ ，则根据针孔相机模型可得到$ {{\boldsymbol{p}}_{{{\rm{c}}}}} $ 到像素坐标系下坐标$ {\boldsymbol{p}} = {[{m}\;\;\;{n}]^{\rm{T}}} $ 的变换关系：$$ \left[ {\begin{array}{*{20}{c}} {{m_1}} \\ {{n_1}} \\ 1 \end{array}} \right] = \frac{1}{{{Z_{\rm{c}}}}}{\boldsymbol{K}}{{\boldsymbol{p}}_{{{\rm{c}}}}} = \frac{1}{{{Z_{\rm{c}}}}}\left[ {\begin{array}{*{20}{c}} {{f_{{m}}}}&0&{{c_{{m}}}} \\ 0&{{f_{{n}}}}&{{c_{{n}}}} \\ 0&0&1 \end{array}} \right]{\boldsymbol{A}}\left[ {\begin{array}{*{20}{c}} {{X_{\rm{w}}}} \\ {{Y_{\rm{w}}}} \\ {{Z_{\rm{w}}}} \\ 1 \end{array}} \right] $$ (8) 式中：
${\boldsymbol{K}}$ 为相机内参数矩阵，一般由标定得到；${f_{{m}}}$ ，${f_{{n}}}$ 分别为像素平面m轴和n轴上的尺度因子；${c_{{m}}}$ ，${c_{{n}}}$ 为相机光心偏移量；A为变换矩阵，$ {\boldsymbol{A}} = [{\boldsymbol{r}}\;\;\;{\boldsymbol{\tau}} ] $ ，$ {\boldsymbol{r}} $ ，${\boldsymbol{\tau}}$ 分别为相机当前位姿的旋转矩阵和平移矩阵。实际安装时相机并不处于探测机器人坐标系中心，故某点p在探测机器人坐标系下的坐标需进行坐标变换得到。
$$ {{\boldsymbol{p}}_{{{{\rm{rq}}}}}} = {}_{{{\rm{c}}}}^{{{\rm{r}}}}{\boldsymbol{T}}{{\boldsymbol{p}}_{{{{\rm{cq}}}}}} = {}_{{{\rm{c}}}}^{{{\rm{r}}}}{\boldsymbol{T}}{}_{{{\rm{w}}}}^{{{\rm{c}}}}{\boldsymbol{T}}{{\boldsymbol{p}}_{{{{\rm{wq}}}}}} = {}_{{{\rm{c}}}}^{{{\rm{r}}}}{\boldsymbol{T}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}} \\ {\begin{array}{*{20}{c}} {{{\bf{0}}^{{{\rm{T}}}}}}&{{1}} \end{array}} \end{array}} \right]{{\boldsymbol{p}}_{{{{\rm{wq}}}}}} $$ (9) 式中：
${{\boldsymbol{p}}_{{\rm{rq}}}}$ ，${{\boldsymbol{p}}_{{\rm{cq}}}}$ 分别为某点p在探测机器人坐标系、相机坐标系下的齐次坐标；${}_{\rm{c}}^{\rm{r}}{\boldsymbol{T}}$ 为探测机器人坐标系和相机坐标系之间的变换矩阵，由于相机固定在探测机器人车身上，所以其中参数均为常数；${}_{\rm{w}}^{\rm{c}}{\boldsymbol{T}}$ 为世界坐标系到相机坐标系的变换矩阵。根据式( 8)和式( 9)可由像素坐标推算得到探测机器人坐标系与世界坐标系下ArUco码之间的位姿变换关系。
1.4 ArUco码识别与解算
由于相机固定在探测机器人车身上，所以求得像素坐标系相对于ArUco码所在的世界坐标系的位姿变换矩阵及像素坐标系与相机坐标系之间的位姿变换矩阵，即可得到相机和ArUco码的位置变换关系。
首先使用张正友标定法对相机进行标定，得到相机的内参矩阵和畸变参数。然后分别通过自适应滤波对图像进行二值化，利用Suzuki算法进行轮廓提取并再次进行滤波，采用Douglas−Peucker算法进行四边形轮廓拟合并且得到4个角点。因ArUco码外侧边界均为黑色，可据此进行筛选得到ArUco码的候选轮廓，从而得到其角点在世界坐标系的位置。最后通过求解单应矩阵即可得到2帧图像间探测机器人的运动估计。
考虑角点位于ArUco所在平面，可设其为
$ {{\boldsymbol{C}}^{\rm{T}}}{{p}} + $ $ {{D}} = 0 $ ，其中${\boldsymbol{C}}$ 为系数矩阵，$ {{D}} $ 为常数。则可推导出2帧图像的像素坐标系中一对匹配的角点对$ {{\boldsymbol{p}}_{{1}}} = $ $ {[{m_1}\;\;\;{n_1}]^{\rm{T}}} $ 和$ {{\boldsymbol{p}}_{{2}}} = {[{m_2}\;\;\;{n_2}]^{\rm{T}}} $ 的转换关系：$$ s{{\boldsymbol{p}}_{{2}}} = {\boldsymbol{K}}({\boldsymbol{r}}{{\boldsymbol{p}}_{{1}}} + {\boldsymbol{\tau }}) = {\boldsymbol{K}}({\boldsymbol{r}} - \frac{{{\boldsymbol{\tau }}{{\boldsymbol{C}}^{\mathbf{T}}}}}{{{D}}}){{\boldsymbol{K}}^{ - 1}}{{\boldsymbol{p}}_{{1}}} = {{\boldsymbol{H}}_{{{3}} \times {{3}}}}{{\boldsymbol{p}}_{{1}}} $$ (10) 式中s为尺度因子。
根据匹配点参数即可求解出
${{\boldsymbol{H}}_{{{3}} \times {{3}}}}$ 矩阵：$$ s\left[ {\begin{array}{*{20}{c}} {{m_2}} \\ {{n_2}} \\ 1 \end{array}} \right] = {{\boldsymbol{H}}_{{{3}} \times {{3}}}}\left[ {\begin{array}{*{20}{c}} {{m_1}} \\ {{n_1}} \\ 1 \end{array}} \right] $$ (11) 展开式（11）得
$$ \left\{ {\begin{array}{*{20}{c}} {{\rm{s}}{m_2} = {m_1}{h_{11}} + {n_1}{h_{12}} + {h_{13}}} \\ {{\rm{s}}{n_2} = {m_1}{h_{21}} + {n_1}{h_{22}} + {h_{23}}} \\ {{\rm{s}}{{_{}^{}}_{}} = {m_1}{h_{31}} + {n_1}{h_{32}} + {h_{33}}} \end{array}} \right. $$ (12) 式中
${h_{\alpha \beta }}$ $(\alpha ,\beta=1,2,3) $ 为${{\boldsymbol{H}}_{{\mathbf{3}} \times {\mathbf{3}}}}$ 矩阵中第$\alpha $ 行第$\beta $ 列元素。进一步消去尺度因子
$ s $ ，且当$ {h_{33}} $ 非0时，设$ {h_{33}} = 1 $ ，将${{\boldsymbol{H}}_{{\mathbf{3}} \times {\mathbf{3}}}}$ 矩阵看作列向量进行化简：$$ \left[ {\begin{array}{*{20}{c}} {{m_1}}&{{n_1}}&1&0&0&0&{ - {m_1}{m_2}}&{ - {n_1}{m_2}}\\ 0&0&0&{{m_1}}&{{n_1}}&1&{ - {m_1}{n_2}}&{ - {n_1}{n_2}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{h_{11}}} \\ {{h_{12}}} \\ {{h_{13}}} \\ {\begin{array}{*{20}{c}} {{h_{21}}} \\ {{h_{22}}} \\ {{h_{23}}} \\ {\begin{array}{*{20}{c}} {{h_{31}}} \\ {{h_{32}}} \end{array}} \end{array}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{m_2}} \\ {{n_2}} \end{array}} \right] $$ (13) 式中有
$ {h_{11}} $ ，$ {h_{12}} $ ，$ {h_{13}} $ ，$ {h_{21}} $ ，$ {h_{22}} $ ，$ {h_{23}} $ ，$ {h_{31}} $ ，$ {h_{32}} $ 8个未知量，故利用4个角点的匹配点对即可建立方程组，完成相邻2帧图像间探测机器人位姿变换的初步计算。此外，根据提取的轮廓可进行二进制编码解算，由于每个ArUco码表示唯一，只有原始角点顺序下的ArUco码可在字典中找到编号，所以对提取的图像进行4次旋转，分别计算其编码值，即可得到一个唯一确定的编号。由于在首次完成ArUco码放置后，每个ArUco码在世界坐标系下的位姿均为固定值，根据其编号可确定ArUco码表示的位姿信息。1.5 姿态和位置信息优化
探测机器人携带单目相机进入隧洞，当检测到ArUco码时，按照1.3和1.4节所述方法计算得到探测机器人和ArUco码之间的相对位姿关系及ArUco码代表的绝对位置信息。考虑到探测机器人在运动过程中传感器存在误差，因此，使用图优化方式进行姿态和位置信息的优化。将全过程待优化的状态变量记为
$\chi $ ，则有$$ {\boldsymbol{\chi}} = [{{\boldsymbol{\eta}} _0}\;\;\;{{\boldsymbol{\eta}} _1}\;\cdots\;\;\;{{\boldsymbol{\eta}} _i}\;\;\;{\lambda _0}\;\;\;{\lambda _1}\;\cdots\;{\lambda _j}] $$ (14) 第k帧图像采集时刻的位姿为
$$ {{\boldsymbol{\eta}} _{{k}}} = [{w_k}\;\;\;{q_k}\;\;\;{v_k}\;\;\;b_k^{\rm{a}}\;\;\;b_k^{\rm{g}}] $$ (15) 式中：
$ {{\boldsymbol{\eta}} _i} $ 为第$i$ 个观测点处的本体位姿；${\lambda _j}$ 为观测到的第$j$ 个路标位置；${w_k}$ 为第k帧图像采集时刻的IMU位置；${q_k}$ 为第k帧图像采集时刻的IMU姿态；${v_k}$ 为第k帧图像采集时刻的IMU速度；$b_k^{\rm{a}}$ ，$b_k^{\rm{g}}$ 分别为陀螺仪和加速度计的偏移。在实际环境下，探测机器人在运动过程中存在噪声，基于高斯分布假设和最大似然估计，可建立如下目标函数：
$$ {J^*} = \arg \max \phi(z,u|{\boldsymbol{\eta}} ,\lambda ) $$ (16) 式中：
$ z,u $ 分别为观测数据和输入数据；$ {\boldsymbol{\eta}} ,\lambda $ 分别为待估计的位姿和路标位置。假设各个时刻的输入数据和观测数据之间相互独立，则联合分布可分解为
$$ \phi(z,u|{\boldsymbol{\eta}} ,\lambda ) = \prod\limits_k {\phi({u_k}|{{\boldsymbol{\eta}} _{k - 1}},{{\boldsymbol{\eta}} _k})\prod\limits_{k,j} {\phi({z_{k,j}}|{{\boldsymbol{\eta}} _k},{\lambda _j})} } $$ (17) 式中：
$ {u_k} $ 为第k帧的输入数据；$ {{\boldsymbol{\eta}} _{k - 1}} $ ，${{\boldsymbol{\eta}} _k}$ 分别为第$k-1 $ 帧和第k帧的位姿估计值；$ {z_{k,j}} $ 为第k帧探测机器人对第$j$ 个路标位置的观测值。设运动模型噪声
$ {{\boldsymbol{e}}_{u,k}} $ 服从正态分布$ N(0,{{\boldsymbol{\sigma}} _{u,k}}) $ ，观测模型噪声${{\boldsymbol{e}}_{z,k,j}}$ 服从正态分布$N(0,{{\boldsymbol{\sigma}} _{z,k,j}})$ ，$ {{\boldsymbol{\sigma}} _{u,k}},{{\boldsymbol{\sigma}} _{z,k,j}} $ 分别为2类噪声的协方差矩阵。那么单独考虑相机观测方程$ {z_{k,j}} = g({{\boldsymbol{\eta}} _k},{\lambda _j}) $ 并联立式（16）、式（17），同时对两侧取负对数并化简得$$ \begin{split} & {\rm{min}}[ - {\rm{ln}}\phi(z|{\boldsymbol{\eta}} ,\lambda )] = G + {\rm{min}}\left[ {\frac{1}{2}{{(z - g({\boldsymbol{\eta}} ,\lambda ))}^{\rm{T}}}{\boldsymbol{\sigma}} _{z,k,j}^{ - 1}(z - g({\boldsymbol{\eta}} ,\lambda ))} \right] =\hfill \\ &\quad\quad\quad\quad\quad\quad\quad\;\; G + \sum\limits_{{{k}} = 1}^{{L}} {\sum\limits_{{{j}} = 1}^{{E}} {\frac{1}{2}{\boldsymbol{e}}_{z,k,j}^{\rm{T}}} } {\boldsymbol{\sigma}} _{z,k,j}^{ - 1}{{\boldsymbol{e}}_{z,k,j}} \hfill \\[-18pt] \end{split} $$ (18) 式中：
$ G $ 为常数；L为进行优化的位姿观测点数；E为观测到的路标数。同理，单独考虑输入情况方程
$l({{\boldsymbol{\eta }}_{k - 1}},{{\boldsymbol{\eta}} _k})$ ，有$$ {J^*} = {\rm{argmin}}(\sum\limits_{{{k}} = 1}^{{L}} {\frac{1}{2}{\boldsymbol{e}}_{u,k}^{\rm{T}}{\boldsymbol{\sigma}} _{u,k}^{ - 1}{{\boldsymbol{e}}_{u,k}}} + \sum\limits_{{{k}} = 1}^{{L}} {\sum\limits_{{{j}} = 1}^{{E}} {\frac{1}{2}{\boldsymbol{e}}_{z,k,j}^{\rm{T}}} } {\boldsymbol{\sigma}} _{z,k,j}^{ - 1}{{\boldsymbol{e}}_{z,k,j}}) $$ (19) 待求解优化函数可表示为
$$ \begin{split} {J^*} =& \mathop {\arg \min }\limits_{\boldsymbol{\chi}} \Bigg\{ \sum\limits_{k \in B} {\mu \rho (||{M_{\rm{IMU}}}(\tilde Z_{{\rm{IMU}}k + 1}^{{\rm{IMU}}k},{\boldsymbol{\chi}} )||^2)} + \Bigg.\\ & \qquad\quad \Bigg. \sum\limits_{(i,j) \in W} {(1 - \mu )\rho (||{M_{\rm{f}}}(\tilde Z_{{f_j}}^{{c_i}},{\boldsymbol{\chi}} )||^2)} \Bigg\} = \hfill \\ &\qquad\quad \mathop {\arg \min }\limits_{\boldsymbol{\chi}} \left\{ {\mu {U_{{\rm{IMU}}}} + (1 - \mu ){U_{{\rm{cam}}}}} \right\} \hfill \end{split} $$ (20) 式中：
$B$ 为第k帧图像采集时刻IMU测量数据集合；$\mu $ 为比例因子，2帧图像之间探测机器人在ArUco码范围内识别的特征点信息越多且运动速度越慢，则$\mu $ 越小，表明该情况下依靠相机的位姿估计准确性更高；$\rho $ 为huber核函数，当匹配偏差超过阈值使其保持线性增长，从而降低错误匹配带来的影响；${M_{\rm{IMU}}}(\tilde Z_{{\rm{IMU}}k + 1}^{{\rm{IMU}}k},{\boldsymbol{\chi}} )$ 为IMU数据的残差项($\tilde Z_{{\rm{IMU}}k + 1}^{{\rm{IMU}}k} $ 为IMU的观测值)，由IMU测量预积分形成的观测值和实际值比较形成，用于约束k和k+1连续2帧；$W$ 为匹配良好的特征点对集合；${M_{\rm{f}}}(\tilde Z_{{f_j}}^{{c_i}},{\boldsymbol{\chi}} )$ 为相机残差项($\tilde Z_{{f_j}}^{{c_i}} $ 为相机观测值)，通过将特征点$ {c_i} $ 投影到第${f_j}$ 帧并与测量值计算重投影残差得到；${U_{{\rm{IMU}}}}$ 和${U_{{\rm{cam}}}}$ 分别为IMU和相机约束。利用里程计信息、IMU信息和相机信息，根据图优化的方式建立如图2所示的约束因子，图中
${t_k}$ 为时间戳对齐至第k帧的时间，v'k为里程计测量速度。在ArUco码附近范围实现探测机器人位姿优化，其中里程计约束提供的是相对约束，转换到IMU坐标系中进行速度和位姿的限制。结合绝对位置信息和局部位姿优化结果更新探测机器人位姿，同时完成对IMU和里程计累计误差的修正。
2. 隧洞探测机器人系统设计
为验证多传感器融合定位算法在尾矿库排洪隧洞实际环境中的运行效果，研究并设计了隧洞探测机器人系统，如图3所示。
尾矿库排洪隧洞探测机器人系统以履带式探测机器人实验平台为本体，基于机器人操作系统（Robot Operating System，ROS）进行软件设计，采用分布式系统，将监视端工控机布置于隧洞入口，而运行多传感器融合定位算法的工控机安装在探测机器人本体上，上下位机之间利用无线网桥进行组网通信。为实现多传感器融合定位算法，隧洞探测机器人实验平台和洞口的监控平台需要进行相应的软硬件设计，且根据系统功能将隧洞探测机器人系统划分为4个子系统：驱动控制系统、主控系统、传感通信系统、人机交互系统。驱动控制系统对接主控系统，进行电动机控制和电动机相关数据传输，保证探测机器人的运动能力。主控系统对接驱动控制系统和传感通信系统，运行多传感器融合定位、定位校正等核心算法，实现探测机器人自主定位和定位校正等功能，并解析执行控制端操作人员指令，通过串口转RS232实现与驱动控制系统的通信，通过USB获取传感器数据，在进行数据分析后将关键数据和音视频通过传感通信系统传输至人机交互系统，并实现人机交互系统的各类指令下传，以实现对探测机器人本体的控制。传感通信系统对接主控系统和人机交互系统，负责探测机器人环境感知和通信，使洞口操作人员可实时查看探测机器人实验平台的运行状态。人机交互系统对接传感通信系统，负责解析操作人员指令并通过传感通信系统将指令发送到主控系统，将主控系统数据实时显示在监视端工控机界面中。
3. 实验验证和分析
为探究搭载多传感器融合定位算法的尾矿库排洪隧洞探测机器人实验平台在实际环境中运行的稳定性，并验证多传感器融合定位算法的有效性，在尾矿库排洪隧洞实际场景中运行和关闭多传感器融合定位算法，分别进行10组行进40 m的重复定位实验，且每轮实验均记录探测机器人在20 m和40 m处的定位误差。运行多传感器融合定位算法需要在20 m处额外记录校正后误差。
在实验开始前需要标定和校正系统误差，探测机器人平均车轮直径为0.21 m，理论轮距为 0.88 m，使用UMBmark算法进行标定，根据实验数据求得校正后左轮直径为0.209 5 m，右轮直径为0.210 5 m，轮间距为0.84 m。
修改里程计参数后，在距起点20 m处的排洪隧洞墙壁上布设ArUco码。在探测机器人系统中利用ROS编写程序，对校正后的里程计和IMU数据进行EKF，利用OpenCV库进行ArUco码识别与相应图像处理操作，并利用Graph SLAM Tutorial实现图优化算法对探测机器人位姿进行优化并校正定位。在尾矿库排洪隧洞中进行重复定位实验，实验环境和实验平台实物如图4所示。
![]()
图 4 尾矿库重复定位实验环境和实验平台Figure 4. Experimental environment and platform for repeated positioning of tailings pond根据实验设计运行和关闭多传感器融合定位算法，在相同环境中各进行10组40 m行进实验，并记录实验数据，各组实验20 m和40 m处定位误差如图5所示。根据图5数据统计定位误差，结果见表1。
表 1 探测机器人定位误差Table 1. Positioning error of detection robot定位位置 平均定位误差/cm 关闭算法实验组 运行算法实验组 行进20 m 22.51 19.77 行进40 m 58.11 21.23 运行多传感器融合定位算法实验组利用行进20 m时的ArUco码进行误差校正。使用米尺测量ArUco码中心和探测机器人车身中心的
$X$ ，$Y$ 方向数据，与位姿估计结果比较并计算距离误差，得误差均值为4.2 mm。从图5、表1可看出，使用多传感器融合定位算法的探测机器人在行进20 m时平均定位误差为19.77 cm，在行进40 m时的平均定位误差为21.23 cm。同样条件下，关闭多传感器融合定位算法实验组在行进20 m时的平均定位误差为22.51 cm，在行进40 m时的平均定位误差为58.11 cm。可见使用多传感器融合定位算法的探测机器人在运行过程中可保持更高定位精度，且在20 m的ArUco码放置处多传感器融合定位算法可有效修正累计误差，使校正后误差均值为4.2 mm。
4. 结论
（1） 针对现有定位算法不适用于网络信号弱、场景重复率高且属于狭长室内环境的尾矿库排洪隧洞探测机器人定位的问题，提出了一种探测机器人多传感器融合定位算法。将排洪隧洞环境的远距离定位以ArUco码为边界点，简化为多段短距离定位。未检测到ArUco码时，利用EKF算法融合里程计和IMU信息，可实现位姿信息的快速更新；检测到ArUco码时，则利用图优化算法将里程计、IMU、相机信息组成约束因子图，并结合ArUco码的绝对位置实现探测机器人局部位姿优化和定位校正，避免误差累计，最终实现狭长空间内探测机器人的远距离高精度定位。
（2） 为验证多传感器融合定位算法的有效性，设计了隧洞探测机器人系统，在尾矿库排洪隧洞实际场景中运行和关闭多传感器融合定位算法，分别进行10组行进40 m的重复定位实验，实验结果表明：多传感器融合定位算法具有较高的稳定性和精确性，能够实现隧洞探测机器人在排洪隧洞环境中的高精度定位，行进20 m的平均定位误差为19.77 cm，40 m的平均定位误差为21.23 cm，且20 m处校正后误差均值为4.2 mm。
致谢：本文得到了北京科技大学、福建马坑矿业股份有限公司、中国安全生产科学研究院和中国计量大学人工智能与机器人团队的支持，对此表示感谢！
-
![]()
图 1 相似材料模拟实验方案及模型
Figure 1. Scheme and model of Similar material simulation experiment
![]()
图 7 含水层失水状态下各观测线的下沉量变化
Figure 7. The variation of subsidence of each observation line under water loss state of aquifer
表 1 研究区地层结构
Table 1 A histogram of stratigraphic structure in the study area
界 系 统 组 厚度(m) 主要岩性 新生界 第四系 全新统 — 40～130 浅黄、灰黄色粘土夹砂层 更新统 第三系上 上新统 — 0～152 灰绿、浅黄，多为粘土夹杂砂层 中新统 第三系下 渐新统 — >205 浅灰、棕色砂泥岩互层，夹杂砂砾岩 始新统 中生界 白垩系 上统 — >647 紫红色粉、细砂岩，砂砾岩 下统 — 844 棕红粉砂岩、泥岩及细中粒砂岩 侏罗系 上统 — >637 凝灰质砂砾岩，凝灰岩和安山岩 三叠系 下统 — 316～446 紫红色砂泥岩 古生界 二叠系 上统 石千峰组 114～260 紫红杂色砂泥岩，夹石英砂岩及砂砾岩 上石盒子组 316～566 灰绿色和浅灰色砂岩，底为石英砂岩且为含煤层 下统 下石盒子组 106～265 灰色砂泥岩及其互层，底含粗砂岩，含煤层 山西组 52～88 上部细砂岩、粗砂岩，下部深灰色泥岩，含煤层 石炭系 上统 太原组 102～148 灰岩为主，夹泥岩及砂岩，含薄煤层 奥陶系 中下统 — 400 中厚层为白云岩及白云质灰岩，夹灰岩 寒武系 上统 土坝组 170～220 硅质结核白云岩，产Heleionellasp.化石 固山组 9～78 白云岩，竹叶状灰岩，鲕状灰岩。 中统 张夏组 146 鲕状灰岩，白云岩产Dameselluasp.化石 徐庄组 190 棕黄砂岩，夹页岩及石灰岩 毛庄组 152 多为砾状灰岩，鲕状灰岩和页岩 下统 馒头组 215 紫色页岩夹灰岩，产Redlichasp.化石 猴家山组 100～150 鲕状灰岩，孔洞灰岩和砂灰岩 凤台组 10～100 页岩，砾岩 上元古界 震旦系 徐淮群 九顶山组 117 白云岩，底部夹竹叶状灰岩 倪园组 92 上部含泥白云岩，夹黄绿色钙质页岩，下部硅质条带白云岩 四顶山组 137 厚层白云岩为主，产蠕形动物化石 九里桥组 119 泥灰岩，砂灰岩 四十里长山组 93 石英岩及钙质砂岩 青白口系 八公
山群刘老碑组 1050 页岩，泥灰岩，石英砂岩，底部铁质砂砾岩，含藻及疑源类化石 伍山组 张店组 下元古界 — — 凤阳群 1 171 千枚岩，白云岩，大理岩，白云质石英片岩，石英岩，含藻化石 上太古界 — — 五河群 >6 422 片麻岩，浅粒岩，变粒岩，斜长角闪岩互层，夹少量大理岩及磁 铁矿层，岩石混合岩化 
下载: 导出CSV
-
[1] 谢和平, 吴立新, 郑德志. 2025 年中国能源消费及煤炭需求预测[J]. 煤炭学报,2019,44(7):1949-1960. XIE Heping, WU Lixin, ZHENG Dezhi. Prediction on the energy consumption and coal demand of China in 2025[J]. Journal of China Coal Society,2019,44(7):1949-1960.
[2] 高永格, 牛矗, 张强, 等. 厚松散层下采煤地表沉陷特征研究[J]. 煤炭科学技术,2019,47(6):192-198. GAO Yongge, NIU Chu, ZHANG Qiang, et al. Study on surface subsidence characteristics of coal mining under thick loose layer[J]. Coal Science and Technology,2019,47(6):192-198.
[3] 张海峰, 李文, 李少刚, 等. 浅埋深厚松散层综放工作面覆岩破坏监测技术[J]. 煤炭科学技术,2014,42(10):24-27. ZHANG Haifeng, LI Wen, LI Shaogang, et al. Research on technology of overlying strata failure monitoring in fully mechanized caving coal face with shallow and thick loose bed[J]. Coal Science and Technology,2014,42(10):24-27.
[4] 陈俊杰, 邹友峰, 郭文兵. 厚松散层下下沉系数与采动程度关系研究[J]. 采矿与安全工程学报,2012,29(2):250-254. DOI: 10.3969/j.issn.1673-3363.2012.02.018 CHEN Junjie, ZOU Youfeng, GUO Wenbing. Study on the relationship between subsidence coefficient and mining degree under a thick alluvium stratum[J]. Journal of Mining & Safety Engineering,2012,29(2):250-254. DOI: 10.3969/j.issn.1673-3363.2012.02.018
[5] 李新岭, 郭文兵, 赵高博. 巨厚松散层土体压缩特性对开采沉陷影响研究[J]. 中国安全科学学报,2018,28(7):135-141. LI Xinling, GUO Wenbing, ZHAO Gaobo. Study on influence of compression characteristics of super-thick alluvium on mining subsidence[J]. China Safety Science Journal,2018,28(7):135-141.
[6] 许国胜, 许胜军, 李德海, 等. 赵固矿区厚松散层下开采地表移动特征实测研究[J]. 中国矿业, 2021, 30(5): 168-172. XU Guosheng, XU Shengjun, LI Dehai, et al. Field research on surface movement characteristic by coal mining under thick alluvium in Zhaogu mine area[J]China Mining Magazine, 2021, 30(5): 168-172.
[7] 许国胜, 李德海, 侯得峰, 等. 厚松散层下开采地表动态移动变形规律实测及预测研究[J]. 岩土力学,2016,37(7):2056-2062. XU Guosheng, LI Dehai, HOU Defeng, et al. Measurement and prediction of the transient surface movement and deformation induced by mining under thick alluvium[J]. Rock and Soil Mechanics,2016,37(7):2056-2062.
[8] 崔希民, 方志海, 左红飞, 等. 开采引起的含水层失水对地表下沉的影响[J]. 煤田地质与勘探,2000,28(5):47-48. DOI: 10.3969/j.issn.1001-1986.2000.05.015 CUI Ximin, FANG Zhihai, ZUO Hongfei, et al. De-watering effect of aquifer to surface subsidence[J]. Coal Geology & Exploration,2000,28(5):47-48. DOI: 10.3969/j.issn.1001-1986.2000.05.015
[9] 崔希民, 邓喀中. 煤矿开采沉陷预计理论与方法研究评述[J]. 煤炭科学技术,2017,45(1):160-169. CUI Ximin, DENG Kazhong. Research review of predicting theory and method for coal mining subsidence[J]. Coal Science and Technology,2017,45(1):160-169.
[10] 李春意, 郭增长. 含水松散层固结对地表沉陷影响的数值分析[J]. 河南理工大学学报(自然科学版),2012,31(1):48-53. LI Chunyi, GUO Zengzhang. Numerical analysis of water-bearing alluvium consolidation influence on surface subsidence[J]. Journal of Henan Polytechnic University(Natural Science),2012,31(1):48-53.
[11] 张丁丁. 兖州矿区第四系厚松散层沉降特性研究[D]. 西安: 西安科技大学, 2015. ZHANG Dingding. Study on settlement characteristics of thick quaternary unconsolidated layer in Yanzhou mining area. [D]. Xi′an: Xi′an University of Science and Technology, 2015.
[12] 张劲满. 厚松散层含水层失水开采条件下覆岩运动及地表移动规律研究[D]. 淮南: 安徽理工大学, 2019. ZHANG Jinman. Study on overburden movement and surface movement under water loss mining in thick loose aquifer[D]. Huainan: Anhui University of Science and Technology, 2019.
[13] 杜锋, 白海波. 厚松散层薄基岩综放开采覆岩破断机理研究[J]. 煤炭学报,2012,37(7):1105-1110. DU Feng, BAI Haibo. Mechanism research of overlying strata activity with fully mechanized caving in thin bedrock with thick alluvium[J]. Journal of China Coal Society,2012,37(7):1105-1110.
[14] 梁庆华, 温兴林, 何刚, 等. 黏土体失水引起的地表下沉计算方法研究[J]. 采矿与安全工程学报,2007,24(1):105-108. DOI: 10.3969/j.issn.1673-3363.2007.01.023 LIANG Qinghua, WEN Xinglin, HE Gang, et al. Study on calculation methods for surface subsidence caused by water loss of clay[J]. Journal of Mining & Safety Engineering,2007,24(1):105-108. DOI: 10.3969/j.issn.1673-3363.2007.01.023
[15] 赵春虎, 靳德武, 王皓, 等. 榆神矿区中深煤层开采覆岩损伤变形与含水层失水模型构建[J]. 煤炭学报,2019,44(7):2227-2235. ZHAO Chunhu, JIN Dewu, WANG Hao, et al. Construction and application of overburden damage and aquifer water loss model in medium-deep buried coal seam mining in Yushen mining area[J]. Journal of China Coal Society,2019,44(7):2227-2235.
[16] 徐良骥, 朱楠, 马荣振, 等. 厚松散承压含水层失水沉降模拟实验研究[J]. 采矿与安全工程学报,2015,32(5):821-826. XU Liangji, ZHU Nan, MA Rongzhen, et al. Water loss settlement simulation of thick unconsolidated confined aquifer layer[J]. Journal of Mining & Safety Engineering,2015,32(5):821-826.
[17] 周大伟. 煤矿开采沉陷中岩土体的协同机理及预测[D]. 徐州: 中国矿业大学, 2014. ZHOU Dawei. The synergy mechanism between rock mass and soil in mining subsidence and its prediction[D]. Xuzhou: China University of Mining and Technology, 2014.
-
期刊类型引用(0)
其他类型引用(5)
计量
- 文章访问数: 108
- HTML全文浏览量: 107
- PDF下载量: 18
- 被引次数: 5
