Damage characteristics and mechanical properties of superhigh-water material consolidated body under triaxial stress
-
摘要:
为研究超高水材料固结体在三向应力状态下的损伤特征及力学特性,进行不同养护时间(1、7、14、21、28 d)下超高水材料固结体的单轴压缩试验。借助PFC3D的平行黏结模型建立不同养护时间下的单轴压缩模型,得出与单轴试验强度相符合的5组PFC3D单轴压缩模拟试验强度。并统计各组平行黏结模型的细观物理力学参数,依据所得参数建立不同养护时间下超高水固结体的三轴压缩模型,施加与轴向应力相等的围压,记录不同养护时间下三轴试验过程中的应力-应变曲线及破坏时的力链分布。分析了超高水固结体在三向应力状态下的损伤特征,研究结果表明:①超高水材料固结体三向应力状态下的极限强度随养护时间的变化规律可由波尔茨曼方程表示。养护时间1~14 d时,极限强度增长最快;养护时间28 d时达到最大极限强度3.1 MPa。②超高水材料固结体三轴压缩模型力链贯通程度随养护时间的变化规律为:养护时间1~28 d,横向接触力链数目分别为4 006、4 561、4 891、5 017、5 062;纵向接触力链数目为4 029、4 439、4 716、4 917、5 123。表明超高水材料固结体的承载能力随养护时间的增加而增强,在养护时间0~14 d增长最快,14~28 d时承载能力趋于稳定。③以受拉力链模拟超高水材料固结体三轴压缩模型破坏过程中的裂隙发育过程,得出养护时间为1 d的超高水材料固结体破坏时的裂隙主要集中试样的上下两端。随着养护时间的增加,裂隙的发育程度增加,破坏时试样的中部开始出现裂隙,与上下两端的裂隙逐步贯通。
Abstract:In order to study the damage characteristics of superhigh-wate consolidated body under triaxial stress state, uniaxial experiments of superhigh-water consolidated body under different curing time (1, 7, 14, 21, 28 d) were carried out. Uniaxial compression models under different curing times were established by the parallel bonding model of PFC3D, and five groups of uniaxial strength of PFC3D simulation experiment were obtained, which were consistent with the uniaxial strength of the uniaxial experiment. Statistics each simulation experiment of uniaxial compression parallel bond, mesoscopic physical and mechanical parameters in the model on the basis of the parameters under different curing time, superhigh-water consolidation triaxial compression model body, the same confining pressure and axial stress is applied, record the triaxial experiments under different curing time in the process of stress-strain curve and the force when the damage distribution chain. Analyzed superhigh-water body of consolidation in the damage characteristics of three to the stress state, the results show that: ① the superhigh-water concretion body three to the stress state of the changing rule of the ultimate strength with curing time can be represented by the Bohr boltzmann equation. When curing time is 1 to 14 days, the ultimate strength increases fastest, and the maximum ultimate strength reaches 3.1 MPa when curing time is 28 days. ② The variation rule of the degree of penetration of force chains in the triaxial compression model of superhigh-water consolidated body with curing time is as follows: within curing time of 1-28 d, the number of transverse contact force chains is 4 006, 4 561, 4 891, 5 017, 5 062, respectively. The number of longitudinal contact force chains is 4 029, 4 439, 4 716, 4 917 and 5 123. The results show that the carrying capacity of the superhigh-water consolidated body increases with the increase of curing time, and increases fastest during the curing time from 0 to 14 days, and tends to be stable during 14 to 28 days. ③ The tensile chain was used to simulate the fracture development in the triaxial compression model of superhigh-water consolidated body. The results show that the cracks concentrate on the upper and lower ends of the specimen when the curing time is 1 d. With the increase of curing time, the cracks in the middle of the specimen begin to increase and finally connect with the cracks at the upper and lower ends of the specimen.
-
0. 引 言
煤矿井下巷道掘进作为煤炭开采的重要环节,相比于综采工作面工艺复杂、智能化水平较低,易造成“采掘失衡”,严重影响煤矿开采效率[1-3]。因此提高巷道掘进装备的智能截割、定向掘进及协同控制等智能化水平,是解决该问题的重要途径,关键在于解决掘进装备的位姿测量及精确定位。
针对复杂工况下掘进机的位姿测量,国内外研究学者做了大量研究,主要包括惯性导航方法、超宽带定位方法、全站仪定位方法、iGPS定位方法、视觉定位方法等方面。惯性定位方法[4-5]具有无源性、自主测量性强、环境适应性高的特点,能够在较短时间为掘进机提供较高精度的位姿数据,但受掘进机运动特性(慢速往复运动、作业时间长、剧烈振动)导致积分运算误差增大,位姿发散严重。超宽带定位方法[6-8]利用站点之间往返的超短脉冲信号实现测距定位,在有效范围内能够提供较高精度位姿信息,但需要多个基站配合,井下空间受限、设备多,信号易受遮挡衰减,测量稳定性需进一步提升。全站仪定位方法[9-10]多用于盾构机的作业定位具有较高定位精度,但煤矿井下巷道环境复杂、高粉尘水雾、空间狭小,全站仪发射的光路很容易被遮挡,无法实时获得掘进机的位置信息。如要获取掘进机的位姿信息必须在机身上安装多个棱镜通过几何关系得到掘进机的位姿信息,无法实现自主测量。iGPS定位方法[11-12]具有较高的实时定位能力,但需要合适的基站部署和精确的信号测量,以确保定位精度。但受粉尘、遮挡等因素,影响信号传播,导致总体定位精度不高,局限性较大。视觉定位方法[13-15]具有系统简单、成本低、非接触无累积误差等特点,多用于掘进设备的局部与全局定位,通常采用人为构造的几何特征标靶作为视觉识别目标,构建基于单目或多目位姿视觉解算模型,实现掘进装备的高精度位姿测量。但易受到巷道顶底板起伏、非结构化受限环境干扰等出现标靶遮挡造成视觉失效。
基于上述分析,单一的测量方式在复杂、恶劣的井下巷道环境下面临累积误差、性能退化、设备干扰、遮挡失效等限制,很难提供实时、稳定的掘进机位姿信息,因此融合多传感器信息的掘进机组合定位方法成为热点。LI等[16]采用紧耦合的方式将惯导和UWB系统组合,试验验证了该方法能够减小系统非视距误差,提高掘进机定位精度。张旭辉等[17]、马宏伟等[18]采用惯导与全站仪组合定位的方法对掘进机位姿进行长距离测量,依靠全站仪精确的位置输出对惯导补偿在一定程度上提高了掘进机的定位精度及自主性。HAN等[19]、王浩然等[20]、刘送永等[21]在惯导定位的基础上结合里程计输出的速度位置信息,采用滤波算法对两者数据进行融合处理实现掘进机的位姿测量,但里程计大都安装在掘进机履带上,受巷道底板环境影响履带容易出现打滑,导致里程计输出非掘进机实际运动状态。毛清华等[22]、张超等[23]、黄东等[24]利用惯导与视觉测量技术融合对掘进机实时定位,区别在于一般采用单/多目视觉、不同的视觉位姿测量模型或基于视觉环境信息(颜色、纹理等)的视觉测量系统位姿输出,以适应井下低照度、粉尘、多干扰等复杂环境。上述煤矿掘进装备组合定位大都采用“惯导+”的组合方式,借助惯导短时定位精度高、自主性强的特点,但受限于长时定位发散的缺点,综合其他传感器优势有效改善了惯性累积误差,提升了系统整体定位性能。
因此,结合团队前期研究基础,提出一种基于视觉与光纤惯导的掘进机组合定位方法,提高复杂环境下掘进机位姿测量可靠性。以三激光指向仪为特征的视觉测量系统、光纤惯导构建掘进机组合定位方案并阐述组合定位原理。介绍了视觉定位与惯导定位方法,通过卡尔曼滤波算法将惯导数据与视觉测量系统位姿数据进行融合实现掘进机的精确定位定向,为实现掘进机智能控制提供数据基础。
1. 基于视觉与惯导的掘进机组合定位方案
1.1 系统组成
掘进机组合定位系统组成如图1所示,主要由悬臂式掘进机、防爆相机、三激光指向仪标靶、光纤惯导和防爆计算机组成。光纤惯导和防爆计算机安装在导航控制箱内;三激光指向仪标靶通过安装架固定连接构成三角形安装在掘进机后方巷道顶板上朝向巷道掘进方向;防爆相机与掘进机机身刚性连接朝后拍摄三激光指向仪标靶;光纤惯导安装在防爆电控箱内;防爆计算机用来图像处理及数据解算。
1.2 组合定位原理
掘进机组合定位原理如图2所示,采用单目视觉测量方法,利用机载防爆相机采集三激光标靶图像,通过图像预处理、形态学变换、特征点提取得到相机坐标系下3个激光光斑的中心坐标,建立基于三点特征的掘进机位姿测量模型,结合矿用全站仪外参标定结果,得到了掘进机机身位姿信息。
光纤惯导通过自身的数学解算平台实时获取陀螺仪和加速度计采集的地球自转角速度分量、重力加速度分量,以惯性坐标系为参考坐标系,在惯性坐标系中测量出载体的角运动和线运动历程,以及重力加速度的方向变化情况,初始标定对准后精确估算出初始姿态矩阵,通过积分运算得到掘进机位置和姿态。融合前对数据进行有效值判断及异常值剔除后,采用卡尔曼滤波算法对这两者数据进行融合,实时输出掘进机的精准位姿信息,并利用融合后的位姿信息对惯导进行校正,实现掘进机的精确定位定向。
2. 掘进机机身位姿测量技术
2.1 系统坐标系定义
为了实现掘进机位姿测量,笔者建立如图3所示坐标系。
1)导航坐标系${o_{\mathrm{g}}}{x_{\mathrm{g}}}{y_{\mathrm{g}}}{z_{\mathrm{g}}}$。以掘进机重心为导航坐标系坐标原点${o_{\mathrm{g}}}$,${x_{\mathrm{g}}}$轴指向地理的东向,${y_{\mathrm{g}}}$轴指向地理北向,${z_{\mathrm{g}}}$轴与$ {x}_{{\mathrm{g}}}、{y}_{{\mathrm{g}}} $轴构成右手坐标系。导航坐标系用来表示惯导定位结果。
2)载体坐标系${o_{\mathrm{ri}}}{x_{\mathrm{ri}}}{y_{\mathrm{ri}}}{z_{\mathrm{ri}}}$。以掘进机重心为坐标原点${o_{\mathrm{ri}}}$,${x_{\mathrm{ri}}}$轴指向掘进机的左方,${y_{\mathrm{ri}}}$轴指向掘进机的前方,${z_{\mathrm{ri}}}$指向天与$ {x}_{{\mathrm{ri}}}、{y}_{{\mathrm{ri}}} $两轴构成右手坐标系。
3)相机坐标系${o}{x}{y}{z}$。相机坐标系以相机光心为坐标原点${o}$,${x}$轴指向相机右方,${y}$轴指向相机下方,${z}$轴指向相机前方。
4)巷道坐标系${o_{\mathrm{la}}}{x_{\mathrm{la}}}{y_{\mathrm{la}}}{z_{\mathrm{la}}}$。以巷道中线起点为坐标原点${o_{\mathrm{la}}}$,${x_{\mathrm{la}}}$轴指向巷道掘进方向的右方,${y_{\mathrm{la}}}$轴与巷道中线重合指向巷道掘进方向,${z_{\mathrm{la}}}$轴与$ {x}_{{\mathrm{la}}}、{y}_{{\mathrm{la}}} $垂直指向巷道顶板方向。为了更好的表征掘进机的运动状态,本文将掘进机的位姿信息均转化到了巷道坐标系。
5)地理坐标系${o_{\mathrm{ge}}}{x_{\mathrm{ge}}}{y_{\mathrm{ge}}}{z_{\mathrm{ge}}}$。地理坐标系以地球中心为坐标原点${o_{\mathrm{ge}}}$。${o_{\mathrm{ge}}}{x_{\mathrm{ge}}}$轴指向本初子午线与赤道的交点,${o_{\mathrm{ge}}}{y_{\mathrm{ge}}}$轴沿地球自转轴指向北极,${o_{\mathrm{ge}}}{z_{\mathrm{ge}}}$轴与其余两轴构成右手坐标系。
2.2 基于三激光指向仪标靶的掘进机位姿视觉测量方法
2.2.1 激光光斑中心定位
传统的光斑中心提取算法有灰度质心法[25]、Hough变换法[26]、高斯曲面拟合算法[27]等,这些算法在特定环境下能够有效的提取光斑中心,但煤矿井下掘进工作面存在粉尘、水雾、非平稳噪声,伴有矿灯、巷道照明灯、地测用激光指向仪等干扰光源,环境背景复杂,难以稳定、高精度的提取激光光斑中心。
考虑井下复杂干扰因素影响及激光指向仪的发光颜色,选用彩色防爆相机进行图像采集,对采集到的彩色图像经过图像处理后利用高斯曲面拟合确定激光光斑中心。根据激光指向仪光斑的HSV颜色空间即色度H、饱和度S、明度V作为颜色空间分量,通过设置参数范围可有效地滤除矿灯、照明灯等其他颜色杂光,区分出激光光斑的像素信息。闭运算操作将图像上光斑的区域边界向光斑中心填充,增大边界区域使光斑边界连续并标记出光斑位置。利用高斯曲面拟合方法得到图像上各区域的光斑中心,再通过求取各光斑之间的距离及比例,可以去除剩余与激光光斑特征不一致的杂光。激光光斑中心定位流程如图4所示,具体的光斑中心定位步骤:① 相机采集彩色光斑图像通过畸变矫正、高斯滤波等预处理后,遍历图像所有像素信息通过设置图像HSV参数,H的范围为[0,180],S的范围为[0,255],V的范围为[0,255],可以得到满足色度H在[60,120],饱和度S在[150,255],明度V在[30,180]的所有光斑$\sum {h\left( {i,j} \right)} $。② 对所有光斑$\sum {h\left( {i,j} \right)} $进行轮廓提取,去除大于激光光斑轮廓的杂光轮廓,剩余包括激光光斑在内的区域$\sum {g\left( {i,j} \right)} $。③ 对所有剩余光斑$\sum {g\left( {i,j} \right)} $进行闭运算使光斑区域连续平滑。④ 再利用高斯曲面拟合方法得到图像区域所有光斑的中心位置。⑤ 求取各光斑之间的距离D及比例S,保存满足给定的距离范围$\left( {{d_{\min }} \sim {d_{\max }}} \right)$及比例范围$\left( {{s_{\min }} \sim {s_{\max }}} \right)$光斑坐标,即激光指向仪光斑中心坐标${O_{p_1}}\left( {i,j} \right)$、${O_{p_2}}\left( {i,j} \right)$、${O_{p_3}}\left( {i,j} \right)$。在杂光背景下提取的激光光斑效果图如图5所示。
2.2.2 基于三激光点的掘进机视觉定位模型
如图6所示建立悬臂式掘进机位姿解算模型,已知空间3个特征点$ {P}_{1}、{P}_{2}、{P}_{3} $构成空间三角形,在相机像平面的投影分别为$ {p}_{1}、{p}_{2}、{p}_{3} $,O为相机光轴中心。
图中三角形存在对应关系:$\Delta O{P_1}{P_2} 与 \Delta O{p_1}{p_2}$、$\Delta O{P_1}{P_3} 与 \Delta O{p_1}{p_3}$、$\Delta O{P_2}{P_3} 与 \Delta O{p_2}{p_3}$,利用余弦定理有:
$$ \left\{\begin{array}{l} l_{{{{OP}}}_{1}}^{2}+l_{{{{OP}}}_{2}}^{2}-2l_{{{{OP}}}_{1}}l_{{{{OP}}}_{2}}\cos\;\alpha =l_{{{{P}}}_{1}{{{P}}}_{2}}^{2}\\ l_{{{{OP}}}_{1}}^{2}+l_{{{{OP}}}_{3}}^{2}-2l_{{{{OP}}}_{1}}l_{{{{OP}}}_{3}}\cos\;\beta =l_{{{{P}}}_{1}{{{P}}}_{3}}^{2}\\ l_{{{{OP}}}_{2}}^{2}+l_{{{{OP}}}_{3}}^{2}-2l_{{{{OP}}}_{2}}l_{{{{OP}}}_{3}}\cos\;\gamma =l_{{{{P}}}_{2}{{{P}}}_{3}}^{2}\end{array} \right. $$ (1) 为了简化计算在式(1)两边除以$O{P_3}^2$,并假设:
$$ \left\{ \begin{gathered} x = \frac{{{l_{{OP_1}}}}}{{l_{{{OP_3}}}}} \\ y = \frac{{{l_{{OP_2}}}}}{{{l_{{OP_3}}}}} \\ \end{gathered} \right. $$ (2) 将式(2)代入式(1),得:
$$ \left\{\begin{array}{l}{x}^{2}+{y}^{2}-2xy{\cos}\;\alpha =l_{{{{P}}}_{1}{{{P}}}_{2}}^{2}/l_{{{{OP}}}_{3}}^{2}\\ {x}_{1}{}^{2}+{y}^{2}-2y{\cos}\;\beta =l_{{{{P}}}_{1}{{{P}}}_{3}}^{2}/l_{{{{OP}}}_{3}}^{2}\\ {x}^{2}+1-2x{\cos}\;\gamma =l_{{{{P}}}_{2}{{{P}}}_{3}}^{2}/l_{{{{OP}}}_{3}}^{2}\end{array} \right. $$ (3) 对式(3),令
$$ \left\{ \begin{gathered} {{a = l_{{{{P}}_1}{{{P}}_2}}^2} \mathord{\left/ {\vphantom {{a = l_{{{{P}}_1}{{{P}}_2}}^2} {l_{{{OP}}_3}^2}}} \right. } {l_{{{OP}}_3}^2}} \\ ab = {{l_{{{{P}}_1}{{{P}}_3}}^2} \mathord{\left/ {\vphantom {{l_{{{{P}}_1}{{{P}}_3}}^2} {l_{{{OP}}_3}^2}}} \right. } {l_{{{OP}}_3}^2}} \\ ac = {{l_{{{{P}}_2}{{{P}}_3}}^2} \mathord{\left/ {\vphantom {{l_{{{{P}}_2}{{{P}}_3}}^2} {l_{{{OP}}_3}^2}}} \right. } {l_{{{OP}}_3}^2}} \\ \end{gathered} \right. $$ (4) 将式(4)带入式(3)中,得:
$$ \left\{\begin{array}{l}{x}^{2}+{y}^{2}-2xy\mathrm{cos}\;\alpha =a\\ {x}_{1}{}^{2}+{y}^{2}-2y\mathrm{cos}\;\beta =ab\\ {x}^{2}+1-2x\mathrm{cos}\;\gamma =ac\end{array} \right. $$ (5) 对式(5)进行化简,得:
$$ \left \{\begin{array}{l}\left(1-b\right){y}^{2}-b{x}^{2}-y\mathrm{cos}\;\beta +2bxy\mathrm{cos}\;\alpha +1=0\\ \left(1-c\right){x}^{2}-c{y}^{2}-x\mathrm{cos}\;\gamma +2cxy\mathrm{cos}\;\alpha +1=0\end{array}\right. $$ (6) 式中:$ \mathrm{cos}\;\alpha 、\mathrm{cos}\;\beta 、\mathrm{cos}\;\lambda $可以通过激光点的图像坐标计算得到,结合矿用全站仪标定获得的$ {P}_{1}、 {P}_{2}、{P}_{3} $空间点在巷道坐标系下的坐标$ {{\mathbf{P}}_{{\mathrm{w}}i}} = [ {x_{{\mathrm{w}}i}}, {y_{{\mathrm{w}}i}}, {z_{{\mathrm{w}}i}} ]^{\text{T}} \left( {i = 1,2,3} \right) $计算得到b、c。x、y作为未知数随着相机的位姿改变而发生变化。
因此对式(6)采用吴消元法求解,该方程最多可以得到4个解$ ({x}_{g},{y}_{g}) $,$ g = 1,2,3,4 $,即可以得到$ {P}_{1} $、${P}_{2} $、${P}_{3} $在相机坐标系下4组坐标$ {{\mathbf{P}}_{{\mathrm{c}}i}} = {\left[ {x_i^m,y_i^m,z_i^m} \right]^{\text{T}}} \left( {i = 1,2,3;m = 1,2,3,4} \right) $。
根据相机透视投影模型,相机坐标系与巷道坐标系间的旋转矩阵R和平移向量t可以表示为
$$ {{\boldsymbol{P}}_{{\mathrm{ci}}}} = {\boldsymbol{R}}{{\boldsymbol{P}}_{{\mathrm{w}}i}} + {\boldsymbol{t}} $$ (7) 根据旋转矩阵R和平移向量t,得到防爆相机在巷道坐标系下的位置和航向角、俯仰角、横滚角信息。但其中只有一组姿态解是最优解,寻求最优解将其转化为最小化重投影误差的问题,建立最小化重投影误差函数为
$$ \mathop \sum \limits_{i = 1}^n \mathop \sum \limits_{j = 1}^m {\left\| {{p_i} - \frac{1}{s}{\boldsymbol{K}}\left( {{{\boldsymbol{R}}_j}{P_i} + {{\boldsymbol{t}}_i}} \right)} \right\|^2} $$ (8) 式中:${p_i}$为第i个激光点在图像坐标系的图像坐标;s为投影像素点坐标转化为齐次坐标的系数;K为相机的内参矩阵;${{\boldsymbol{R}}_j},{{\boldsymbol{t}}_j}$分别为式(7)得到的第j个旋转矩阵和平移矩阵。最后再结合系统初始标定的掘进机机身坐标系与相机坐标系之间的转换关系,得到掘进机机身在巷道坐标系下的位姿信息。
2.3 基于光纤惯导的掘进机定位方法
2.3.1 姿态更新
为了便于计算本文采用四元数法[28]进行掘进机姿态角计算。四元数可定义为
$$ {\boldsymbol{Q}} = {q_0} + {q_1}{\boldsymbol{i}} + {q_2}{\boldsymbol{j}} + {q_3}{\boldsymbol{k}} $$ (9) 式中:${q_0},{q_1},{q_2},{q_3}$为实数;i, j, k为相互正交的单位向量。机身坐标系到地理坐标系的姿态转化矩阵${\boldsymbol{C}}_{\mathrm{o}}^n$用四元数可表示为
$$ {\boldsymbol{C}}_{\mathrm{o}}^n = \left[ {\begin{array}{*{20}{c}} {q_0^2 + q_1^2 - q_2^2 - q_3^2} & {2\left( {{q_1}{q_2} - {q_0}{q_3}} \right)} & {2\left( {{q_1}{q_3} + {q_0}{q_2}} \right)} \\ {2\left( {{q_1}{q_2} - {q_0}{q_3}} \right)} & {q_0^2 - q_1^2 + q_2^2 - q_3^2} & {2\left( {{q_2}{q_3} + {q_0}{q_1}} \right)} \\ {2\left( {{q_1}{q_3} + {q_0}{q_2}} \right)} & {2\left( {{q_2}{q_3} - {q_0}{q_1}} \right)} & {q_0^2 - q_1^2 - q_2^2 + q_3^2} \end{array}} \right] $$ (10) 姿态更新周期为${T_k} = {t_{k + 1}} - {t_k}$ 时,姿态更新可表示为
$$ {\boldsymbol{Q}}_{\mathrm{o}}^b\left( {{{{t}}_{k + 1}}} \right) = \frac{1}{2}{\boldsymbol{Q}}_{\mathrm{o}}^b\left( {{{{t}}_k}} \right) \otimes {\boldsymbol{Q}}_{\mathrm{o}}^b\left( {{{{T}}_k}} \right) $$ (11) 得到姿态转化矩阵后即可求解掘进机在地理坐标系下的姿态角:
$$ \begin{array}{*{20}{l}} {\alpha = \arcsin \left[ {2\left( {{q_2}{q_3} - {q_0}{q_1}} \right)} \right]} \\ {\beta = \arctan \left[ { - \dfrac{{2\left( {{q_2}{q_3} - {q_0}{q_1}} \right)}}{{q_0^2 - q_1^2 - q_2^2 + q_3^2}}} \right]} \\ {\gamma = \arctan \left[ { - \dfrac{{2\left( {{q_1}{q_3} + {q_0}{q_2}} \right)}}{{q_0^2 - q_1^2 + q_2^2 - q_3^2}}} \right]} \end{array} $$ (12) 2.3.2 速度更新
将掘进机加速度投影至导航坐标系下,可得到掘进机速度更新方程:
$$ {\dot V^n} = C_{\mathrm{o}}^n{f^{\mathrm{o}}} - \left( {2\omega _{{\mathrm{ie}}}^n + \omega _{{\mathrm{en}}}^n} \right) \times {V^n} + {g^n} $$ (13) 其中
$$ \begin{array}{*{20}{c}} {\omega _{{\mathrm{ie}}}^n = {{\left[ {\begin{array}{*{20}{c}} 0&{{\omega _{{\mathrm{ie}}}}\cos\; L}&{{\omega _{{\mathrm{ie}}}}\sin \;L} \end{array}} \right]}^{\text{T}}}} \\ {{\mathbf{\omega }}_{{\mathrm{en}}}^n = {{\left[ { - \dfrac{{{V_{\text{N}}}}}{{{R_{\mathrm{M}}} + h}}\quad \dfrac{{{V_{\text{E}}}}}{{{R_{\mathrm{N}}} + h}}\quad \dfrac{{{V_{\text{E}}}}}{{{R_{\mathrm{N}}} + h}}\tan\; L} \right]}^{\text{T}}}} \end{array} $$ (14) 式中:${f^{\mathrm{o}}}$为加速度计测量值,m/s2;$ \omega _{{\mathrm{ie}}}^n $为地球坐标系相对于地心坐标系的旋转角速度,rad/s;$\omega _{{\mathrm{en}}}^n$为导航坐标系相对于地球坐标系的旋转角速度,rad/s;${g^n}$为导航坐标系下重力加速度,m/s2;${{\boldsymbol{V}}^{n}} = \left[ {{V_{\mathrm{E}}},{V_{\mathrm{N}}},{V_{\mathrm{Z}}}} \right]$, ${V_{\mathrm{E}}},{V_{\mathrm{N}}},{V_{\mathrm{Z}}}$为东向、北向、天向的速度,m/s;$L$为地理纬度,(°);h为高度,m;再去除有害加速度$ \left( {2\omega _{{\mathrm{ie}}}^n + \omega _{{\mathrm{en}}}^n} \right) \times {{\boldsymbol{V}}^n} + {g^n} $才能获得掘进机的加速度。
2.3.3 位置更新
根据速度方程,获得惯导的位置微分方程后再积分得到掘进机的位置:
$$ \left\{ {\begin{array}{*{20}{l}} {{L_{{{s}} + 1}} = {L_{{s}}} + \displaystyle\int_{{t_{{s}}}}^{{t_{{{s}} + 1}}} {\dfrac{{{V_{\mathrm{N}}}}}{{{R_{\mathrm{M}}} + h}}{\mathrm{d}}t} } \\ {{\lambda _{{{s}} + 1}} = {\lambda _{{s}}} + \displaystyle\int_{{t_{{s}}}}^{{t_{{{s}} + 1}}} {\dfrac{{{V_{\mathrm{E}}}}}{{\left( {{R_{\mathrm{N}}} + h} \right)\cos\; L}}{\mathrm{d}}t} } \\ {{h_{{{s}} + 1}} = {h_{{s}}} + \displaystyle\int_{{t_{{s}}}}^{{t_{{{s}} + 1}}} {{V_{\mathrm{Z}}}{\mathrm{d}}t} } \end{array}} \right. $$ (15) 式中:${L_{{{s}} + 1}}$为${t_{{{s}} + 1}}$时刻掘进机在地理坐标系下的经度,(°);${\lambda _{{{s}} + 1}} $为${t_{{{s}} + 1}} $时刻掘进机在地理坐标系下的纬度,(°);${h_{{{s}} + 1}} $为${t_{{{s}} + 1}} $时刻掘进机在地理坐标系下的高度,m;${L_{{s}}}$为掘进机在${t_{{s}}}$时刻掘进机在地理坐标系下的经度,(°);${\lambda _{{s}}} $为${t_{{s}}} $时刻掘进机在地理坐标系下的纬度,(°);${h_{{s}}} $为${t_{{s}}} $时刻掘进机在地理坐标系下的高度,m;${R_{\mathrm{M}}}、{R_{\mathrm{N}}}$分别为地球卯酉圈和子午圈的曲率半径,m。
通过惯导解算得到的掘进机位姿均是在导航坐标系下的位姿,掘进机位姿表示通常是在巷道坐标系下,掘进机在导航坐标系的位姿表示为${{\boldsymbol{M}}_{\mathrm{n}}}$,假设导航坐标系与巷道坐标系的转换矩阵为${{\boldsymbol{M}}_{{\mathrm{nh}}}}$,即掘进机在巷道坐标系的位姿可表示为
$$ {{\boldsymbol{M}}_{\mathrm{h}}} = {{\boldsymbol{M}}_{{\mathrm{nh}}}} {{\boldsymbol{M}}_{\mathrm{n}}} $$ (16) 3. 基于卡尔曼滤波的掘进机位姿数据融合
采用卡尔曼滤波算法对光纤惯导和视觉测量位姿数据进行融合,在融合前对各传感器数据进行有效性判断,并在时间实现配准,运用卡尔曼滤波融合算法融合视觉与惯导数据,最终得到掘进机的位姿信息,同时采用融合后的位姿数据对惯导进行补偿校正。
3.1 数据时间配准
视觉测量与惯导测量两种测量方式的采样频率及计算速度各不相同,在数据融合时要考虑两种测量结果在时间上同步。通常采用插值方法进行时间配准,为使插值误差尽可能小,根据插值原则插值节点应尽可能位于插值区间中部,由高采样率传感器测得的数据向低采样率传感器测得数据进行时间配准。本文采用拉格朗日三点插值[29]实现惯导位姿数据向视觉位姿数据进行时间配准。根据插值原理,假设$ {t_{n - 1}},{t_n},{t_{n + 1}} $时刻传感器测量数据分别为${x_{n - 1}}, {x_n},{x_{n + 1}}$由于掘进机的移动速度很小,因此可以将$ {t_{n - 1}},{t_n},{t_{n + 1}} $看作是等间距的,若在时间$t = {t_n} + \Delta t$处插值,运用拉格朗日三点插值可以计算出t时刻的测量值为
$$ \begin{split} x_t= & \frac{\left(t-t_n\right)\left(t-t_{n-1}\right)}{\left(t_{n-1}-t_n\right)\left(t_{n-1}-t_{n+1}\right)} x_{n-1} +\frac{\left(t-t_{n-1}\right)\left(t-t_{n+1}\right)}{\left(t_n-t_{n-1}\right)\left(t_n-t_{n+1}\right)} x_n +\\& \frac{\left(t-t_{n-1}\right)\left(t-t_n\right)}{\left(t_{n+1}-t_{n-1}\right)\left(t_{n+1}-t_n\right)} x_{n+1} \end{split} $$ (17) 惯导的采样周期为0.05 s,视觉系统的采样周期为0.1 s,所有位姿数据都带有时间戳,如图7所示,假设利用拉格朗日三点插值方法将惯导数据向视觉系统的550 ms点处进行时间配准,惯导在该点附近的数据为550、600、650 ms点,提取时间$ t = 550{\text{ ms}} $,$ t_{n-1}=550 \mathrm{\;ms} $,$t_{n}=600 \mathrm{\;ms} $,$ t_{n+1}=650 \mathrm{\;ms} $,根据式(17)完成两者数据的时间配准。
3.2 误差方程
3.2.1 视觉测量系统
视觉测量给出的掘进机位姿误差多由激光光斑提取误差及位姿测量误差引起,可以在图像处理算法、视觉硬件等方面进一步优化,减少视觉测量系统误差。虽然上述方法在一定程度上可以改善误差对测量系统的影响,但仍不能完全消除误差。因此,建立光斑中心定位与特征点三维坐标之间的误差传递函数,根据视觉测量误差理论,对光斑中心坐标求导,按泰勒级数展开并忽略高次项,得到特征点图像坐标与三维坐标测量误差为
$$ \delta {X_{{\mathrm{w}}i}} = \frac{{\partial {X_{{\mathrm{w}}i}}}}{{\partial {u_i}}}\delta {u_i} + \frac{{\partial {X_{{\mathrm{w}}i}}}}{{\partial {v_i}}}\delta {v_i} $$ (18) $$ \delta {Y_{{\mathrm{w}}i}} = \frac{{\partial {Y_{{\mathrm{w}}i}}}}{{\partial {u_i}}}\delta {u_i} + \frac{{\partial {Y_{{\mathrm{w}}i}}}}{{\partial {v_i}}}\delta {v_i} $$ (19) 式中:$ \delta u_{1} $,$ \delta v_{i} $为光斑图像坐标提取误差;i为特征点个数,i=1,2,3。由式(18)、式(19)得到特征点图像与三维坐标之间的误差传递关系为
$$ \delta {\boldsymbol{W}} = {{\boldsymbol{E}}_{{\mathrm{w}}P}}\delta {\boldsymbol{P}} $$ (20) 式中:$ {{\boldsymbol{E}}_{{\mathrm{w}}P}} = {\left[ \begin{gathered} \frac{{\partial {X_{{\mathrm{w}}1}}}}{{\partial {u_1}}}\quad \frac{{\partial {X_{{\mathrm{w}}1}}}}{{\partial {v_1}}}\quad \cdots \quad \frac{{\partial {X_{{\mathrm{w}}1}}}}{{\partial {u_3}}}\quad \frac{{\partial {X_{{\mathrm{w}}1}}}}{{\partial {v_3}}} \\[-2pt] \quad \vdots \quad \quad \;\; \vdots \quad \quad \vdots \quad \quad \;\; \vdots \quad \quad \;\; \vdots \\[-2pt] \frac{{\partial {Y_{{\mathrm{w}}3}}}}{{\partial {u_1}}}\quad \frac{{\partial {Y_{{\mathrm{w}}3}}}}{{\partial {v_1}}}\quad \; \cdots \quad \frac{{\partial {Y_{{\mathrm{w}}3}}}}{{\partial {u_3}}}\quad \frac{{\partial {Y_{{\mathrm{w}}3}}}}{{\partial {v_3}}} \\ \end{gathered} \right)_{6 \times 6}} $为光斑图像坐标系与世界坐标系之间的误差传递模型;$\delta {\boldsymbol{P}} = {\left( {\delta {u_1},\delta {v_1},\delta {u_2},\delta {v_2},\delta {u_3},\delta {v_3}} \right)^{\text{T}}}$为特征点图像坐标误差,px;$\delta {\boldsymbol{W}} = {\left( {\delta {X_{{\mathrm{w}}1}},\delta {Y_{{\mathrm{w}}1}},\delta {X_{{\mathrm{w}}2}},\delta {Y_{{\mathrm{w}}2}},\delta {X_{{\mathrm{w}}3}},\delta {Y_{{\mathrm{w}}3}}} \right)^{\text{T}}}$为世界坐标测量误差,m。
根据视觉位姿解算模型,位姿参数与光斑的三维坐标有关,可以表示为
$$ \begin{split} f( {t_x},{t_y},{t_z},{\theta _x}, & {\theta _y},{\theta _z},{X_{{\mathrm{w}}1}},{Y_{{\mathrm{w}}1}},{Z_{{\mathrm{w}}1}},{X_{{\mathrm{w}}2}},{Y_{{\mathrm{w}}2}},{Z_{{\mathrm{w}}2}},\\&{X_{{\mathrm{w}}3}}, {Y_{{\mathrm{w}}3}},{Z_{{\mathrm{w}}3}} ) = 0 \end{split} $$ (21) 对式(21)进行泰勒展开,利用最小二乘法得到三维坐标与位姿解算的误差传递模型:
$$ \delta {\boldsymbol{q}} = {\left( {{\boldsymbol{J}}_{\mathrm{f}}^{\text{T}}{{\boldsymbol{J}}_{\mathrm{f}}}} \right)^{ - 1}}{\boldsymbol{J}}_{\mathrm{f}}^{\text{T}}\delta {\boldsymbol{W}} $$ (22) 式中:J为式(21)的雅可比矩阵,令$\delta {{\boldsymbol{E}}_P} = {\left( {{\boldsymbol{J}}_{\mathrm{f}}^{\text{T}}{{\boldsymbol{J}}_{\mathrm{f}}}} \right)^{ - 1}}{\boldsymbol{J}}_{\mathrm{f}}^{\text{T}}$,结合式(20)得到视觉测量系统误差为
$$ \delta {\boldsymbol{E}} = {{\boldsymbol{E}}_{\mathrm{P}}}{{\boldsymbol{E}}_{{\mathrm{wP}}}}\delta {\boldsymbol{W}} $$ (23) 3.2.2 惯导系统
受惯导本身测量原理及外界因素干扰,导致系统出现测量误差或偏差,致使系统性能下降。通常由外界因素如振动、电磁干扰等影响,可采取物理减震、隔离、屏蔽等措施;由原理性导致的可采用滤波技术、定时校准、校正等措施,进一步降低误差,但仍不能完全消除。因此,根据文献[30]建立惯导的误差方程为:
$$ \dot \phi = \phi \times \omega _{{\mathrm{in}}}^n + \delta \omega _{{\mathrm{in}}}^n - {\varepsilon ^n} $$ (24) $$ \begin{split} \delta {\dot V^n} = & - \phi {f^n} - \left( {2\omega _{{\mathrm{ie}}}^n + \omega _{{\mathrm{en}}}^n} \right) \delta {V^n} - \\& \left( {2\delta \omega _{{\mathrm{ie}}}^n + \delta \omega _{{\mathrm{en}}}^n} \right) {V^n} + {\nabla ^n} \end{split} $$ (25) $$ \left\{ \begin{gathered} \delta \dot L = \frac{{\delta {V_{\mathrm{N}}}}}{{{R_{\mathrm{M}}} + h}} - \delta h\frac{{{V_{\mathrm{N}}}}}{{{{\left( {{R_{\mathrm{M}}} + h} \right)}^2}}} \\ \delta \dot \lambda = \frac{{\delta {V_{\mathrm{E}}}}}{{{R_{\mathrm{N}}} + h}}\sec \;L + \delta L\frac{{{V_{\mathrm{E}}}}}{{\left( {{R_{\mathrm{N}}} + h} \right)}}\tan\; L\sec\; L - \\\qquad \delta h\frac{{{V_{\mathrm{E}}}\sec\; L}}{{{{\left( {{R_{\mathrm{N}}} + h} \right)}^2}}} \\ \delta \dot h = \delta {V_{\mathrm{Z}}} \\ \end{gathered} \right. $$ (26) 式中:$ \phi $为姿态误差,(°);$ \omega _{{\mathrm{in}}}^n $为导航坐标系相对于地心坐标系的旋转角速度,rad/s;$ \delta \omega _{{\mathrm{in}}}^n $为$ \omega _{{\mathrm{in}}}^n $的误差,rad/s;$ {\varepsilon ^n} $为陀螺仪零偏误差,rad/s;$ \delta {V^n} $为速度误差,m/s;$ {f^n} $为导航坐标系下加速度计测得的比力,m/s2;$ {{\boldsymbol{V}}^n} $为惯导得到的导航坐标系下掘进机速度,m/s;$ {\nabla ^n} $为加速度计的零偏误差,m/s2;$ \delta L $为纬度误差,(°);$ \delta \lambda $为经度误差,(°);$ \delta {{h}} $高度误差,m;$ \delta {V_{\mathrm{N}}} $,$ \delta {V_{\mathrm{E}}} $,$ {V_{\mathrm{Z}}} $别为北向、东向、天向速度误差,m/s。
3.3 卡尔曼滤波器设计
采用卡尔曼滤波器对惯导数据和视觉数据进行融合,系统的状态向量表示为
$$ \begin{split} {\boldsymbol{X}} = & [ {\phi _{\mathrm{N}}},{\phi _{\mathrm{E}}},{\phi _{\mathrm{Z}}},\delta {v_{\mathrm{N}}},\delta {v_{\mathrm{E}}},\delta {v_{\mathrm{Z}}},\delta L,\delta \lambda ,\delta h,\\& {\varepsilon _{\mathrm{N}}},{\varepsilon _{\mathrm{E}}},{\varepsilon _{\mathrm{Z}}},{\nabla _{\mathrm{N}}},{\nabla _{\mathrm{E}}},{\nabla _{\mathrm{Z}}},\delta {\boldsymbol{E}} ]^{\text{T}} \end{split} $$ (27) 式中:$ {\phi _{\mathrm{N}}},{\phi _{\mathrm{E}}},{\phi _{\mathrm{Z}}} $为姿态角误差,(°);$ \delta {v_{\mathrm{N}}},\delta {v_{\mathrm{E}}},\delta {v_{\mathrm{Z}}} $为速度误差,m/s;$ {\varepsilon _{\mathrm{N}}},{\varepsilon _{\mathrm{E}}},{\varepsilon _{\mathrm{Z}}} $为3个轴向的陀螺仪漂移,rad/s;$ {\nabla _{\mathrm{N}}},{\nabla _{\mathrm{E}}},{\nabla _{\mathrm{Z}}} $为3个轴向的加速度漂移,m/s2;$\delta {\boldsymbol{E}}$为视觉测量系统误差,m。
以视觉测量和惯导解算的掘进机位置与姿态数据之差作为系统的观测量为
$$ {\boldsymbol{Z}} = \left[ \begin{gathered} {\phi _{{\mathrm{INS}}}} - {\phi _{\mathrm{C}}} \\ {{\boldsymbol{P}}_{{\mathrm{INS}}}} - {{\boldsymbol{P}}_{\mathrm{C}}} \\ \end{gathered} \right] $$ (28) 系统卡尔曼滤波状态方程和观测方程为
$$ \left\{ \begin{gathered} {{\boldsymbol{X}}_k} = {\boldsymbol{A}}{{\boldsymbol{X}}_{k - 1}} + {\boldsymbol{\varGamma }}{{\boldsymbol{W}}_{k - 1}} \\ {{\boldsymbol{Z}}_k} = {\boldsymbol{H}}{{\boldsymbol{X}}_k} + {{\boldsymbol{V}}_k} \\ \end{gathered} \right. $$ (29) 式中:${{\boldsymbol{X}}_k}$为$k$时刻的状态量;${\boldsymbol{A}}$为状态转移矩阵;${{\boldsymbol{X}}_{k - 1}}$为$k - 1$时刻的状态量;${\boldsymbol{\varGamma }}$为噪声驱动矩阵;${{\boldsymbol{W}}_k}$为$k - 1$时刻的输入噪声;${\boldsymbol{H}}$为量测矩阵;${{\boldsymbol{Z}}_k}$为$k$时刻的观测值;${{\boldsymbol{V}}_k}$为$k$时刻的观测噪声。
已知卡尔曼滤波状态方程和观测方程后,就可以对下一时刻状态进行预测
$$ \begin{gathered} {{{\boldsymbol{\hat X}}}^ - }_k = {\boldsymbol{A}}{{{\boldsymbol{\hat X}}}^{}}_{k - 1} \\ {\boldsymbol{P}}_k^ - = {\boldsymbol{A}}{{\boldsymbol{P}}_{k - 1}}{{\boldsymbol{A}}^T} + {\boldsymbol{D}} \\ \end{gathered} $$ (30) 式中:$ {{\boldsymbol{\hat X}}^ - }_k $为$k$时刻的先验估计值(预测值);$ {{\boldsymbol{\hat X}}^ - }_{k - 1} $为$k - 1$时刻的后验估计值(最优估计);$ {\boldsymbol{P}}_k^ - $为先验估计的协方差矩阵;$ {{\boldsymbol{P}}_{k - 1}} $为$k - 1$时刻的后验估计协方差矩阵;${\boldsymbol{D}}$为噪声协方差矩阵。
对卡尔曼滤波器状态进行更新
$$ \begin{gathered} {{\boldsymbol{{ K}}}_k} = {\boldsymbol{P}}_k^ - {{\boldsymbol{H}}^{\rm T}}{\left( {{\boldsymbol{HP}}_k^ - {{\boldsymbol{H}}^{\text{T}}} + {{\boldsymbol{R}}_{\mathrm{e}}}} \right)^{ - 1}} \\ {{{\boldsymbol{\hat X}}}_k} = {\boldsymbol{\hat X}}_k^ - + {{\boldsymbol{K}}_k}\left( {{{\boldsymbol{Z}}_k} - {\boldsymbol{H\hat X}}_k^ - } \right) \\ {{\boldsymbol{P}}_k} = \left( {{\boldsymbol{I}} - {{\boldsymbol{K}}_k}{\boldsymbol{H}}} \right){\boldsymbol{P}}_k^ - \\ \end{gathered} $$ (31) 式中:$ {{\boldsymbol{K}}_k} $为卡尔曼滤波增益;$ {\boldsymbol{P}}_k^ - $为$k$时刻预测误差协方差矩阵;$ {{\boldsymbol{R}}_{\mathrm{e}}} $为观测噪声协方差矩阵;$ {{\boldsymbol{\hat X}}_k} $为$k$时刻的后验估计值;${\boldsymbol{I}}$为单位矩阵。
3.4 传感器输出有效性判断及异常值剔除
为解决视觉测量系统因视觉标靶短时遮挡造成的数据丢失及各传感器出现的异常值(以下均成称为无效数据)对系统精度、稳定性的影响,在两者数据融合前对系统中每个传感器i分别进行预测估计$ {}^i{{\boldsymbol{\hat X}}^ - }_k = {\boldsymbol{A}} {}^i{{\boldsymbol{\hat X}}^{}}_{k - 1} $,可得传感器i的预测观测值
$$ {}^i{{\boldsymbol{\hat Z}}_k} = {\boldsymbol{H}} {}^i{{\boldsymbol{\hat X}}_k} $$ (32) 式中:${}^i{{\boldsymbol{\hat Z}}_k}$为$k$时刻传感器i的预测观测值。给定阈值${\varepsilon _1}$,${\varepsilon _2}$$\left( {{\varepsilon _1} < {\varepsilon _2}} \right)$,若新的观测${}^i{{\boldsymbol{Z}}_k}$满足:
$$ {\varepsilon _1} < \left| {{}^i{{\boldsymbol{Z}}_k} - {}^i{{{\boldsymbol{\hat Z}}}_k}} \right| < {\varepsilon _2} $$ (33) 则认为观测值${}^i{{\boldsymbol{Z}}_k}$数据有效。不在该范围内则认为观测值${}^i{{\boldsymbol{Z}}_k}$数据无效,在判断确定无效数据的传感器后,系统只对有效数据的传感器计算局部估计及其观测量进行状态更新,即将无效传感器对应的卡尔曼滤波增益对应元素设置为0,以此剔除系统中的无效数据,降低该传感器对融合结果的影响,确保系统的稳定性和鲁棒性。
4. 试验验证
4.1 试验平台搭建与试验设计
为了验证笔者提出的掘进机组合定位方法的有效性及组合定位系统功能,搭建悬臂式掘进机位姿测量实验平台如图8所示,试验平台包括移动履带机器人、三激光指向仪标靶、防爆相机、光纤惯导、防爆计算机、数字全站仪。用烟雾制造机模拟井下粉尘环境,以第三方测量仪器数字全站仪测量的掘进机机身位姿为真实值,构建全站仪统一全局坐标系(X轴指向右、Y轴指向上、Z轴指向履带机器人前进方向)。试验验证设置如下:① 基于三激光指向仪标靶的掘进机位姿视觉测量精度评价试验,由于组合定位系统位置基准由视觉测量系统提供,精确的视觉定位数据是组合定位系统数据有效融合的保证。② 视觉定位、惯导定位、组合定位系统定位结果对比试验,进一步验证本文提出的视觉短时失效时有效值判断及处理方法的有效性。
4.2 基于三激光指向仪标靶的掘进机位姿视觉测量
试验过程中首先对三激光指向仪标靶完成测量,将视觉测量系统得到的位姿转化到全站仪坐标系下,以全站仪直接测量的履带机器人位姿为真实值,对比视觉测量系统得到的测量值。在楼道移动履带机器人模拟掘进机在巷道内运动,利用防爆相机采集三激光指向仪标靶图像,规定履带机器人距离三激光标靶初始距离为10,随后每次移动10距离,测量距离分别为20、30、40、50、60 m。并在每个位置采集不少于30张激光标靶图像,在计算机中完成图像处理、激光光斑中心定位,计算得到掘进机机身位姿,同时对比全站仪测量的履带机器人真实位姿,对掘进机视觉测量系统精度进行验证。
履带机器人在楼道移动过程中机身位姿测量值(视觉测量结果)和真实值(全站仪测量结果)见表1。结果表明:本文算法能够在60 m测试范围内实现机身位姿测量。不同距离下测量误差结果如图9所示,结果表明在60 m测量范围内本文所述的悬臂式掘进机位置视觉测量方法计算得到的掘进机机身X,Y,Z 3个轴方向上的误差均在35 mm以内,掘进机机身姿态角误差均在0.5°以内。
表 1 不同距离时视觉位姿测量结果Table 1. Visual pose measurement results at different distances项目 距离/m x/mm y/mm z/mm 俯仰角/
(°)横滚角/
(°)航向角/
(°)测量值 10 −457 567 9997 0.56 2.91 4.93 20 −471 544 20022 1.80 2.80 4.84 30 −564 561 29983 −1.19 2.38 5.03 40 −516 563 40054 0.07 2.96 4.17 50 −535 541 50071 0.81 3.06 5.09 60 −488 546 59975 −2.09 3.01 4.16 真实值 10 −452 556 10005 0.69 2.76 4.76 20 −480 559 20009 1.66 2.59 4.62 30 −551 542 30005 −1.36 2.61 4.70 40 −537 540 40027 0.26 2.58 4.54 50 −559 512 50041 0.54 2.61 4.68 60 −516 514 60009 −1.74 2.54 4.58 4.3 对比试验
在位姿测量试验开始前,惯导根据全站仪给出的起始点坐标进行对准,由于试验场地受限无法长距离运行,因此将惯导静置20 min后开始试验,模拟掘进机长时间作业运动,视觉系统实时采集履带机器人不同位姿时的标靶图像,通过位姿解算模型计算获得履带机器人实时位姿信息。通过卡尔曼滤波器对视觉与惯导数据融合,利用融合后的数据对惯导进行校正与补偿。
试验中操作履带机器人从规定的起点沿着楼道模拟掘进机向前行进,并伴有转向动作,通过全站仪跟踪测量(X轴指向右、Y轴指向上、Z轴指向履带机器人前进方向)向前行进约50 m,期间采用惯导、视觉定位、组合定位方法分别对履带机器人进行定位测量,并在履带机器人移动过程中某一时刻短时遮挡标靶,模拟实际掘进机运行过程中视觉标靶遮挡的情况。机器人运动轨迹及不同方法定位结果如图10所示,可以看出,随着履带机器人的前进,惯导定位由于积分累计原因定位偏差逐渐增大,视觉系统定位结果与全站仪测量趋势一致,组合定位结果更接近全站仪得到的实际位移,位置最大偏差在30 mm以内,并且视觉标靶短时遮挡时组合定位结果出现波动但偏差仍在允许范围内。
4.4 组合定位系统现场应用
基于本方法研发的掘进机组合定位系统在山西某矿胶带运输巷掘进工作面进行了工业性试验。试验设备包括EBZ260型悬臂式掘进机、防爆相机、光纤惯导系统、计算机、三激光指向仪标靶,系统硬件及安装如图11所示。
井下应用时,为保证组合定位系统的功能,根据巷道条件及各设备的安装位置,将三激光指向仪标靶通过固定架安装在巷道顶板上并指向巷道掘进方向,根据井下防爆要求将光纤惯导固定在特制减震架上与计算机共同安装在导航防爆控制箱内。防爆相机安装在掘进机机尾部靠右侧保证激光指向仪在相机视角范围内,并将实时采集的标靶图像通过网线传输至导航控制箱计算机完成图像处理及位姿解算。截取掘进机作业过程中部分位姿数据如图12所示,分析截取的
2700 个位姿数据,可以看出,组合定位系统输出稳定,能够实时测量掘进的位姿信息。在井下近2个月测试运行期间系统稳定运行,机身定位误差20 m范围内在30 mm以内,姿态角误差均在0.5°以内,实现了掘进机位姿实时测量,满足巷道掘进施工精度要求[31],并能为掘进机智能控制提供准确可靠的位姿数据。
5. 结 论
1)为解决煤矿复杂巷道环境下掘进装备单一位姿测量方法存在适应性、稳定性差的问题,提出一种基于视觉与光纤惯导融合的悬臂式掘进机组合定位方法。
2)构建基于三激光指向仪标靶的悬臂式掘进机位姿视觉测量模型,推导了视觉系统与光纤惯导的误差方程,设计基于卡尔曼滤波的视觉与光纤惯导数据融合算法,利用融合后的位姿数据对惯导进行补偿校正,消除了惯导长时间积分运算导致的累积误差。
3)提出基于预测传感器观测向量的有效值判断及异常值剔除方法,在数据融合前对传感器输出进行有效性判断,解决了因视觉标靶短时遮挡导致数据丢失对融合结果的影响,提高了定位系统的稳定性。
4)在某矿井下进行了组合定位现场应用,结果表明组合定位方法得到的掘进机机身位置误差总体在30 mm以内,态角误差均在0.5°以内,且能持续稳定输出,满足巷掘进施工及掘进机智能截割数据精度要求。
-
表 1 颗粒细观参数
Table 1 Particle microparameter
养护
时间/d最小颗粒
半径/mm粒径比 孔隙率 阻尼 密度/
(kg∙m−3)1 0.1 1.5 0.35 0.7 2 600 7 0.1 1.5 0.35 0.7 2 650 14 0.1 1.5 0.35 0.7 2 700 21 0.1 1.5 0.35 0.7 2 750 28 0.1 1.5 0.35 0.7 2 800 表 2 平行黏结模型细观参数
Table 2 Parallel bonding model microscopic parameters
养护
时间/d有效
模量/GPa刚度比 摩擦
因数平行黏结
有效模量/
GPa拉强度/
MPa黏结
强度/MPa1 0.110 1.5 0.23 0.110 0.179 1.45 7 0.181 1.6 0.25 0.181 0.189 1.55 14 0.215 1.7 0.27 0.215 0.198 1.61 21 0.247 1.8 0.28 0.247 0.217 1.73 28 0.275 1.9 0.30 0.275 0.259 1.95 表 3 不同养护时间下的力链分布
Table 3 Distribution of force chain under different curing time
养护时间/d 1 7 14 21 28 横向力链分布 纵向力链分布 -
[1] 丁 玉,冯光明,王成真. 超高水充填材料基本性能试验研究[J]. 煤炭学报,2011,36(7):1087−1092. doi: 10.13225/j.cnki.jccs.2011.07.027 DING Yu,FENG Guangming,WANG Chengzhen. Experimental research on basic properties of superhigh-water packing material[J]. Journal of China Coal Society,2011,36(7):1087−1092. doi: 10.13225/j.cnki.jccs.2011.07.027
[2] 孙春东,刘树轮,李继升. 超高水材料在煤矿的系列应用技术[J]. 煤炭科学技术,2017,45(8):42−47. doi: 10.13199/j.cnki.cst.2017.08.008 SUN Chundong,LIU Shulun,LI Jisheng. Application technology of ultra high-water material in coal mine[J]. Coal Science and Technology,2017,45(8):42−47. doi: 10.13199/j.cnki.cst.2017.08.008
[3] 冯光明,丁 玉,朱红菊,等. 矿用超高水充填材料及其结构的试验研究[J]. 中国矿业大学学报,2010,39(6):813−819. FENG Guangming,DING Yu,ZHU Hongju,et al. Experimental research on a superhigh-water packing material for mining and its micromorphology[J]. Journal of China University of Mining & Technology,2010,39(6):813−819.
[4] 郭文兵,马志宝,白二虎. 我国煤矿“三下一上”采煤技术现状与展望[J]. 煤炭科学技术,2020,48(9):16−26. doi: 10.13199/j.cnki.cst.2020.09.002 GUO Wenbing,MA Zhibao,BAI Erhu. Current status and prospect of coal mining technology under building, waterbodies and railways, and above confined water in China[J]. Coal Science and Technology,2020,48(9):16−26. doi: 10.13199/j.cnki.cst.2020.09.002
[5] 王方田,李 岗,班建光,等. 深部开采充填体与煤柱协同承载效应研究[J]. 采矿与安全工程学报,2020,37(2):311−318. WANG Fangtian,LI Gang,BAN Jianguang,et al. Synergistic bearing effect of backfilling body and coal pillar in deep mining[J]. Journal of Mining & Safety Engineering,2020,37(2):311−318.
[6] 杜兆文,陈绍杰,尹大伟,等. 氯盐侵蚀环境下膏体充填体稳定性试验研究[J]. 中国矿业大学学报,2021,50(3):532−538,547. DU Zhaowen,CHEN Shaojie,YIN Dawei,et al. Experimental study of the stability of paste backfill under chlorine erosion environment[J]. Journal of China University of Mining & Technology,2021,50(3):532−538,547.
[7] 任 帅,王方田,李少涛,等. 深井超高水充填工作面小煤柱稳定性规律及控制技术[J]. 煤矿安全,2021,52(6):243−249,254. REN Shuai,WANG Fangtian,LI Shaotao,et al. Stability law and control technology of small coal pillar in super high water backfilling face withdeep mining depth[J]. Safety in Coal Mines,2021,52(6):243−249,254.
[8] 李召峰,张 晨,张 健,等. 不同水饱和度充填体力学性能及损伤机制研究[J]. 采矿与安全工程学报,2021,38(5):1063−1069, 107. LI Zhaofeng,ZHANG Chen,ZHANG Jian,et al. Study on mechanical properties and damage mechanism of backfill with different water saturation[J]. Journal of Mining & Safety Engineering,2021,38(5):1063−1069, 107.
[9] 王旭锋,孙春东,张东升,等. 超高水材料充填胶结体工程特性试验研究[J]. 采矿与安全工程学报,2014,31(6):852−856. doi: 10.13545/j.issn1673-3363.2014.06.004 WANG Xufeng,SUN Chundong,ZHANG Dongsheng,et al. Experimental study on engineering characteristics of super-high water filling body[J]. Journal of Mining & Safety Engineering,2014,31(6):852−856. doi: 10.13545/j.issn1673-3363.2014.06.004
[10] 贾凯军,冯光明,王誉钦,等. 俯斜开采超高水材料袋式充填体失稳机理及防治[J]. 中国矿业大学学报,2015,44(3):409−415. doi: 10.13247/j.cnki.jcumt.000323 JIA Kaijun,FENG Guangming,WANG Yuqin,et al. Instability mechanism and prevention of bag filling body constructed by super-high-water material under the condition of down-dip mining[J]. Journal of China University of Mining & Technology,2015,44(3):409−415. doi: 10.13247/j.cnki.jcumt.000323
[11] 赵家巍,苏 腾,荣腾龙,等. 超高水充填材料侧限压缩应力应变本构关系研究[J]. 采矿与安全工程学报,2020,37(2):394−400. doi: 10.13545/j.cnki.jmse.2020.02.020 ZHAO Jiawei,SU Teng,RONG Tenglong,et al. A constitutive model of superhigh-water filling material under confined compression[J]. Journal of Mining & Safety Engineering,2020,37(2):394−400. doi: 10.13545/j.cnki.jmse.2020.02.020
[12] 苏承东,张振华. 大理岩三轴压缩的塑性变形与能量特征分析[J]. 岩石力学与工程学报,2008,27(2):273−280. doi: 10.3321/j.issn:1000-6915.2008.02.007 SU Chengdong,ZHANG Zhenhua. Analysis of plastic deformation and energy property of marble under pseudo-triaxial compression[J]. Chinese Journal of Rock Mechanics and Engineering,2008,27(2):273−280. doi: 10.3321/j.issn:1000-6915.2008.02.007
[13] 刁兆丰,刘长武,张连卫,等. 掺污泥改性高水材料典型加载方式下的变形特征和力学响应[J]. 科学技术与工程,2018,18(10):262−266. doi: 10.3969/j.issn.1671-1815.2018.10.045 DIAO Zhaofeng,LIU Changwu,ZHANG Lianwei,et al. Deformation characteristics and mechanical response of modified high water material mixed with sludge under typical loading methods[J]. Science Technology and Engineering,2018,18(10):262−266. doi: 10.3969/j.issn.1671-1815.2018.10.045
[14] 冯 波,刘长武,谢 辉,等. 改性高水材料尺寸与形状效应研究[J]. 工程科学与技术,2017,49(S2):121−127. FENG Bo,LIU Changwu,XIE Hui,et al. Experimental study on the size and the shape of high-water-content material that modified[J]. Advanced Engineering Sciences,2017,49(S2):121−127.
[15] 张 钊,刘长武,王一冰,等. 改性高水材料抗压、抗剪强度特征及对比分析[J]. 工程科学学报,2021,43(4):552−560. ZHANG Zhao,LIU Changwu,WANG Yibing,et al. Characteristic and comparative analysis of compressive and shear strengths of modified high-water materials[J]. Chinese Journal of Engineering,2021,43(4):552−560.
[16] 冯 波,刘长武,谢 辉,等. 粉煤灰改性高水材料力学性能试验研究及机理分析[J]. 工程科学学报,2018,40(10):1187−1195. FENG Bo,LIU Changwu,XIE Hui,et al. Experimental study and analysis of the mechanical properties of high-water-content materials modified with fly ash[J]. Chinese Journal of Engineering,2018,40(10):1187−1195.
[17] 赵永辉,冉洪宇,冯国瑞,等. 单轴压缩下不同高宽比矸石胶结充填体损伤演化及破坏特征研究[J]. 采矿与安全工程学报,2022,39(4):674−682. ZHAO Yonghui,RAN Hongyu,FENG Guorui,et al. Damage evolution and failure characteristics of cemented gangue backfill body with different height-width ratios under uniaxial compression[J]. Journal of Mining & Safety Engineering,2022,39(4):674−682.
[18] 郭育霞,赵永辉,冯国瑞,等. 矸石胶结充填体单轴压缩损伤破坏尺寸效应研究[J]. 岩石力学与工程学报,2021,40(12):2434−2444. doi: 10.13722/j.cnki.jrme.2021.0527 GUO Yuxia,ZHAO Yonghui,FENG Guorui,et al. Study on damage size effect of cemented gangue backfill body under uniaxial compression[J]. Chinese Journal of Rock Mechanics and Engineering,2021,40(12):2434−2444. doi: 10.13722/j.cnki.jrme.2021.0527
[19] 侯永强,尹升华,杨世兴,等. 动态荷载下胶结充填体力学响应及能量损伤演化过程研究[J]. 岩土力学,2022,43(S1):145−156. HOU Yongqiang,YIN Shenghua,YANG Shixing,et al. Mechanical response and energy damage evolution process of cemented backfill under impact loading[J]. Rock and Soil Mechanics,2022,43(S1):145−156.
[20] 张学朋,王 刚,蒋宇静,等. 基于颗粒离散元模型的花岗岩压缩试验模拟研究[J]. 岩土力学,2014,35(S1):99−105. doi: 10.16285/j.rsm.2014.s1.009 ZHANG Xuepeng,WANG Gang,JIANG Yujing,et al. Simulation research on granite compression test based on particle discrete element model[J]. Rock and Soil Mechanics,2014,35(S1):99−105. doi: 10.16285/j.rsm.2014.s1.009
[21] 石崇颗粒流(PFC5.0)数值模拟技术及应用[M]. 北京: 中国建筑工业出版社, 2018. [22] Itasca Consulting Group Inc. Manual of particle flow code in 2-dimension (Version 3.10)[M]. Minneapolis: Itasca Consulting Group Inc, 2004.
[23] 冯光明,李乃梁,丁 玉. 煤矿超高水材料充填工艺系统与实践[J]. 矿业工程研究,2019,34(1):13−22. FENG Guangming,LI Nailiang,DING Yu. Process system of filling with supper high-water materials and practice in coal mine[J]. Mineral Engineering Research,2019,34(1):13−22.
[24] 胡炳南,刘鹏亮,崔 锋,等. 我国充填采煤技术回顾及发展现状[J]. 煤炭科学技术,2020,48(9):39−47. doi: 10.13199/j.cnki.cst.2020.09.004 HU Bingnan,LIU Pengliang,CUI Feng,et al. Review and development status of backfill coal mining technology in China[J]. Coal Science and Technology,2020,48(9):39−47. doi: 10.13199/j.cnki.cst.2020.09.004
[25] 冯光明. 超高水充填材料及其充填开采技术研究与应用[D]. 徐州: 中国矿业大学, 2009. FENG Gangming. Research on the superhigh-water packing material and filling mining technology and their application[D]. Xuzhou: China University of Mining and Technology, 2009.
[26] 李 博. 超高水材料开放式充填固结体承载特性试验研究[D]. 徐州: 中国矿业大学, 2016. LI Bo. Experimental study on bearing capacity of open filling consolidation of super high water materials[D]. Xuzhou: China University of Mining and Technology, 2016.
[27] 罗怀廷,李俊孟,黄艳利,等. 露天矿排放矸石三轴压缩宏-细观力学特性研究[J]. 采矿与安全工程学报,2018,35(1):170−178. doi: 10.13545/j.cnki.jmse.2018.01.024 LUO Huaiting,LI Junmeng,HUANG Yanli,et al. Study on macro and microcosmic mechanical properties of crushed gangue from pit mine under triaxial compression[J]. Journal of Mining & Safety Engineering,2018,35(1):170−178. doi: 10.13545/j.cnki.jmse.2018.01.024