基于小波包和样本熵的水泵机组振动特征提取

宋 媛，彭利鸿，赖冠文，张嘉勋，肖志怀，宋丽波

(1.武汉大学动力与机械学院，武汉 430072；2.广东省水科院，广州 510000；3.黄河水利水电开发总公司，河南 济源 454681)

随着我国经济水平快速稳定发展，南水北调工程的开工，一大批大、中型泵站将投入运行，水泵一旦发生故障将会造成巨大的经济损失和严重的社会影响。泵站运行可靠性受到越来越多的重视，因此，有必要对水泵故障进行及时诊断,以提高水泵机组运行的可靠性[1]。

水泵机组振动信号中承载着丰富的机械设备运行状态信息，80%的机组故障在振动信号中有所反映[2]。但泵站的工作环境恶劣，由于采集的振动信号受到影响，其中的有用信息和噪声掺杂在一起，有用信息被抑制。如何从海量状态信息中提取出有用的机组故障特征是机组故障诊断的一大难点。因此，从振动信号中提取故障特征在机组故障诊断领域获得了广泛应用，如：彭文季[3]等利用快速傅里叶变换把振动信号从时域转换到频域，然后提取其幅值特征进行机组故障诊断；赵道利[4]等对机组振动信号经小波分解后的单尺度信号进行傅里叶变换，提取信号的故障特征。但是，傅里叶变换实质上是信号的整体变换，其应用是建立在信号平稳性假设上的，而振动信号一般是非平稳信号，其波动性比较大，信号中的干扰比较多，其中的有用信息在信号中的成分一般是低频部分或平稳部分，傅里叶变换无法实现信号在时间域上的特征，这使得利用傅里叶变换处理非平稳信号受到限制。

为了克服傅里叶变换的使用缺陷，近年来，一种基于时频分析的方法——小波分析开始在机械设备的状态监测与故障诊断领域有了较快的发展和进步。小波分析相对于傅里叶分析的优越之处在于它对复杂、非平稳信号的处理时，具有局部化分析能力，可以根据需要对分析对象进行任意细节的分析，是对传统傅里叶变换的巨大改进。但小波分析同样存在缺陷，利用小波变换处理信号时，只是进一步分解信号的低频部分，而不再分解信号的高频部分，它不能很好地分解和表示包含大量细节信息的中、高频信号，如机械非平稳振动信号、地震信号等[5]。因此在小波分析的基础上，后来的学者又提出了小波包分析。小波包分析是小波分析的进一步发展，克服了小波分析在分析信号高频部分时的不足，可以更加精细分解信号的高频部分，且自动匹配信号特征，因此小波包分析对非平稳信号局部的细节分析精度更好，能更加精细的展现信号的特征。样本熵[6]对数据长度的依赖性很小，同时具有很好的抗干扰能力。因此将小波包变换和样本熵相结合的分析方法引入到水泵机组的振动信号故障特征提取中，为信号故障诊断提供一个更为准确可靠的特征提取算法。

1 小波包分析

1.1 小波包理论

小波包变换是基于小波变换的进一步发展，能够提供比小波变换更高的分辨率。设φ(t)正交尺度函数，ψ(t)是相应的正交小波，{Vk}是由φ(t)生成的多分辨分析空间序列，Wk-1是Vk中关于Vk-1正交补空间。这时L2(R)能写成正交和分解，即：

L2(R)=…⊕W-1⊕W0⊕W1⊕…

(1)

而{2k/2ψ(2kt-l):l∈Z}是Wk的一个规范正交基。

为了对分解后的各个频带进一步进行细分，即把时频进一步局部化，设：

ψ0(t)=φ(t)，ψ1(t)=ψ(t)

则φ和ψ的两尺度关系可写成：

(2)

式中：g(k)=(-1)kh(1-k)。由上式确定的函数ψn(t)，n=2l或2l+1，l=0，1，…，称为关于尺度函数φ(t)的小波包。

1.2 小波包算法

设gnj(t)∈ψnj，则gnj可表示为：

(3)

小波包分解算法：由{dj+1,nl}求{dj,2nl}与{dj,2n+1l}。

(4)

小波包重构算法：由{dj+1,nl}以及{dj,2nl}求{dj,2n+1l}。

(5)

式中：dj,2nl和dj,2n+1l为下一层小波包分解结果；dj+1,nl为上层分解的结果；j为尺度系数；l为位置系数；n为频率；k为变量；h和g分别为正交共轭低通和高通滤波器。

小波包对高频小波系数做了进一步的分解，将随分解尺度j增大而变宽的频谱窗口进一步分割细化，克服了正交小波基随着分解尺度j的增大，空间分辨率增大而频率分辨率降低的缺陷，从而提高了频率分辨率，更能适应信号的变化，对信号局部特征的刻画更加清楚，时频特性更好。

2 样本熵

样本熵(SampEn，Sample Entropy)是由Richman和Moornan[6]提出的一种新的时间序列复杂性测度方法，可以用SampEn(m，r，N)来表示。其中，N为长度，r为相似容限，维数为m及m+1。样本熵旨在降低近似熵的误差，与已知随机部分有更加紧密的一致性。样本熵是一种与近似熵类似但精度更高的方法[6]。样本熵算法[7]如下：

Step1：设原始数据为u(1)，u(2)，… ，u(N)共N个点。

Step2：按序号连续顺序组成一组m维矢量：从Xm(1),到Xm(N-m),其中:Xm(i)={u(i),u(i+1),…,u(i+m-1) },1≤i≤N-m+1。这些矢量代表着从第i个点开始连续的m个u的值。

Step3：定义矢量Xm(i)和Xm(j)间的距离d[Xm(i),Xm(j)]为两者对应元素中差值最大的绝对值, 即:

d[Xm(i),Xm(j)]=max(|u(i+k)-u(j+k)|)

(6)

其中k=0～m-1;i,j=1～N-m,j≠i。

Step4：给定阈值r，对每个i≤N-m的值，统计d[Xm(i),Xm(j)]小于r的数目(称为模板匹配数，记为Nmi)及此数目与距离总数N-m-1 的比值，记作：

(7)

对所有i对应的Bmi(r)求平均值，记为：

(8)

Step5：增加维数为m+1，按序号连续顺序组成一组m+1维矢量：

从Xm+1(1)到Xm+1(N-m),其中:Xm+1(i)=[u(i),u(i+1),…,u(i+m)]。这些矢量代表着从第i个点开始连续的m+1个u的值。

Step6：定义矢量Xm+1(i)和Xm+1(j)间的距离d[Xm+1(i),Xm+1(j)]为两者对应元素中差值最大的绝对值，即：

d[Xm+1(i),Xm+1(j)]=max(|u(i+k)-u(j+k)|)

(9)

其中k=0～m;i,j=1～N-m,j≠i。

Step7：给定阈值r，对每个i≤N-m的值，统计d[Xm+1(i),Xm+1(j)]小于r的数目(称为模板匹配数，记为Nm+1i及此数目与距离总数N-m-1的比值，记作：

(10)

对所有i对应的Bm+1i(r)求平均值，记为：

(11)

理论上，此序列的样本熵为：

(12)

样本熵的计算过程很明显地反映出m，r取值对于样本熵值的巨大影响。根据Pincus[8]的研究结果，当m=1或2，r=0.1Std～0.25Std(Std是原始数据的标准差)时计算得到的样本熵值具有较为合理的统计特性。

样本熵具有以下优点：

(1)样本熵不包含自身数据段的比较，因此它是条件概率的负平均自然对数的精确值，因此样本熵的计算不依赖数据长度；

(2)样本熵具有更好的一致性。即如一时间序列比另一时间序列有较高的值的话，那对于其他m和r值，也具有较高的值；

(3)样本熵对于丢失数据不敏感。即使数据丢失多达1/3，对计算值影响依然很小。

综合上述优点可知，样本熵算法可以应用于随机过程的研究。

3 基于小波包变换和样本熵的特征提取

时频域特征量可以同时反映分析对象的时域和频域特性。和小波分析相比，小波包分解可以对信号的低频部分进行处理，同时对信号的高频部分进行处理，具有更好的时频分辨率。当机械设备处于异常运行状态或故障状态时，对应信号的不规则性和不平稳性可以反映设备的状态，对故障的发生和发展也会有所体现。首先使用小波包对信号进行分解，得到n个频带系数。对小波包频带系数进行重构，计算各频带的样本熵值。以各频带信号的样本熵值为元素，构造特征向量V=[SamEn1,SamEn2,SamEn3,…,SamEnn]。

3.1 试验数据来源及处理

图1为本文采用的信号采集系统。本采集系统主要由转子振动试验台、转子台控制器、前置器以及计算机软件系统构成。转子振动试验台配备有一台直流电机，配有DH5600转轴控制器，通过改变转子转速、刚度、质量不平衡、轴的摩擦或冲击条件以及联轴器的型式来模拟旋转机械振动情况。该设备可以模拟水泵机组的常见故障(如不平衡、不对中、碰磨等)。最终的振动信号通过传感器(测量振动和转速)传递到前置器，进行相关处理后，输送到计算机软件系统供存储、分析使用。

图1 转子振动试验台

利用转子试验台模拟了水泵机组运行过程中转子的4种典型工况：①正常工况；②不平衡工况；③不对中工况；④碰磨工况。在转子试验台上，通过对转子的主轴或者转盘进行设置来模拟不同工况。让转子正常运转得到的就是正常工况下的振动信号，要模拟不平衡工况，就需要在转盘上的螺纹孔里嵌入一个两克重的质量块；模拟不对中工况需要将法兰处两个轴的位置错开放置；模拟碰磨工况则需要将碰磨螺栓旋入碰磨螺纹箱，采集到的信号即为所需的故障信号。信号转速和采样频率都设置为相同大小，分别为1 200 r/min和2 048 Hz，共采集100组数据，每组数据都包含2 048个采样点。通过小波改进阈值降噪[9]对采集信号进行3层分解，进行降噪处理。图2为不同运行状态下降噪后的信号时域图。

图2 不同运行状态下降噪后的信号

从图2可以看出，在经过了降噪处理后，4种运行状态下的振动信号中的大量噪声成分被去除，留下的是真实有用的适合后续故障诊断特征提取及模式识别的信号。并且在四种运行状态下，去噪后可以直观地看到不同运行状态采集到的信号形状各不相同，在分辨率允许的情况下，肉眼可以根据信号波形对水泵机组运行状态进行粗略的分类。但是，为了能精确区分出水泵机组运行状态，需要将这些信号波形转化为可以量化的特征量。将这些特征量输入到模式识别模型中可以进行水泵机组的智能诊断。

3.2 振动信号特征提取

本文使用小波包分析和样本熵相结合的方法对降噪信号进行特征提取，借此区分水泵机组几种不同的运行状态。

首先使用DB8小波[10]对降噪后的信号进行小波包分解，分3层，其小波包分解树结构如图3所示。经小波包3层分解后，包含8个节点，每个节点对应于不同的频带。信号的采样频率为2 048 Hz，由香农采样定理可知，信号的奈奎斯特频率为1 024 Hz。可知，经小波包分解后对应节点对应的频率范围如表1所示。

图3 小波包3层分解树形结构

表1 小波包节点对应的频率范围

图4～图7是不同运行状态下的水泵机组振动信号经小波包分解后得到的不同节点的小波包系数。从图中可以看出，不同的运行状态下，信号的能量集中在不同的频带范围内，这为信号的特征提取提供了可能。计算各频带的样本熵，作为特征向量，记作V=[SamEn1,SamEn2,SamEn3,SamEn4,SamEn5,SamEn6,SamEn7,SamEn8]。

从表2中的样本熵值可以看出，不同运行状态下，相同频带的样本熵值有很大不同，可以很明显地区分不同运行状态，对水泵机组的不同运行状态具有较高的区分度。

4 基于小波包变换和样本熵的特征提取方法验证

为了定性定量地验证小波包变换与样本熵结合的振动特征提取方法的实际效果。本文分别采集了水泵机组四种不同运行状态下的振动信号若干组，使用学习向量量化(learning vector quantization，LVQ)神经网络[11]对样本进行分类。

LVQ算法的基本思想是：计算距离输入向量最近的竞争层神经元，从而找到与之相连接的线性输出层神经元，若输入向量的类别与线性输出层神经元所对应的类别一致，则对应的竞争层神经元权值沿着输入向量的方向移动；反之，则对应的竞争层神经元权值沿着输入向量的反方向移动。

LVQ神经网络的结构如图8所示，它由三层组成：即输入层、隐含层和输出层。该网络在输入层与隐含层间为完全连接，而在隐含层与输出层间为部分连接，每个输出神经元与隐含神经元的不同组相连接。

图4 正常运行状态下的振动信号小波包分解

图5 不平衡运行状态下的振动信号的小波包分解

图6 不对中运行状态下的振动信号的小波包分解

图7 碰磨运行状态下的振动信号的小波包分解

状态SamEn1SamEn2SamEn3SamEn4SamEn5SamEn6SamEn7SamEn8正常不平衡不对中碰磨0.41060.35120.52080.45330.47560.37300.63240.58980.49000.38520.66110.60290.48610.38130.62820.57790.73420.74600.76090.73751.03531.01230.85850.83901.08120.99960.67520.81011.46081.42211.13731.2664

LVQ神经网络有以下优点：

(1)结构简单，只需要三层网络就可以实现模式识别；

(2)不存在BP网络有可能陷入局部最小问题；

(3)收敛速度比BP网络更快；

(4)不需要将输入向量进行归一化、正交化，只需要直接计算输入向量与竞争层之间的距离，从而实现模式识别,因此简单易行，识别效率更高。

图8 LVQ神经网络结构图

使用LVQ神经网络验证基于小波包分解与样本熵的水泵机组振动信号特征提取方法的有效性，随机选取样本中60组用于训练样本，40组作为测试样本。利用神经网络工具箱[12]中的newlvq( )函数创建LVQ神经网络，其中隐含层神经元个数为20。LVQ神经网络测试结果如图9所示。

图9 测试集实际分类和预测分类图

从图9可以看出，40组测试样本中仅有一组数据分类出现偏差，由此得出测试样本的分类效果识别率达到97.5%，分类效果较好。由此说明基于小波包变换和样本熵相结合的特征提取方法对水泵机组不同运行状态具有较好的识别率。

5 结 语

本文在前人研究的基础上，在信号特征提取方面提出了新的观点和方法，该方法将小波包变换和样本熵两种工具结合起来，既考虑了信号的细节成分，同时样本熵算法不依赖数据长度，因而使故障特征提取更精确。通过对转子试验台采集的振动信号的处理与分类验证表明，所提出的方法具有一定的优越性，可以作为信号处理方法的借鉴。

□

[1] 赵康德.基于小波变换和Lipschitz指数的水泵故障诊断研究[D]. 江苏扬州：扬州大学，2010.

[2] 梁武科,张彦宁,罗兴锜.水电机组故障诊断系统信号特征的提取[J].大电机技术,2003,(4):53-56.

[3] 彭文季,罗兴锜,赵道利.基于频谱法于径向基函数网络的水电机组振动故障诊断[J].中国电机工程学报，2006,(9):155-158.

[4] 赵道利,梁武科,罗兴锜,等.水电机组振动信号的子带能量特征提取方法研究[J].水力发电学报，2004,(6):116-119,115.

[5] 单立场.小波变换在吸毒者心率变异性信号分析中的应用[D].重庆:重庆大学,2007.

[6] Richman JS, Moorman JR. Physiological time-series analysis using approximate entropy and sample entropy[J]. American Journal of Physiology-Heart and Circulatory Physiology, 2002,278(6):2 039-2 049.

[7] R Alcaraz, J J Rieta. A review on sample entropy applications for the non-invasive analysis of atrial fibrillation electrocardiograms[J]. Biomedical Signal Processing and Control, 2010

[8] Pincus S M. Assessing serial irregularity and its implications for health [J]. Ann. N. Y. Acad. Sci, 2002,954:245-267.

[9] 安周鹏,肖志怀,孙召辉,等.改进小波阈值降噪算法在水电机组信号处理中的应用[J].中国农村水利水电，2014,(12):165-168,172.

[10] 费佩燕,刘曙光.几种常见小波应用性能分析[C]∥中国电子学会第七届学术年会论文集，2001.

[11] 刘文轩.基于LVQ神经网络的配电网故障定位方法[J].电力系统保护与控制,2012,40(5):90-95.

[12] 史 峰.MATLAB神经网络30个案例分析[M]. 北京:北京航空航天大学出版社,2010.