Research on relationship between tailing slurry concentration and ultrasonic attenuation coefficient
-
摘要: 针对目前基于超声波法的尾矿浆浓度测量方法没有直接表达浓度与散射衰减系数、黏滞衰减系数的相关关系的问题，采用数值计算方法，研究了超声波在尾矿浆中的衰减特性，得到了声波衰减系数与浓度的相关关系，并通过实验数据的分析论证了数值计算的准确性。研究结果表明：浓度在10%以内时，浓度越大，声速越小；浓度在5%以内时，超声波的衰减系数与浓度近似成正比关系；当超声波频率、粒径一定，浓度在10%以内时，浓度越大，超声波的衰减系数越大。研究结果为超声波在尾矿浆浓度在线监测中的实际应用提供了依据。Abstract: In view of problem that the current tailing slurry concentration measurement method based on ultrasonic method does not directly express correlation between concentration and scattering attenuation coefficients and viscous attenuation coefficients, the attenuation characteristics of ultrasonic wave in tailing slurry were studied by numerical calculation method, and the correlation between acoustic attenuation coefficient and concentration was obtained, and the accuracy of numerical calculation is demonstrated through experimental data. The research results show that within the concentration of 10%, the greater the concentration, the lower the speed of sound; within the concentration of 5%, the attenuation coefficient of ultrasonic wave is approximately proportional to the concentration; when the ultrasonic frequency and particle size are constant, within the concentration of 10%,the greater the concentration, the greater the attenuation coefficient of the ultrasound. The research results provide a basis for the practical application of ultrasound in online monitoring of tailings slurry concentration.
-
0. 引言
移动机器人是实现煤矿智能化和无人化的关键，对提升煤矿企业安全性、降低矿工劳动强度和优化管理效率具有划时代意义[1-2]。基于同时定位与建图（Simultaneous Localization and Mapping, SLAM）技术的移动机器人能够快速、准确、自动化地采集空间数据，进行空间智能感知和环境地图构建，以满足煤矿井下自主作业的要求[3]。煤矿井下无全球导航卫星系统（Global Navigation Satellite System, GNSS）信号且地形复杂，低光照和弱纹理场景易使单一传感器的SLAM精度不足，机器人的前端位姿估计退化，甚至失效。具体表现：相机在低光照、弱纹理和高粉尘的煤矿巷道中无法稳健使用[4]；激光雷达（Light Laser Detection and Ranging, LiDAR）可提供准确的深度测量，不受场景照明和纹理变化的影响，但在狭长巷道中易导致位姿估计退化[5]；惯性测量单元（Inertial Measurement Unit，IMU）响应速度快，基本不受外界环境影响，但其测量误差随时间累计[6]。
多传感器数据融合的SLAM较单一传感器的视觉SLAM、激光雷达SLAM稳定性更好，且具有更高的定位精度和更丰富的环境信息[7-8]。目前，多传感器数据融合的SLAM系统主要分为视觉惯性系统（Visual Inertial System, VIS）、激光惯性系统（LiDAR Inertial System, LIS）及激光视觉惯性系统（LiDAR Visual Inertial System, LVIS）。VIS方面，文献[9]通过单独处理图像和IMU的测量数据，再融合各自的位姿状态估计，得到里程计的位姿，但整体精度较低。文献[10]将视觉与IMU紧耦合，联合优化估计状态，即在视觉良好的环境下可实现高精度的位姿估计与建图，但在低照度、弱纹理场景表现欠佳。文献[11]提出的点线特征融合视觉惯导里程计（Point and Line Features-Visual Inertial Odometry，PL−VIO）和文献[12]提出的线面特征结合的SLAM对点线特征进行联合优化，克服了弱纹理环境下无法进行特征跟踪的问题，但未充分考虑低照度对点线特征跟踪的影响。LIS方面，文献[13]提出一种松耦合的LiDAR惯性里程计及建图系统，根据曲率大小提取点线特征，实现准确的位姿估计，但没有后端闭环优化模块，容易累计漂移。文献[14]提出一种紧耦合的激光雷达惯性里程计（Lidar Inertial Odometry via Smoothing and Mapping, LIO−SAM）方法，该方法将LiDAR里程计因子、IMU预积分因子、全球定位系统（Global Positioning System，GPS）因子及回环检测因子引入因子图中，实现高精度、实时的移动机器人轨迹估计和地图构建。该系统通过ISAM（Incremental Smoothing and Mapping）实现局部和全局的优化，表现出较高的鲁棒性和精度，但尚未充分考虑退化场景对建图的影响。为了将LiDAR、视觉、IMU有效融合用于定位与建图，研究者提出了LVIS。文献[15]提出了V−LOAM（Visual-Lidar Odometry and Mapping）算法，它融合了视觉里程计与LiDAR里程计，用视觉里程计进行位姿估计，但该算法执行过程按逐帧顺序进行，全局整体一致性不强。文献[16]提出了VIL−SLAM（Visual-Inertial SLAM with Points and Lines）算法，该算法是基于立体视觉、IMU、LiDAR融合的SLAM系统，将视觉惯性里程计与激光里程计相结合，提高了定位和地图构建的精度，但该算法计算量大、硬件成本高，在光照和遮挡等复杂环境下易出现定位和地图构建错误。文献[17]提出了一种LiDAR与视觉融合的方法，该方法的LiDAR里程计位姿是由视觉惯性里程计给予先验信息初始化，同时可优化两者的位姿，提高位姿估计的精度，但在低光照时易发生初始化错误，导致LiDAR里程计位姿误差较大。文献[18]提出了紧耦合的激光视觉惯导SLAM（Lidar-Visual-Inertial Odometry via Smoothing and Mapping，LVI−SAM）系统，该系统被分成LIS子系统和VIS子系统，视觉为LiDAR提供初始位姿估计值，LiDAR给视觉提供深度信息，在后端利用视觉信息进行闭环检测，利用LiDAR进行闭环优化。由于视觉部分依赖于图像特征，所以可能会受到光照和纹理的影响。
针对煤矿井下复杂环境导致多传感器融合SLAM存在机器人前端位姿估计退化失效和后端融合精度不足的问题，本文提出一种LiDAR−视觉−惯性自适应融合SLAM方法。首先，在视觉图像数据预处理阶段引入图像增强算法，以减少低光照场景对视觉特征提取的影响。然后，动态设置检测阈值来自适应选取当前最优位姿，提高位姿估计精度。最后，基于滑动窗口实现多传感器数据联合非线性优化，实现一种高精度、低漂移的多传感器数据融合SLAM系统。
1. 方法概述
LiDAR−视觉−惯性自适应融合SLAM方法可分为多源传感器数据预处理、LiDAR/视觉/IMU里程计、后端优化3个部分，如图1所示。在多源传感器数据预处理部分，首先对LiDAR非地面点云数据进行聚类分割，提取线面特征，然后利用IMU预积分状态进行畸变校正，最后对低照度图像进行增强，再提取视觉点线特征。用IMU预积分状态为LiDAR特征匹配与视觉特征跟踪提供位姿初值。根据LiDAR相邻帧的线面特征匹配得到机器人位姿，之后进行视觉点线特征跟踪，分别计算LiDAR、视觉、IMU位姿变化值，通过设定动态阈值来检测前端里程计的稳定性，自适应选取最优位姿。对不同传感器构建残差项，包括点云匹配残差、IMU预积分残差、视觉点线残差、边缘化残差。为了兼顾精度与实时性，基于滑动窗口实现激光点云特征、视觉特征、IMU测量的多源数据联合非线性优化，实现煤矿井下连续可用、精确可靠的SLAM。
![]()
图 1 LiDAR−视觉−惯性融合SLAM方法原理Figure 1. Principle of LiDAR-visual -inertial measurement unit (IMU) fusion simultaneous localization and mapping (SLAM) method2. 方法关键技术
2.1 多源传感器数据预处理
2.1.1 IMU预积分
通过IMU预积分可获得机器人状态。为避免每次IMU测量过程中增加新的状态量，通常在2帧之间增加1个重新参数化过程来实现运动约束，避免重复积分。
k时刻机器人的状态量为
$$ {{\boldsymbol{I}}_k} = \left[ {{{\boldsymbol{R}}_k},{p_k},{v_k},{b_k}} \right] $$ (1) 式中：$ {{\boldsymbol{R}}_k} $为k时刻的旋转矩阵；$ {p_k},{v_k},{b_k} $分别为IMU预积分获得的位置、速度、偏置。
基于文献[16]中的IMU预积分方法，将关键帧的状态加入到滑动窗口中，执行光束法平差（Bundle Adjustment, BA）优化，以实现IMU坐标系到世界坐标系的变换。可得到连续IMU关键帧$ m_{i} $和$ m_{i+1} $之间的相对运动测量值：
$$ \begin{split} & \boldsymbol{r}_{\mathcal{B}}(\hat{\text{z}}_{m_i+1}^{m_i},\boldsymbol{I})= \\ &\quad\left[ \begin{array}{l}\delta\varphi_{m_{i+1}}^{m_i} \\ \delta\mu_{m_{i+1}}^{m_i} \\ \delta\theta_{m_{i+1}}^{m_i} \\ \delta b_a \\ \delta b_g\end{array} \right]=\left[ \begin{array}{c}\boldsymbol{R}_{\mathrm{w}}^{m_i}\left(p_{m_{i+1}}^{\mathrm{w}}-p_{m_i}^{\mathrm{w}}+\dfrac{1}{2}g^{\mathrm{w}}\Delta t^2+v_{m_i}^{\mathrm{w}}\Delta t\right)-\hat{\varphi}_{m_{i+1}}^{m_i} \\ \boldsymbol{R}_{\mathrm{w}}^{m_i}(v_{m_i+1}^{\mathrm{w}}+g^{\mathrm{w}}\Delta t-v_{m_i}^{\mathrm{w}})-\mu_{m_{i+1}}^{m_i} \\ 2[\left(q_{m_i}^{\mathrm{w}}\right)^{-1}\otimes q_{m_{i+1}}^{\mathrm{w}}\otimes(\hat{\theta}_{m_{i+1}}^{m_i})^{-1}]_{xy\text{z}} \\ b_a^{m_{i+1}}-b_a^{m_i} \\ b_g^{m_{i+1}}-b_g^{m_i}\end{array} \right] \end{split} $$ (2) 式中：${\boldsymbol{r}}_{\mathcal{B}}(\hat{{\textit{z}}}_{m_{i+1}}^{m_i}, {{\boldsymbol{I}}}) $为IMU预积分残差；$ \mathcal{B} $为滑动窗口中预积分的IMU测量集；$ \hat{{\textit{z}}}_{m_{i+1}}^{m_i} $为第i帧关键帧$ m_{i} $与第i+1帧关键帧$ m_{i+1} $的残差值；I为机器人的状态量；$ \delta $为预积分；$ \varphi_{m_{i+1}}^{m_{i}} $，$ \mu_{m_{i+1}}^{m_{i}} $，$ \theta_{m_{i+1}}^{m_{i}} $为相邻2个关键帧之间的IMU观测值；$ b_{a} $和$ b_{g}$分别为加速度a偏置和重力加速度g偏置；$ {\boldsymbol{R}}_{{\mathrm{w}}}^{m_{i}} $为关键帧$ m_{i} $在世界坐标系下的旋转矩阵；$ p_{m_{i}}^{{{\mathrm{w}}}} $，$ p_{m_{i}+1}^{{\mathrm{w}}} $，$ v_{m_{i}}^{{\mathrm{w}}} $，$ v_{m_{i+1}}^{{\mathrm{w}}} $，$ q_{m_{i}}^{{\mathrm{w}}} $，$ q_{m_{i+1}}^{{\mathrm{w}}} $分别为关键帧$ m_{i} $与关键帧$ m_{i+1} $的IMU坐标系到世界坐标系的平移、速度和旋转；$ g^{{\mathrm{w}}} $为世界坐标系下的重力加速度；$ \Delta t $为关键帧$ m_{i} $到关键帧$ m_{i+1} $所需时间；$ \hat{\varphi}_{m_{i+1}}^{m_i} $，$ \hat{\theta}_{m_{i+1}}^{m_i}$为含有噪声的关键帧$ m_{i} $与关键帧$ m_{i+1} $之间的IMU观测值；$ x, y, {\textit{z}} $为三轴方向；$ b_{a}^{m_{i}} $，$ b_{a}^{m_{i+1}} $分别为关键帧$ m_{i} $与关键帧$ m_{i+1} $的加速度a偏置；$ b_{g}^{m_{i}} $，$ b_{g}^{m_{i+1}} $分别为关键帧$ m_{i} $与关键帧$ m_{i+1} $的重力加速度g偏置。
2.1.2 激光点云数据预处理
LiDAR点云数据预处理主要有点云数据分割、线面特征提取和运动畸变校正。首先，基于随机样本一致性（Random Sample Consensus, RANSAC）[19]的快速点云分割算法对LiDAR点云数据进行分割处理。然后，利用深度图计算出分割后的LiDAR点云数据的曲率，将曲率较大的非地面点标记为边缘特征点，将曲率较小的标记为平面特征点，从而提取煤矿巷道线特征和平面特征。最后，为了消除因传感器移动导致的点云数据畸变问题，使用IMU的预积分数据对LiDAR点云数据进行运动补偿。
2.1.3 视觉图像预处理
煤矿井下高粉尘、弱纹理、灯光区域与背光区域亮度反差大，图像中物体表面容易出现反光、高光现象，使得图像整体或局部光照不均，严重影响视觉特征提取与跟踪，易使视觉SLAM发生漂移甚至失效。因此，需在视觉特征提取前对图像质量进行增强，以突出图像中的特征信息。在HSV空间下，本文将自适应Gamma校正[20-21]和对比度受限的自适应直方图均衡化（Contrast Limited Adaptive Histogram Equalization, CLAHE）算法[22]相融合，以对图像进行增强，从而提高光照不均环境下视觉特征提取数量。图像增强算法流程如图2所示。首先，将图像从RGB空间变换到HSV空间，避免在RGB颜色空间中处理造成图像失真。其次，只对亮度分量进行自适应Gamma校正处理。然后，对Gamma校正结果进行CLAHE算法处理，将图像中分布不均匀的直方图变换为均匀分布的直方图，有效拉伸灰度值的动态范围，以提高图像对比度。最后，将未处理的色调分量、饱和度分量和处理后的亮度分量进行融合，并逆变换回RGB颜色空间，得到最终增强图像。
对增强图像进行视觉点线特征提取。首先，使用Kanade-Lucas-Tomasi（KLT）特征跟踪稀疏光流算法[23]提取和跟踪每个新图像帧中的初始特征点，保持每帧图像中点特征数量为200～400个。然后，使用RANSAC和基本矩阵模型进行异常值抑制[19]，剔除离群值大的特征点。最后，采用线段检测器（Line Segment Detector, LSD）[24]提取线段，并用线带描述符（Line Band Descriptor, LBD）[25]描述线特征。
2.2 LiDAR、视觉、IMU里程计
2.2.1 点云数据线面特征匹配
考虑到煤矿巷道环境中平面特征远多于线特征，采用平面特征进行点到面距离最小约束，估计出Z方向、翻滚角、俯仰角的变化量；在此基础上进行点到线距离最小约束，估计出X方向、Y方向、偏航角的变化量，联合地面点和角点优化估计连续帧的相对位姿；采用麦夸特（Levenberg-Marquardt，LM）算法[13]迭代求解最优位姿，迭代计算的初值为当前帧IMU预积分的位姿，可有效减少迭代次数，且消除零初值或匀速运动假设初值带来的误匹配。
$$ \begin{cases}e_{A}^{ {{\mathrm{line}} }}=\varGamma_{A}\left(Q_{A}-\left({\boldsymbol{R}} P_{A}+t\right)\right) & Q_{A}, P_{A} \in F_{A}^{ {{\mathrm{line}} }} \\ e_{B}^{ {{\mathrm{plane}} }}=\varLambda_{B}^{\rm{T}}\left(Q_{B}-\left({\boldsymbol{R}} P_{B}+t\right)\right) & Q_{B}, P_{B} \in F_{B+1}^{ {{\mathrm{plane}} }}\end{cases} $$ (3) $$ \min _{{\boldsymbol{I}}}\left\{\underset{P \in F_A}{\varSigma}\left\|e_{A}^{ {{\mathrm{point}} } \rightarrow { {\mathrm{line}}, {\mathrm{plane}} }}\right\|^2+\underset{(A, A+1) \in \mathcal{B}}{\varSigma}\left\|e^{{\mathrm{i m u}}}\right\|^2+E_{ {{\mathrm{imuodom}} }}^{ {{\mathrm{prior}} }}\right\} $$ (4) 式中：$ e_{A}^{ {{\mathrm{line}} }} $和$ e_{{B}}^{ {{\mathrm{plane}} }} $分别为第A个点到线、第B条线到面的距离；$ \varGamma_{A} $和$ \varLambda_{B} $分别为特征对应的主方向和法线方向；QA和PA为在线特征点集$ F_{A+1}^{{{\mathrm{line}}}} $集合中的线特征点；R为旋转矩阵；t为时间；QB和PB为在面特征点集$ F_{B+ 1}^{{\mathrm{plane}}} $集合中的面特征点；$ \underset{P \in F_{A}} {\varSigma}e_{A}^{ {{\mathrm{point}} } \rightarrow { {\mathrm{line}} }, { {\mathrm{plane}} }}$为点到线、点到面的残差；$\underset{P \in F_{{A}}}{\varSigma}$为特征点$ P $在线特征点集合$ F_{A} $上的点；$ \underset{(A,A+1) \in \mathcal{B}}{\varSigma}\left\|e^{{\mathrm{i m u}}}\right\|^2$和$ E_{ {{\mathrm{imuodom}} }}^{ {{\mathrm{prior}} }} $分别为IMU预积分因子和IMU里程计提供的预测位姿先验值。
2.2.2 视觉点线特征跟踪
首先选择相似的线段进行跟踪，并添加线段长度相近的图像帧上的特征点集，确保在新图像至少有2个特征点的线段。然后，根据文献[26]中的跟踪线段策略，对2个连续图像帧线进行跟踪。最后，为提高视觉位姿的精度，采用点特征投影残差、线特征投影残差和IMU的测量结果构造局部残差约束。当插入1个新的关键帧时，对当前建图执行1个光束法平差调整。在滑动窗口中，通过最小化所有测量残差的先验项和代价项，对所有状态变量进行优化。
$$ \begin{split} \mathcal{C}=&\min _{{{\boldsymbol{I}}}}\Bigg\{ \underset{i \in \mathcal{B}}{\varSigma}\left\|{\boldsymbol{r}}_{\mathcal{B}}(\hat{{\textit{z}}}_{m_{i+1}}^{m_i}, {{\boldsymbol{I}}})\right\|_{\varSigma_{m_{i+1}}^{m_i}}^2+\underset{(i, j) \in \mathcal{F}_p}{\varSigma}\left\|{\boldsymbol{r}}_{\mathcal{F}_{\mathcal{P}}}(\hat{{\textit{z}}}_j^{c_i}, {{\boldsymbol{I}}})\right\|_{\varSigma_j^{c_i}}^2+\\& \underset{(i, k) \in \mathcal{F}_\mathcal{L}}{\varSigma}\left\|{\boldsymbol{r}}_{\mathcal{F}_\mathcal{L}}(\hat{{\textit{z}}}_k^{c_i}, {{\boldsymbol{I}}})\right\|_{\varSigma_k^{c_k}}^2+\left\|{\boldsymbol{r}}_\iota-{\boldsymbol{J}}_\iota {{\boldsymbol{I}}}\right\|_{\varSigma_\iota}^2\Bigg\} \end{split} $$ (5) 式中：$\varSigma_{m_{i+1}}^{m_i}$为相邻关键帧的集合；$ {\boldsymbol{r}}_{\mathcal{F}_{\mathcal{P}}}(\hat{{\textit{z}}}_j^{c_i}, {{\boldsymbol{I}}})$和$ {\boldsymbol{r}}_{\mathcal{F}_{\mathcal{L}}}(\hat{{\textit{z}}}_k^{c_i}, {{\boldsymbol{I}}})$分别为点特征投影残差和线特征投影残差；$ \mathcal{F}_ \mathcal{P} $和$ \mathcal{F}_ \mathcal{L} $分别为在滑动窗口中至少观察到2次的点特征集和线特征集；$ \hat{{\textit{z}}}_{j}^{c_{i}} $和$ \hat{{\textit{z}}}_{k}^{c_{i}} $分别为第$ i $帧图像帧$ c_{{i}} $中的第$ j $个视觉点特征残差值和第$ k $个视觉线特征残差值；$ \varSigma_{j}^{c_{i}} $，$ \varSigma_{k}^{c_{i}} $分别为在第$ i $帧图像帧$ c_{{i}} $中点特征集合、线特征集合；$\left\|{\boldsymbol{r}}_\iota-{\boldsymbol{J}}_\iota {{\boldsymbol{I}}}\right\|_{\varSigma_\iota}^2 $为边缘化残差的先验信息；$ {\boldsymbol{r}}_\iota$为边缘化残差信息；$ {\boldsymbol{J}}_\iota$为边缘化雅可比矩阵；$ {\varSigma_\iota}$为边缘化集合。
对于点特征，重投影误差可定义为观测点与重投影点在同一像素平面上的距离。假设第$ {j} $个视觉点特征在第$ h $帧图像帧$ c_{h} $中首先被观察到，则第$ i$帧图像帧$ c_{i} $中的重投影误差为
$$ \boldsymbol{r}_{\mathcal{F}_P}(\hat{\text{z}}_j^{c_i},\boldsymbol{I})=\left[ \begin{array}{l}\dfrac{x^{c_i}}{\text{z}^{c_i}}-\mathcal{U}_j^{c_i} \\ \dfrac{y^{c_i}}{\text{z}^{c_i}}-\mathcal{V}_j^{c_i}\end{array} \right] $$ (6) $$ \begin{split} {\boldsymbol{D}}_j^{c_i}=&\left[ \begin{array}{l} x^{c_i} \\ y_i^{c_i} \\ {\textit{z}}^{c_i} \end{array} \right]={\boldsymbol{R}}_{\varrho c}\Bigg({\boldsymbol{R}}_{ \mathbf{w}_{\varrho_i}}\Bigg({\boldsymbol{R}}_{\varrho_i \mathbf{w}}\Bigg({\boldsymbol{R}}_{c \varrho} \frac{1}{\lambda_j}\left[ \begin{array}{c} \mathcal{U}_j^{c_h} \\ \mathcal{V}_j^{c_h} \\ 1 \end{array} \right]+\mathcal{P}_{c \varrho}\Bigg)+\\&\mathcal{P}_{\varrho_h \mathbf{w}}-\mathcal{P}_{\varrho_i \mathbf{w}}\Bigg)-\mathcal{P}_{c \varrho}\Bigg) \end{split}$$ (7) 式中：$ {({x^{{c_i}}},{y^{{c_i}}},{{\textit{z}}^{{c_i}}})^{\text{T}}} $为图像帧ci的坐标；$ [\mathcal{U}_j^{c_i}, \mathcal{V}_j^{c_i}]^{\rm{T}} $为图像帧$ c_{i} $的第j个视觉点特征的像素坐标；$ {{\boldsymbol{D}}}_j^{c_i}$为图像帧$ c_{{i}} $的第$ j $个视觉点特征的重投影误差；$ {\boldsymbol{R}}_{\varrho c},\;{\boldsymbol{R}}_{\mathrm{w}_{ \varrho_i}},\; {\boldsymbol{R}}_{\varrho_i \mathrm{w}},\;{\boldsymbol{R}}_{c \varrho}$分别为IMU到相机、世界坐标到IMU、IMU到世界坐标、相机到IMU的旋转矩阵；$ \dfrac{1}{\lambda_{j}} $为第$ j $个视觉点特征的缩放值；$ [\mathcal{U}_j^{c_h}, \mathcal{V}_j^{c_h}]^{\rm{T}}$为图像帧ch的第$ j$个视觉点特征的像素坐标；$ \mathcal{P}_{c \varrho}, \mathcal{P}_{\varrho_h {\rm{w}}}, \mathcal{P}_{\varrho_\iota {\rm{w}}}$分别为IMU与相机、世界坐标与IMU相互转换的平移量。
对于线特征，重投影误差可定义为观测线端点到同一像素平面内重投影线的距离。将一条三维直线$ \mathcal{L}_k^w$投影到二维图像平面上，即可得到图像帧$ c_i$的二维投影线段$ {L}_k^{ c_i}$。
$$ L_k^{c_i}=\boldsymbol{K}{\boldsymbol{N}_c}=\left[ \begin{array}{ccc} f_{\mathcal{V}} & 0 & 0 \\ 0 & f_{\mathcal{U}} & 0 \\ -f_{\mathcal{V}} c_{\mathcal{U}} & f_{\mathcal{U}} c_{\mathcal{V}} & f_{\mathcal{U}} f_{\mathcal{V}} \end{array} \right] \boldsymbol{N}_c=\left[ \begin{array}{l} L_1 \\ L_2 \\ L_3 \end{array} \right] $$ (8) $$ {\boldsymbol{\gamma}}_{\mathcal{F}_{\mathcal{c}}}\left(\hat{{\textit{z}}}_k^{c_i}, {{\boldsymbol{I}}}\right)=d\left(I_k, L_k^{c_i}\right)=\left[\frac{s_k^{\rm{T}}}{\sqrt{L_1^2+L_2^2}} \frac{l_k^{\rm{T}} L_k^{c_i}}{\sqrt{L_1^2+L_2^2}}\right]^{\rm{T}} $$ (9) 式中：$L_{k}^{c_{i}} $为图像帧ci中的第k个视觉线特征；$ \boldsymbol{K} $为三维直线的内参矩阵；$ \boldsymbol{N}_{c} $为图像帧$ c $的旋转矩阵；$f_{\mathcal{V}}$，$ f_\mathcal{U} $分别为二维像素坐标系下的像素值；$ c_\mathcal{U} $，$ c_\mathcal{V} $分别为二维像素坐标系下图像帧c的像素转换值；$ L_{1} $，$ L_{2} $，$ L_{3} $分别为三轴方向的线段分量；${\boldsymbol{\gamma}}_{\mathcal{F}_{\mathcal{c}}} $为视觉线特征重投影误差；$ d\left(I_{k}, L_{k}^{c_{i}}\right) $为图像帧$ c_{i} $中的第$ k $个视觉线特征与投影到二维像素平面的第$ k $个视觉线特征的距离；$ s_{k} $和$ l_{k} $为二维像素平面观测到的端点。
2.2.3 自适应位姿选取
首先，设定一个动态阈值，根据IMU预积分对三轴方向的加速度积分，得到三轴方向的速度，计算出IMU在世界坐标系下的位移，将一定周期内的位移设置为动态阈值。然后，将LiDAR−IMU里程计位姿、视觉−IMU里程计位姿、LiDAR−视觉−IMU里程计位姿转换到世界坐标系下，分别计算出一定周期内的位姿偏差。最后，将动态阈值与LiDAR−IMU里程计位姿偏差、视觉−IMU里程计位姿偏差、LiDAR−视觉−IMU里程计位姿偏差进行检验，以此判断上述3个里程计输出的位姿是否准确，从而根据动态阈值自适应选取最优位姿。设定动态阈值前，需先对三轴加速度积分，并计算三轴在世界坐标系下的位移。
$$ S = \frac{1}{\eta }\sum\limits_{K = 1}^n {V({t_K})} = \frac{1}{\eta }\sum\limits_{K = 1}^n {\left({V_0} + \sum\limits_{i = 1}^K {a({t_i})} \frac{1}{\eta }\right)} $$ (10) 式中：$S$为动态阈值；$ \eta $为IMU的频率；K为帧数；n为总帧数；$ V $为IMU预积分后的速度；$ {t_K} $为第$K$帧的时间；V0为IMU预积分后的初速度；$a({t_i})$为${t_i}$时刻的加速度。
$S$是通过IMU预积分后的速度计算得来的，IMU误差会随着时间增加而累计。为此，在IMU预积分周期内，将上一周期内的平均速度作为IMU初速度来校准当前时刻的IMU预积分速度，以防止累计误差。
$$ \left\{ \begin{gathered} {({x_{{\mathrm{lio}}}},{y_{{\mathrm{lio}}}},{{\textit{z}}_{{\mathrm{lio}}}})^{\rm T}} = {\boldsymbol{U}} _n^{{\mathrm{lio}}{\mathrm{w}}} - {\boldsymbol{U}} _{n - 1}^{{\mathrm{lio}}{\mathrm{w}}} \\ {({x_{{\mathrm{cio}}}},{y_{{\mathrm{cio}}}},{{\textit{z}}_{{\mathrm{cio}}}})^{\rm T}} = {\boldsymbol{U}} _n^{{\mathrm{cio}}{\mathrm{w}}} - {\boldsymbol{U}} _{n - 1}^{{\mathrm{cio}}{\mathrm{w}}} \\ {({x_{{\mathrm{lci}}}},{y_{{\mathrm{lci}}}},{{\textit{z}}_{{\mathrm{lci}}}})^{\rm T}} = {\boldsymbol{U}} _n^{{\mathrm{lci}}{\mathrm{w}}} - {\boldsymbol{U}} _{n - 1}^{{\mathrm{lci}}{\mathrm{w}}} \\ \end{gathered} \right. $$ (11) 式中：$ \left( {{x_{{\mathrm{lio}}}},{y_{{\mathrm{lio}}}},{{\textit{z}}_{{\mathrm{lio}}}}} \right) $为LiDAR−IMU里程计位姿坐标；$ \left( {{x_{{\mathrm{cio}}}},{y_{{\mathrm{cio}}}},{{\textit{z}}_{{\mathrm{cio}}}}} \right) $为视觉−IMU里程计位姿坐标；$ ( {x_{{\mathrm{lci}}}}, {y_{{\mathrm{lci}}}},{{\textit{z}}_{{\mathrm{lci}}}} ) $为LiDAR−视觉−IMU里程计位姿坐标；${\boldsymbol{U}} $为平移矩阵。
在$\Delta t$周期内，根据IMU对加速度积分计算得到机器人位姿$({x_{{\mathrm{imu}}}},{y_{{\mathrm{imu}}}},{{\textit{z}}_{{\mathrm{imu}}}})$，根据LiDAR−视觉−IMU里程计位姿得到机器人位置差${H_{{\mathrm{lci}}}}$、LiDAR−IMU里程计位姿得到机器人位置差值${H_{{\mathrm{lio}}}}$、视觉−IMU里程计位姿得到机器人位置差值${H_{{\mathrm{cio}}}}$。
$$ \left\{ \begin{gathered} {H_{{\mathrm{lci}}}} = \sqrt {{{({x_{{\mathrm{lci}}}} - {x_{{\mathrm{imu}}}})}^2} + {{({y_{{\mathrm{lci}}}} - {y_{{\mathrm{imu}}}})}^2} + {{({{\textit{z}}_{{\mathrm{lci}}}} - {{\textit{z}}_{{\mathrm{imu}}}})}^2}} \\ {H_{{\mathrm{lio}}}} = \sqrt {{{({x_{{\mathrm{lio}}}} - {x_{{\mathrm{imu}}}})}^2} + {{({y_{{\mathrm{lio}}}} - {y_{{\mathrm{imu}}}})}^2} + {{({{\textit{z}}_{{\mathrm{lio}}}} - {{\textit{z}}_{{\mathrm{imu}}}})}^2}} \\ {H_{{\mathrm{cio}}}} = \sqrt {{{({x_{{\mathrm{cio}}}} - {x_{{\mathrm{imu}}}})}^2} + {{({y_{{\mathrm{cio}}}} - {y_{{\mathrm{imu}}}})}^2} + {{({{\textit{z}}_{{\mathrm{cio}}}} - {{\textit{z}}_{{\mathrm{imu}}}})}^2}} \\ \end{gathered} \right. $$ (12) 自适应位姿选取方法通过输入动态阈值$S$和位置误差（${H_{{\mathrm{lci}}}}$,${H_{{\mathrm{lio}}}}$,${H_{{\mathrm{cio}}}}$），对LiDAR−视觉−IMU里程计的位置误差${H_{{\mathrm{lci}}}}$与动态阈值$S$进行判断。在$ \Delta t $周期内：当${H_{{\mathrm{lci}}}} ＜ S$时，选择LiDAR−视觉−IMU里程计的位姿为后端优化的初始值。当${H_{{\mathrm{lci}}}} ＞ S$时，需比较${H_{{\mathrm{lio}}}}$与${H_{{\mathrm{cio}}}}$的大小，若$ {\mathrm{{\rm H}}}_{\mathrm{l}\mathrm{i}\mathrm{o}} ＜ {\mathrm{{\rm H}}}_{\mathrm{c}\mathrm{i}\mathrm{o}} $且${H_{{\mathrm{lio}}}} ＜ {H_{{\mathrm{lci}}}}$，选择LiDAR−IMU里程计的位姿；若${H_{{\mathrm{lio}}}} ＜ {H_{{\mathrm{cio}}}}$且${H_{{\mathrm{lio}}}} ＞ {H_{{\mathrm{lci}}}}$，选择LiDAR−视觉−IMU里程计的位姿；若${H_{{\mathrm{cio}}}} ＜ {H_{{\mathrm{lio}}}}$且${H_{{\mathrm{cio}}}} ＜ {H_{{\mathrm{lci}}}}$，选择视觉−IMU里程计的位姿；若${H_{{\mathrm{cio}}}} ＜ {H_{{\mathrm{lio}}}}$且${H_{{\mathrm{cio}}}} ＞ {H_{{\mathrm{lci}}}}$时，选择LiDAR−视觉−IMU里程计位姿。经过动态阈值$S$检验与位姿自适应筛选，可选取最优位姿作为后端优化的初始值。
2.3 后端优化
在位姿估计和IMU状态预测后，使用滑动窗口紧耦合优化方法估计系统的最优状态。为了确保较高的精度与实时性，本文选择滑动窗口大小为10，当窗口满时，需对窗口内的状态节点进行滑动，从而实现局部光束法平差。
采用非线性优化方法对LiDAR、视觉、IMU传感器的观测值进行联合优化。根据非线性优化理论[27]，先构建不同类型观测值相对于待估状态的量测方程，再构建残差项，进一步将不同类型的残差加起来，得到整个优化问题的代价函数。
$$ \begin{split}& \min _{{{\boldsymbol{I}}}}\left\{\underset{(\varPsi \varUpsilon) \in \mathcal{B}}{\varSigma}\left\|{\boldsymbol{r}}_{\mathcal{B}}(\hat{{\textit{z}}}_{m_{i+1}}^{m_i}, {{\boldsymbol{I}}})\right\|+\underset{(\varPsi \varUpsilon) \in \mathcal{B}}{\varSigma} e_{\varPsi \varUpsilon}^{{\mathrm{L I S}}} {\boldsymbol{W}}_{\varPsi \varUpsilon}^{-1} e_{\varPsi \varUpsilon}^{{\mathrm{L I S}}}+\right.\\&\quad \left.\underset{(\varPsi \varUpsilon) \in \mathcal{B}}{\varSigma} e_{\varPsi \varUpsilon}^{{\mathrm{C I S}}^{\rm{T}}} {\boldsymbol{W}}_{\varPsi \varUpsilon}^{-1} e_{\varPsi \varUpsilon}^{{\mathrm{C I S}}}+E_m\right\} \end{split} $$ (13) 式中：$ \underset{(\varPsi \varUpsilon) \in \mathcal{B}}{\varSigma}\left\|{\boldsymbol{r}}_{\mathcal{B}}(\hat{z}_{m_{i+1}}^{m_i}, {{\boldsymbol{I}}})\right\|$为IMU预积分因子；$ \underset{(\varPsi \varUpsilon) \in {\mathcal{B}}} {\varSigma} $为关键帧$ \varPsi$和关键帧$ \varUpsilon$的IMU数据集合；$ \underset{(\varPsi \varUpsilon) \in \mathcal{B}}{\varSigma} e_{\varPsi \varUpsilon}^{{\mathrm{LI S}}^{\rm{T}}} {\boldsymbol{W}}_{\varPsi \varUpsilon}^{-1} e_{\varPsi \varUpsilon}^{{\mathrm{LI S}}}$为加权后的LiDAR线面残差因子；$ e_{\varPsi \varUpsilon}^{{\mathrm{L I S}}}$为LiDAR关键帧$ \varPsi $和关键帧$ \varUpsilon$之间的位姿残差项；$ {\boldsymbol{W}}_{\varPsi \varUpsilon}^{-1} $为协方差矩阵权重；$\underset{(\varPsi, \varUpsilon){\in \mathcal{B}}}{\varSigma} e_{\varPsi \varUpsilon}^{{\mathrm{C I S}}^{\rm{T}}} {\boldsymbol{W}}_{\varPsi \varUpsilon}^{-1} e_{\varPsi \varUpsilon}^{{\mathrm{C I S}}}$为加权后的视觉点线残差因子；$ e_{\varPsi \varUpsilon}^{{\mathrm{C I S}}} $为视觉关键帧$ \varPsi $和关键帧$ \varUpsilon$之间的位姿残差项；$ E_{m} $为边缘化残差项。
残差项构建后，根据因子图理论将上述涉及到的多源数据抽象为测量因子，采用基于最大后验估计的滑动窗口因子图优化算法构建因子图优化模型，包括LiDAR里程计（LiDAR Odometry，LO）因子、IMU预积分因子、视觉里程计（Visual Odometry，VO）因子、激光视觉惯导里程计（LiDAR Visual IMU Odometry，LVIO）因子。状态节点变量用橘黄色圆圈表示，包括移动机器人的姿态、位置、速度等信息，如图3所示。首先，通过关键帧选取策略对实时采集的LiDAR扫描帧和序列图像提取相应的关键帧，并将关键帧对应的LiDAR和相机状态节点加入因子图，同时将对应的LO因子和VO因子加入因子图。然后，对LiDAR扫描帧和图像间的IMU数据进行预积分，得到IMU预积分因子，添加IMU预积分因子。最后，构建该时刻的LVIO因子，并加入因子图进行优化。
3. 实验及结果分析
3.1 试验介绍
为了验证本文方法的性能，自主搭建了LiDAR、相机、IMU多源传感器的硬件采集平台，在狭长走廊和煤矿井下环境进行试验，如图4所示。本文系统基于ROS进行时间同步和坐标变换，非线性优化问题使用Ceres库实现。提前使用文献[28]中方法标定LiDAR和相机，用Kalibr[29]标定IMU和相机。上述硬件设备型号和详细参数见表1。
表 1 多源传感器设备型号与信息Table 1. Sensors model and information设备 设备型号 设备信息 LiDAR VLP−16 频率：10 Hz ；最大测量距离：150 m IMU Ellipse2−N 频率：200 Hz ；横滚/俯仰误差为±0.1°；航向误差为0.5° 相机 Zed−2i 频率：30 Hz； 图像分辨率（H×V）：1 280×720 控制台 Autolabor−PC CPU：AMD Ryzen3 3 200G、DDR4 8 G 机器人 Autolabor−Pro1 驱动模式：四驱；位移速度：0.5～1.5 m/s；旋转角速度：0.56 rad/s 狭长走廊中布设控制点L1—L6，光照不均匀，结构对称；煤矿井下布设控制点M1—M9，存在低光照巷道，需翻越轨道，存在大量的喷浆表面和对称巷道，如图5所示。在移动机器人上粘贴测量标志，使用全站仪测量标志点坐标为参考值，以移动机器人在对应控制点位置停止时间段内的平均位置估值为测量值。
![]()
图 5 狭长走廊与煤矿巷道环境及控制点布设情况Figure 5. Control points layout in narrow corridor and coal mine roadway environment3.2 图像增强试验
为验证本文图像增强算法的有效性，分别与常用的Retinex算法、自适应Gamma校正算法、CLAHE算法和文献[30]所提算法进行对比。选取煤矿巷道低照度区域的（M7，M9）数据进行图像增强和特征点匹配，定性分析图像增强效果和定量评估特征点匹配情况。
5种图像增强算法处理结果如图6所示。可看出Retinex算法调节了图像的整体亮度，但同时给图像带来了一定噪声，使整张图像出现泛白现象；自适应Gamma校正算法对灰度级进行拉伸，使得灰度级更加平滑，对暗区域增强效果更显著，会降低光照区域的增强效果；CLAHE算法虽然提高了图像的对比度，保留了丰富的细节信息，但对暗区域的亮度增强效果不明显；文献[30]引入熵指数来控制增强程度，再根据图像局部特征进行自适应增强，但增强的亮度不如CLAHE算法明显；本文算法利用自适应Gamma校正算法和CLAHE算法优势互补，在增强图像亮度的同时提高了图像对比度，图像细节体现更丰富，说明本文算法对煤矿巷道图像亮度有明显改善，更符合人眼视觉感知，充分展现出环境中的特征信息。
通过对比图像增强前后特征点提取和跟踪情况，对图像质量进行定量评估，结果见表2。可看出原始图像提取的特征点很少，且匹配成功率不超过39.6%，无法满足视觉SLAM要求。经过Retinex算法增强后提取的特征点数量有所增加，由于增强过程中引入了噪声，导致后续RANSAC去除误匹配点的匹配成功率较低。自适应Gamma校正算法、CLAHE算法和文献[30]所提算法虽能够提高匹配成功率，增强视觉里程计的稳定性，但视觉位姿估计结果不够准确。本文算法提升了背光区和光照区的亮度和对比度，增加了图像中的特征信息，在特征点提取数量和匹配质量方面都有大幅度提升，匹配成功率达90.7%，提高了视觉里程计的鲁棒性和精度。
表 2 特征提取与跟踪对比Table 2. Features extraction and tracking comparison算法 特征点数量/个 匹配数量/个 匹配成功率% M7 M9 M7 M9 M7 M9 原图 273 308 89 122 32.6 39.6 Retinex算法 368 390 242 269 65.8 69.0 自适应Gamma校正算法 330 367 250 286 75.8 77.9 CLAHE算法 329 356 276 303 83.9 85.1 文献[30]算法 347 374 297 331 85.6 88.5 本文算法 365 392 331 365 90.7 93.1 3.3 多源传感器数据融合SLAM试验
为了验证激光雷达−相机−惯性自适应融合SLAM方法在低照度、弱纹理、退化场景中的表现。将本文方法与图像增强LVI−SAM（Image Enhancement LVI−SAM, IE−LVI−SAM）、激光雷达惯性测深和测绘（Lidar Inertial Odometry and Mapping，LIO−Mapping）、面向快速和旋转（Oriented Fast and Rotated Brief，ORB）的SLAM2（ORB−SLAM2）、图像增强ORB−SLAM2（IE−ORB−SLAM2）等方法从定位轨迹和建图效果2个方面进行对比。不同方法的平面轨迹对比如图7所示。从图7（a）可看出，在狭长走廊中，ORB−SLAM2方法面对低照度的走廊白墙缺乏足够的视觉特征用于约束位姿，尺度估计失效，定位轨迹明显缩小；IE−ORB−SLAM2方法提高了图像亮度和对比度，提取特征点和匹配成功率有所增加，较ORB−SLAM2方法的定位轨迹更接近控制点；LIO−Mapping方法由于走廊纵向末端提取的点云几何特征稀疏，其特征点云出现误匹配，导致LiDAR在前进方向发生退化，定位轨迹在Y方向缩短；IE−LVI−SAM方法在图像增强的基础上，将LiDAR里程计、视觉里程计、IMU预积分约束和回环约束进行联合优化，全局定位轨迹相对较好，但在转弯处视觉和IMU位姿估计误差较大，导致融合位姿有一定漂移。本文方法定位轨迹较IE−LVI−SAM方法更平滑。从图7（b）可看出，在煤矿巷道场景中，ORB−SLAM2方法在低照度喷浆墙面上提取特征较少，加上翻越轨道时剧烈晃动，其定位轨迹偏差最大。IE−ORB−SLAM2方法虽引入图像增强算法增强图像视觉效果，但面对喷浆墙面只能提取一些不稳定的特征点，在M9控制点翻越轨道时，轨迹漂移较大。LIO−Mapping方法在巷道弱几何结构环境下，LiDAR前端线面特征约束较少，且在转弯处，以IMU状态为迭代初值的位姿估计误差较大，定位轨迹发生倾斜。IE−LVI−SAM方法在M7控制点附近增强图像的亮度和对比度，避免过暗图像对视觉里程计的影响，但后续长直喷浆墙面对视觉和LiDAR里程计仍有较大影响，联合优化位姿有一定漂移。本文方法根据位姿差值，自适应选取最优位姿，减少单一传感器失效对联合优化的影响，估计轨迹更靠近控制点。
进一步对上述定位结果进行量化分析，绝对定位误差如图8所示， 定位轨迹的标准差（Standard Deviation，STD）和均方根误差（Root Mean Squared Error，RMSE）见表3。可看出，在弱纹理墙面和低照度环境中，采用ORB−SLAM2方法时视觉里程计无法对狭长走廊与煤矿巷道进行精准的位姿估计，定位轨迹RMSE较大；IE−ORB−SLAM2方法提高了特征提取数量及匹配正确率，较好地约束视觉里程计，狭长走廊与煤矿巷道中RMSE小于3 m，但缺乏尺度信息；LIO−Mapping方法虽不受低照度的影响，但在特征相似的环境中LiDAR易发生退化，导致狭长走廊与煤矿巷道中X，Y轴发生退化；IE−LVI−SAM方法提高了视觉里程计的稳定性，但是狭长走廊与煤矿巷道中LiDAR发生退化时，不能较好地约束X，Y轴的位移量，在狭长走廊与煤矿巷道中RMSE范围为0.28～0.39 m；本文方法的RMSE小于0.19 m。
表 3 定位轨迹的标准差和均方根误差Table 3. Root mean squared error and standard deviation of positioning trajectorym 误差 狭长走廊 煤矿巷道 本文方法 IE−LVI−SAM LIO−Mapping IE−ORB−SLAM2 ORB−SLAM2 本文方法 IE−LVI−SAM LIO−Mapping IE−ORB−SLAM2 ORB−SLAM2 STD 0.06 0.27 0.78 0.89 3.39 0.08 0.10 0.24 0.57 4.05 RMSE 0.15 0.39 1.25 2.58 8.12 0.19 0.28 0.75 0.60 6.82 3.4 建图效果分析
对本文方法与IE−LVI−SAM方法在狭长走廊与煤矿巷道得到的全局地图（图9、图10）定性分析建图效果。从图9可看出，由于狭长走廊口字型场景中存在相似结构，视觉尺度估计误差较大，LiDAR在长走廊中发生退化，IE−LVI−SAM方法估计的局部轨迹发生漂移，狭长走廊转弯处点云地图出现重影现象；在煤矿巷道场景下，弱纹理的喷浆巷道导致IE−LVI−SAM方法定位轨迹向Y轴漂移，运输巷道点云也向Y轴偏移，点云地图有明显错位现象。
从图10可看出，在狭长走廊中，本文方法可较好地提取几何特征，融合多源传感器数据进行高精度定位，构建的点云地图轮廓与真实场景相似度较高；在煤矿巷道中，移动机器人进入狭长巷道时，本文方法可有效抑制LiDAR退化和提升视觉特征匹配，构建的采掘巷道、综掘巷道、运输巷道点云地图与煤矿井下巷道几何特征一致，说明本文方法在面对煤矿井下复杂环境时，可以实现较高精度的定位与建图。
4. 结论
1） 提出了一种激光雷达−视觉−惯性自适应融合SLAM方法。首先，采用基于自适应Gamma校正和CLAHE的图像增强算法提升低照度图像的亮度与对比度，提高视觉里程计的鲁棒性。其次，通过设定动态阈值来检测前端里程计的稳定性，在此基础上自适应选取移动机器人在煤矿井下的最优位姿，为后端优化迭代提供位姿初始值。最后，在后端基于滑动窗口实现多源传感器数据联合非线性优化，从而提高定位与建图的精度和鲁棒性，并有效抑制轨迹发散与地图扭曲现象。
2） 对原始图像增强前后的特征提取数量及匹配正确率进行试验验证，结果表明，自适应Gamma校正和CLAHE相融合的算法能显著提升背光区和光照区的亮度和对比度，增加图像中的特征信息，大幅提升特征点提取和匹配质量，匹配成功率达90.7%。
3） 为验证本文方法性能，在狭长走廊和煤矿巷道场景下进行试验验证，结果表明，激光雷达−视觉−惯性自适应融合SLAM方法在狭长走廊场景下可很好地提取走廊内几何特征，较好地约束三轴方向的位姿，定位RMSE为0.15 m，构建的点云地图一致性较高；在煤矿巷道场景下，激光雷达−视觉−惯性自适应融合SLAM方法动态筛选出最优位姿，可较好地补偿三轴位置偏差，定位RMSE为0.19 m，构建的点云地图可以真实地反映煤矿井下环境。
-
期刊类型引用(4)
1. 陈旭升，代勇，温承超，秦冠一，周智晨，吴佳欣. 基于多传感器融合的无人车定位方法. 沈阳理工大学学报. 2025(02): 34-40+47 . 
百度学术
2. 牟琦，梁鑫，郭媛婕，王煜豪，李占利. 边缘感知增强的煤矿井下视觉SLAM方法. 煤田地质与勘探. 2025(03): 231-242 . 百度学术
3. 崔邵云，鲍久圣，胡德平，袁晓明，张可琨，阴妍，王茂森，朱晨钟. SLAM技术及其在矿山无人驾驶领域的研究现状与发展趋势. 工矿自动化. 2024(10): 38-52 . 本站查看
4. 杨林，赵春雷. 煤矿智能化巡检机器人系统设计与现场应用. 山西焦煤科技. 2024(12): 42-46 . 百度学术
其他类型引用(0)
下载:
计量
- 文章访问数: 104
- HTML全文浏览量: 10
- PDF下载量: 16
- 被引次数: 4
