基于惯性测量单元的刮板输送机形态监测
2023-04-29魏东李祖旭司垒谭超王忠宾梁斌肖俊鹏
魏东 李祖旭 司垒 谭超 王忠宾 梁斌 肖俊鹏
摘要：刮板输送机作为综采工作面的核心运输装备，准确感知其形态是提升其带载能力、缓解传动冲击、改善综采工作面直线度的重要前提。目前常用的刮板输送机形态间接测量方法难以准确表征其形态，导致测量模型误差较大。针对该问题，采用惯性测量单元直接测量刮板输送机中部槽原始位姿信息，实现刮板输送机形态数据的准确获取。采用融合Heursure阈值规则和新阈值函数的小波阈值去噪方法滤除中部槽运动加速度信号中的噪声干扰，在此基础上分析了中部槽运动特征，设计了基于随机森林的中部槽运动状态识别模型，根据运动状态识别结果采用不同的策略更新中部槽位置，减小了随时间累计的 IMU 数据误差，提升了 IMU 位置解算精度。设计了改进哈里斯鹰优化（HHO）算法优化无迹卡尔曼滤波（UKF）进行中部槽姿态解算，通过实验验证了该方法解算的姿态角满足中部槽姿态测量要求。搭建了刮板输送机形态监测实验平台，对基于运动状态识别和改进 HHO 优化 UKF 的刮板输送机形态解算方法进行实验验证，结果表明：刮板输送机进行单次推溜且步距为250 mm 时，由10节中部槽组成的刮板输送机在底板水平工况下，X，Y 轴方向上位移的最大累计误差分别为6.4，8.4 mm，Z 轴方向上位移始终保持不变，俯仰角、横滚角和航向角的最大累计误差分别为?0.148，?0.035，0.457?;在底板起伏工况下，X，Y，Z 轴方向上位移的最大累计误差分别为6.6，11.5，6.9 mm，俯仰角、横滚角和航向角的最大累计误差分别为?0.540，?0.157，0.817?。该方法可有效抑制累计误差，降低测量误差，实现刮板输送机形态的准确感知。
关键词：刮板输送机；形态监测；惯性测量单元；位姿解算；运动状态识别
中图分类号： TD634.2 文献标志码： A
Shape monitoring of scraper conveyor based on inertial measurement unit
WEI Dong1， LI Zuxu1， SI Lei1， TAN Chao1， WANG Zhongbin1， LIANG Bin1， XIAO Junpeng2
（1. School of Mechatronic Engineering， China University of Mining and Technology， Xuzhou 221116， China；
2. Inner Mongolia Jarud Zhahanaoer Coal Industry Co.， Ltd.， Tongliao 029114， China）
Abstract： Scraper conveyor is the core transportation equipment of the fully mechanized working face. Accurately perceiving its form is an important prerequisite to enhance its carrying capacity， alleviate the transmission impact， and improve the straightness of fully mechanized working face. The commonly used indirect measurement methods for the shape of scraper conveyors are difficult to accurately characterize their shape， resulting insignificant measurement model errors. To address this issue， an inertial measurement unit is used to directly measure the original pose information of the middle trough of scraper conveyor， achieving accurate acquisition of the shape data of scraper conveyor. A wavelet thresholding denoising method that combines Heursure threshold rules and a new threshold function is used to filter out noise interference in the acceleration signal of the middle trough. Based on this， the motion features of the middle trough are analyzed， and a middle trough motion state recognition model based on random forest algorithm is designed. Based on the motion staterecognition results， different strategies are used to update the position of the middle trough. It reduces the accumulated IMU data error over time and improves the precision of IMU position calculation. The improved Harris hawk optimization （HHO） algorithm unscented Kalman filter （UKF） is designed for middle trough attitude calculation. It is verified through experiments that the attitude angle calculated by this method meets the requirements of middle trough attitude measurement. The experimental platform for shape monitoring of scraper conveyors is constructed. It conducts experimental verification on the shape calculation method of scraper conveyors based on motion state recognition and improved HHO optimized UKF. The results show that when the scraper conveyor performs a single sliding with a step distance of 250 mm， the maximum cumulative errors of displacement in the X and Y directions of the scraper conveyor composed of 10 middle troughs are 6.4 mm and 8.4 mm respectively under the horizontal working condition of bottom plate. It remains unchanged in the Z direction. The maximum cumulative errors of pitch angle， roll angle， and heading angle are ?0.148°， ?0.035°， and 0.457° respectively. Under the working condition of floor undulation， the maximum cumulative errors of displacement in the X， Y， and Z directions are 6.6 mm， 11.5 mm， and 6.9 mm respectively. The maximum cumulative errors of pitch angle， roll angle， and heading angle are ?0.540°， ?0.157°， and 0.817° respectively. This method can effectively suppress cumulative errors， reduce measurement errors， and achieve accurate perception of the shape of the scraper conveyor.
Key words： scraper conveyor; shape monitoring; inertial measurement unit; pose calculation; motion state recognition
0 引言
在煤炭开采中，可弯曲重型刮板输送机是综采工作面的核心运输装备[1]。随着我国大型煤矿数量增多，对井下运输设备的各项性能提出了更高要求[2]。刮板输送机逐渐向着长距离、大运量和大功率方向发展，其稳定性是影响煤矿安全高效生产的关键因素。若刮板输送机相邻中部槽间的水平夹角过大，不仅会减少哑铃销的使用寿命及刮板输送机机头电动机的负载，还会加大采煤机的运行阻力，且若中部槽起伏过大，在采煤作业过程中也会导致采煤机截割误差不断累计[3]。良好的刮板输送机形态既可为采煤机提供合理的运行轨迹，又可为液压支架的移架动作提供精确的基准，有利于减小采煤机的运行阻力和截割误差，提高液压支架的调直精度[4-6]。
目前，刮板输送机形态监测方法主要有间接测量和直接测量2种。传统间接测量方法主要通过监测采煤机的运行轨迹或推移千斤顶的位移来实现。文献[7-9]将惯导系统解算出的采煤机位姿信息映射到刮板输送机坐标系下，并计算出推移误差反馈给液压支架控制器，对刮板输送机的调直过程进行补偿；文献[10]采用起始点处第4—6节中部槽的俯仰角作为先验信息，滚动计算出采煤机运行过程中所经过的各节中部槽的俯仰角，从而获得刮板输送机竖直方向上的形态。间接测量方法使用的传感器数量少且测量方式简单，但存在监测参数不完整、测量结果误差大、实时性差等缺点[11]。
在综采工作面中，刮板输送机一般由80～200节中部槽组成，刮板输送机形态实则是各节不同位姿状态的中部槽的整体表征[12]。因此，直接测量各节中部槽的位置和姿态是实现刮板输送机形态监测更有效的方式。文献[13]基于中部槽结构尺寸和装配约束，建立了刮板输送机形态测量模型，利用振动传感器捕捉采煤机经过相邻中部槽连接间隙时产生的振动信号，分割出每节中部槽的位姿数据。文献[14]提出在相邻中部槽间安装 FBG 张拉裂隙传感器，根据 FBG 波长与中部槽角度的线性关系及 FBG 波长的蠕变性能，实现中部槽直线度实时监测。文献[15]提出一种基于视觉测量的综采工作面直线度测量方法，利用梯形窗口匹配技术，将复杂的综采工作面直线度监测简化为刮板输送机电缆槽边缘的护板位置标定，从而实现刮板输送机直线度监测。
尽管国内外学者对刮板输送机形态监测进行了大量研究，但仍存在以下问题：
1）由于装配关系的限制，难以建立如采煤机位姿轨迹等间接测量信息与刮板输送机形态间的直接映射关系，所以，采用间接测量方法无法准确表征刮板输送机形态，所建立的刮板输送机形态测量模型误差较大。
2）刮板输送机中部槽的位姿变化主要发生在强冲击和噪声干扰大的推溜过程中，其他时间段处于基本静止状态，具有静止时间长、运动时间短和噪声干扰大的特点，若直接利用惯性测量单元（InertialMeasurement Unit，IMU）的输出数据进行位置解算，会在中部槽处于基本静止状态时，因噪声干扰产生较大的累计误差。
针对上述问题，本文以刮板输送机中部槽为研究对象，利用 IMU 直接获取中部槽原始位姿信息，提出了一种基于运动状态识别的中部槽位置解算方法，设计了一种基于改进哈里斯鹰优化（Harris Hawks Optimization，HHO）算法优化无迹卡尔曼滤波（Unscented Kalman Filter，UKF）的中部槽姿态解算方法，并通过实验验证了所提方法可提高刮板输送机中部槽的姿态解算精度，进而实现刮板输送机形态的有效监测。
1 基于 IMU 的刮板输送机形态测量
1.1 刮板输送机形态测量坐标系建立
为准确监测刮板输送机形态，需建立合适的测量坐标系。将 IMU 固定在中部槽销耳附近，IMU 和中部槽可视为刚体，其坐标系为载体坐标系OXbYbZb，如图1所示。 Yb 垂直指向工作面且平行于第1节中部槽宽度方向，Xb平行于第1节中部槽铺设方向，Zb满足右手规则，向上为正，此处OXbYbZb为“前右上”坐标系。导航坐标系OXnYnZn用于确定中部槽在基准坐标系下的导航参数，即位置和姿態，通常为“东北天”坐标系，其中Xn，Yn，Zn 分别指向东、北和天方向。
为将各节中部槽的位姿数据表述在同一坐标系下，建立基准坐标系。惯性导航中往往使用惯性坐标系（原点在地心）作为基准坐标系，但由于综采工作面长度为100～200 m，而地球半径为6400 km，刮板输送机所在地面曲率很小，不需像飞行器一样以地心为原点建立惯性坐标系，可在工作面某一固定位置建立适用于刮板输送机形态监测的固定坐标系。以第1节中部槽初始位置左下顶点为原点，以其所在“东北天”坐标系建立坐标系，并将其定义为基准坐标系OwXwYwZw，进而建立刮板输送机形态测量坐标系，如图2所示。
1.2 中部槽位姿实时解算矩阵建立
为了实时获取各节中部槽的位姿信息，需要实现姿态旋转矩阵的实时解算。根据刮板输送机的工作环境特点，选用四元数法进行姿态旋转矩阵的实时更新。四元数由4个元 q0?q3构成，其定义为
四元数可看作是复数域的扩展，q0和 qv 分别为其实部和虚部，i，j，k 为基向量。如果 q2和 q3均为0，则四元数退化为二维空间下的纯复数。
在实时解算过程中，四元数随姿态变化不断更新，通过四元数微分方程（式（2））可解算出每一时刻下的四元数。
式中：M（Q）为由四元数组成的系数矩阵；ω为载体坐标系相对于导航坐标系的角速度向量。
式中ωX，ωY，ωZ 为三轴角速度。
四元数微分方程的矩阵形式可表示为
通过式（4）可求得四元数的实时数据。采用四元数表示载体坐标系到导航坐标系的姿态旋转矩阵：
l'q0（2）+ q1（2）- q2（2）- q3（2） 2（q1q2- q0q3） 2（q1q3+ q0q2）」I
' I
利用姿态旋转矩阵可实现加速度的去重处理及姿态角的实时更新。
1.3 刮板输送机形态表征
为实现刮板输送机中部槽在基准坐标系的表征，需先确定导航坐标系下各节中部槽的初始位姿信息，将其映射到基准坐标系下，再结合位置和姿态的累计特性即可求解。定义 t 時刻第i节中部槽测量点在基准坐标系下的位置为[Xi（t） Yi（t） Zi（t）]姿态角为[θi（t）γi（t）ψi（t）]，θi（t），γi（t），ψi（t）分别为 t 时刻第i节中部槽的俯仰角、横滚角、航向角。通过测量手段确定各节中部槽在基准坐标系下的初始位置信息，结合初始姿态角解算方法可获得中部槽在导航坐标系下的初始姿态角。本文建立的基准坐标系基于“东北天”坐标系，可将导航坐标系下的初始姿态角视为基准坐标系下的初始姿态角。设第i节中部槽的初始位置为[Xi（0） Yi（0）Zi（0）]T，初始姿态角为[θi（0）γi（0）ψi（0）]T。在初始位姿的基础上，结合提出的运动状态识别方法和改进 HHO 优化 UKF 对 IMU 原始位姿信息进行坐标转换和累计计算，便可不断更新各节中部槽在基准坐标系下的位置和姿态，实现刮板输送机的形态表征。
2 基于运动状态识别的中部槽位置解算
针对 IMU 输出的加速度信号中的噪声干扰，采用小波阈值去噪方法减小噪声影响，结合运动信号的特征信息，建立基于随机森林（Random Forest，RF）算法的中部槽运动状态识别模型。基于运动状态识别结果，采用不同的位置更新策略，实现中部槽位置的精确解算。
2.1 中部槽惯性传感数据去噪
IMU 安装于刮板输送机中部槽上，加速度计的高频敏感特性导致其输出信号中含有大量干扰信号，如刮板与槽体间的碰撞摩擦、槽体与底板间的摩擦、落煤冲击等引起的噪声干扰。为降低噪声干扰对中部槽位置解算精度的影响，选用小波阈值去噪方法处理采集到的加速度信号，降低噪声干扰，为中部槽运动状态识别和位置解算提供准确的数据参考。小波阈值去噪方法通过小波变换对小波函数进行时间平移和尺度伸缩，实现对非平稳信号的多尺度时频特性分析，可分辨信号中更多的细节信息，适用于非平稳信号的去噪处理[16]。
在小波阈值去噪中，选取合适的阈值估计方法和阈值函数是决定去噪效果的关键。Heursure阈值规则以其去噪力度大、灵活性强、不易发生过保留的特点，更加适合于复杂运动环境下物体运动信号的去噪处理。对于小波基函数和分解层数的选取，需要结合信号的具体特性。在实验室环境下，多次采集中部槽的运动信号进行分析，得出小波基函数为 sym4、分解层数为5时，分解出的小波系数最符合信号的原始特性。因此，选择 sym4为小波阈值去噪的基函数，分解层数为5，结合新阈值函数和Heursure阈值规则进行小波阈值去噪。
2.2 中部槽加速度信号特征提取
误差随时间累计是 IMU 存在测量误差的主要因素。本文通过分析推溜过程中中部槽的运动状态，研究加速度信号的多种时频特征，结合中部槽的运动状态特点，提取出能够有效区分中部槽运动状态的信号特征，为后续中部槽运动状态识别提供适配的特征输入。
2.2.1 中部槽运动状态分析
以单推溜过程中的 Y 轴加速度信号为代表，进行中部槽运动状态分析。利用实验室测试平台模拟中部槽静止?推溜?静止的运动过程，采集整个过程中 IMU 加速度计数据，采样频率为100 Hz。设定中部槽静止时间为10 s，推溜时间为5 s，模拟中部槽真实推溜步距0.8 m。对采集的 IMU 加速度信号进行小波阈值去噪，利用原始积分算法进行位置解算，结果如图3所示（以 Y 轴方向为例， g 为重力加速度）。
从图3可看出，中部槽在0～10 s处于静止状态；10 s 后开始执行推溜动作，10～10.15 s加速度信号出现正向突变，中部槽从静止状态变为加速运动状态；10.15～15.15 s 中部槽处于匀速运动状态；15.15～15.28 s 加速度出现反向突变，表示中部槽由匀速运动状态进入减速运动状态；15.28 s 后中部槽再次进入静止状态，推溜动作结束。上述运动过程中，中部槽的运动状态依次为静止、加速、匀速和减速，分别定义为 S1?S4。实际中部槽的运动过程是这4种运动状态的顺序转换。由于中部槽被一直推向煤壁，仅凭突变信号的正负便可判断出中部槽在 Y 轴方向的 S2和 S4状态，但考虑到底板起伏工况下中部槽在 X 轴和 Z 轴的运动方向不确定，本文定义 X 轴和 Z 轴方向上的 S2状态包括正向加速和反向加速， S4状态包括正向减速和反向减速。
分析图3中位移曲线可知，在 S1状态下，由于去噪后的加速度信号并非真实信号，依然含有残余噪声，对加速度积分解算出的速度不为0，导致解算结果中出现位移变化，与中部槽实际运动状态不符。经历 S2，S3，S4状态后，中部槽再次进入静止状态，由于加速度信号并非理想信号，导致速度难以回归至0，即使中部槽已进入 S1状态，也无法保证解算出的位移量为0，即位置仍在变化。
实际工况中，S1状态对应的时间段在中部槽4种运动状态中占比最大，通常情况下占单次截割周期的98%以上，该状态下直接通过积分运算求解中部槽位置会造成较大的累计误差。因此，有必要对中部槽的运动状态进行识别，根据识别结果采用不同的位置解算策略，重点在于抑制 S1状态下的位置累计误差，提高整个刮板输送机形态监测过程中对中部槽位置的解算精度。
2.2.2 加速度信号特征提取
根据中部槽静止时间长、运动时间短及各阶段信号变化明显的特点，可提取能够表征中部槽运动状态的加速度信号特征，对中部槽的4种运动状态进行识别。假设采集的加速度信号為 n 维时间序列，对其可提取以下特征：均值特征 F1、高度特征F2、脉冲因子特征 F3、偏度特征 F4、峰度特征 F5、标准差特征 F6、平均能量特征 F7、均方根频率特征F8、频率标准差特征 F9和能量熵特征 F10。
为验证上述10种特征在中部槽运动状态识别中的适用性，并筛选出能够有效区分中部槽运动状态的加速度特征，利用实验测试平台，采集100个推溜周期的加速度数据（以 Y 轴加速度数据为例），每个周期内的加速度序列长度为2500（100 Hz×25 s），对各阶段特征值进行归一化处理，结果如图4所示。
由图4可知，特征 F1—F10均能区分 S2和 S3这2个状态，但特征 F5，F6不能有效区分 S3和 S42个连续状态，特征 F9不能有效区分 S1和 S2这2个连续状态。特征 F1—F4，F7，F8，F10均能在连续时间序列上对4种运动状态进行区分，但特征 F1，F2，F10在跨时间序列上存在2个状态特征值近似的情况，F3，F4，F7，F8特征值在4种状态间的区分度明显，且 F3，F4，F7，F84种特征包含了信号的时域特征和频域特征分析，可信度较高。因此，选用 F3，F4，F7，F84种特征构成特征向量，作为后续状态识别算法的输入值。
2.3 中部槽运动状态识别
由2.2节分析可知，中部槽的运动过程是 S1，S2，S3，S44种运动状态间的顺序转换，可通过特定的运动信号特征进行区分，将中部槽的运动状态识别问题转换为分类问题。 RF 算法是一种结合多个弱分类器组成决策树的集成分类算法[17]，属于引导聚集（Bootstrap aggregating，Bagging）算法类型[18]，通过两重随机过程提升样本分类的随机性能，且多决策树决策提高了分类准确性。因此，提出了一种基于RF 算法的中部槽运动状态识别方法，根据识别出的运动状态，采用相应的位置解算策略，以提高中部槽位置解算精度。基于分类回归树（Classification andRegression Tree，CART）算法运算速度快、适用于大样本分类的特点，采用 CART 算法作为 RF 的决策树分裂算法，以特征向量[F3 F4 F7 F8]T 作为 RF 算法的输入，以 [S1 S2 S3 S4]T 为目标状态向量，作为RF 算法的输出，建立中部槽运动状态识别模型。由于中部槽的运动状态不断处于 S1?S2?S3?S4?S1 循环转换中，为提高运动状态识别的准确性，本文在 4 种运动状态的识别中添加时序约束，即 S1 仅出现在S2 之前及 S4 之后，S3 仅出现在 S2 之后及 S4 之前。
为验证中部槽运动状态识别模型的有效性，在实验室环境下进行中部槽推溜实验。由于 Y 轴方向为中部槽主要运动方向，以 Y 轴加速度信号为代表，给出运动状态识别过程。对加速度信号进行去噪和归一化处理后，采集2000组数据（采样频率为100 Hz），每种运动状态各采集500组，其特征数据见表1。1—500组对应 S1状态，501—1000组对应 S2状态，1001—1500组对应 S3状态，1501—2000组对应 S4 状态。
将上述数据的60%作为训练集，其余40%数据作为测试集，将其特征向量代入 RF 算法，运动状态识别结果见表2。可看出中部槽运动状态识别模型准确率整体较高，引入时序约束后，S1和 S3状态识别准确率达100%，S2和 S4状态因运动时间较短，数据突变程度不一致，导致识别准确率有所下降，但也可达96.4%和97.9%。
2.4 基于运动状态的中部槽位置解算
根据识别出的中部槽运动状态，采用不同的位置解算策略：中部槽处于 S1静止状态时，不进行位移累计计算，即保持上一时刻的位置不变；中部槽处于 S2，S3，S4运动状态时，对小波阈值去噪处理后的加速度信号进行积分解算。该策略抑制了长时间 S1状态下产生的位置累计误差，整体上提高了中部槽位置解算精度。
IMU 测得的加速度为载体坐标系下的比力信息，在进行位置求解时需要利用姿态旋转矩阵转换至导航坐标系下，以去除重力影响。定义载体坐标系和导航坐标系下的加速度分别为fb =[ fX（b） fY（b） fZ（b）]T和fn =[ fX（n） fY（n） fZ（n）]T，则有
式中ln 为平移向量。
IMU 实际输出为离散数据，定义 k 时刻第i节中部槽的加速度为fiw（k），位置和速度信号为xiw（k）， xi（w）（k）=[Xi （k） Yi （k） Zi （k） viX （k） viY （k） viZ （k）]T ，[Xi （k） Yi （k） Zi （k）]T，[viX （k） viY （k） viZ （k）]T 分别为 k 时刻第i节中部槽的位置、速度，则 k+1时刻的位置和速度信号为
xi（w）（k+1）=Ψxi（w）（k）+Ωfiw （k+1） （7）
式中Ψ， Ω为系数矩阵。
3 基于改进 HHO 优化 UKF 的中部槽姿态解算
为提高姿态解算精度，利用 UKF 算法作为陀螺仪、加速度计和磁力计的数据融合算法[19-20]，结合非线性收敛因子和动态惯性权重，提出一种基于改进 HHO 优化 UKF 的中部槽姿态解算方法，通过优化 UKF 的过程噪声协方差，提高中部槽姿态解算精度。
3.1 改进 HHO 算法
HHO 算法具有原理简单、参数少、搜索方式多样化等优点[21]，但在参数分布规律和寻找全局最优解方面还存在较大的改善空间，主要表现为以下2个方面。
1）HHO 算法中，猎物的逃逸能量 E 决定了算法在全局搜索和局部开发间的状态切换。E越小， HHO 算法的局部开发能力越强；E 越大，HHO 算法的全局搜索能力越强。可见 E 对 HHO 算法搜索能力的平衡至关重要。
E = E0E1 （8）
式中：E0为[?1，1]内的随机数； E1为收敛因子。
E1=2（1-d=dmax ） （9）
式中 d，dmax分别为当前迭代次数和最大迭代次数。
逃逸能量 E 主要受收敛因子 E1的影响，而 E1从2线性递减至0，使得在 HHO 算法寻优的中后期仅满足|E|≤1，导致算法在中后期仅能进行局部开发，破坏了算法局部开发和全局搜索的平衡，易出现早熟现象。
2）HHO 算法中猎物位置为当前全局最优解，影响哈里斯鹰个体位置的变化。若下次迭代完成后猎物位置并未更新，算法可能陷入局部最优，导致无法收敛到全局最优解。
针对以上2个问题，在传统 HHO 算法的基础上，引入卡方分布的逆累积分布函数改进收敛因子 E1，使其呈非线性递减，以保证 HHO 算法局部开发与全局搜索的平衡；利用自适应权重改进猎物位置的更新策略，在算法陷入局部最优时，可跳至其邻域范围内进行更新，有助于提高算法收敛的速度及跳出局部最优的能力。
1）收敛因子的非线性化。采用卡方分布的逆累积分布函数实现收敛因子的非线性变化。卡方分布是数理统计中的一种分布方式[22]，累积分布函数Pc（p|N）是对变量 p 的概率密度函数的积分，表达式为
式中：N 为自由度；Γ（·）为 Gamma 函数。
改进后的非线性收敛因子为
式中 Pc?1（y|N）为逆累积分布函数，可借助于Matlab中的 chi2inv（p，N）函数进行求解。
将改进前后的收敛因子 E1 和逃逸能量 E 变化曲线进行对比，结果如图 5 所示。
由图 5（a）可知，改进后的收敛因子 E1 呈非线性递减趋势，前期和中期收敛缓慢，算法有充足的机会寻找到全局最优解，后期收敛迅速，使算法在寻找到全局最优解后迅速收敛。逃逸能量 E 可在整个迭代过程中实现|E|＞1，使改进 HHO 算法在后期主要进行局部开發的情况下，依然存在全局搜索的可能，从而在整个收敛过程中提高了算法局部开发和全局搜索的平衡能力，有利于提高算法的收敛精度和收敛速度。
2）动态惯性权重改进更新策略。为了提高 HHO 算法跳出局部最优的能力，引入动态惯性权重wd 对猎物位置Lrabbit进行调整[23]。
wd = w（w）m（m）ax（in）-（wmax s（-） n-smin ） s（s） savg（savg） （12）
式中：wmin和wmax分别为最小权重和最大权重；s 为当前适应度；smin和savg分别为最近10次的猎物位置对应的适应度最小值和平均值。
根据上述动态惯性权重，对猎物位置进行更新：
Lrabbit （d+1）= Lrabbit （d）+ wdLrabbit （d） （13）
式中Lrabbit（d）为第 d 次迭代后猎物位置。
如果Lrabbit对应的适应度超过3次未被更新，则使用较大的权重扰动Lrabbit，使其跳出当前位置，在邻域内重新搜索猎物位置，提高了算法跳出局部最优的能力，反之使用较小的权重更新Lrabbit。使用动态惯性权重可加速算法收敛，提高算法跳出局部最优的能力。
3.2 基于改进 HHO 优化 UKF 的中部槽姿态求解
中部槽姿态解算系统为非线性系统，利用 UKF 求解中部槽姿态具有较高的滤波精度，但 UKF 采用传统卡尔曼滤波算法的滤波框架，使用过程噪声协方差表示数学模型与真实系统间的近似程度，对滤波精度具有重要影响。在实际应用中，过程噪声协方差常取经验值，但中部槽姿态解算系统的外部环境可能发生变化，过程噪声协方差也应随之改变，以保证系统模型的可靠性。对此，采用改进 HHO 算法优化 UKF 中的过程噪声协方差，提高中部槽姿态解算精度。
3.2.1 中部槽初始姿态角确定
中部槽的三轴姿态向量为[θγψ]T，初始姿态角[θ0 γ0 ψ0]T 可由 IMU 的加速度计和磁力计获得。
1）初始俯仰角和横滚角的确定。静止状态或匀速状态下的加速度fb =[ fX（b） fY（b） fZ（b）]T 与导航坐标系下的重力向量f0n =[00 ?g]之间的关系为
式中Cn（b）为导航坐标系到载体坐标系的姿态转换矩阵。
将[θ0 γ0 ψ0]T 代入Cn（b），得
由此得出初始俯仰角θ0和初始横滚角γ0：
上述加速度计输出数据需进行归一化处理。
2）初始航向角的确定。九轴 IMU 包含地磁计，可用来求解初始航向角ψ0。设归一化后的磁力计在载体坐标系下的输出为[mX（b） mY（b） mZ（b）]T，为得到航向角，可保持ψ=0并将载体沿俯仰角和横滚角方向分别旋转θ和γ， 使磁力计平面与导航坐标系平行。该过程所对应的姿态旋转矩阵为
设导航坐标系下的磁场强度为[mX（n） mY（n） mZ（n）]T，可得
经过处理可得水平方向上的2个磁场分量：
由于 X 轴和 Y 轴方向的向量和指向磁北，可得磁航角：
将求得的初始俯仰角θ0和横滚角γ0代入式（20），可得初始航向角：
式中Δψ为磁北与真北之间的磁偏角。
3）初始四元数的确定。使用四元数更新姿态旋转矩阵，因此需确定初始四元数。已知欧拉角与四元数间的转换关系为
将θ0，γ0，ψ0代入式（22），即可获取初始四元数[q0，0 q1，0 q2，0 q3，0]T。
3.2.2 中部槽姿态角确定
1）状态方程的建立。中部槽姿态解算系统中的状态信息为姿态角，选用四元数为系统状态信息。通过求解四元数微分方程（式（23））计算四元数。
式中ωbX，ωbY，ωbZ分别为陀螺仪的三轴过程偏差。
定义系统的状态向量为 Ak=[q0，k q1，k q2，k q3，k]T，建立系统的离散化状态方程：
式中：G（·）为非线性状态方程；Wk?1为过程噪声，满足 e（Wk?1）=0，e（·）为期望。
式中：T 为两状态的时间间隔；ωbX，k，ωbY，k，ωbZ，k为 k 时刻陀螺仪的三轴过程偏差。
将式（24）展开整理后可得
2）觀测方程的建立。为融合加速度计、磁力计与陀螺仪的姿态信息，需分别以加速度计和磁力计的测量值为观测量，建立双通道的观测方程。
将归一化后的 k 时刻三轴加速度计数据[fX（b）;k fY;（b）k fZ;（b）k ]T作为 UKF 的观测量，单位为 g。此时
观测方程为
式中：h（·）为非线性观测方程；Vk为观测噪声，满足 e（Vk）=0及 e（VkVkT）=Rk，Rk为观测噪声协方差，与测量误差及传感器特性有关，本文通过统计长时间采集的加速度计和磁力计静态数据得出。
将式（27）展开可得
已知载体坐标系和导航坐标系下的磁场强度分别为[mX（b） mY（b） mZ（b）]T和[mX（n） mY（n） mZ（n）]T 。一般认为导航坐标系下仅有北向和天向存在磁场数据，东向无磁场分布，即水平方向上的磁场强度为 X 轴和 Y 轴方向的向量和，所以导航坐标系下的磁场强度可表示为
因此，磁力计数据为观测量时的观测方程为
3）过程噪声协方差的定义。过程噪声协方差Uk?1=e（Wk?1 WkT-1），将式（25）代入可得
式中σωX，k，σωY，k，σωZ，k为过程偏差的标准差。
3.2.3 中部槽姿态解算
定义改进 HHO 算法的适应度函数为
式中：A 为中部槽的真实状态；A（?）k 为中部槽的估计状态。
基于改进 HHO 优化 UKF 的中部槽姿态解算过程如下：
步骤1：分3步完成参数的初始化设置。
1）利用安装在中部槽上的 IMU，采集其在上电2 s 内的数据进行均值处理，代入式（16）和式（21）求得初始姿态角，根据式（22）求得初始四元数状态值A（?）0，初始协方差矩阵 P（?）0= e[（A0- A（?）0）（A0- A（?）0）T ]，计算初始标准差，并设置陀螺仪的初始过程偏差。
2）确定无迹变换的相关参数κ， λ， α， β。
3）计算 Sigma 点的均值权值和方差权值。
步骤2：对 k?1时刻的后验最优状态估计A（?）k-1进行 Sigma 点采样，得到 Sigma 点集。
步骤3：计算 Sigma 点集的预测值。
步骤4：进行 UKF，获得先验状态估计和先验协方差估计。
步骤5：对先验状态估计再进行一次 Sigma 点采样，得到新的 Sigma 点集。
步骤6：将新的 Sigma 点集代入式（28），得到观测预测值。
步骤7：计算观测预测值均值、方差矩阵和协方差矩阵。
步骤8：更新卡尔曼滤波增益、四元数状态和协方差矩阵。
步骤9：初始化 HHO 算法相关参数。
步骤10：根据式（8）和式（11）计算猎物的逃逸能量。
步骤11：根据猎物的逃逸能量，进入不同的搜索阶段进行位置更新。
步骤12：根据动态惯性权重更新策略，更新当前全局最优位置。
步骤13：若迭代次数达到最大值，则结束优化过程，并输出全局最优位置，否则跳转至步骤11继续执行。
步骤14：根据全局最优位置更新过程噪声协方差。
步骤15：将更新后的四元数状态转换成欧拉角，得到俯仰角、横滚角、航向角。
步骤16：如果 UKF 迭代次数达到终止条件，则结束循环，否则跳转至步骤2继续执行。
以上步骤对应的流程如图6所示。
3.2.4 中部槽姿态解算验证
在实验室环境下，采用旋转平台和中部槽模型对中部槽姿态解算方法进行实验验证，实验装置如图7所示。采集安装于中部槽上的 IMU 原始数据，分别采用改进 HHO 优化 UKF、扩展卡尔曼滤波（Extended Kalman Filter，EKF）、UKF、HHO 优化 UKF 4种算法进行姿态解算，比较4种算法的解算效果。
根据刮板输送机的工作原理可知，形成 S 弯区域的过程中，中部槽的航向角在三轴姿态角中变化最为明显。为达到更好的验证效果，采用小角度转动进行实验设计，通过指令控制旋转平台绕轴往复转动3?。中部槽模型的具体转动过程：0～50 s模型保持静止，50～55 s模型逆时针水平旋转3?， 55～108 s 模型保持静止，108～113 s模型顺时针水平转动3?， 113～168 s模型保持静止，168～173 s模型顺时针水平转动3?， 173～227 s模型保持静止，227～232 s模型逆时针水平旋转3?， 232～316 s模型保持静止，316～322 s模型逆时针水平旋转3?， 322～351 s模型保持静止，351～356 s模型顺时针水平旋转3?， 此后模型保持静止。
实验中，UKF 参数设置：κ=0，α=0.01，β=2；改进 HHO 优化 UKF 参数设置：最大迭代次数为250，种群规模为30，搜索范围为[?1，1]，最小、最大权重分别为0.95，0.45。此外，为了对比姿态的动态变化过程，以传感器内置算法解算出的姿态角作为参考值。实验结果如图8所示。
根据实验结果，统计4种算法解算结果的最大误差和平均绝对值误差，以此评价4种姿态解算算法的精度，结果见表3。
结合图8和表3对4种算法的解算结果进行以下分析。
1）從图8（a）可看出，4种算法的航向角解算结果均与中部槽实际运动趋势一致。在0～113 s 4种算法的解算结果均与参考值接近，但113 s之后，EKF、UKF、HHO 优化 UKF 3种算法的解算结果与参考值存在较大偏差，而改进 HHO 优化 UKF 解算结果与参考值之间的偏差较小。整个航向角测量中，改进 HHO 优化 UKF 解算结果的最大误差为0.193?， 平均绝对值误差为0.057?， 在4种算法中精度最高。
2）从图8（c）和图8（e）可看出，在水平旋转过程中，传感器输出的俯仰角和横滚角波动较大，但均接近0，这是因为传感器自身解算算法存在误差所致。结合所用实验平台水平旋转的稳定性及算法输出结果均接近0，可推导出俯仰角和横滚角的理论值应保持在0附近。因此，4种算法的解算结果均优于传感器自身的解算结果，这是因为算法中对多传感器数据进行了融合。
3）若均以0作为俯仰角和横滚角参考值，在横滚角测量中，改进 HHO 优化 UKF 解算结果的最大误差为?0.003?， 平均绝对值误差为5.188×10?4 ?。在俯仰角测量中，改进 HHO 优化 UKF 解算结果的最大误差为0.010?， 平均绝对值误差为5.805×10?4 ?， 在4种算法中误差最小。
可见，改进 HHO 优化 UKF 在中部槽姿态解算中精度最高，符合中部槽的姿态测量要求。
4 实验验证
4.1 实验平台搭建
为进行可重复性实验，参照 SGZ?800型刮板输送机进行中部槽模型的简化设计。模型尺寸为300 mm×200 mm×56 mm（长×宽×高），可实现相邻中部槽在水平和垂直方向上的最大允许夹角为3?。采用10节中部槽模型模拟刮板输送机的不同工作形态，利用电动推杆模拟液压支架推移千斤顶的推溜动作，电动推杆行程为0～800 mm，推移速度在0～180 mm/s 之间可调，通过 PLC 实现中部槽顺序推溜控制。安装于中部槽模型上的数据采集单元实时获取各节中部槽的角速度、加速度和磁力计数据，传送至对应的 Powerlink 从站中进行位姿解算，然后将解算结果传送至 Powerlink 主站中进行整合，最终将数据通过以太网 UDP 协议传送至上位机的刮板输送机形态监测系统进行数据更新。刮板输送机形态监测实验方案及平台分别如图9和图10所示。
4.2 实验验证
4.2.1 底板水平工况下刮板输送机形态监测实验
通过底板水平工况下的刮板输送机形态监测实验，验证本文刮板输送机形态监测技术的可行性。真实工况下中部槽的推移步距为800 mm，但考虑中部槽模型的整体推溜效果，设定最大推移步距为250 mm。起始时刻各节中部槽静止，然后顺序推动中部槽直至形成 S 弯区域，如图11所示。
基准坐标系下各节中部槽测量点位置如图12所示。黑色曲线为基准坐标系下各节中部槽测量点的初始位置连线，可利用激光雷达测得；蓝色和红色曲线分别为激光雷达和本文方法得到的测量点在终止位置的连线，分别作为参考值和解算值。可看出中部槽推溜前后由直线形态变为 S 弯形态。
为便于对各方向上的位置进行对比分析，将终止状态下的中部槽测量点位置投影到 XOY 和 XOZ 平面内，并分析其误差，如图13所示。
从图13（a）可看出，一次推溜之后，S 弯区域主要由第2?9节中部槽构成，且第2?5节中部槽和第6?9节中部槽弯曲方向大致相反。从图13（b）可看出，Z轴方向上测量点位置的参考值在±0.2 mm 内波动，反映出中部槽所处底板较为水平，符合设定的实验场景，因此，Z 轴方向上中部槽运动不明显。本文所提基于运动状态识别的中部槽位置解算方法将 Z 轴方向的运动状态识别为静止（S1状态），解算出的10节中部槽的 Z 轴位置均保持在19.0 mm处。从图13（c）和图13（d）可看出，10节中部槽在 X 轴方向位置的最大累计误差为6.4 mm，Y 轴方向位置的最大累计误差为8.4 mm，满足实际工况下的位置精度要求（刮板输送机拉线误差不超过±100 mm）。
在中部槽推溜前后，采用高精度倾角传感器测量中部槽的初始姿态角和终止姿态角，将其作为参考值，与改进 HHO 优化 UKF 解算出的中部槽姿态角进行对比，结果如图14、图15所示。
图14（a）、图14（b）和图14（c）分别为由 HHO 优化 UKF 解算的三轴姿态角和初始参考值的对比，图14（c）中由于引入磁力计修正航向角，可测出中部槽的方位，所以初始航向角在20?附近，而俯仰角和横滚角在0附近。图14（d）为10节中部槽的初始姿态角误差，俯仰角、横滚角、航向角最大误差分别为0.012，0.006，0.145?， 俯仰角和横滚角误差较小，原因是改进 HHO 优化 UKF 引入了加速度计进行角度补偿。
从图15（a）和图15（b）可看出，各节中部槽的俯仰角和横滚角均小于0.1?， 且相邻中部槽间相对转角较小，可得出 S 弯形态下的刮板输送机在垂直方向上基本对齐，符合底板水平工况。从图15（c）可看出，改进 HHO 优化 UKF 解算的航向角符合实际 S 弯区域的变化趋势。值得注意的是，该算法解算的第1节中部槽的初始姿态和终止姿态未发生变化，原因是第2节中部槽的转动角度全部位于中部槽模型允许的转动范围内，相邻中部槽间没有出现联动现象，第1节中部槽相对静止，图12也可验证这一现象。由图15（d）可知，10节中部槽终止状态下俯仰角、横滚角、航向角的最大累计误差分别为 ?0.148，?0.035，0.457?。实验结果表明，改进 HHO 优化 UKF 的解算精度高，满足中部槽姿态解算的精度要求。
4.2.2 底板起伏工况下刮板输送机形态监测实验
为了更加全面地验证本文刮板输送机形态监测方法的有效性，设计了底板起伏模拟实验平台，利用表面起伏的硬质胶皮模拟实际工作面的截割底板，如图16所示。各节中部槽在该底板上做推溜运动，三轴姿态角均会发生明显变化，其余实验条件与底板水平工况实验一致。
在基准坐标系下，推溜前后的各节中部槽测量点位置如图17所示，蓝色曲线为激光雷达测得的测量点终止位置的连线（参考值），红色曲线为本文方法解算的测量点位置的连线（解算值）。
为便于对各方向上的位置进行对比分析，将终止状态下的中部槽测量点位置投影到 XOY 和 XOZ 平面内，如图18所示。
从图18（a）可看出，一次推溜运动后，10节中部槽呈现 S 弯形态，符合中部槽的实际形态分布，其中第 2?9 节中部槽均处于弯曲段。从图 18（b）可看出，各节中部槽测量点处存在高度差，各节中部槽处于不同程度的起伏状态，最大高度差为 84 mm。
10 节中部槽的位置累计误差如图 19 所示。可看出三轴方向上的位置累计误差整体呈递增趋势，这是由于2?10节中部槽依次完成推溜，随着中部槽编号增大，推溜距离增大，推溜时间变长，尽管加速度信号经过小波阈值去噪，但中部槽在运动过程中残存的噪声仍会带来一定误差，该误差随时间逐渐累计增大。 Y 轴方向是中部槽的主要运动方向，因此整体累计误差较大，且呈明显上升趋势，而 X 轴和 Z 轴方向由于底板起伏，误差出现波动现象，但整体呈上升趋势。X，Y，Z 三轴方向的位置最大累计误差均出现在第10节中部槽处，分别为6.6，11.5，6.9 mm，满足刮板输送机形态监测的精度要求。
将高精度倾角传感器测得的各节中部槽初始姿态（参考值）与本文方法的解算值进行对比，结果如图20所示。图20（a）、图20（b）和图20（c）中的角度波动均可反映出实验所用模拟底板呈起伏状态。由图20（d）可知，俯仰角、横滚角和航向角的最大累计误差分别为0.055，0.006，0.301?。
推溜结束后，各节中部槽的终止姿态角如图21所示。
图21（a）和图21（b）中的角度波动可反映刮板输送机模型呈起伏状态。由图20（c）可知，改进 HHO优化 UKF 解算出的航向角符合实际 S 弯区域的变化趋势。由图20（d）可知，本次实验中10节中部槽终止俯仰角、横滚角、航向角的最大累计误差分别为?0.540，?0.157，0.817?。由此可知改进 HHO 优化 UKF 的求解精度较高，满足中部槽姿态解算的精度要求。
5结论
1）采用间接测量方法无法准确表征刮板输送机形态，所建立的刮板输送机形态测量模型误差较大，因此，建立了基于 IMU 的刮板输送机形态测量模型，采用 IMU 直接测量各节中部槽的位姿信息，实现了刮板输送机形態数据的准确获取。
2）针对中部槽运动加速度信号中的噪声干扰问题，提出了融合Heursure阈值规则和新阈值函数的小波阈值去噪方法。在此基础上，分析了中部槽的运动特征，设计了基于 RF 算法的中部槽运动状态识别模型，根据运动状态识别结果，采用不同的策略更新中部槽位置，有效降低了随时间累计的 IMU 数据误差，提升了 IMU 位置信息的解算精度。
3）设计了基于改进 HHO 优化 UKF 的中部槽姿态解算方法，提升了刮板输送机中部槽的姿态解算精度。实验结果表明：采用该方法解算出的中部槽俯仰角、横滚角、航向角的最大误差分别为0.010，?0.003，0.193?， 平均绝对值误差分别为5.805×10?4，5.188×10?4，0.057?，满足中部槽姿态测量要求。
4）搭建了刮板输送机形态监测实验平台，进行刮板输送机形态监测实验数据采集与分析，结果表明：在底板水平工况下，X 轴和 Y 轴方向上位置的最大累计误差分别为6.4，8.4 mm，Z 轴方向上始终保持位置不变，俯仰角、横滚角和航向角的最大累计误差分别为?0.148，?0.035，0.457?;在底板起伏工况下，X， Y，Z 轴方向上位置的最大累计误差分别为6.6，11.5，6.9 mm，俯仰角、横滚角和航向角的最大累计误差分别为?0.540，?0.157，0.817?。实验结果验证了所提刮板输送机形态监测方法的正确性和有效性，具有较高的推广应用价值。
参考文献（References）：
[1] 司垒，李嘉豪，谭超，等.矿用刮板输送机垂直冲击下负载电流特性研究[J].煤炭科学技术，2023，51（2）：400-411.
SI Lei，LI Jiahao，TAN Chao，et al. Study on load current characteristics of scraper conveyor under vertical impact[J]. Coal Science and Technology，2023，51（2）：400-411.
[2] 徐亮.我国煤炭开发建设现状与“十四五”展望[J].中国煤炭，2021，47（3）：44-48.
XU Liang. The present situation and expectation of coal exploitation and construction in China's 14th Five-Year Plan period[J]. China Coal，2021，47（3）：44-48.
[3] 毛君.刮板输送机动力学行为分析与控制理论研究[D].阜新：辽宁工程技术大学，2006.
MAO Jun. Dynamical behavior analysis and control theory research of scraper conveyor[D]. Fuxin： Liaoning Technical University，2006.
[4] 王学文，李素华，谢嘉成，等.机器人运动学与时序预测融合驱动的刮板输送机调直方法[J].煤炭学报，2021，46（2）：652-666.
WANG Xuewen， LI Suhua， XIE Jiacheng， et al. Straightening method of scraper conveyor driven by robot kinematics and time series prediction[J]. Journalof China Coal Society，2021，46（2）：652-666.
[5] 刘婷.刮板输送机S弯区域水平弯曲角优化关键技术研究[D].徐州：中国矿业大学，2019.
LIU Ting. Research on key technologies for optimizing horizontal bending angle in S-bending area of scraper conveyor[D]. Xuzhou：China University of Mining and Technology，2019.
[6] 乔春光，王学文，谢嘉成，等.刮板输送机水平面形态检测方法[J].工矿自动化，2018，44（8）：52-57.
QIAO Chunguang，WANG Xuewen，XIE Jiacheng， et al. Horizontal shape detection method for scraper conveyor[J]. Industry and Mine Automation，2018，44（8）：52-57.
[7] KELLY M S，HAINSWORTH D W. The landmark longwall automation project[J]. Acarp Project，2005，1：1-7.
[8] REID D C，HAINSWORTH D W，RALSTON J C， et al. Shearer guidance：a major advance in longwall mining[C]. The 4th International Conference on Field and Service Robotics，2003：469-476.
[9] RALSTON J C，REID D C，DUNN M T，et al. Longwall automation：delivering enabling technology to achieve safer and more productive underground mining[J]. International Journal of Mining Science and Technology，2015，25（6）：865-876.
[10] 喬春光，王学文，谢嘉成，等.基于采煤机运行轨迹的刮板输送机竖直面形态解算方法[J].工程设计学报，2018，25（5）：495-502.
QIAO Chunguang，WANG Xuewen，XIE Jiacheng， et al. Vertical plane shape calculation method of scraper conveyor based on running track of shearer[J]. Chinese Journal of Engineering Design，2018，25（5）：495-502.
[11] 刘旭，王世博，韩子晨，等.刮板输送机调直方法研究[J].煤矿机械，2020，41（6）：54-57.
LIU Xu，WANG Shibo，HAN Zichen，et al. Study on face alignment method of scraper conveyor[J]. Coal Mine Machinery，2020，41（6）：54-57.
[12] 李宁.智能工作面中刮板输送机直线度感知技术[J].机电工程技术，2021，50（3）：186-188.
LI Ning. Straightness sensing technology of scraper conveyor of intelligent face[J]. Mechanical & Electrical Engineering Technology，2021，50（3）：186-188.
[13] 王超，李威，杨海，等.基于航位推测的刮板输送机形态检测研究[J].煤炭学报，2017，42（8）：2173-2180.
WANG Chao，LI Wei，YANG Hai，et al. Scraper conveyor shape detection based on dead reckoning[J]. Journal of China Coal Society，2017，42（8）：2173-2180.
[14] 权志桥.智能工作面刮板输送机直线度FBG感知研究[D].徐州：中国矿业大学，2021.
QUAN Zhiqiao. Research on straightness FBG perception of scraper conveyor in intelligent working face[D]. Xuzhou：China University of Mining and Technology，2021.
[15] 刘鹏坤.基于视觉测量的综采工作面直线度控制研究[D].北京：中国矿业大学（北京），2020.
LIU Pengkun. Research on straightness control of fully mechanized coal face based on vision measurement [D]. Beijing：China University of Mining and Technology- Beijing，2020.
[16] DONOHO D L， JOHNSTONE I M. Adapting to unknown smoothness via wavelet shrinkage[J]. Journal of the American Statistical Association，1995，90（432）：1200-1224.
[17] BREIMAN L. Random forests[J]. Machine Learning，2001，45（1）：5-32.
[18] BREIMAN L. Bagging predictors[J]. MachineLearning，1996，24（2）：123-140.
[19] JULIER S，UHLMANN J K. Unscented filtering and nonlinear estimation[J]. Proceedings of theIEEE，2004，92（3）：401-422.
[20] 闫小龙，陈国光，田晓丽.两步快速可重构无迹卡尔曼滤波算法测量导弹滚转角[J].仪器仪表学报，2018，39（6）：140-147.
YAN Xiaolong，CHEN Guoguang，TIAN Xiaoli. Two- step fast reconfigurable unscented Kalman filter algorithm for measurement of missile rolling angle[J]. Chinese Journal of Scientific Instrument，2018，39（6）：140-147.
[21] HEIDARI A A，MIRJALILI S，FARIS H，et al. Harris hawks optimization： algorithm and applications[J]. Future Generation Computer Systems，2019，97：849-872.
[22] 翟光，王妍欣，孙一勇.基于低轨星网的多目标协同跟踪滤波技术[J].系统工程与电子技术，2022，44（6）：1957-1967.
ZHAI Guang， WANG Yanxin， SUN Yiyong. Cooperative tracking filtering technology of multi-target based on low orbit satellite constellation[J]. Systems Engineering and Electronics，2022，44（6）：1957-1967.
[23] 湯安迪，韩统，徐登武，等.混沌精英哈里斯鹰优化算法[J].计算机应用，2021，41（8）：2265-2272.
TANG Andi，HAN Tong，XU Dengwu，et al. Chaotic elite Harris hawks optimization algorithm[J]. Journal of Computer Applications，2021，41（8）：2265-2272.
