中科大机械设备非平稳信号的故障诊断原理及应用讲义05谐波小波原理及其工程应用

第五章谐波小波原理及其工程应用小波分析中被广泛使用的Daubechies类小波与样条小波都是实小波，它们没有明确的解析式，对信号的小波分解是通过构造相应的正交滤波器{hk}和{gk}运用Mallat算法实现的。除了这两类小波，其它类型的一些小波基函数也被陆续构造出来并且得到了深入研究和工程运用[1-4]。本章介绍的谐波小波是由剑桥大学D.E.Newland教授提出的。谐波小波是一种复小波，在频域紧支，有明确的函数表达式，其伸缩与平移构成了L2(R)空间的规范正交基[5]。谐波小波分解算法是通过信号的快速傅立叶变换（FFT）及其逆变换（IFFT）实现的，算法速度快，精度高，因而具有很好的工程应用价值。5．1谐波小波的定义及正交性小波是满足允许条件的函数，如果一个小波具有完全“盒形”的频谱将是非常理想的。从这一考虑出发，设有实偶函数we(t)和实奇函数wo(t),它们的傅利叶变换分别为其它i/4πωi/4πωw(t)={exp(i4πt)–exp(i2πt)}/i2πt称(5.1.4)式定义的函数为谐波小波(harmonicwavelet)，它是复小波，在频域紧支，且具有完全“盒形”的频谱。其实部与虚部如图5.1.2所示：Re(wRe(w(t))tt图5.1.2谐波小波的实部和虚部根据小波理论对谐波小波进行伸缩、平移就生成谐波小波函数族（j,k∈Zw(2jt-k)={exp(i4π(2jt-k))-exp(i2π(2jt-k))}/i2π(2它在时间尺度上是（5.1.4）被拉伸或压缩的结果，而位置会沿着时间轴运动k个新尺度可以证明谐波小波构成一个正交系。设w(t)伸缩平移得到函数族为v(t)，即1/4π,j=11/8π,j=21/16π,j=30图5.1.3不同层谐波小波的频谱令p=2jt–k，则t=(p+k)/2j,dt=2-jdp，于是说明随着小波层（即j）的变大，谐波小波的频谱宽度倍增而幅值降低，如图5.1.3。对于谐波小波w(t)及其伸缩族w(2jt–k)(j,k∈Z)，计算它们的内积：且由于傅利叶变换是L2保范的[7]，得当j≠0，W(ω)与V(ω)在频域中总处于不同的频段，因而总有说明处于不同层的谐波小波总是正交的。说明处于第零层的谐波小波也是正交的。对其他层，以上结论可以类似得到。这样，就证明了w(t)及其伸缩平移函数族（5.1.5）构成信号的正交基。因而，以谐波小波作为基函数系就可以将信号既不交迭，又无遗漏地分解到相互独立的空间，实现将信号成分分解到不同频段。在故障诊断过程中通过谐波小波分解，可以使得故障信息从强烈的信号背景中分离出来，有利于机器故障特征的提取。从频谱图5.1.3可以看出，谐波小波对信号的分析频宽从高频到低频是以1/2关系逐渐减小的，对信号的低频部分划分比较细，而高频部分划分比较粗，这说明谐波小波分解是一种小波分解。5．2Newland快速算法及时频剖面图谐波小波构成了L2(R)空间的规范正交基，则任何信号x(t)∈L2(R)都可以表示为谐波小波的线性和，即这就是信号的谐波小波展开。aj,k为函数x(t)的小波展开系数，aj,k=<x(t),w(2jt_k)>(j,k∈Ζ)用求内积的方法计算小波展开系数运算量太大，是很不实用的。因此谐波小波的提出者Newland给出了一种快速算法[5]，可以快速而精确地求得谐波小波分解，对谐波小波运用于工程实践有很大好处。Newland快速算法是谐波小波分解算法，算法通过信号的FFT和IFFT实现[5]。设有离散信号x(r)，r=0,…,N–1,其中N=2n，其谐波小波分解为as,s=0,…,N–1。令as由Fs经分段、对每一段作IFFT得到，下两式为其表达式：aN-s=as(5.2.1)图5.2.1为一长度16的实序列数据的谐波小波分解示意图。FFTFFTFFFFFFFFFFFFIFFTx0x1x2x3x4x5x6x7x8x9x10x11x12x13x14x15a0a1a2a图5.2.1谐波小波算法图示系数as为一复序列，它包含信号时频两方面的信息。as某一段对应谐波小波分解的每层2j项。as关于aN/2两两互为复共轭[5]，则|as|=|aN-s|，所以在谐波小波时频图中通j=-1,0j=1,j=2,j=log2(N/4)图5.2.2as的分段与j的关系谐波小波分解结果一般用小波时频图（WaveletTime-FrequencyMap）和时频图的等高线图（ContourPlot）直观表示[5]。小波时频图构成如图5.2.3所示：水平面称基平面，上面一轴表示时间，与之垂直的另一轴表示小波层。运用上面的分段结果，将属于不同j的as沿频率轴排列，属相同j的as沿时间轴排列。这样基平面就被划分成网格。在各网格以as模的平方为高作柱体就构成了谐波小波时频图。小波时频图是随|as|2起伏的面。这里高度取lg|as|2。tj图5.2.3谐波小波时频图构成小波时频图物理意义由Parseval公式得到。序列xr与谐波小波分解系数as有[5，6]：式左边为信号的平均能量，右边为|as|2的加权和，说明谐波小波分解结果表明不同频率和时间的谐波小波能量对整个信号能量贡献的大小。谐波小波时频图是分解结果的直观表示，其起伏对应不同谐波小波能量的相对大小。通过谐波小波时频图，我们可以知道在什么时间什么频率成分对信号组成有重要影响。作等高线图同样可以直观表示谐波小波的分解结果，如下节例所示。及谐波小波分解时频图。该信号是单一频率的，所以谐波小波分解只有一个层有值，在小波时频图上表现为对应的层有峰值。/图5.2.4信号x(t)的时域波形与小波时频图谐波小波分解对信号中存在的奇异点极其敏感，可以用来检测信号中存在的微小的奇异点。例如有模拟信号r=r/320在正弦的8个周期设有奇异点x1(i×64+16)=x1(i×64+16)+0.04(i=0,…,7)(5.2.3)jjt/st/s(a)x1的时域局部细化(b)x1的谐波小波等高线图jtjt/s(c)x1的谐波小波时频图(d)x2的谐波小波等高线图(e)x2的谐波小波时频图tt/s(f)x2的谐波小波时频剖面图（j=7）即在特定位置有一点值与前一点相同。从图5.2.5（a）信号的时域波形细化可以看到这些微小的平台（第一、四波谷处）；谐波小波分解得到的等高线图5.2.5（b）和小波时频图5.2.5（c）上明显地有它们相对的峰。说明信号中微小的奇异成分从原信号中被分离出来，或者说，信号的微弱成分被“放大”了。这在故障诊断的工程实际中是很有意义的。系统故障最初一般很微弱，被完全淹没在整个信号背景中。通常的信号处理方法难以将它分离出来因而不能予以发现和排除，使“小”故障逐渐发展为“大”故障，最终影响系统运行。谐波小波对微小奇异信号的分离、放大功能为这一类潜在故障的诊断提供了一种工具。（5.2.2）式为一纯正弦信号，实际应用中，这种不被噪声污染的完全“干净”的信号是不存在的。可以想象，微小的噪声干扰对谐波小波的分解结果影响都会很大。这时，谐波小波分解是否能提取出信号中包含的较弱的奇异信息就成为问题。例如对（5.2.2）式引入随机干扰x2(r)=sin(2π×15×tr)+rand/10r=r/320同时在正弦的某些位置修正x2(i×64+16)=x2(i×64+16)+0.2(i=0,…,7)这是噪声背景下有周期性奇异点的信号。其谐波小波分解结果如图5.2.5（d,e,f）所示。这里尽管奇异信号的幅值较（5.2.3）式增大至5倍，小波时频图（e）和等高线图（d）还是显得非常杂乱，几乎没有提供什么信号的特征信息。说明通常所用的小波时频图和等高线图在实际分析中已显无力。为了弥补小波时频图和等高线图的不足，我们提出用一组通过j垂直于基平面并平行于时间轴的平面来剖时频图，称平面与谐波小波时频图的交线为谐波小波时频剖面图（Time-FrequenceProfilePlot）。时频剖面图将三维的时频图分离为一组平面图形，因而可以细致地研究小波分解的每一层。作图5.2.5（e）时频图的谐波小波时频剖面图（图5.2.5（f发现在最高层（j=7）的时频剖面图有8个明显较高的、等间隔的峰。这些峰反映了式（5.2.5）中的突变信号，说明谐波小波时频剖面图可以有效检测出噪声背景下的信号突变成分。由于信号中迭加的突跳是宽带信号，所以在谐波小波分解的最高层，即信号高半频有分量。谐波小波时频剖面图表明了谐波小波分解同一层内小波幅值（能量）的大小，通过它我们可以分频段研究信号，又可以将一组时频图综合考虑研究整个信号的性质。同时，剖面图将小波分解结果变成一组一维信号，因而可用一维信号的处理方法来更充分地研5．3谐波小波滤波大型旋转机械状态监测与故障诊断中，往往利用机组同一截面两路相互垂直振动信号的合成轴心轨迹来监测其运行状态和识别故障类型。当设备出现故障时，信号表现出非平稳特性，因而传统的基于保相滤波的FFT轴心轨迹已不太适用。小波对处理非平稳信号是非常有效的，我们是否可以用相互垂直的X方向与Y方向的小波分解结果来合成轴心轨迹呢？众所周知，大部分小波，比如样条小波与Daubechies类小波都没有函数表达式，它们的小波变换是通过Mallat算法和相应的小波包Wickerhause树形算法实现的。Mallat算法分解时要隔二抽一，从而使得小波分解各层的数据点数和采样频率随尺度逐渐减小。直接将转子垂直、水平方向振动信号小波分解同一尺度同一频段的数据合成轴心轨迹，将使轴心轨迹不具有可比性，且由于数据点数减少、采样频率降低会使合成的轴心轨迹失真，这种直接的合成是不合适的。本节将述及应用谐波小波方法将信号相同尺度，相同频段的成分从原信号中分离出来，且保持数据点数与采样频率不变，进而实现旋转机械振动信号不同尺度不同频段轴心轨迹的合成与分析。最后，利用分形方法研究了机组不同层不同频段轴心轨迹的不规则度。作为轴心轨迹的不规则度指标，它本质上描述了机组运行的非线性、非平稳特性。谐波小波实际上是一个完全理想的带通滤波器，这由它具有的“盒形”（Box-like）频谱可以推知。因此谐波小波具有良好的滤波特性。实际上，同样可以用下面的方法定义谐波小波[4]：其它其中m,n决定了谐波小波变换的尺度（j且n=2m，当m=0时，n=1。可以看出谐波小波在时域中紧支特性很差，这是因为它的“盒形”频特性导致的，且光滑性远比Daubechies类小波高。同时，由于谐波小波具有相互成90°相差的实部偶小波和虚部奇小波，由数字信号处理基本知识可知[8]，实部偶小波和虚部奇小波都是零相移滤波器，这就是说谐波小波具有“锁定”信号相位的能力。我们知道振动信号的相位信息在诊断中是十分重要的，转子不平衡类型识别，转子裂纹等故障对相位信息都十分敏感，特别是在建立振动信号轴心轨迹时，相位无疑更重要的。另一个方面，由于Daubechies类小波本身的不规则性，若应用该类小波建立轴心轨迹，将由于小波本身的特性导致轴心轨迹看似能量突变点多而复杂，容易引起误诊。Lemare-Meyer小波相对于Daubechies类小波，其光滑性虽大大提高，但它和Daubechies类小波一样，不具有明显的时、频域表达式因而不得不采用Mallat塔形算法。我们知道，Mallat塔形算法是一个隔二抽取算法，这种算法结果使得不同尺度下各频段序列的采样频率不一致，随着尺度增加，采频减半，数据点数减半，那么将导致不同尺度下的轴心轨迹突变点增多，很易导致误诊。当然，我们可通过剔除某些频道序列后，用重构算法进行重构之后，再利用X方向、Y方向信号进行合成来进行轴心轨迹分析，但对厂方技术人员应用来说过于繁由于谐波小波的光滑性，“盒形”谱特性，零相移特性以及明显的数学表达式，使得我们可构造出不同尺度下各频段序列数据点数不变、采样频率不变的算法，最终成功一般地，对一定尺度的小波ψ(t)，信号x(t)的小波变换可简单表示为：则x(t)相对于尺度j的谐波小波ψm,n(t)的小波变换为[4]：根据傅立叶变换的性质，m,n由于谐波小波有明确的频域函数表达式，而x(t)的傅立叶变换可通过FFT得到，Λ则上式是很容易计算的。得到了Wx(⑴)，作它的傅立叶逆变换（IFFT）就得到了由m,nΛ决定的尺度j下的信号谐波小波变换，该过程实现如图5.3.1所示：FFTΛ图5.3.1尺度j下的信号谐波小波变换设离散信号x(t)的采样频率为Js，则分析频率JN=Js/2，此时对不同的m,n，信号频率划分如图5.3.2所示：JJN/2JN/4JN/2 JN/4m=00JNmmnn图5.3.2不同m,n对信号频率的划分如何理解上述算法与谐波小波原理的关系呢？仔细研究图5.2.1所示的谐波小波分解算法，考虑到傅立叶变换与傅立叶逆变换之间的互逆性会很自然地得到谐波小波分解的重构算法。例如，图5.2.1所示的16点数据分解的谐波小波重构算法可如下图示：IFFTIFFTFsx(r)谐波小波各层系数FFT图5.3.3谐波小波重构算法为了对信号的某一特定频段的成分进行研究，在对信号的谐波小波分解进行重构时可将其它频段的谐波小波系数置为“0”，只保留该段的小波系数，由于谐波小波的正交性，如此重构的结果只包含信号该频段的成分，其他成分都被剔除了。这个算法与本节开始所给出的算法是一致的，实际是谐波小波重构算法的延伸，是对信号进行了滤波,我们称这一过程为谐波小波滤波。由该算法描述可以看到，计算过程并未采用基于隔二抽取的Mallat算法，因此保证了信号各频段成分点数不变，采样频率不变，这样就可以实现机组同一截面互相垂直两个方向振动信号的轴心轨迹合成。另一方面，如前所述，谐波小波是一个理想带通滤波器，谐波小波滤波可以将任何2(R)正交、无冗余、无泄露地分解到相互独立的频段上。各频段是相互正交的，因而互不干扰。特定频段的成分与信号的其它频率成分通过谐波小波分解被分离了，从而消除了其它频段成分对该频率段的影响，使一些被淹没的较弱的信号得以突显出来，等于提高了信噪比，无疑对研究机组运行中出现的微弱故障，提取故障特征量很有帮助。一个旋转机组如果轴承系统是各向同性的话，其各谐波分量的轴心轨迹将是一个圆。当系统发生故障或出现异常，轴心轨迹将变得十分复杂且不规则，研究其各谐波分量的轴心轨迹会得到许多故障信息。谐波小波滤波算法实现了不同频段（层）信号的分离，又保持了各段数据点数不变，采样频率不变，克服了小波分解Mallat算法引起的合成轴心轨迹的失真和信息丢失的缺点，使机组振动信号各谐波分量的轴心轨迹分析更加有效和精确，对工程中运用轴心轨迹更好地研究机组故障有很大的帮助。与小波包分解类似，谐波小波滤波也可以处理信号高半频成分。小波包分解第j层有2j段，频段宽fN/2j，m=i×fN/2j，n=(i+1)×fN/2j，对信号的频率划分如图5.3.4所示:mn(m)0fN/2fNmn(m)n(m)n(m)nL___________上_____L___________L___________0fN/4fN/23fN/4fN图5.3.4谐波小波包变换图示5．4谐波小波应用举例5．4．1谐波小波分解及时频剖面图在旋转机组振动特征提取中的运用某炼油厂重催装置由烟机、风机、齿轮箱、电机组成。在运行中用涡流传感器采集了各机组轴瓦处径向位移振动信号。数据参数为：机组转速：5850r/min采样频率：2000Hz机组中风机1瓦和相邻的烟机2瓦振动较大。用谐波小波分解分析风机1瓦径向位移振动的一组数据（数据长1024得到结果如图5.4.1所示。t/s（a）时频图（b）时频剖面图（j=8）JJ/Hz0.20tt/s（c）图（b）的傅立叶变换（d）时域波形细化图5.4.1某炼油厂重催振动数据的谐波小波分解图5.4.1（a）谐波小波时频图上两条显著的峰带反映了机组运行的工频和一个较低频率。时频图最高层（j=8）有较低起伏的峰与谷，表示信号的高频段有一些微弱的分量，但时频图很杂乱，难以得到更多的信息。这时作时频图最高层的时频剖面图（图5.4.1（b看到时频剖面图上有许多等间隔的峰。谐波小波分解的最高层包含信号的高半频成分，说明信号中有等间隔出现的高频分量。为从剖面图推知这些干扰的周期，作图5.4.1（b）的傅立叶频谱如图5.4.1（c）所示。谱图表明（b）中峰出现的周期为1/97.656Hz=10.24ms,而97.656Hz正是重催系统运行的工频，说明系统运行中每个周期都有一个高频的扰动产生。图5.4.1（d）为该组振动数据的波形细化，研究表明这些扰动就是其中箭头所示的一系列细微“毛刺”。分析风机1瓦相邻测点（烟机2瓦）的数据得到了相同的结果，表明它们是由同一原因造成的。剖面图频谱5.4.1（c）中二倍工频较大。综合考虑这些因素，推测信号中的高频扰动是由于烟机、风机轴系的不对中引起的，它们之间的膜片式联轴器为了补偿不对中而出现每周期一次错动，造成微弱的振动。这些分析为预知维修提供了依据。5．4．2谐波小波轴心轨迹阵列的实现及其不规则度描述图5.4.2是某大型化肥厂CO2压缩机发生喘振时，高压缸水平方向（X方向）和垂直方向（Y方向）由涡流式位移传感器拾取的振动信号，其中图5.4.2(a)为X方向，图5.4.2(b)为Y方向。测量参数为：转子转速：6530rmp采样频率：2000Hz图5.4.2(c)是原始轴心轨迹，从图中看轴心轨迹较为复杂且不规则，加之较小的高倍工频分量影响使得轴心轨迹有一些局部能量突变点，且其分形盒维数也比较大（分形内容详见第10章）。令d,k表示信号小波分解第j层第k频段的分形盒维数，Oj,k表示相应的轴心轨迹。如图5.4.2(d)所示d,0=1.(a)X方向振动信号的时域波形(b)Y方向振动信号的时域波形logNrlogNr(c)原始轴心轨迹O0,0(d)原始轴心轨迹分形盒维数d,0图5.4.2原始信号的时域波形、轴心轨迹和轴心轨迹分形盒维数图5.4.3(a)、(b)是X方向、Y方向信号的第2层谐波小波包分解，从中看出主要小波成分集中在第0频段，为抓住主要矛盾，对第0频段小波进行合成得到的轴心轨迹如图5.4.3(c)所示，第0频段小波对应的是低频喘振、工频振动和二倍频振动的特征，从图5.4.3(c)中看高倍工频分量影响已剔除，轴心轨迹光滑度提高，不规则度减少，其分形盒维数相对原始轴心轨迹也有所减少，如图5.4.3(d)所示为d,0=1.3536。(a)X方向信号的第2层谐波小波包分解(b)Y方向信号的第2层谐波小波包分解logNrlogNr(c)第0频段合成轴心轨迹O2,0(d)第0频段合成轴心轨迹的分形盒维数d,0图5.4.3X方向、Y方向信号的第2层谐波小波包分解，第0频段合成轴心轨迹及其分形盒维数图5.4.4(a)、(b)是X方向、Y方向信号的