煤炭开采不断在向智能化和无人化的方向发展，想要实现这一目标就必须使综采工作面采煤机、刮板输送机和液压支架协同工作，提高生产的自动化水平，其中实现对综采装备运行姿态和位置的实时监测是非常关键的环节^[1-3]。刮板输送机不仅承担运送落煤的工作，还与底板耦合，兼作采煤机的运行轨道，使采煤机可依靠牵引机构沿着中部槽行走从而截割煤壁，它还与液压支架相连接，在支架推移装置的推动下不断向前运动^[4-5]。因此，刮板输送机作为煤炭开采的关键装备，其工作姿态的准确监测非常重要。

采煤机的定位方法有很多^[6]，如齿轮计数法、红外对射法、超声波反射法和无线传感器网络(wireless sensor networks, WSN)定位法等^[7]，但上述方法的定位精度都较低，且易受外界环境影响，而捷联惯导系统在工作时能够利用载体自身惯性信息进行定位, 不依靠外界信息，还可结合其他定位方式提高定位精度，具有很大的应用潜力。Reid等^[8]提出利用捷联惯导系统来检测采煤机运行轨迹。张智喆等^[9]分析了采煤机在刮板输送机行走时的运动关系, 构建了以采煤机行走路径来检测刮板输送机直线度的数学模型。石勇等^[10]对刮板输送机因推溜而形成的弯曲段进行了计算，实现了液压支架的自动化跟机。冯帅^[11]用实时校正方法减小了底板轨迹误差，使得采煤机定位更加准确。武培林等^[12]在刮板输送机每节中部槽侧板上都安装了倾角传感器和电位器以获取准确的工作面底板状态的变化趋势。尹力等^[13]提出采煤机俯仰角是前后两个支撑滑靴连线的夹角，俩滑靴之间包含4~6节中部槽。

现有的刮板输送机形态检测方法是在每节中部槽上均布置传感器，但这种方法需要安装大量传感器，且井下对传感器防爆性要求高，布置大量矿用传感器会导致采集信息困难、成本高、采煤工作不可靠，因此该方法在实际应用中的局限性很大。而已提出的采用捷联惯导定位技术来检测刮板输送机形态的方法^[14]是在竖直平面直接用测得的采煤机的机身俯仰角来表示刮板输送机各中部槽的角度, 在底板不平整的工况条件下误差很大。

综上，本文采用捷联惯导定位技术和滚动计算方法来检测刮板输送机各节中部槽的角度，通过测量采煤机左滑靴所在中部槽和俩支撑滑靴机身之间对应中部槽的俯仰角，根据机身长度不同，测量4~6节中部槽的俯仰角作为先验信息，结合采煤机运行轨迹中的位置和姿态信息，根据采煤机在刮板输送机上运行时的耦合关系建立刮板输送机竖直平面弯曲形态解算模型，利用该模型计算采煤机右滑靴所在中部槽的俯仰角，并将它作为已知信息随采煤机运行滚动计算，解算出刮板输送机各节中部槽俯仰角，得到刮板输送机实际形态。由于综采工作面环境恶劣，地质条件复杂，进行井下实验非常困难，因此设计了采煤机与刮板输送机联合实验平台来模拟井下底板不平整、地形起伏等实际工况。通过搭建该测试平台来进行实验验证，以测试在实际工况下基于采煤机运行轨迹解算刮板输送机形态的精度。

1 采煤机运行轨迹解算 1.1 采煤机坐标系的建立

对采煤机进行姿态解算是实现其定位的基础环节，设采煤机与捷联惯导系统固联的坐标系为载体坐标系O-X_pY_pZ_p，建立地理坐标系即导航坐标系O-X_nY_nZ_n时，将载体的重心定义为原点^[15-16]，如图 1所示。捷联惯导系统安装在采煤机上，该系统直接检测出的信息均是在载体坐标系下的，因此需要利用余弦矩阵将载体坐标系下的位姿信息转换为地理坐标系下的值。

1.2 采煤机位姿解算

捷联惯导系统中的三轴陀螺仪可以测量三轴角速度，加速度计可测得载体运行的加速度。两者输出的都是比力，经过一次积分可以得到载体的姿态角和运动速度，再经过一次积分就可以获得载体的位置。

采煤机位姿解算的方法很多，如四元数法、等效旋转矢量法和Rodrigues参数法等。四元数法通过求解包含4个变量的微分方程进行姿态更新，操作简单，便于计算，但在角振动剧烈的条件下并不适用，而等效旋转矢量法可补偿误差，在姿态更新的周期内提高算法阶数从而提高计算的准确度^[17]。因此，本文采用三子样等效旋转矢量算法对采煤机运行过程中的位姿进行求解。

1) 等效旋转矢量计算。

由文献[18]可知计算等效旋转矢量Φ的方程为：

$ \mathit{\boldsymbol{ \boldsymbol{\dot \varPhi} }} = \mathit{\boldsymbol{\omega }} + \frac{1}{2}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \times \mathit{\boldsymbol{\omega }} + \frac{1}{{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^2}}}\left[ {1 - \frac{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\sin \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}}{{2\left( {1 - \cos \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \right)}}} \right]\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \times \left( {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \times \mathit{\boldsymbol{\omega }}} \right) $

(1)

式中：ω为圆锥运动角速度。

Φ的计算周期一般比较短且其值较小，因此可以忽略Φ的高阶项，得到计算中常用的简化方程为：

$ \mathit{\boldsymbol{ \boldsymbol{\dot \varPhi} }} = \mathit{\boldsymbol{\omega }} + \frac{1}{2}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \times \mathit{\boldsymbol{\omega }} + \frac{1}{{12}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \times \left( {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \times \mathit{\boldsymbol{\omega }}} \right) $

(2)

定义α为陀螺仪的半锥角，T为姿态更新周期，则陀螺仪在采样周期内的角度变化量Δθ_i为:

$ \Delta {\mathit{\boldsymbol{\theta }}_i} = \int\limits_{t + \frac{{i - 1}}{3}}^{t + \frac{i}{3}} {\omega \left( \tau \right){\rm{d}}\tau } = \left[ \begin{array}{l} - \frac{2}{3}\omega {\sin ^2}\frac{\alpha }{2}\\ - \sin \alpha \sin \frac{{\omega T}}{6}\sin \omega \left( {t + \frac{{2i - 1}}{6}T} \right)\\ 2\sin \alpha \sin \frac{{\omega T}}{6}\cos \omega \left( {t + \frac{{2i - 1}}{6}T} \right) \end{array} \right],i = 1,2,3, \cdots $

(3)

式中：t为某一时刻，τ为采样时间。

则三子样等效旋转矢量为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \Delta {\mathit{\boldsymbol{\theta }}_1} + \Delta {\mathit{\boldsymbol{\theta }}_2} + \Delta {\mathit{\boldsymbol{\theta }}_3} + \frac{9}{{20}}\left( {\Delta {\mathit{\boldsymbol{\theta }}_1} \times \Delta {\mathit{\boldsymbol{\theta }}_2}} \right) + }\\ {\frac{{27}}{{40}}\Delta {\mathit{\boldsymbol{\theta }}_2} \times \left( {\Delta {\mathit{\boldsymbol{\theta }}_3} - \Delta {\mathit{\boldsymbol{\theta }}_1}} \right)} \end{array} $

(4)

2) 四元数更新方程建立。

设Φ为T到T+t的等效旋转矢量，姿态四元数可表示为：

$ \mathit{\boldsymbol{q}}\left( t \right) = \cos \frac{{\left\| \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \right\|}}{2} + \frac{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}{{\left\| \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \right\|}}\sin \frac{{\left\| \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \right\|}}{2} $

(5)

T到T+t的姿态更新方程为：

$ \mathit{\boldsymbol{Q}}\left( {T + t} \right) = \mathit{\boldsymbol{Q}}\left( T \right) \otimes \mathit{\boldsymbol{q}}\left( t \right) $

(6)

3) 姿态矩阵和姿态角求解。

设采煤机俯仰角为θ, 航向角为ψ, 横滚角为γ，则从载体坐标系到地理坐标系(导航坐标系)的余弦矩阵C_pⁿ为：

$ \mathit{\boldsymbol{C}}_{\rm{p}}^{\rm{n}} = \left[ {\begin{array}{*{20}{c}} {\cos \gamma \cos \psi - \sin \gamma \sin \theta \sin \psi }&{\sin \psi \cos \theta }&{\sin \gamma \cos \psi - \cos \gamma \sin \psi \sin \theta }\\ { - \cos \gamma \sin \psi + \sin \gamma \sin \theta \cos \psi }&{\cos \psi \cos \theta }&{ - \sin \gamma \sin \psi - \cos \gamma \sin \psi \cos \theta }\\ { - \sin \gamma \cos \theta }&{\sin \theta }&{\cos \gamma \cos \theta } \end{array}} \right] $

(7)

用四元数表示的姿态矩阵T_q为：

$ {\mathit{\boldsymbol{T}}_{\rm{q}}} = \left[ {\begin{array}{*{20}{c}} {q_0^2 + q_1^2 - q_2^2 - q_3^2}&{2\left( {{q_1}{q_3} - {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] $

(8)

余弦矩阵与用四元数表示的姿态矩阵的关系为：

$ \mathit{\boldsymbol{C}}_{\rm{p}}^{\rm{n}} = {\mathit{\boldsymbol{T}}_{\rm{q}}} $

(9)

联立式(7)至式(9)可解得采煤机俯仰角θ、航向角ψ和横滚角γ分别为：

$ \left\{ \begin{array}{l} \theta = \arcsin \left[ {2\left( {{q_2}{q_3} + {q_0}{q_1}} \right)} \right]\\ \psi = \arctan \left[ {\frac{{2\left( {{q_1}{q_2} + {q_0}{q_3}} \right)}}{{2q_0^2 - q_1^2 + q_2^2 - q_3^2}}} \right]\\ \gamma = \arctan \left[ { - \frac{{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} \right. $

(10)

由式(10)求取俯仰角θ、航向角ψ和横滚角γ，其值可根据文献[19]中的真值换算表确定。

4) 位置解算。

采煤机速度更新微分方程为^[20]：

$ {{\mathit{\boldsymbol{\dot V}}}_{\rm{n}}} = \mathit{\boldsymbol{C}}_{\rm{p}}^{\rm{n}}{\mathit{\boldsymbol{f}}_{\rm{p}}} - \left( {2\mathit{\boldsymbol{\omega }}_{{\rm{ic}}}^{\rm{n}} + \mathit{\boldsymbol{\omega }}_{{\rm{en}}}^{\rm{n}}} \right) \times {\mathit{\boldsymbol{V}}_{\rm{n}}} + {\mathit{\boldsymbol{g}}_{\rm{n}}} $

(11)

式中:f_p和g_n分别为加速度计测得的比力和重力加速度矢量, ω_ieⁿ和ω_enⁿ分别为地球自转角速度和采煤机位置角速度。对加速度计测得的加速度进行2次积分，可求得采煤机在各方向上的位置，结合3个姿态角便可得到采煤机的运行轨迹。

2 刮板输送机竖直面形态解算 2.1 刮板输送机竖直面形态解算方法

刮板输送机是由哑铃销连接各节中部槽形成的，因此它可以在一定范围内自适应弯曲，从而与底板耦合。采煤机俩支撑滑靴中心点连线与连线在水平面上的投影线的夹角为机身俯仰角，由于俩支撑滑靴与中部槽实时接触，所以两者之间空间位置关系会直接影响采煤机机身的俯仰角。

在采煤机运行过程中，其左右支撑滑靴分别与中部槽的铲煤板接触。在刮板输送机竖直面形态解算过程中，为了准确识别采煤机支撑滑靴所经过的每节中部槽的俯仰角，以第1节中部槽为起点，建立如图 2所示的坐标系。测量初始位置处采煤机左支撑滑靴和俩支撑滑靴之间对应的多节中部槽的俯仰角θ_i，并将它作为已知量，结合利用采煤机与刮板输送机的耦合关系建立的解算模型，可滚动计算出在采煤机运行过程中其右支撑滑靴所接触中部槽的俯仰角，根据获得的各中部槽俯仰角可解算出刮板输送机的实际形态。

2.2 刮板输送机竖直面形态解算模型

已知θ_i为中部槽的俯仰角，S为由捷联惯导系统解算出的采煤机行程，L_z为每节中部槽长度，L_h为左支撑滑靴特征点O₁到B的距离，L_n为左支撑滑靴特征点O₁到滑靴底面的距离，β为BO₁与左支撑滑靴底面的夹角，如图 2(b)所示，γ为CO₂与右支撑滑靴底面的夹角，如图 2(c)所示，因为俩滑靴结构相同，所以β=γ，根据滑靴结构可得出：

$ \beta = \gamma = \arcsin \frac{{{L_{\rm{n}}}}}{{{L_{\rm{h}}}}} $

(12)

采煤机左支撑滑靴走完n节中部槽后，其特征点O₁位于第n+1节中部槽上的p(p为采煤机左支撑滑靴特征点O₁在所处中部槽上投影点与中部槽左端点的距离)处，则此时采煤机左支撑滑靴特征点O₁的坐标为：

$ \left\{ \begin{array}{l} {X_{{O_1}}} = {L_{\rm{z}}}\sum\limits_{i = 1}^n {\cos {\theta _i}} + p\cos {\theta _{n + 1}} + {L_{\rm{h}}}\cos \beta \cos {\theta _{n + 1}}\\ {Y_{{O_1}}} = {L_{\rm{z}}}\sum\limits_{i = 1}^n {\sin {\theta _i}} + p\sin {\theta _{n + 1}} + {L_{\rm{h}}}\sin \beta \sin {\theta _{n + 1}} \end{array} \right. $

(13)

则采煤机机身俯仰角θ_JS为：

$ \tan {\theta _{{\rm{JS}}}} = \frac{{{Y_{{O_2}}} - {Y_{{O_1}}}}}{{{X_{{O_2}}} - {X_{{O_1}}}}} $

(14)

此时将俩支撑滑靴之间多节中部槽的俯仰角作为已知信息，采煤机俩支撑滑靴之间机身长度为L_JS，采煤机右支撑滑靴特征点O₂位于第n+m+2节(其中m为采煤机机身之间中部槽的节数，m=4~6)中部槽的q(q为采煤机右支撑滑靴特征点O₂在所处中部槽上投影点与中部槽左端点的距离)处，则右支撑滑靴特征点O₂的坐标为：

$ \left\{ \begin{array}{l} {X_{{O_2}}} = {L_{\rm{z}}}\sum\limits_{i = 1}^{n + m + 1} {\cos {\theta _i} + q\cos {\theta _{n + m + 2}}} + {L_{\rm{h}}}\cos \gamma cos{\theta _{n + m + 2}}\\ {Y_{{O_2}}} = {L_{\rm{z}}}\sum\limits_{i = 1}^{n + m + 1} {\sin {\theta _i} + q\sin {\theta _{n + m + 2}}} + {L_{\rm{h}}}\sin \gamma \sin {\theta _{n + m + 2}}\\ {L_{{\rm{JS}}}}{\theta _{{\rm{JS}}}} = {L_{\rm{z}}}\sum\limits_{i = n + 2}^{n + m + 2} {\cos {\theta _i} + \left( {{L_{\rm{z}}} - p} \right)\cos {\theta _{n + 1}} + q\cos {\theta _{n + m + 2}}} \end{array} \right. $

(15)

当采煤机左支撑滑靴特征点O₁位于第n+1节中部槽时，此中部槽俯仰角已知，右支撑滑靴特征点O₂位于第n+m+2节中部槽，利用式(12)至式(15)就可求出右支撑滑靴特征点O₂所在中部槽的俯仰角。以此类推，采煤机沿刮板输送机行走时，利用求解出的左支撑滑靴特征点O₁的坐标就可计算出与采煤机右支撑滑靴所接触的中部槽的俯仰角，并将计算出中部槽的俯仰角作为下一节中部槽俯仰角计算时的已知信息进行滚动求解，即可求出采煤机运动过程中所经过的所有中部槽的俯仰角。

刮板输送机的初始坐标为(X₀, Y₀)，每节中部槽长度一定，利用上述方法计算出各中部槽的俯仰角，代入式(16)就可解算出刮板输送机竖直平面的形态。

$ \left\{ \begin{array}{l} {X_g} = {X_0} + {L_{\rm{z}}}\sum\limits_{i = 1}^{n + m + 2} {\cos {\theta _i}} \\ {Y_g} = {Y_0} + {L_{\rm{z}}}\sum\limits_{i = 1}^{n + m + 2} {\sin {\theta _i}} \end{array} \right. $

(16)

2.3 刮板输送机竖直面形态解算流程

综上所述，基于采煤机运行轨迹的刮板输送机竖直面形态解算流程如图 3所示。

3 实验测试 3.1 实验平台搭建

根据实验室现有的综采三机设备等比例缩小设计了样机实验平台，如图 4所示。实验设备包括：MTi-300捷联惯导系统、SVT626T双轴倾角传感器、等比例缩小的3D打印采煤机、可弯曲刮板输送机模型以及定位主机。刮板输送机模型可以纵向弯曲1°~3°，可模拟刮板输送机处于不同弯曲形状时的采煤机支撑滑靴与中部槽的接触情况。采煤机内部装有轮胎，可模拟采煤机行走。

3.2 采煤机运行轨迹检测实验

基于样机实验平台，在采煤机左支撑滑靴上方安装捷联惯导系统和倾角传感器，开启采煤机运行开关，实时采集倾角传感器数据和捷联惯导系统的加速度、角速度信息并传输至上位机，如图 4(c)所示。选取13个特征点，各特征点的数据见表 1。

3.3 刮板输送机竖直面形态测量实验

利用搭建好的实验台，模拟底板不平的实际工况，将SVT626T双轴倾角传感器布置在已设置好弯曲形状的刮板输送机模型上，如图 4(b)所示。刮板输送机模型各节中部槽的长、宽、高分别为200，11，85 mm。每节中部槽上均标记4个点，依次将倾角传感器放置在每一段中部槽的边界点和标记点上，根据倾角传感器所测得的信息得到各节中部槽的俯仰角，通过计算得到刮板输送机竖直平面的形态。

3.4 实验结果及分析

采用等效旋转矢量法解算捷联惯导系统所测得的信息，得到采煤机沿轨道运行过程中的姿态，获得行驶过程中其机身俯仰角的变化曲线。将计算出的采煤机机身俯仰角与使用倾角传感器测得的实际角度作比较，如图 5所示，结果表明：解算角度与实际角度之间的最大误差为1.9°，计算精度较高，可以满足刮板输送机中部槽俯仰角的计算要求。

如图 6所示为刮板输送机各中部槽俯仰角与运行过程中采煤机机身俯仰角的变化曲线，从图中可发现在实际工作中采煤机机身俯仰角并不能准确反映刮板输送机中部槽的起伏变化，采煤机机身俯仰角的波动要小于各中部槽俯仰角的变化，这说明采煤机对底板起伏变化不敏感。

利用建立的实验平台进行4组实验，并利用所构建的刮板输送机竖直面形态解算模型根据测得的采煤机运行轨迹计算出刮板输送机竖直平面的形态。图 7和表 2为刮板输送机竖直面的布置形态和解算形态的比较。表 2所示的13个特征监测点涵盖了刮板输送机全部弯曲段，而且各弯曲段坡度不同，在每一个上坡、下坡行程中均已取点。在底板水平的工作面上，采煤机机身俯仰角与刮板输送机中部槽俯仰角相等；在底板弯曲的工作面上，采煤机机身俯仰角与刮板输送机中部槽俯仰角不相等。若所有弯曲段和水平段的特征监测点解算值均与实际设置值一致，则可说明刮板输送机竖直面形态解算模型具有较高的准确性。

由图 7和表 2可知，根据解算模型计算出的刮板输送机竖直面形态与设置模型基本一致，最大误差小于15 mm，可以满足实际工作的要求。

4 结论与展望

1) 本方法只需测量4~6节中部槽的俯仰角作为先验信息，根据采煤机运行中与刮板输送机的耦合关系建立数学模型，基于采煤机运行轨迹不断滚动解算，便可获得刮板输送机竖直面形态，通过实验验证，对刮板输送机竖直面形态的检测精度可满足工作要求。

2) 可实时准确地掌握刮板输送机竖直面的形态变化，避免因刮板输送机弯曲角度过大而使零部件损坏或电动机过载，使刮板输送机保持正常工作，延长其使用寿命。

3) 本方法可为采煤机下滚筒的调高提供信息，保证综采工作面的平整性。

4) 结合下滚筒截割高度可以实现下一截割循环底板曲线预测，对后续的刮板输送机和液压支架动作提供信息支持。

下一步会考虑利用捷联惯导系统测得的采煤机航向角，采用形态解析法建立解算刮板输送机水平面形态的模型，与本文竖直平面的形态解算模型相结合，得到刮板输送机三维空间的形态。