基于改进EWT-多尺度熵和KELM的球磨机负荷识别方法

2020-05-15罗小燕戴聪聪程铁栋蔡改贫刘鑫刘吉顺

化工学报 2020年3期

罗小燕，戴聪聪，程铁栋，蔡改贫，刘鑫，刘吉顺

（江西理工大学机电工程学院，江西赣州341000）

引言

球磨机安全稳定工作是保障磨矿系统高效运行的关键，其负荷状态检测技术日益受到重视[1-2]。当前磨矿作业中球磨机内部的工作环境非常复杂多变，所以想要对球磨机负荷进行识别是非常困难的。相比球磨机轴承振动信号，磨矿过程筒体产生的振动信号中蕴藏着丰富的负荷信息，因而分析筒体振动信号与负荷之间的关系，实现球磨机负荷状态的准确识别，不仅能高效降低磨矿过程产生的能耗，提高生产效益，而且对磨机进行优化控制具有很好的指导意义[3-5]。由于球磨机磨矿原理就是利用钢球对物料的不断冲击来达到研磨的目的，所以此过程具有很大的随机性，导致其筒体振动信号有较强的非线性、非平稳性。汤健等[6]采用EMD、EEMD 和HVD 算法对磨机筒体振动信号分解，并求取多个子信号的频谱特征进行融合，建立频谱特征与磨机负荷参数的映射模型，通过提出IBBSEN 方法验证了模型具有泛化能力强等优势；赵立杰等[7]构建了基于振动信号EMD 分解的磨机负荷参数软测量模型，通过互信息法选择IMF 频谱特征，采用KPLS 选择性集成方法实现多尺度频谱特征的信息融合，并实验验证了方法的有效性。尽管EMD 和小波变换在处理非线性信号具有一定的优势[8-9]，但也存在模态混叠及小波基函数选取等问题，所以此种方法应用在分析球磨机振动信号上还有很多的不足之处。

2013 年，Gills[10]提出了EWT 理论，以小波分析原理为基础，通过对信号频谱自适应分割，利用正交小波滤波器将信号分解为多个模态分量，每个模态分量都能含有不同频率特征信息的，从而实现信号特征提取[11-13]。由于EWT算法是在成熟的小波和EMD 分解理论基础上得到的启发，相较于小波技术和EMD，具有充足的理论基础和较高的计算效率[14]；另外EWT 根据信号内容构造经验小波，并获得多个窄带频率分量，即具有EMD 的自适应性，能够有效地分解信号，所以适用于球磨机振动信号这类非平稳信号处理，近年来EWT在其他领域也得到广泛应用[15-18]。然而，EWT 存在频谱分割不准确等问题，导致分解效果不好，还需进一步优化。

多尺度熵于2002年由Costa等[19]提出，它可以有效地反映信号的复杂度与时间尺度的相关性，多尺度特征分析具有精度高、可行性强等特点，能够很好地描述信号在不同尺度上的相关特征[20-23]，但多尺度熵特征在磨机负荷检测中还很少被应用。然而，如何提高由多尺度熵值构成的特征向量进行负荷分类的准确率又是个关键。Huang 等[24]将核方法引入到ELM 中，提出核极限学习机(KELM)，解决了ELM算法随机初始化的问题，但KELM算法中，核参数的选取一般是人为设置，具有随机性，导致分类精度较低。目前，核参数选取方法主要通过交叉验证法，该方法需要反复训练测试，当被优化的核参数较多时，计算量大且算法难以实现[25-27]。核排列通过调节不同的核参数值，找到判断矩阵与理想矩阵的最小距离，来获取核排列值最大，此时的核参数最优[28-29]。因此，本文尝试利用EWT 对磨机振动信号进行分解，考虑到信号具有瞬态性和多样性的特点，在EWT 基础上提出了一种新的频谱分割方法，对各个IMF 分量进行相关性分析提取有效IMF分量进行重构，再求出重构信号的多尺度熵特征最后通过构建的KTA-KELM 分类识别模型实现磨机不同负荷振动信号的准确识别。

1 改进EWT-多尺度熵特征提取算法

1.1 改进经验小波变换

EWT 是将小波变换和EMD 结合起来的一种提取信号显著模态的时频分析方法。EWT 方法最主要的是如何分割傅里叶频谱，但球磨机不同负荷振动信号频域特性较为复杂，依靠基于局部极大值的频谱分割方法可能出现过分割或欠分割的问题。鉴于此，提出了一种可靠的频谱分割方法，它的主要流程如下：

(1) 获取时域离散信号f(n)，由于在截取记录时，样本选择长度不当或者其他的外界原因（传感器的零点漂移，基础运动等引起的信号波形偏移）导致趋势项的产生，因此首先利用最小二乘法对信号进行去趋势项处理;

(2) 利用傅里叶变换得到输入信号的频谱F(w)，并归一化到[0,π]；

(3)然后利用三次样条插值对信号的频谱求取包络线，搜寻信号频谱中的所有局部极大值，将极大值进行降序处理A1≥A2≥…≥AM，然后以r =AM+ α(A1- AM)为阈值对信号进行裁剪，以阈值与包络线的交点以及包络线的极小值点为频谱分割点wn；

(4)小波函数ψn(w)和尺度函数φn(w)利用分割的频谱进行构造；

(5) 应用傅里叶反变换计算F(w) × ψn(w)和F(w) × φn(w)，从而可得到各个分量的时域表示。

1.2 多尺度熵及多尺度熵偏均值算法

1.2.1 多尺度熵样本熵的不足之处是仅仅只能够反映单一尺度上的原始时间序列的复杂程度，因此，为反映原始时间序列在不同尺度下熵的自相似性，多尺度熵(multi-scale entropy，MSE)在原始时间序列下进行粗粒化处理，得到原始时间序列在多个尺度下的变化，算法步骤如下。

(1) 若某一时间序列为 {x(i) =x(1),x(2),…,x(N)}，长度为N，通过设置嵌入维数m和相似容限r，得到一组粗粒序列：

式中，τ = 1,2,… 为尺度因子，一般τmax= 15～20。

当τ = 1 时，yj(1)就对应原始序列，对于非零整数τ，时间序列Xi被割分成N/τ 个长为τ 的粗粒序列yj(τ)。

(2)计算出每个粗粒化后的时间序列样本熵，然后将其画成尺度因子τ 的解析函数，因此定义多尺度熵，即

1.2.2 多尺度熵偏均值多尺度熵偏均值(PMMSE)是建立在多尺度熵值和样本熵基础上的一个变化趋势的综合性指标，可以更加全面和完整地反映出信号的复杂度。具体计算步骤如下：

(1) 假设某一原始时间序列X ={x(n),n =1,2,…,N}，计算该组数据的多尺度熵记为：

(2)计算偏斜度Ske，计算公式为：

式中，MSEa、MSEb和MSEc为多尺度熵的均值、中位数和标准差。

(3) 原始序列的多尺度熵偏均值可根据式(5)求得

1.2.3 基于改进EWT-多尺度熵的球磨机振动信号特征提取根据对改进的EWT 算法和多尺度熵理论的研究，同时结合球磨机振动信号的特点，得到基于改进EWT和多尺度熵偏均值的特征提取算法。具体步骤如下。

(1)对采集到的磨机振动信号经改进EWT 算法分解后得到IMFi(i = 1,2,…,n)；

(2) 通过式(6)求出每个IMF 分量与原始信号的相关系数，再结合式(7)求出的阈值来筛选敏感的IMF分量。

相关系数阈值计算公式为[30]：

当前，随着企业的发展，项目全过程的预算在房产公司中的应用越来越多，应用十分普遍，通过房产公司的不断总结经验，累积与沉淀，已经探索了一条以房产项目费用为中心，对资金流动的管控为主要措施的全过程预算形式。即项目的成本费用是主要的管控目标，成本费用的支出计划为基点，对成本的进出作为重点管控为主要措施，费用的考核评价为考核标准的项目管理过程与成本费用管控方式。

其中，uh为阈值，ui为第i 个IMF 分量与原始信号的相关系数。选取相关系数值大于阈值uh的IMF分量，作为有效的IMF分量。

(3)根据选择出的有效IMF 分量进行重构，得到磨机不同负荷状态下的重构振动信号；

(4)计算重构振动信号的MSE和PMMSE。

2 基于KLEM的磨机负荷识别模型

2.1 核排列

KTA 的核心思想是：通过对核函数与学习目标之间相似性的度量，选取具有相关度最大的核函数，从而达到核参数优化的目的。假设存在某个数据集T，包含了l 个样本，定义核函数的核矩阵[K]i,j= K(xi,xj)，目标矩阵为Y，则核排列A(K,Y)为：

通过替换不同的核函数来获取不同的核矩阵，选择使A(K,Y)最大化的核参数值，此时，核矩阵与理想目标矩阵间距离最短，如式(10)所示。

2.2 核极限学习机

KELM是在ELM基础上通过用核函数代替隐含层的权值矩阵来进行计算，然而核参数选取直接决定了KELM 模型的分类精度，例如Gauss 核函数，如式(11)所示。

因此，通过KTA 算法计算所选择的每个核参数所对应的核排列的值，来寻找最优核参数σ，从而提高KELM模型分类精度。

2.3 KTA-KELM 识别模型建立

基于KTA-KELM 球磨机负荷状态识别模型的建立，具体步骤如下：

(1)将采集的球磨机不同负荷状态下的振动信号经改进EWT 分解，通过相关系数法选取敏感IMF分量重构信号，并计算其多尺度熵，构建20 维的特征向量集。

(2)初始化KTA算法参数：核排列值Ai和高斯核函数的径向宽度σi，利用式(9)计算Ai对应的σi，进而求σi的核矩阵及核排列值Aj。当Aj＞Ai时，则更新Ai和σi并重新计算，当Aj≤Ai时，则σi为最佳核参数。

(3) 利用最佳核参数构建KELM 的状态识别模型，流程图如图1所示。

图1 KTA-KELM 磨机负荷状态识别模型建立流程图Fig.1 Flow chart of KTA-KELM mill load state identification model

3 改进EWT的模拟仿真信号分析

为了验证改进EWT 方法的提取信号特征分量的能力，构造如式(12)所示仿真信号x(t)进行考察，n(t)为白噪声，信噪比为3，t ∈[0,1]，并与传统EWT和EMD 方法进行对比研究，它的时域波形如图2所示。

图2 仿真信号x(t)的时域波形Fig.2 Time domain waveform of simulation signal x(t)

经验小波分解是通过一个尺度函数和N-1 个小波函数分别滤波后得到，在频谱中寻找所有的极大值，将极大值进行降序处理，依据设定的模态个数来确定保留的峰值个数，以两个相邻极大值的中点作为分界点对频谱进行划分，然后利用尺度函数和小波函数进行滤波。图3 为改进EWT 的频谱分割图及相应的小波滤波器组，其中α = 0.1，图4为传统EWT的频谱分割图及相应的小波滤波器组。图5为仿真信号x(t)经改进EWT 和传统EWT 分解结果，图6为EMD分解结果。

图3 改进EWT频谱分割图及相应的小波滤波器组Fig.3 Improvement of EWT spectrum segmentation map and corresponding wavelet filter banks

图4 传统EWT的频谱分割图及相应的小波滤波器组Fig.4 Spectrum segmentation of traditional EWT and corresponding wavelet filter banks

图5 改进EWT和传统EWT分解结果(红色虚线代表原始信号，蓝色实线代表分解结果)Fig.5 Improved EWT and traditional EWT decomposition results

图6 EMD分解结果Fig.6 EMD decomposition results

图5(a)中分量f2～f5分别对应于信号x3(t)～x1(t)，由图5(a)可看出，信号包含的噪声被很好地分解出来了，且各个分量的吻合度都非常高，并且x3(t)下的两个模态也被独立地分解出来，克服了模态混叠现象。而图5(b)传统经验小波变换能将噪声分解出来，但是分量x1(t)、x2(t)、x3(t)都出现畸形。这是由于传统的EWT分割方法过于简单，在分析局部噪声或非平稳信号时，由噪声和非平稳分量产生的一些局部极大值可能出现并错误的保持在峰值序列中，而一些有用的极大值可能不保持在峰值序列中，导致了不当的分割。而改进的EWT 利用包络谱表示光谱的变化趋势，这可以减少噪声和非平稳分量的影响，大大增加了频谱分割的可靠性。结合图6 可以看出，EMD 分解中几个分量出现了严重的模态混叠现象，另外分解出了多个低频分量，这些低频分量原本属于同一分量的部分信息，但因EMD 终止条件的不合理性，导致过分解。因此造成混叠现象的原因一方面是终止条件的不合理性，另一方面是当信号的时间尺度存在跳跃性变化时，也会产生模态混叠现象，而混叠现象会影响后期的特征提取。

通过比较两种方法的仿真信号可知，EMD 分量虽然可以根据信号自动估算分解的层数，但分解的分量较多，而且这些分量不具有明确的物理意义，这样既消耗计算时间又影响算法性能。而改进EWT 方法能够有效地分解出原始信号的不同分量，更好地抑制了模态混叠的产生，分解效果要优于EMD和传统EWT方法。

4 磨机负荷识别实验分析

4.1 数据采集

本次实验采用型号为φ330 mm×330 mm 的Bond指数干式球磨机作为实验对象，DH5922N动态数据采集仪、DH131 振动加速度传感器、DH5857-1电荷适配器、变频器、电能表、圆柱直筒型万向轴承以及PC机等组成筒体振动信号采集系统，采样频率设置为20 kHz，均匀地改变负荷参数充填率和料球比，采集三种不同负荷下筒体振动信号的多个样本。实验装置如图7 所示。其中磨机负荷划分为：充填率10%～20%为欠负荷，充填率20%～40%为正常负荷，充填率40%～60%为过负荷。

图7 实验装置Fig.7 Experimental device diagram

4.2 振动信号的分解与重构

限于篇幅，选取球磨机不同负荷状态下筒体振动信号分别为欠负荷(料球比0.3，充填率15%)、正常负荷(料球比0.5，充填率30%)、过负荷(料球比0.7，充填率50%)，如图8所示。

图8 原始信号Fig.8 Original signal

从图8 可以看出，三种负荷状态下的筒体振动信号存在着大量的噪声，且在时域内观察信号幅值虽然存在一定的差异，但变化规律不明显，因此采用改进EWT 算法对不同负荷各10 组筒体振动信号进行预处理，每组信号得到10 个IMF 分量，计算各IMF 分量与原始信号的相关系数值，结果如图9所示。

图9 IMF分量与原始信号的相关系数Fig.9 Correlation coefficient between IMF component and original signal

图9中各负荷状态下样本相关系数的平均阈值分别为0.2742(欠负荷)、0.2431(正常负荷)、0.2399(过负荷)。由误差棒图可以看出，不同负荷状态下振动信号经改进EWT 分解得到的IMF 分量相关系数值变化趋势不同，各样本IMF 分量相关系数值误差上下限都较小，将欠负荷状态下IMF2、IMF7、IMF8、IMF10确定为敏感模态分量，正常负荷状态下IMF2、IMF3、IMF6、IMF7、IMF10确定为敏感模态分量，过负荷状态下IMF1、IMF2、IMF3、IMF10确定为有效模态分量，其余模态分量则视为虚假分量。因此，选取以上包含原始信号特征最多的IMF 分量进行重构，重构信号如图10所示。

比较图8 和图10 可以看出，不同负荷重构信号都有效地去除了高频噪声，且信号的幅值及趋势变化得到了完整的保留，所以原始信号的丰富特征信息可以体现在重构信号中。以信噪比为评价指标，对改进EWT算法的去噪效果进行分析，计算结果如表1所示。

由表1 可以看出，欠负荷的重构振动信号的信噪比提高了17.45 dB，正常负荷的重构振动信号的信噪比提高了13.59 dB，过负荷的重构振动信号的信噪比提高了18.30 dB，表明改进EWT 算法对筒体振动信号的去噪效果很好。

表1 不同负荷信号去噪效果比较Table 1 Comparison of denoising effects of different load signals

图10 重构信号Fig.10 Reconstructed signal

为了进一步验证改进EWT算法的优越性，分别利用EMD、EWT算法对欠负荷状态下原始振动信号分解，然后利用相关系数法选取敏感IMF 分量进行重构，重构信号如图11所示。

通过对比图8(a)、图10(a)和图11 可知，通过改进EWT 算法处理得到的筒体振动信号明显比EMD、EWT 算法去除了更多的高频噪声，信号更光滑，且信号的幅值得到了很好的保留，为了进一步量化EMD 和EWT 算法的去噪效果，分别计算信噪比，结果如表2所示。

表2 不同算法去噪后的信噪比Table 2 Signal to noise ratio after denoising with different algorithms

由表2 分析可知，EMD 去噪后筒体振动信号的信噪比为12.65 dB，EWT 去噪后信噪比为18.47 dB，改进EWT 去噪后信噪比远远大于EMD 和EWT 算法，表明改进EWT算法对磨机筒体振动信号的去噪效果最好。

4.3 重构信号的特征提取

每种球磨机负荷状态取5 组样本信号，然后计算各个样本的样本熵值，如表3所示。

表3 3种负荷状态振动信号的样本熵值Table 3 Sample entropy values of vibration signals in three load states

由表3 分析可得，磨机不同负荷状态筒体振动信号的样本熵的平均值差异较大，同种负荷状态下样本熵值在平均值上下波动，且较为稳定。对比三种不同负荷筒体振动信号样本熵值发现，欠负荷状态的样本熵值相对较大，这是由于欠负荷状态下筒体内的钢球和矿料相对较少，筒体运转时，矿料和钢球在下落过程中与其他矿料、钢球以及筒壁发生碰撞，能量主要被用于钢球与筒壁和钢球与钢球之间的碰撞，所以产生的筒体振动信号复杂程度较高。而在过负荷状态下样本熵值相对较小，是由于在过负荷状态下磨机筒体内的填充率过高，使磨机在运转时筒体内的物料和钢球只能进行蠕动，所以采集到的振动信号随机性小，即复杂程度低。正常负荷状态下，磨机的能量主要进行研磨，所以产生的振动信号复杂程度比较适中。

图11 欠负荷筒体振动信号重构Fig.11 Vibration signal reconstruction of under-loaded cylinder

在正常负荷与过负荷两种状态下样本熵值存在相差不大、个别出现交叉重叠的问题，区分效果不明显。因此，对于磨机振动信号的分析引入多尺度熵，根据多尺度定义可以看出，计算过程中需要考虑三个重要参数:数据长度N、尺度因子τ 以及相似容限r。现分别讨论各参数对球磨机振动信号多尺度熵计算的影响。本文研究选取τmax= 20。

在尺度因子τ = 20以及相似容限r=0.2的条件下，研究不同负荷下振动信号多尺度熵随振动信号不同数据长度N的计算结果，如图12所示。

从图12 中可以看出，数据长度N=100000 时，多尺度熵值变化趋于稳定，且欠负荷在数据长度为100000 时，多尺度熵值在0.8～0.845 范围波动，正常负荷在0.53～0.58 范围下变化，过负荷在0.35～0.42范围下波动，说明三种负荷状态存在一定的区分度，但当数据长度大于或小于100000 时，多尺度熵值波动变化较为复杂，过负荷与正常负荷会出现在尺度为3 下的多尺度熵值重叠，因此选取振动信号数据长度N=100000较为合适。

在数据长度N = 100000 以及m=20 的条件下，研究不同负荷下振动信号多尺度熵随相似容限r 的计算结果，如图13所示。

由图13 分析可以看出，随着相似容限r 逐渐增大，各负荷状态的多尺度熵值先逐渐减小后趋于稳定，当相似容限r=0.2 时，三种负荷状态多尺度熵值变化范围较小为0～0.05，且三种负荷状态区分明显，对噪声的敏感程度较高，因此，选取相似容限r=0.2。

本文最终取N = 100000，τ = 20，r = 0.2 进行相应的多尺度熵和多尺度熵偏均值计算。图14 为不同算法预处理后三种负荷状态重构信号的多尺度熵值变化。图15 为不同算法下三种负荷状态重构信号的多尺度熵偏均值变化。

比较分析图14可知：直接进行多尺度熵值计算时，三种负荷状态在多个尺度上出现重合，几乎无法分类；EMD+MSE 的负荷分类效果得到改善，过负荷和欠负荷在尺度为1、2、3、4、5、6、8、9、10、11、12、14、15、16、17、18、19、20 上能够较好地区分，欠负荷与正常负荷在尺度为2、3、5、8、9、11、16、17 和20 上可以很好地区分，但过负荷与正常负荷也基本无法分类；EWT+MSE 分类效果进一步得到提升，过负荷和欠负荷可以完全区分开，过负荷与正常负荷在尺度为1、2、3、4、5、6、7、8、9、10、11、12、16、17、20上可以完全区分，在尺度为13、14、15、18、19上也无法区分，欠负荷与正常负荷在尺度为1、2、4、5、7、8、9、10、11、13、14、15、16、17、18、19、20 上有很好的区分度；改进EWT-MSE 在任何尺度下，三种不同负荷状态都能很好的区分开，分类效果明显由于其他3 种方法，且多尺度熵值关系为：欠负荷＞正常负荷＞过负荷。

通过图15的对比分析，可知多尺度熵偏均值可以很好的解决多尺度上的混叠现象，使磨机三种负荷状态能够清晰地区分开来。改进EWT-PMMSE在样本数为20 时，可以很好地区别开；而EWTPMMSE 在样本为第9、10、14、16、17 个时存在重叠，不能很好地区分，但总体来看，三种不同负荷筒体振动信号的多尺度熵偏均值大小关系也表现为：欠负荷＞正常负荷＞过负荷。

4.3 负荷状态识别与分类

图12 不同负荷下振动信号多尺度熵随不同数据长度N变化Fig.12 Variation of multi-scale entropy of vibration signal with different data length N under different loads

每种负荷状态训练样本选取60 组，共计180组；测试样本各30 组，共计90 组。将三种不同负荷状态振动信号提取出的多尺度熵值特征向量进行归一化处理，并将欠负荷、正常负荷、过负荷的类别标签分别设置为1、2、3，通过KTA-KELM 算法的学习建立球磨机负荷状态识别模型。为了更好地证明KTA-KELM 算法的稳定性，求取总体准确识别率，将其与KELM 算法进行对比，如图16 所示。其中KELM 的核参数为高斯核函数，核半径设置为0.5，两种算法对不同负荷状态样本执行30 次。

图13 不同负荷下振动信号多尺度熵随相似容限r变化Fig.13 Multi-scale entropy of vibration signal varies with similar tolerance r under different loads

由图16 可以看出，KTA-KELM 算法的稳定性高，随算法执行次数增加，总体准确识别率的变化波动范围较小，为95.3%～97.6%；KELM 算法因核参数固定选择为0.5，每次算法执行时会造成欠拟合和过拟合现象，导致识别准确率不稳定，波动范围较大为91%～96.3%。

图14 不同算法下三种负荷状态重构信号的多尺度熵值变化Fig.14 Multiscale entropy variation of three load state reconstruction signals with different algorithms

图15 不同算法下三种负荷状态重构信号的多尺度熵偏均值Fig.15 Multiscale entropy partial means of three load state reconstruction signals with different algorithms

最后利用KTA-KELM 进行球磨机负荷状态诊断，并与SVM 算法的诊断结果进行对比。其中SVM也采用高斯核函数。图17为负荷状态识别结果。

图中红色*表示测试样本预测类别，蓝色○表示测试样本实际类别，若蓝色和红色重合时，表示负荷识别正确，否则错误。从图17 中可以看出，KTA-KELM 算法在正常负荷状态下出现2 个错判，在过负荷状态下出现1 个错判；SVM 算法在正常负荷状态下出现3 个错判，在过负荷状态下出现2 个错判，但两者都在欠负荷状态下完全识别。因此，KTA-KELM 的负荷识别准确率较高，相比SVM提高了3.4%。