多尺度熵算法（精选7篇）
多尺度熵算法 篇1
摘要：提出了一种多尺度最大信息熵(Max Information Entropy,MIE)及梯度的图像融合算法,该算法在对源图像多尺度分解的基础上,根据低频小波系数及高频小波系数的特点,把信息熵引用到小波低频系数的选择中,根据局部信息熵的大小确定小波系数的选择;而高频采用基于最大梯度值的方法,最后对所选小波系数进行重构,即可得到融合图像。二者的结合,对图像的细节处理更加细致,又有效地消除冗余信息。通过实验分析,结果表明该算法与其他基于区域的方法相比,提高了融合效果。
关键词：多尺度分解,信息熵,梯度,融合算法
图像融合是数字图像处理研究的重要领域之一。随着数字化设备、图像处理技术的发展,图像融合技术越来越得到人们的重视。在图像融合中,采用的融合规则直接影响到图像的融合质量,这也是图像融合中至今没有得到很好解决的问题。
一般情况下,图像某一区域的局部特征并不能完全由单个像素表达,各个像素之间往往有较强的相关性。但目前已有的一些像素级图像融合算法通常根据各个图像分解层上对应元素的大小来确定融合图像相应分解层上的元素值,其融合规则一般为对应元素值选大、选小或对应元素值的加权平均。正是由于忽略了像素之间的相关性,融合效果常常不甚理想,很难达到融合要求。本文提出了一种多尺度最大信息熵和梯度的图像融合算法,在多尺度分解的基础上,把局部信息熵引用到低频系数的选择中,用局部信息熵来代替单个像素或者相邻像素均值,充分利用图像的特征信息,因此,在最后的融合结果中能更好地保留细节信息。
1 图像多尺度的分解与重构
1.1mallat算法
mallat提出了小波的快速分解与重构算法,利用两个一维滤波器对二维图像实现快速小波分解 , 低通(H)和高通(G)为两个一维镜像滤波算子 ,其下标r、c分别对应图像的行和列 ,按照mallat算法 , 则在尺度j-1有如下分解公式:
与之相应的小波重构公式为:
Cj-1=H*rH*cCj+H*rG*cCj+G*rH*cCj+G*rG*cCj(j=J,J-1,…,1) (2)
式中 ,H*、G*分别为H、G的共轭转置矩阵。
在Matlab中,可以选择小波工具箱中的Wavedec2()及Waverec2(),只需给出相应的小波基及变换级数,函数即可完成图像的小波分解与重构。
1.2 小波基的选择
小波基满足以下条件:正交性、支撑集 、对称性 、规则性 、消失矩阶数 。尽管满足上述条件的任何小波基都可以实现图像的小波分解,但并不是任何分解均能满足我们的要求,同一幅图像,用不同的小波基进行分解所得到的效果是不同的。我们希望经小波分解后得到的三个方向的细节分量具有高度的局部相关性,而整体相关性被大部分甚至完全解除。本文选用具有对称性的双正交小波bior6.8,构造具有线性相位的双正交滤波器组,进行小波的分解与重构。
1.3 小波变换级数的确定
图像融合的结果与分解和重构的变换级数有一定的关系, 级数太大会引起严重的失真,且程序的时间复杂性大幅增加;级数太小,高于所能达到的分解分辨率的细节信息将会丢失,达不到分解效果。通常,可通过下面方法确定变换级数的大小。
先定义两个变量:
(1)小波分解子图面积比
设小波分解中,第k层分解子图的大小,即该子图的面积为Sk,上一层即第k-1层分解子图像的面积为Sk-1,则面积比r为:
(2)均方根误差
设两幅源图像第k层分解的低频子带分别为Ak和Bk,则它们的均方根误差Ek为:
式中I、J为低频子带的维数。
基于上面两个变量的最佳小波分解层数选取规则如下:
①计算出小波分解子图像的面积比大于0.5的最小分解层数m和面积比小于0.5的最大层数n,m和n满足:m=n+1
②计算出对应层次的均方根误差Em和En
③确定最佳分解层数K
同时考虑到计算的复杂性,分解层数为4时效果最好。根据不同的图像,通常选取3～5层比较合适。
2 一种基于信息熵和梯度的多尺度图像融合算法
根据小波分解系数的特点,即低频分量是对其父图像的粗劣逼近,而高频分量反应了其父图像在三个不同方向的细节表示,所以,针对小波子图的不同分量,选用不同的融合准则,对小波变换后的低频分量,计算以某点为中心大小为M×N的窗口的局部信息熵,该窗口局部信息熵越大,则说明该窗口图像信息越丰富,对三个高频分量进行梯度计算,根据梯度值大小进行高频分量的系数融合,最后通过小波逆变换重构图像。具体过程如图1所示。
2.1 图像融合算法
定义AlX(x,y)为源图像X的第l级小波分解后在(x,y)的低频分量小波系数,DklX(x,y)为源图像X的第l级分解后的第k方向在(x,y)的高频分量小波系数。同理定义Y的第l级小波分解后在(x,y)的低频系数和高频系数为AlY(x,y)和DklY(x,y)。定义融合后的高频小波系数为DklF(x,y),低频小波系数为AlF(x,y)
2.1.1 基于信息熵的低频小波系数选择
根据前面的分析,低频反映了源图像的大致轮廓,经过小波分解以后,仍能看出源图像的基本信息,在一个M×N大小的窗口内,该窗口的信息熵基本反映了该窗口的信息的多少,亦即反应了以点(x1,y1)为中心图像信息量的大小。所以在低频系数的融合上,用基于信息熵的融合算法。
Setp 1: 计算以(x1,y1)为中心窗口大小为M×N的局部信息熵HX、HY:
其中ph为灰度为h的像素在窗口M×N内出现的概率
同理
Step 2: 比较两信息熵的大小,确定融合小波系数AlF(x1,y1):
Step 3: 循环扫描X、Y第l级低频子图像,求出融合图像第l级低频分量小波系数。
2.1.2 基于梯度的高频小波系数选择
小波系数的高频分量包含了图像的细节部分,根据它的特点,采用了与低频分量不同的融合准则。考虑到小波系数变化较大的值,往往表征了图像的细节,而细节正是该点梯度变化较大的地方。因而根据各个小波系数的梯度进行融合。Grad代表梯度幅度,此处用Roberts算子。
2.2 图像融合算法步骤
融合的基本步骤如下:
Step 1:对已配准的两幅源图像进行小波变换,建立图像的小波金字塔序列;
Step 2:对各分解层分别进行融合处理,各分解层上的不同频率分量采用不同的融合算子进行融合处理,最终得到融合后的小波分解系数;
Step 3:对融合后所得的小波系数进行小波逆变换 (即进行图像重构),所得到的重构图像即为融合图像。
在图像融合中,考虑到相邻像素的相关性,当组合后某图像变换系数来自于图像X变换后的系数而其相邻像素组合后的系数大多来自于图像Y变换后的系数,反之亦然,这时就要进行一次性校验,亦即对该点小波系数重新进行赋值,以确保某像素组合后的图像变换系数预期相邻像素组合后的系数大多来自于同一幅图像变换后的系数。
3 实验结果与分析
本文使用Lena作为测试图像,分别对Lena图左右部分进行高斯模糊模拟两幅多聚焦图像,为了验证本算法的性能,实验还将本文提出的融合算法与小波系数最大值融合法以及小波系数均值融合法进行比较,高频融合原则均采用梯度最大值法。结果如图2所示。
(a)左半部分模糊 (b)右半部分模糊 (c)本章算法 (d)系数最大值法 (e)小波系数均值法
从主观效果来看,小波系数最大值法及小波系数均值法的融合图像都存在轻微马赛克,本算法能够较好的保留细节部分,同时融合了两幅图像的信息,使融合后的图像得到了增强,相对系数最大值及系数均值融合算法有明显的优势。
为了更好地比较几种算法的优劣,文中对MIE算法、最大值融合算法、均值融合算法进行了实验比较,结果如表1所示。
从表中各评价参量数据来看,本文融合算法的结果包含的信息量、图像平均梯度和信噪比这些数据都有了较大的提升,说明算法较好地提高了图像的质量。
4 结束语
本文提出了基于多尺度的最大信息熵和梯度的图像融合算法,在整个过程中,利用小波变换的多分辨率的特征,对图像信息由粗到精,从整体到局部,逐层进行分析,因此对于图像的细节处理得更加细致,有效地消除冗余信息,保留了图像的细节信息,得到了全局清晰的图像。从实验的结果来看,本算法在一定程度上提高了融合的效果,和其他基于区域的融合算法相比,融合效果有明显的提高。
参考文献
[1]敬忠良,肖钢,李振华.图像融合——理论与应用[M].北京:高等教育出版社,2007:9-10
[2]LIU Gang.Liu Xue-qin.Performance measure for image fusion considering region information[J]浙江大学学报A(英文版),2007,4:168-171
[3]Youcef Chibani.Selective synthetic aperture radar and panchromatic image fusion by using the劋trous wavelet decomposition[J].EURASIP Journal on Applied Signal Processing,2005,14:2207-2214
[4]Lewis J J.Region-based image fusion using complex wavelets[C].Processing of the 7th International Conference on InformationFusion,Stockholm,2004:555-562
[5]Piella G.A general framework for multiresolution image fusion:from pixels to regions[J].Information Fusion,2003,4(4):259-280
[6]WANG Hai-hui.A new multiwavelet-based approach to image fusion[J].Journal of Mathem-atical Imaging and Vision,2004,21(3):177-192
[7]Ioannidou S,Karathanassi V.Investigation of the dual-tree complex and shift-invariant discrete wavelet transforms on quickbirdimagefusion.J-.IEEE Geoscience and Remote Sensing Letters,2007,4(1):166-170
多尺度熵算法 篇2
关键词：多尺度排列熵,偏均值,液压泵,故障特征
0 引言
分析壳体振动信号是液压泵故障诊断与预测的主要方法之一，基于振动信号分析结果可以判定液压泵故障的类型和程度[1]。对于轴向柱塞泵而言，由于液压油的压缩性、泵源与液压回路的流固耦合作用以及工作过程中泵体本身固有的机械冲击，导致液压泵的振动信号表现出很强的非平稳、非线性特性[2]。传统的小波分析、经验模态分解等方法能够很好地解决振动信号的非平稳性问题，但采用以上方法处理具有非线性特性的液压泵振动信号具有一定的局限性。
随着非线性理论的发展，许多非线性方法被应用到机械信号处理中，如分形维数、近似熵和样本熵等[3,4]。排列熵（permutation entropy,PE）是Bandt等[5]提出的一种新的时间序列复杂度指标，与Lyapunov指数、分形维数以及样本熵等相比，它在概念上更容易理解，且具有计算简单、运行速度快和对噪声鲁棒性强的优点，排列熵已被广泛应用于脑电信号、心音信号、地磁信号以及机械信号[6]处理中。多尺度排列熵（multi-scale permutation entropy,MPE）是Aziz等在排列熵的基础上提出的[7]，相关研究表明，多尺度排列熵具有比排列熵更好的鲁棒性[8]。在国内，刘永斌等[9]将排列熵用于旋转机械的状态监测，并分析了不同的延迟时间和嵌入维数对排列熵的影响。冯辅周团队对排列熵进行了更进一步的研究，将排列熵用于机械设备的状态监测和故障预测，并对排列熵的参数优化方法进行了研究[10,11]。郑近德等[12]将多尺度排列熵用于轴承的故障特征提取，并与支持向量机相结合，有效地实现了轴承的故障诊断。但现有关于排列熵的研究大都是基于排列熵良好的突变检测性能展开的，多数是采用排列熵检测某一系统的异常或者预测某一故障的发展变化趋势，将排列熵用作故障特征指标实现系统多故障识别的研究还比较少，将排列熵用于液压泵故障识别的研究更是很少见到报道。
本文将排列熵引入液压泵的故障识别中，在对液压泵振动信号排列熵和多尺度排列熵进行研究的基础上，提出了一种综合液压泵振动信号在多个尺度上排列熵值和排列熵值随尺度变化趋势的指标———多尺度排列熵偏均值（partial mean of multi-scale permutation entropy,PMMPE）作为液压泵的故障特征指标，对液压泵实测信号的分析结果验证了所提指标的有效性和优越性。
1 排列熵相关理论
1.1 排列熵算法及其参数优化
排列熵是一种衡量一维时间序列复杂度和随机性的指标，它可以很好地检测出系统的动力学突变[5]。对于一维时间序列{x(j),j=1,2，…，n}，以嵌入维数为m、延迟时间为r对其进行相空间重构，可以得到如下形式的矩阵：
重构的矩阵共包括K行，每一行是矩阵的一个重构分量。对每个重构分量中的元素按其数值大小进行升序排列，然后提取每个元素在排序前重构分量中所在列的索引组成一个符号序列。对于重构的m维列向量，可能出现的符号序列共有m！种，计算第k种排列形式的符号序列出现的概率，记为Pk，则该时间序列的排列熵可以由下式求得：
对Hp进行归一化可得
Hp的大小反映了时间序列的复杂程度和随机性。机械设备发生某种故障时，故障越严重，其振动信号的随机性越小、复杂度越低，此时振动信号的排列熵越小；反之，机械设备处于正常状态时其振动信号的随机性最大，排列熵值也最大[9]。Hp值的变化能够很好地反映机械设备故障程度的变化，常被用作机械设备状态监测和异常检测的指标。对于液压泵而言，当其处于不同的故障状态时，其振动信号的复杂度和随机性也各不相同，因此，排列熵应该可以作为液压泵的故障特征指标用于液压泵的故障识别。
根据排列熵的计算步骤可知，相空间重构时的延迟时间和嵌入维数是影响排列熵算法的两个主要参数，人为地确定这两个参数具有一定的主观性和随机性，针对该问题，饶国强等[11]对比分析了采用互信息法和伪近邻法独立确定两参数与采用关联积分法（C-C算法）联合确定两参数的效果，发现独立确定的参数求得的排列熵具有更好的突变检测效果，也即参数优化能够提高排列熵区分突变前后两状态的能力。
1.2 多尺度排列熵
多尺度排列熵是在多个尺度上计算时间序列的排列熵，求时间序列的多尺度排列熵首先要将时间序列多尺度化即粗粒化[12]。振动信号序列X={x(j),j=1,2，…，n}可以根据下式进行粗粒化：
yi(l）表示尺度为l的粗粒化序列，尺度因子l决定了时间序列的粗粒化程度，在粗粒化过程中时间序列的长度也相应缩短，尺度因子为l时，粗粒化序列长度为ent(n/l），当l=1时粗粒化序列就是原时间序列。
在对时间序列粗粒化处理后，计算每个粗粒化序列的排列熵即可得到多尺度排列熵。为了使每个尺度下的排列熵具有更好的故障识别效果，在计算每个尺度上的排列熵时有必要采用互信息法和伪近邻法优选延迟时间和嵌入维数，基于优化的参数计算各尺度排列熵。
1.3 多尺度排列熵偏均值
在关于多尺度排列熵的研究中，大部分文献都没有提出一种综合多个尺度上排列熵值的指标。文献[13]在对多尺度熵研究的基础上提出了综合时间序列在多个尺度上的非线性信息的指标———多尺度熵偏均值。该指标在计算轴承振动信号多尺度熵的基础上，结合偏均值的概念，综合了多个尺度样本熵值和熵值的变化趋势两方面的信息，能更加全面地反映振动信号所包含的信息。本文在对排列熵和多尺度排列熵研究的基础上，提出了多尺度排列熵偏均值的概念，并采用多尺度排列熵偏均值作为液压泵故障识别的特征参量，以期得到更好的故障识别效果。
某一振动信号序列X={x(j),j=1,2，…，n}的多尺度排列熵偏均值计算步骤如下：
(1）确定多尺度排列熵的最大尺度因子s。
(2）在某一尺度l(l=1,2，…，s）下对振动信号进行粗粒化，采用互信息法和伪近邻法确定该粗粒化序列的最佳延迟时间和嵌入维数，然后计算该粗粒化序列的排列熵Hp(l）；计算所有尺度下粗粒化序列的排列熵可得到X的多尺度排列熵Hmp(X）={Hp(1),Hp(2），…，Hp(s）}。
(3）计算多尺度排列熵的偏斜度Ske，即该序列的偏态绝对值与其标准差的比值，其计算公式如下：
式中，Hmmp、Hcmp和Hdmp分别表示多尺度排列熵Hmp的均值、中位数和标准差。
(4）该振动信号的多尺度排列熵偏均值可按下式求得：
2 液压泵故障特征提取
2.1 振动信号采集
实测液压泵振动信号采自液压泵试验台，液压泵型号为SY-10MCY14-1EL，驱动电机型号为Y132M-4，其额定转速为1480r/min。选用CA-YD-139型压电式加速度传感器与液压泵端盖进行刚性连接，见图1，使用DH-5920动态信号测试分析系统进行数据采集。采集正常、滑靴磨损、松靴、双松靴以及严重松靴5种故障状态下的液压泵振动信号，试验中液压泵故障采用装备检修时换下的出现故障的柱塞代替正常柱塞的方式进行模拟，试验所用部分柱塞见图2。试验中振动信号采样频率为10kHz，每种故障采集10组数据，每组数据采样点数为2048，采样间隔为1min。试验过程中试验台主溢流阀压力为10MPa，电机转速为其额定转速。
采集到的5种状态下液压泵振动信号波形如图3所示。从振动信号的波形图可以看出，不同故障模式下液压泵振动信号的幅值不同，正常液压泵振动信号的幅值最小，但根据时域波形图无法判断液压泵的故障。
2.2 基于排列熵的液压泵故障特征提取
首先分析排列熵作为特征指标区分液压泵不同故障的能力。根据文献[9]的经验，本文取排列熵计算过程中相空间重构的嵌入维数m=5，延迟时间r=3，根据排列熵计算步骤求得5种故障模式下各组样本的排列熵，如图4所示。由图4可以看出，不同故障模式下液压泵振动信号的排列熵具有不同的波动区间，且波动强度也不相同。正常信号的排列熵值最大，波动性最小；严重松靴状态下振动信号的排列熵值最小，波动性最大。排列熵能够较好地衡量不同故障模式下液压泵振动信号的复杂度和随机性。但从图4也可以看出，不同故障模式下的排列熵波动区间有一定的重叠和交叉，直接采用排列熵作为液压泵的故障特征可能会引起误判。
采用互信息法和伪近邻法对排列熵计算过程中的延迟时间和嵌入维数进行优选，基于优选的参数对信号序列进行相空间重构，然后计算其排列熵，以期得到更好的故障识别效果。限于篇幅，此处取正常信号中的第三组数据介绍其参数优选过程。首先采用互信息法确定延迟时间r，求得互信息（mutual information,MI）随延迟时间变化的曲线，如图5所示，根据互信息法确定延迟时间的规则，选定延迟时间r=4。在确定延迟时间的基础上[11]，采用伪近邻法优选嵌入维数，其中最大嵌入维数设置为8，判据一设置为20，判据二设置为2，伪近邻率（ratio of false neighbor,RFN）随着嵌入维数变化的曲线如图6所示，在嵌入维数为4处伪近邻率不再随着嵌入维数的增加而减小，则取嵌入维数m=4。
采用互信息法和伪近邻法求得的5种故障模式下各组样本的延迟时间和嵌入维数如表1所示。基于优选的延迟时间和嵌入维数对各组样本进行相空间重构，并计算其排列熵，得到优化参数下的排列熵如图7所示。
可以看出，与图4相比，图7中不同故障模式下排列熵间的区分度更好，同一故障模式下不同样本的排列熵间的差异更小，说明参数优选能够有效提高排列熵区分液压泵不同故障的能力。但是，不同故障模式下排列熵的波动区间仍有重叠现象。另外，在参数优选条件下，滑靴磨损信号的排列熵值明显变大，造成该变化的原因在于参数优选确定的延迟时间和嵌入维数明显区别于前文排列熵计算时的延迟时间和嵌入维数。以上分析表明，在计算排列熵过程中有必要对延迟时间和嵌入维数进行优选。
2.3 基于多尺度排列熵偏均值的液压泵故障特征提取
排列熵只能反映振动信号在单个尺度上的复杂度和随机性，为了衡量振动信号在多个尺度上的复杂度和随机性，并将多个尺度上的复杂度用一个指标反映出来，本文计算了液压泵振动信号的多尺度排列熵偏均值Hmppc，并对其作为液压泵故障特征指标的可行性和有效性进行分析。
按照多尺度熵偏均值的计算步骤，首先计算5种故障模式下每组样本的多尺度排列熵Hmp，此处从每种故障模式的样本中各选一组进行分析，取最大尺度因子为12，求得5组样本的多尺度排列熵曲线，如图8所示。由图8知，液压泵振动信号的多尺度排列熵熵值Hmp随着尺度因子l的增大呈现递减的趋势，这说明随着尺度的增大，粗粒化序列的复杂度和随机性降低。另外，不同故障类型振动信号的多尺度排列熵曲线具有不同的下降速率，说明不同故障振动信号随着尺度的增大其复杂度降低的速率不同。为了更好地区分液压泵故障，采用多尺度排列熵偏均值作为液压泵的故障特征指标。
计算5种故障模式下各组样本的多尺度排列熵Hmp，并计算它们的多尺度排列熵偏均值Hmppc，结果如图9所示。对比图9和图7可知，多尺度排列熵偏均值能够更好地区分液压泵故障，且具有更好的稳定性。实测液压泵振动信号的分析结果证明了本文提出的多尺度排列熵偏均值作为液压泵故障特征指标的有效性和优越性。
3 结论
多尺度熵算法 篇3
目前，基于混沌、分形、耗散结构、时域分析、频域分析及状态空间等方法对流化床内脉动行为的研究日趋增多[1,2,3,4]。赵明阳等通过高阶统计量对气-固流化床压力脉动信号进行研究，结果表明流化床压力脉动信号是非高斯性，并提出特征值T，获得了不同流态下信号偏离高斯性的程度[5];秦伟刚等利用希尔伯特-黄变换(HHT)中的经验模态分解(EMD)方法和盲信号分离-三阶累积量方法对压差信号进行分析，结果表明EMD和三阶累积量结合能有效揭示气-固两相流的压差特性[6];王春华等采用分形标度计算方法对气-固流化床压力脉动信号进行分析，研究了分形标度值与最大Lyapunov指数间的关系，得出分形标度值对于混沌特性具有表征作用的结论[7];张少峰等对压力脉动信号进行了时域、频域和自相关性分析，表明流体流动和颗粒运动所引发的压力脉动能量频带分别集中在0～10Hz和30～40Hz，压力脉动的概率密度近似呈正态分布[8];文献[9,10]对垂直上升气-液两相流电导脉动信号采用多尺度熵进行分析，根据不同尺度的多尺度熵脉动特征揭示了不同流型的动力学特性。
笔者将针对气-固两相流化床内的压力脉动信号进行多尺度熵分析，以进一步了解流化床内的流动特性，以指导气-固流化床的实际生产。
1 理论基础
1.1 样本熵
给定一个长度为N的时间序列x={x(1)，…，x(i)，…，x(N)}，其样本熵的计算步骤如下:
a．构造一个m维的模板向量xm(i)={x(i),x(i+1)，…，x(i+1-m)}，1≤i≤(N-m+1);
b．两个向量之间的距离
e．增加向量的维数从m到m+1，重复步骤b～d，得到Bm+1(r)，计算样本熵
1.2 多尺度熵
多尺度熵算法的原理与步骤如下:
a．给定一个长度为n的一维时间序列，即{u(i),i=1,2，…，n}，用公式构建连续粗粒化的时间序列;
e．增加向量的维数从m到m+1，重复步骤b～d，得到Bim+1(r,n);
g．计算粗粒化后各个尺度τ所对应的时间序列的样本熵值，即得到多尺度熵。
为便于分析，上述计算过程中序列的匹配长度m取2，阈值r为原始时间序列标准差σ的0.10～0.25倍，数据长度8 000点，最大粗粒化尺度20。
2 压力脉动信号的获取
实验装置由动力系统、循环流化床与压力检测系统组成，如图1所示。循环流化床主体部分立管的横截面积为0.12m×0.70m，床高2.50m。床底材料为直径0.32mm的石英砂，密度为2 600kg/m3。压力脉动信号在距离布风板200mm处采集，采用Kistler7261型传感器，量程-5～5kPa，响应频率1Hz，采样频率400Hz，采样时间60s。
实验过程中石英砂的总量保持不变，气相速度的变化范围在0.6～5.0m/s，随着气流速度的增加，在流化床依次观察到固定床、鼓泡流化床、湍流流化床和气力输送床4种流型，其对应的压力脉动信号如图2所示。在每种流化型态下，分别采样17段压力脉动时间序列，长度75 000，经过奈奎斯特频率低通滤波由计算机记录。
3 流化床流型的多尺度熵分析
不同工况下，固定床压力脉动信号多尺度熵的计算结果如图3所示，在尺度1～10，随着尺度的增加样本熵值近似呈线性增长，在第10个尺度后样本熵值基本保持不变;样本熵值在小尺度上对流动工况的变化不太敏感。分析其原因，在固定床流型下，气体主要从固体颗粒间的缝隙中通过流化床，粒子的搅动由气体射流，粒子搅动的运动规则毫无规律或者是混沌的，压力脉动信号反映的是整体行为，脉动幅度相对较小，夹杂着较多的周期性成分。
不同工况下，鼓泡流化床压力脉动信号多尺度熵的计算结果如图4所示，在尺度1～6，随着尺度的增加样本熵值近似呈线性增大趋势，且样本熵值对工况的变化不敏感;尺度8后，每种工况的样本熵值增长缓慢，且样本熵值对工况的变化十分敏感。这是因为在鼓泡床流化型态下，由于在流化床内产生明显的气泡，气泡引起的颗粒运动在床层运动中占主导作用，气泡的运动相对颗粒的搅动有规律得多，而随着气相速度的逐渐增大，气泡团聚并产生的大气泡在乳化相的剧烈扰动下破碎，大气泡量明显减少，床内的总体趋势表现为由鼓泡态的大尺寸、少量气泡的状态向小尺寸、多数量气泡的状态演变，床层慢慢进入湍流状态，压力脉动信号中随机分量迅速增强。
不同工况下，湍流流化床压力脉动信号多尺度熵的计算结果如图5所示，在尺度1～12，随着尺度的增加样本熵值近似呈线性增长;在尺度10后，样本熵值的增长十分缓慢，而且在大尺度上样本熵值对流动工况的变化比较敏感。这是由于在流化床湍流流态下，气泡破裂，床面不存在或很难区分;颗粒浓度随高度连续下降，需要一定的颗粒循环量来维持颗粒总量，床密度不依赖颗粒循环倍率。
不同工况下，气力输送床压力脉动信号多尺度熵的计算结果如图6所示，在尺度1～10，随着尺度的增加样本熵值增长趋势近似呈线性;在尺度10后，样本熵值对流动工况的变化较敏感，而且随着气相速度的增加，在相同尺度上样本熵值增大。这是由于颗粒在循环装置中相互碰撞或颗粒与壁面碰撞，颗粒随着气体的运动十分复杂，气相和颗粒处于混沌状态，压力脉动信号接近于固定床的压力脉动信号，脉动幅度相对最小，其内部还夹杂着较多的周期性成分。
从图3～6可以看出，利用固定床压力脉动信号的多尺度熵特征分析可知，多尺度熵可以在不同尺度上很好地揭示固定床、鼓泡流化床、湍流流化床和气力输送床的动力学特性。样本熵值在小尺度上近似呈线性增长，不同流型样本熵值的增长速率有很大差异，同一流型的样本熵值增长速率差别不大。可以利用这一显著特征作为区分流型的新标准。样本熵值的增长速率在熵值曲线上就是曲线的斜率，可以通过最小二乘法拟合得到。笔者将其定义为多尺度熵率(Rate of MSE)。
不同流动工况下的多尺度熵率分布情况如图7所示，4种流型的多尺度熵率差别显著，其中固定床的多尺度熵率0.100～0.150，鼓泡流化床的多尺度熵基本位于0.045以下，湍流流化床的多尺度熵率在0.050～0.070，气力输送床的多尺度熵基本在0.070～0.100。
4结束语
多尺度熵作为一种非线性分析方法可以在多个尺度上表征信号的复杂性，并且对序列长度有较好的鲁棒性。笔者将多尺度熵分析应用到气-固两相流压力脉动信号分析中，揭示了流化床不同流型内部的动力学特性，并进一步揭示了不同流型之间的动力学差异。此外，在小尺度下4种流化型态的多尺度熵值都是随着尺度的增加近似呈线性增长，在大尺度下熵值的增长十分缓慢。在流化床气-固两相流不同流型多尺度熵表现出的不同变化趋势基础上，提出了一个新的参数———多尺度熵率，能较好地区分4种典型流型，为流化床气-固两相流流型辨识提供了新方法。
摘要：采用多尺度熵算法对气-固流化床内4种典型流型的压力脉动信号进行研究,分析了其内部多尺度动力学特性。结果表明:不同流型由于其动力学特性不同,导致压力脉动信号的多尺度排列熵值存在差异。在小尺度下4种流化型态的多尺度熵值都是随着尺度的增加近似呈线性增长,在大尺度下熵值的增长十分缓慢。此外,利用多尺度熵值变化速率特征识别气-固流化床内的流型,是一种新的流型辨识方法。
关键词：气-固两相流化床,多尺度熵算法,动力学特性
参考文献
[1]Li J H,Mooson K.Particle-Fluid Two-phase Flow-the Energy-Minimization Multi-Scale Method[M].Beijing:Metallurgical Industry Press,1994.
[2]Zhao G B,Chen J Z,Yang Y R.Predictive Model and Deterministic Mechanism in a Bubbling Fluidized Bed[J].AICh E J,2001,47(7):1524~1532.
[3]Fan L T,Neogi D,Yashima M,et al.Stochastic Analysis of a Three-phase Fluidized Bed:Fractal Approach[J].AICh E J,1990,36(10):1529~1535.
[4]Franca F,Acikgoz M,Lahey R T,et al.The Use of Fractal Techniques for Flow Regime Identification[J].Int J Multiphase Flow,1991,17(4):545~552.
[5]赵明阳,黄志尧,王保良,等.气固流化床压力波动信号的高阶统计量特性研究[J].仪器仪表学报,2001,22(z2):11~12.
[6]秦伟刚,王超,张文彪,等.气-固两相流差压信号的EMD和三阶累积量分析[J].仪器仪表学报,2014,35(4):762~767.
[7]王春华,仲兆平,鄂加强.气-固流化床压力脉动信号分形标度分析[J].中国电机工程学报,2010,30(8):45~49.
[8]张少峰,王琦,高川博,等.液固两相外循环流化床压力波动信号的统计及频谱分析[J].过程工程学报,2006,6(6):878~883.
[9]Zhu L,Jin N D,Gao Z K,et al.Multi-scale Cross Entropy Analysis for Inclined Oil-water Two-phase Countercurrent Flow Patterns[J].Chemical Engineering Science,2011,66(23):6099~6108.
多尺度熵算法 篇4
时频分析方法由于能够提供振动信号时域和频域局部信息而在故障诊断领域得到了广泛应用, 是目前故障诊断的重要手段[1,2]。该方法通过对信号进行小波分析[3]或经验模态分解[4,5], 将非平稳信号分解为若干个简单的平稳信号之和, 然后对每个分量进行处理, 提取时频域信息, 进而得到原始信号的完整时频信息。然而, 由于机械运转过程中的摩擦、振动以及负载等因素, 机械系统振动信号往往表现出非线性行为, 采用时频分析的方法, 将信号分解为平稳信号, 难免有一定的局限性。非线性分析的方法可以不经过对原始信号分解, 而能够直接提取隐藏在机械系统振动信号中其他方法无法提取的故障信息[6]。
近年来, 越来越多的非线性分析方法被用于机械故障诊断, 如关联维数、近似熵、样本熵、多尺度熵等。徐玉秀等[7]研究了旋转机械的分形特征及故障诊断, Yan等[8]将近似熵应用于机械系统健康状态监测, Zhang等[9]将多尺度熵应用于滚动轴承的故障诊断。然而, 分形维数的计算依赖数据长度, 且较耗时, 不适合在线监测;近似熵相对一致性较差[9];多尺度熵[10,11]是基于样本熵而定义的, 计算较耗时, 且受时间序列的非平稳性和异常值的影响。
排列熵[12]是一种新的随机性和动力学突变的检测方法, 它具有计算简单, 抗噪能力强, 且得到较稳定的系统特征值所需时间序列短, 适合在线监测等优点, 在肌电信号处理[13]、心率信号处理[14]、气温复杂度[15]等方面都取得了良好的效果。Yan等[6]将其应用于旋转机械振动信号的特征提取, 并将其与近似熵和Lempel-Ziv复杂度进行了对比, 结果表明, 排列熵能够有效地检测和放大振动信号的动态变化, 并且能够表征滚动轴承在不同状态下的工况特征。然而, 与传统的基于单一尺度分析的非线性参数类似, 排列熵只是检测时间序列在单一尺度上的随机性和动力学突变。Aziz等[16]提出了多尺度排列熵 (multiscale permutation entropy, MPE) 的概念, 用于衡量时间序列在不同尺度下的复杂性和随机性, 并通过分析生理信号, 将其与多尺度熵进行了对比, 结果表明, 相对于多尺度熵, MPE更具有鲁棒性。由于机械系统比较复杂, 振动信号不仅在单一尺度上包含有重要信息, 而且在其他尺度上也包含有重要信息, 因此, 对振动信号进行多尺度分析是一种有效的方法。
本文将MPE引入到机械故障诊断领域, 应用于滚动轴承的故障特征的提取。由于正常滚动轴承的振动信号是随机振动信号, 而当滚动轴承发生故障时, 振动信号随机性和动力学行为会发生变化。因此, 本文考虑用MPE来衡量振动信号的随机性变化和动力学突变, 并将MPE值作为特征参数, 提取滚动轴承的故障特征。在此基础上, 结合支持向量机 (support vector mechine, SVM) [17,18]作为模式识别分类, 提出一种基于MPE和SVM的滚动轴承故障诊断方法, 并将其应用于滚动轴承实验数据的分析, 结果表明, 新提出的方法能够有效地诊断滚动轴承的故障类型。
1 排列熵及多尺度排列熵原理和算法
1.1 排列熵算法
排列熵的原理在于不考虑数据具体值, 而是基于相邻数据的对比。下面说明其计算方法。
考虑长度为N的时间序列{x (i) , i=1, 2, …, N}, 对其进行相空间重构, 得到如下的时间序列:
式中, m为嵌入维数;λ为时延。
将X (i) 的m个数据按照升序重新排列, 即
如果存在x (i+ (ji1-1) λ) =x (i+ (ji2-1) λ) , 此时按j值的大小来进行排序, 即当jk1<jk2, 有x (i+ (ji1-1) λ) ≤x (i+ (ji2-1) λ) , 所以, 任意一个数据X (i) 都可以得到一组符号序列:
其中, g=1, 2, …, k, k≤m!, m个不同的符号{j1, j2, …, jm}共有m!种不同的排列, 对应地, 共有m!种不同的符号序列, S (g) 是m!种符号序列中的一种。计算每一种符号序列出现的概率Pg (g=1, 2, …, k) , , 此时, 时间序列{x (i) , i=1, 2, …, N}的排列熵就可以按照Shannon熵的形式定义为
注意到, 当Pg=1/m!时, Hp (m) 达到最大值ln (m!) , 因此, 可以通过ln (m!) 将排列熵Hp (m) 进行标准化处理, 即
显然, Hp的取值范围是0≤Hp≤1。Hp值的大小表示时间序列的复杂和随机程度。Hp越大, 说明时间序列越随机, 反之, 则说明时间序列越规则。Hp值的变化反映和放大了时间序列的局部细微变化。
1.2 排列熵参数的选取
在排列熵的计算中, 有3个参数值需要考虑和设定, 即时间序列长度N、嵌入维数m和时延λ。Bandt等[12]建议嵌入维数m取3~7, 因为如果m等于1或2, 此时重构的序列中包含太少的状态, 算法失去意义和有效性, 不能检测时间序列的动力学突变。但是, 如果m取值过大, 也不合适, 因为相空间的重构将会均匀化时间序列, 此时不仅计算比较耗时, 而且也无法反映序列的细微变化[8,15]。
为研究数据长度N对排列熵值的影响, 以长度分别为128、256、512、1024和2048的高斯白噪声信号为例, 求得对应排列熵值, 分别记为PE1~PE5, 如图1所示, 它们在不同嵌入维数下差值如表1所示。
由图1和表1可以发现, 以嵌入维数m=6为例, 数据长度分别为1024和512时, 熵值相差0.0659, 而数据长度分别为2048和1024时, 则熵值仅相差0.0309, 因此, 此时选数据长度为1024较合适。而对m=5而言, 数据长度分别为1024和256的信号的熵值仅相差0.0502, 此时, 数据长度为256已经可以估计合理的排列熵值。一般, 嵌入维数较小时, 数据长度则要求越小。
时延λ对时间序列的计算影响较小, 以长度为512的高斯白噪声信号为例, 在不同λ下的排列熵值随嵌入维数的变化关系如图2所示, 由图可以看出, 时延对信号熵值的影响较小, 因此, 本文取λ=1。
1.3 多尺度排列熵定义
多尺度排列熵定义为不同尺度下的排列熵, 计算方法如下:
(1) 考虑时间序列{x (i) , i=1, 2, …, N}, 对其进行粗粒化处理, 得到粗粒化序列{yj (τ) }。yj (τ) 的表达式为
其中, [N/τ]表示对N/τ取整;τ为尺度因子, τ=1, 2, …。显然τ=1, 粗粒化序列即为原始序列;τ>1时, 原始序列被粗粒化为长度为[N/τ]的粗粒序列。
(2) 计算每个粗粒序列的排列熵, 并画成尺度因子的函数, 上述过程即称为多尺度排列熵分析。
尺度因子的最大值一般取大于10即可, 但要保证粗粒化序列长度[N/τ]不影响熵值的计算。
为了选取合适的计算MPE的嵌入维数, 仍以高斯白噪声为例, 数据长度为2048, 尺度因子最大值为12, λ=1, 在m分别为4、5、6和7时, 求得它们的MPE, 相对耗时分别为0.1880s、0.6710s、3.8290s和27.6710s, 将MPE画成尺度因子的函数, 如图3所示。
由图3可以看出, 若m取值太小, 则PE值随尺度因子的增大而减小, 但m越大, 计算越耗时, 因此, 本文选取m=6。此外, 由图3可以看出, 高斯白噪声的MPE随着尺度因子的增大而单调递减, 这说明白噪声只在最小尺度上包含有主要信息。
2 基于MFE的滚动轴承故障诊断方法及应用
在上述理论的基础上, 本文提出基于MPE和SVM的滚动轴承故障诊断方法。首先, 从滚动轴承的原始振动信号从提取MPE;其次, 依据MPE提取合适的故障特征向量;第三, 采用SVM进行故障分类, 从而实现滚动轴承故障类别的诊断。
为了说明本文方法的有效性, 本文将该方法应用于实验数据分析。本文实验数据采用Case Western Reserve University (CWRU) 轴承数据中心提供的滚动轴承试验数据[9]。测试轴承为6205-2RS JEM SKF深沟球轴承, 电机功率约为2206.4963W, 转速为1730r/min, 采用电火花加工技术在轴承上布置单点故障, 故障直径为0.5334mm, 深度为0.2794mm。在此情况下采集到正常 (normal, 简称NORM) 、内圈故障 (inner race fault, IRF) 、外圈故障 (outer race fault, ORF) 和滚动体故障 (rolling element fault, REF) 4种状态的振动信号, 各30组数据, 数据长度为2048, 采样频率为12kHz, 4种状态轴承的振动信号时域波形如图4所示。
从图4不易发现正常和故障轴承振动信号的明显区别, 尤其是正常和滚动体故障, 以及内圈故障和外圈故障。因此, 本文首先对振动信号进行MPE分析, 取嵌入维数m=6, 时延λ=1, 最大尺度因子为12, 4种状态滚动轴承的多尺度排列熵画成尺度因子的函数, 如图5所示。
在尺度因子等于1时, 即为原始振动信号的排列熵, 由于熵值比较接近, 无法明显地区别三种故障和正常轴承的类型, 因此有必要对振动信号进行多尺度分析。为此, 以多尺度排列熵值为特征参数, 同时建立基于SVM的分类器, 进行训练和测试。如果采用全部的12个特征值进行训练, 会造成信息的冗余, 且训练比较耗时, 也需要较多的训练样本, 且由图5也可以看出, 前几个尺度的熵值表征了振动信号的主要信息, 因此, 采用前4个尺度的排列熵值作为特征向量, 即T= (PE1, PE2, PE3, PE4) 。因此, 本文的方法如下:
首先, 提取特征参数, 即对振动信号进行MPE分析, 提取特征参数T。正常、滚动体故障、内圈故障和外圈故障4种状态, 每种状态取30个样本, 故每种状态可得到30个表征故障特征的特征向量, 共得到120个特征向量。
其次, 训练分类器。由于有3种故障状态和正常状态, 因此, 需建立3个SVM, 其中SVM1为正常对3种故障分类器, SVM2为内圈故障对滚动体和外圈故障分类器, SVM3为滚动体故障和外圈故障分类器。每种状态随机抽取10个样本进行训练, 并将每组30个样本用来测试。经过训练, SVM1和SVM3采用径向基核函数, SVM2采用多项式核函数。基于SVM的多故障分类器如图6所示。
最后, 测试分类器。对已训练的SVM1、SVM2和SVM3, 用全部样本进行测试, 详细测试样本输出结果如表2所示。
由表2可以看出, 本文提出的方法有很好的效果, 在全部样本用来测试中, 只有一组外圈故障的样本被错分为滚动体故障, 其他都得到了正确的分类, 正确识别率为99.17%。为了比较, 下面建立以BP神经网络为基础的多分类器[9,19,20], BP分类器除输入层外, 第一层隐含层有8个节点, 第二层输出层有4个节点。为表述方便, 标记正常为1类, 内圈故障为2类, 滚动体为3类, 外圈故障为4类。BP分类器的训练和测试样本与支持向量机相同, 其分类结果如图7所示。
图7中, 训练样本作测试时, 全部分类正确, 而测试样本作测试时有3组分类错误, 准确率为97.5%。这说明支持向量机的分类效果要优于BP神经网络。而且, 在训练时间上, 支持向量机也比BP神经网络短得多。
为了说明进行多尺度分析的必要性, 下面选取尺度因子等于1时 (即原始信号) 的排列熵值作为特征参数, 分别通过SVM和BP神经网络分类器进行训练和测试。其中, SVM分类器中, 有6组样本分类错误, 如表3所示;而BP神经网络分类器也有6组样本分类错误, 如图8所示。
从上文可以看出, 原始信号单一尺度的排列熵作为特征参数, 分类效果不理想, 这说明了单一尺度上原始振动信号的排列熵并不能反映故障的本质, 而多个尺度上的排列熵值则能够更好地实现分类。另外易发现, 特征值的选取对分类结果的影响尤为关键。本文选取特征值为前4个尺度上的特征值, 主要基于以下原因考虑:如果特征值过少, 不能完全反映故障的特征信息, 而特征值过多会造成信息冗余, 且需要增加训练样本和训练时间, 因此, 本文选取了前4个尺度上的特征值。文献[9]中选择多尺度熵值的统计量时, 将最大值、最小值、代数平均、几何平均和标准差作为特征向量, 但统计量方法忽略了特征值之间的内在关系, 因此, 采用前4个尺度因子的排列熵值作为特征参数。
3 结束语
机械系统发生故障时, 振动信号会在不同尺度上表现出不同程度的随机性和动力学突变, 基于此, 本文提出了一种新的基于多尺度排列熵和支持向量机的滚动轴承故障诊断方法, 并将支持向量机与BP神经网络分类效果进行了对比。结果表明, 支持向量机在训练时间和准确率方面, 都优于BP神经网络。此外, 本文还将特征向量包含多个尺度上熵值与特征向量仅包含原始信号单一尺度上排列熵值进行了对比, 结果表明, 原始信号单一尺度上排列熵值的不能全部反映故障的本质, 而多尺度的排列熵值则有很好的诊断效果。本文提出的方法为故障诊断提供了一种新的思路和手段。
多尺度熵算法 篇5
关键词：小波函数,经验模态分解,多尺度熵,故障诊断
0引言
在常见的机械设备中,滚动轴承是整个设备的重要部件,其在运作中的好坏直接影响到设备的生产效率, 因此,准确地识别出运作过程中的状态至关重要。滚动轴承出现故障是因为设备内部各部件周而复始的运作外加载荷转速等其他原因,致使滚动轴承出现一些内外圈以及滚动体的故障,而且在故障出现的同时,伴随的振动信号也表现出剧烈的非线性和非平稳性的特征。传统对非平稳和非线性信号的处理方法,比如小波分析、分形维数、支持向量机等方法已大量的被用在滚动轴承故障的诊断领域,这些方法的渗透极大丰富了非平稳信号
特征识别的研究体系[1～3]。
HHT是由N.E.Huang提出的一种更具有适应性的处理复杂非线性和非平稳信号的分析方法[4],此方法通过经验模态分解(EMD),使信号本身分解出一组互异的基底,且分解结果具有高度的自适应性。在分解过程中, 传统的均值包络建立是采用三次样条插值伴随边界效应、过冲和欠冲等问题。S. R. Qin提出采用分段平均幂函数的方法求取信号的上下包络线抑制这种过冲现象[5]。 E.Delechelle提出一种利用四阶抛物型偏微分方程解决均值包络问题的方法[6]。
Richman提出了改进的复杂尺度测量新方法即样本熵,样本熵具有得到稳定估计值所需的干扰能力强、抗噪声和参数取值范围内一致性好等显著特点[7]。由于样本熵的尺度单一导致只是衡量振动信号本身尺度下的复杂性。Costa在样本熵的研究基础上,提出了一种新的时间序列复杂度的衡量方法—多尺度熵[8～10],这一次研究使熵的内容得到了全面的扩充。
上述两大类方法,国内研究人员都分别在故障诊断识别方面应用过,并取得了一定的成效,但对识别的效果及精确度仍有待改善。针对HHT,本文提出基于改进的经验模态分解,即利用小波函数作为插值函数求原始振动信号的包络,经过多次实验模拟,发现可以良好解决欠冲、过冲现象,基于这种研究下,再与多尺度熵结合起来应用到滚动轴承的振动信号当中,给出了一种全新的滚动轴承故障诊断识别方法。
1基于小波函数的经验模态分解
EMD分解是HHT方法的核心,EMD方法将任意复合信号分解成有限固有模态函数(IMF)之和。EMD分解对一个x(t)信号进行分解采用步骤如下:
1)首先获得数据的最大值和最小值所有点,再利用小波函数作为基函数进行插值,并把最大值最小值同时连接,形成上下包络线,原始数据处于两线之间。
2)计算出上下包络均值m1,用原始数据减去均值:
判断h1是否符文IMF两个条件(信号零点的个数与信号极值点的个数相等或相差1和整个时间序列关于时间轴局部对称),如果符合,则得到第一个分量。
3)若不符合,把h1看做成原始数据,并把以上2步骤再次运行,得到新的均值m11,判断h11=h1-m1是否符合条件,不符合需再做k次处理,h1(k-1)-m1k=h1k,直至h1k符合为止,为C1=h1k,则C1是信号x(t)的第一个分量。
4)将C1从x(t)中分离出来,得到:
将r1看成原始数据,重复上述过程获得x(t)的第二个分量C2,依次做n次处理,即可获得n个分量,即:
5)最终得到一个单调的函数即循环结束。
式中,rn为残余项。
基于小波函数EMD的优点:
1)一层层将数据基于其本质特征筛分的过程,作用类似一组滤波器。
2)模态波形对称,时间尺度Ci从小到大依次分离。
3)有效的抑制边界效应、过冲和欠冲等问题。
2多尺度熵
多尺度熵的计算是在样本熵的基础上,将原始数据粗粒化得到的在各个尺度上的样本熵值组成的一组数列,即时间序列在不同尺度下的样本熵。如果一个序列和另一个序列在同种尺度下,前者的熵值比后者高,这说明前者的时间序列的复杂性要高于后者,这就体现出了前者和后者两者之间具有强烈的差异性。多尺度熵的计算过程如下:
1)对给定原始信号数据X={x1,x2,x3,L,xN},长度是N,即N为序列长度,给出嵌入维数m,并定义,。
2)定义X(i)与X(j),它们之间距离D(i,j)为两个样本对应元素差值绝对值的最大值。
其中,i,j=1,2,3L,n-m,k=1,2,3Lm-1.。
3)计算所有的X(i)与其他在子样本X(j)的距离, 给出相似容限r,统计并计算出D(i,j)<r的个数,标记成count D(i,j)最后算出平均值:
4)计算Bim(r)的均值定义为Bm(r)。
5)重构数据,令嵌入维数等于m+1,重复过程1-4步骤计算出Bm+1(r)。
6)单一尺度上时间序列的样本熵计算:
一般N均为有限数,即:
7)对于长度为N原始数据X根据尺度因子 τ 的取值不同把原始数据分割成长度为τ的数据单元,再将这N/ τ组单元分别组内求均值,利用所求的均值点形成一个新的时间序列,这个过程也称将原始数据按τ的粗粒化过程。而单元组内求均值的相应点为:
8)对于τ 的不同取值,也就形成了对原始X的多尺度化,对每个尺度化后新的时间序列依次求样本熵,即得到了多尺度熵。
多尺度熵的计算跟嵌入维数m、相似容限r以及尺度因子τ 都有关联。其中,m的取值一般为2,相似容限rr==00..11~~00..2255**δ (其中 δ为原始数据的标准差),对于尺度因子一般取值不超过20。
多尺度熵的优点:
1)特征明确,其熵值的大小直接反映出时间序列产生新模式概率的大小。
2)具有较强的抗干扰性。
3)对原始信号进行多尺度分析,使得分析更具完备性和系统性。
3基于经验模态分解和多尺度熵的滚动轴承故障诊断
数据是从美国华盛顿凯斯西储大学实验室的滚动轴承数据中截断获取,被测轴承是SKF6205,轴承的损坏是由人为的加工制作完成,之后通过振动加速度传感器获得各种工作状态的振动数据。
本文数据是轴承转速在1797r/min的情况下采集,采样频率为12k HZ,主要分析滚动轴承的四种工作状态, 分别是:滚动轴承正常、滚动轴承内圈故障、滚动轴承外圈故障、滚动体故障,每种状态下的数据共6个小样本, 每个小样本数据长度为6000,且每个小样本采集用时30s。所有状态共计24个小样本,耗时12min。在进行多尺度熵计算过程中,m取值为2,相似容限r=0.15*δ ,尺度因子τ==1155。第一次进行样本采集得到的四种工作状态原始数据样本对比以及利用改进的EMD对四种工作状态分解得到的IMF分别如图1~图4所示。
从图1~图4中可以直观地判断出四种工作状态所描绘的振动信号有差异性,正常的工作状态与滚动体出现故障对比内圈出现故障和外圈出现故障的折射的波动性差异明显,但正常工作状态和滚动体出现故障仍然无法准确进行区别,当把振幅作为识别特征,内圈故障和外圈故障也出现了相互无法明确分离的情况。
将改进的经验模态分解应用到原始振动信号数据得到的IMF1和IMF2,利用多尺度熵进行计算,结果如图5和图6所示。
从图5中可以看出,第一,从四种状态的多尺度熵值聚合程度来看,滚动体故障和外圈故障相比原始数据多尺度熵值相对聚集。第二,图5中的四种工作状态的IMF1多尺度熵,外圈故障可以准确的与其他三种情况分离开来。其他工作状态仍然出现了熵值近似情况。 出现这些情况,总的来说原始信号经过改进的EMD分解得到的IMF1虽然振动信号得到了平稳化,而且为原始信号的高频成份,但对于轴承振动而言,这并不能代表原始信号作为特征进行故障识别。而从图6看出, IMF2的多尺度熵值更加不理想,原始信号经过经验模态分解得到的各自IMF是从高频到低频,IMF1的频率高于IMF2,得到的IMF2是原始信号减去IMF1再进行改进的EMD分解而得,对原始信号而言,特征体现次于IMF1。
对得到的固有模态函数IMF1和IMF2取和进行多尺度熵计算,结果如图7所示。
分析图7可以明显发现四种状态基于多尺度熵值可以很好地分离开来,这说明了基于改进的EMD对信号进行分解后,得到固有模态函数IMF1与IMF2两者的组合具有代表原始数据特征的特性,也就意味着要从数据的整体性把握对振动信号的故障分析,加上后续利用多尺度熵的方法,对重构后的新数据进行计算后,得出结论,尤其是在尺度因子选为12即τ==1212时,可以很好地进行滚动轴承振动信号的故障分离与识别。
基于上述研究,利用MATLAB的LIBSVM包对滚动轴承故障状态进行分类,进而实现诊断识别。首先以上述24组四种工作状态的IMF1与IMF2加和的多尺度熵作为训练集,并做出标记,令正常、内圈故障、滚动体故障以及外圈故障分别为1、2、3、4。然后在实验室的滚动轴承数据中,在不含以上数据的数据段中,把随机采取四种状态的6000点作为测试集,测试结果显示,利用改进的经验模态分解得到的IMF1与IMF2加和,其再进行多尺度熵形成的数组,利用LIBSVM之后可以精确的将新数据的四种情况进行分类,分别是正常为1、内圈为2、滚动体为3和外圈为4,从而达到了诊断目的。
4结论
多尺度熵算法 篇6
滚动轴承是旋转机械中应用最广泛的机械部件,也是最容易出现故障的机械部件,其运行状态正常与否将直接影响到整个机组的性能,因此对滚动轴承的故障诊断方法和监测技术的研究具有重要意义[1]。滚动轴承振动信号包含丰富的故障特征信息,当滚动轴承发生故障时,振动信号呈现出非线性非平稳特性。 局部均值分解(local mean decomposition,LMD)是一种自适应时频分析方法,可以将复杂的非平稳信号分解为若干个乘积函数(product function,PF)分量之和。每一个PF分量都是由一个纯调频信号和一个包络信号相乘得到的,将所有PF分量的瞬时幅值和瞬时频率组合即可得到原始信号完整的时频分布[2]。目前,LMD方法在机械故障诊断领域得到广泛应用[3,4,5,6,7]。
时域统计量(峭度、偏度、能量等)能反映振动信号的分布特性,因此根据时域统计量指标可以实现故障的判断[8,9]。程军圣等[10]提出了一种基于LMD时域统计量和神经网络的滚动轴承故障诊断方法,将时域统计量作为神经网络的输入特征量来进行故障诊断。时域统计量作为故障特征参数时,容易受到噪声等因素的影响,从而影响故障诊断的准确率。
熵是表征系统规则度及复杂度的物理量,近年来被应用于机械故障诊断领域,多尺度熵(multi-scale entropy,MSE)能衡量时间序列的复杂性,是时间序列在维数变化时产生新模式概率大小的量度,具有很强的抗噪声和抗干扰能力,是量化非线性时间序列复杂度的良好工具[11,12,13,14,15]。神经网络作为一种智能化的数据处理方法,具有很强的处理非线性关系数据的能力,其中,概率神经网络(probabilistic neural network,PNN)能够以任意精度逼近任何连续非线性函数,同时具有自组织、自学习和并行处理能力[16,17]。基于此,本文将局部均值分解、多尺度熵和概率神经网络相结合,提出了一种基于LMD多尺度熵和概率神经网络的滚动轴承故障诊断方法,并将该方法与基于LMD时域统计量和神经网络的滚动轴承故障诊断方法进行了对比。实验结果证明,该方法是有效可行的。
1 故障诊断方法
LMD将振动信号分解为一系列的PF分量PF1,PF2,…,PFn和一个残余量R,每个PF分量由一个包络信号(PF分量的瞬时幅值)和一个纯调频信号相乘得到,PF分量的瞬时频率则由纯调频信号直接求出。将所有PF分量的瞬时幅值和瞬时频率组合便可以得到原始信号完整的时频分布,其分解过程详见文献[2]。
多尺度熵基于样本熵,其计算方法如下:
(1)设一离散原始时间序列为{x1,x2,…,xN},对原始时间序列进行粗粒化变换,得到新的时间序列:
其中,N为离散时间序列长度;τ为尺度因子。原始序列被分割成τ段且每段长为N/τ 的粗粒序列。τ=1时,新的时间序列就是原始序列。
(2)给定模式维数m和相似容限r(r>0),构造时间序列的m维向量:
(3)计算向量Xm(i)与Xm(j)之间的距离:
(4)对每个i,计算Xm(i)与Xm(j)的距离,统计其距离小于r的数目,记为L,将此数目与距离总数N -m+1的比值记作Cmi(r),即
(5)将Cmi(r)的平均值记作
(6)把维数加1,变成m+1,重复步骤(1)~步骤(5),计算Φm+1(r)。
当N为有限值时,按上述步骤得出的是序列长度为N时的样本熵估计值:
重复上述过程,得到不同尺度下的样本熵值,即为多尺度熵。
概率神经网络是一种前馈型神经网络,包括输入层、模式层、求和层和竞争层,其结构如图1所示。输入层将输入样本传给模式层的各节点。模式层与输入层之间通过连接权wij相连,进行加权求和,并通过非线性算子g(zj):
运算后,传递给求和层。该层各个模式单元的输出为
其中,X为输入样本向量,Wi为输入层到模式层的连接权值;δ为平滑系数,它对分类起着关键性的作用。求和层则简单地将由对应样本中同一类的模式层传来的输出(属于某类的概率)进行累加,即
式中,m为训练样本向量个数。
竞争层(输出层)接收从求和层输出的各类概率密度函数,概率密度函数最大的那个神经元输出为1,其所对应的类即为样本模式识别结果,其余神经元输出为0。
2 基于LMD和神经网络的故障诊断
基于LMD多尺度熵和概率神经网络的故障诊断流程如图2所示,实验中,分别选择故障直径为0.1778mm、0.3556mm、0.5334mm的滚动体故障信号、内圈故障信号和外圈故障信号以及1个正常信号,并定义直径为0.1778mm的轴承故障为轻度损伤,直径为0.3556mm的轴承故障为中度损伤、直径为0.5334mm的轴承故障为重度损伤。因此,滚动轴承的工作状态共有10种,将这10种工作状态用序号1~10表示,如表1所示。
基于LMD多尺度熵和概率神经网络的滚动轴承故障诊断方法具体步骤如下:
(1)对滚动轴承10种状态分别进行10次采样,将得到的100个振动信号作为样本。
(2)对每一种状态下的原始振动信号进行LMD分解,得到各个PF分量。
(3)选择前n个PF分量作为研究对象,计算前n个PF分量的多尺度熵。
(4)构造特征向量
(5)将构造的特征向量作为概率神经网络的输入向量,对网络进行训练。
(6)将待测样本输入到训练后的概率神经网络,以神经网络的输出确定轴承的工作状态。
3 实验研究
实验数据为美国西储大学电气工程实验室的滚动轴承实验数据。采样频率为12kHz,采样点数为5000。对不同状态的振动信号进行LMD分解,其中,故障直径为0.1778mm(轻度损伤)内圈故障信号的LMD分解结果如图3所示。
由于主要故障信息集中分布在前几个PF分量,为了避免多尺度熵特征向量的样本数据冗长,本文选取前3个PF分量进行实验分析。分别对滚动轴承的10种状态进行采样,每种状态选取10组数据作为样本进行LMD分解,将分解后前3个PF分量的多尺度熵(尺度因子分别选取5、10、15)作为特征向量,输入到概率神经网络进行训练。内圈轻度损伤和滚动体中度损伤故障的多尺度熵如表2、表3所示(注:PF1(5)中的5为尺度因子)。
从表2、表3可以看出,滚动轴承在相同工作状态下,各PF分量同一尺度因子的熵值比较接近。不同工作状态下,各PF分量在同一尺度因子下的熵值相差较大。
作为对比,选取前3个PF分量,利用时域统计量峭度Q、偏度P和能量比N作为神经网络的输入特征量,其中,滚动体中度损伤故障的时域统计量如表4所示。从表4可以看出,同一分量下的同种时域统计量的数值存在较大波动,不利于滚动轴承的故障诊断。因此,本文利用多尺度熵值作为故障诊断的依据。将多尺度熵作为神经网络的输入特征向量进行模式识别。
注:PF1(Q)中,Q为时域统计量峭度。
再次对滚动轴承的10种状态进行采样,每种状态采样10次,共得到100组数据(作为待测样本),将待测样本进行LMD分解后,计算出各分量的多尺度熵和时域统计量,分别将各分量的多尺度熵和时域统计量作为输入向量输入到概率神经网络进行滚动轴承工作状态的识别。表5所示为10组待测样本的诊断结果。
注:表中轴承状态“滚-轻”代表滚动体轻度损伤,其他类似。
从表4中可以看出,测试样本中,基于多尺度熵方法的诊断结果与真实工作状态一致,基于时域统计量方法出现2处错误诊断。进一步将基于多尺度熵和基于时域统计量方法的识别率进行对比。对滚动轴承10种状态各进行40次采样,每种状态选取10组(共100组)数据作为训练样本,将每种状态剩下的30组(共300组)数据作为待检测样本进行测试,各诊断方法的识别率如表6所示。
从表6中可以看出,基于LMD分解的时域统计量概率神经网络诊断方法的识别率为89.3%,基于LMD分解的多尺度熵概率神经网络诊断方法的识别率为98.7%。与将时域统计量作为模式识别的特征向量相比,将多尺度熵作为模式识别的特征向量具有较高的诊断识别率。实验结果表明,LMD多尺度熵和PNN结合的滚动轴承故障诊断方法是可行的。因此本文研究的将LMD分解、多尺度熵、概率神经网络相结合的方法能够有效地进行滚动轴承故障诊断。
4 结语
本文提出了一种基于LMD多尺度熵和概率神经网络的故障诊断方法。通过局域均值分解对滚动轴承故障振动信号进行分解,得到了若干PF分量。将多尺度熵作为表征故障特征的参量,求取前三个PF分量的多尺度熵并将其作为概率神经网络的输入向量,进行滚动轴承工作状态的模式识别。实验结果表明:该方法能够有效识别出滚动轴承中不同位置、不同损伤程度的故障,提高了诊断精度和识别的能力。
摘要：研究了一种基于LMD多尺度熵和概率神经网络的滚动轴承故障诊断方法。该方法将故障信号自适应地分解为若干乘积函数分量,然后将各分量的多尺度熵作为故障特征向量输入概率神经网络进行模式识别,实现了对损伤位置和损伤程度的诊断。将该方法与基于LMD时域统计量和神经网络的滚动轴承故障诊断方法进行了对比。实验结果表明,基于LMD多尺度熵和概率神经网络的方法能对滚动轴承故障进行有效的识别与诊断。
多尺度熵算法 篇7
1 Contourlet 变换
Contourlet变换先进行多分辨率分解然后在进行多方向分解,首先用拉普拉斯金字塔(LP,Laplacian Pyramid)变换进行多尺度分解,并捕捉图像信号的奇异点,接着使用方向性滤波器组(DFB,Directional Filter Bank)进行多方向分解,将同个方向上的奇异点合成,捕捉二维信号的高频成分。图1给出了Contourlet变换的流程图。LP变换能够有效的避免由方向性滤波器所造成的低频分量的“泄漏”。
Contourlet变换的核心LP和DFB,因此Contourlet变换也称“塔式方向滤波器组”。而Contourlet变换的最终是用类似于轮廓段 (contour segment)的基结构来逼近表示二维图像信息。图2表示一个Contourlet变换可能分解的频率。由于拉普拉斯金字塔和方向性滤波器组具备的完全重构特性,因此其组合PDFB也具备能完全重构。但是由于拉普拉斯金字塔本身具有冗余性,因此Con-tourlet变换具有冗余性。
图像经Contourlet分解之后绝大多数系数幅值接近零,系数表现相当稀疏,而在原图像的边缘细节以及纹理特征部分稀疏幅值比较大,且在他们在各尺度上具有有一定的相关性。图3给出了lena图的Contourlet分解的例子,图像经3级LP分解,以及8个方向的分解,从图我们可以看到最大多系数都接近零,而在图像的轮廓和边缘部分系数的幅值比较大。
2 基于层结构的多阈值Contourlet降噪算法
阈值降噪是图像降噪处理中最常用的方法,其关键是选取一个适当的阈值把图像噪声系数和信号系数区分开来,然后对噪声系数进行置零或者修改其值,随后通过反变换重构图像,达到图像降噪目的。
文献[8]中提出了基于层结构的Contourlet多阈值图像降噪算法。文献中提出“算法使用Contourlet变换来代替小波变换,Con-tourlet的基函数能有效地对分段光滑的线段进行表示,而对奇异点则影响不大”,在降噪处理的过程中,即使将噪声系数误认为图像细节系数给予保留,或者相反,将图像细节系数误认噪声给予删除,也不会使降噪重构的图像出现显著数值较大的孤立点,从而造成图像边缘细节或纹理特征稀疏的丢失。其方法如下:
1)首先,对带噪声的图像进行CT变换,确定对图像进行分解的级数J和各个级数J上所要分解的方向数s。
2)根据下面式子计算出各个尺度和各个方向上噪声方差估计值,
wgi,j(s) 为带噪声图像高频系数。
3)对图像CT变换后的系数方差进行估计
4) 阈值调整:对低分解层的高频系数进行调整。
5) Contourlet反变换:将处理后的Contourlet变换系数进行反变换,可以重构经过降噪之后的图像。
3 基于Contourlet的相关尺度降噪算法
文献[10]提出了一种基于Contourlet变换的相关降噪新算法,其思想主要是将图像的Contourlet系数跟含噪声图像的Contourlet变换系数进行比较,随着分辨率越高,图像的细节部分对应的Contourlet变换系数越稳定,而噪声则衰减越大,即代表图像纹理特征或边缘细节的Contourlet系数相关比较强,而代表噪声的Contourlet变换系数相关比较弱或不相关。采用了利用小波系数相关量来确定图像边缘信息,并且与阈值函数相关结合。其过程如下:
4 基于Contourlet相关尺度多阈值降噪算法
采用硬阈值函数进行降噪处理后的信号虽然可以比较好地保留图像边缘纹理等细节信息,都由于这种理想滤波器会造成图像出现振铃或者伪吉布斯效应等视觉失真,而基于相关尺度的降噪方法,有时候易将过多噪声系数误判为有用信号,为了取得比较好的降噪效果,该文对以上两种方法进行综合,采用了相关尺度多阈值降噪方法,将阈值函数与相关尺度结合起来。在基于层机构的阈值基础上,对于判断为噪声的Contourlet系数进行相关量的计算,具体算法如下:
1)首先对带有噪声的图像进行CT变换,一般情况下确定对图像进行分解的级数J为3级和各个级数J上所要分解的方向数s为8。
2)运用文献[8]的方法,计算3个尺度以及8个方向上噪声方差估计值,并对对图像的Contourlet系数的方差进行估计最后得到阈值函数:
5 结果及分析
本文采用标准的大小为512*512像素测试图象Lena, Peppers为实验图象,进行仿真实验,比较了多阈值降噪、相关降噪和本文提出的多阈值相关降噪效果的比较,分解级数为4,4,8;表1为各个降噪方法后的信噪比(SNR),图4是lena图各种方法降噪后的结果。
从表1可以看出虽然尺度相关多阈值降噪的SNR不是最高,但在局部图上可以看出此方法在边缘细节的保持上效果更好。
6 结论