分数阶对称平移不变过完备小波的构造及其在轴承故障诊断中的应用

工程科学学报，第37卷，第3期：378-384,2015年3月 Chinese Journal of Engineering,Vol.37,No.3:378-384,March 2015 DOI:10.13374/j.issn2095-9389.2015.03.018:http://journals.ustb.edu.cn 分数阶对称平移不变过完备小波的构造及其在轴承 故障诊断中的应用 沈政伟区，史天，申亚男 北京科技大学数理学院，北京100083 ☒通信作者，E-mail:szw2628@sina.com 摘要提出了一种分数阶的对称性近似平移不变过完备小波的构造方法.首先，给出一种构造具有对称性且具有最小长 度的低通滤波器方法.其次，通过拓普利兹矩阵分解法求出对应的具有近似平移不变性的高通滤波器，此方法比其他分解方 法具有更低的计算复杂度.此外，利用此构造方法，也得到具有更高阶消失矩的分数阶过完备小波变换.最后，将构造出的分 数阶对称平移不变过完备小波应用到轴承故障诊断中.实验结果表明，提出的小波变换能有效地提取出轴承的故障特征. 关键词小波变换：对称滤波器；轴承：故障诊断 分类号TP391.4 Construction of a symmetrical shift-invariant fractional overcomplete wavelet and its application in bearing fault diagnosis SHEN Zheng-wei,SHI Tian,SHEN Ya-nan School of Mathematics,University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:szw2628@sina.com ABSTRACT This article introduces the design of symmetrical approximately shift-invariant fractional overcomplete wavelet trans- forms.First,a design scheme for the symmetrical low-pass filter with minimum-ength was proposed,and then the corresponding high- pass filters with approximately shift-invariant properties were constructed via Toeplitz matrix factorization,which has a lower computa- tional complexity than other methods.In addition,fractional overcomplete wavelet transforms could be designed with higher vanishing moments through the method proposed.Subsequently,a bearing fault diagnosis scheme was proposed using the symmetrical shift- invariant fractional overcomplete wavelet transforms.Experimental results show that the bearing faults can be detected effectively using the symmetrical shift-invariant fractional overcomplete wavelet transforms. KEY WORDS wavelet transforms:symmetrical filters:bearings:fault diagnosis 分数阶过完备小波仅在低通通道对信号进行下采 仍然能够准确地监测到信号中的奇异点，因此方便机 样，因此会产生数据冗余，这种适当的冗余为信号的处 械故障检测中对高频故障信号的准确定位 理，特别是为信号去噪提供了更多的自由度.同时，它 目前，分数阶过完备小波的构造方法还有待改进， 可以对待处理信号从时域和频域同时提高采样率，从 计算效率还可以进一步提高，而且线性相位特性在已 而能更有效地对信号进行时频分析”.此外，分数阶 有的滤波器设计中网并未考虑.小波滤波器具有线性 过完备小波还具有近似平移不变性，因此可以增强信 相位特性指的是滤波器应该具有对称性，这种性质对 号分析的准确性.尤其在信号发生平移的情况下，它 于由于信号自身产生的延迟或者平移具有很强的鲁棒 收稿日期：201401-10 基金项目：教育部新世纪优秀人才支持计划资助项目(NCET-110574)

工程科学学报，第 37 卷，第 3 期: 378--384，2015 年 3 月 Chinese Journal of Engineering，Vol． 37，No． 3: 378--384，March 2015 DOI: 10． 13374 /j． issn2095--9389． 2015． 03． 018; http: / /journals． ustb． edu． cn 分数阶对称平移不变过完备小波的构造及其在轴承 故障诊断中的应用 沈政伟，史 天，申亚男 北京科技大学数理学院，北京 100083  通信作者，E-mail: szw2628@ sina． com 摘 要 提出了一种分数阶的对称性近似平移不变过完备小波的构造方法． 首先，给出一种构造具有对称性且具有最小长 度的低通滤波器方法． 其次，通过拓普利兹矩阵分解法求出对应的具有近似平移不变性的高通滤波器，此方法比其他分解方 法具有更低的计算复杂度． 此外，利用此构造方法，也得到具有更高阶消失矩的分数阶过完备小波变换． 最后，将构造出的分 数阶对称平移不变过完备小波应用到轴承故障诊断中． 实验结果表明，提出的小波变换能有效地提取出轴承的故障特征． 关键词 小波变换; 对称滤波器; 轴承; 故障诊断 分类号 TP391. 4 Construction of a symmetrical shift-invariant fractional overcomplete wavelet and its application in bearing fault diagnosis SHEN Zheng-wei ，SHI Tian，SHEN Ya-nan School of Mathematics，University of Science and Technology Beijing，Beijing 100083，China  Corresponding author，E-mail: szw2628@ sina． com ABSTＲACT This article introduces the design of symmetrical approximately shift-invariant fractional overcomplete wavelet trans￾forms． First，a design scheme for the symmetrical low-pass filter with minimum-length was proposed，and then the corresponding high￾pass filters with approximately shift-invariant properties were constructed via Toeplitz matrix factorization，which has a lower computa￾tional complexity than other methods． In addition，fractional overcomplete wavelet transforms could be designed with higher vanishing moments through the method proposed． Subsequently，a bearing fault diagnosis scheme was proposed using the symmetrical shift￾invariant fractional overcomplete wavelet transforms． Experimental results show that the bearing faults can be detected effectively using the symmetrical shift-invariant fractional overcomplete wavelet transforms． KEY WOＲDS wavelet transforms; symmetrical filters; bearings; fault diagnosis 收稿日期: 2014--01--10 基金项目: 教育部新世纪优秀人才支持计划资助项目( NCET--11--0574) 分数阶过完备小波仅在低通通道对信号进行下采 样，因此会产生数据冗余，这种适当的冗余为信号的处 理，特别是为信号去噪提供了更多的自由度． 同时，它 可以对待处理信号从时域和频域同时提高采样率，从 而能更有效地对信号进行时频分析［1］． 此外，分数阶 过完备小波还具有近似平移不变性，因此可以增强信 号分析的准确性． 尤其在信号发生平移的情况下，它 仍然能够准确地监测到信号中的奇异点，因此方便机 械故障检测中对高频故障信号的准确定位． 目前，分数阶过完备小波的构造方法还有待改进， 计算效率还可以进一步提高，而且线性相位特性在已 有的滤波器设计中［2］并未考虑． 小波滤波器具有线性 相位特性指的是滤波器应该具有对称性，这种性质对 于由于信号自身产生的延迟或者平移具有很强的鲁棒

沈政伟等：分数阶对称平移不变过完备小波的构造及其在轴承故障诊断中的应用 ·379· 性.在信号处理中，系统具有线性相位特性是必要的， 非常繁琐，除了需要用计算复杂度高的矩阵谱因式分 因为这样的系统非常适合处理已经发生平移或者变形 解方法，还需要在此基础上对仿酉矩阵(paraunitary 的信号.本文在对分数阶过完备小波变换理论研究的 matrix)进行频率响应修正.本文在此基础上，提出了 基础上，提出了一种分数阶对称过完备小波低通滤波 一种改进的分数阶对称平移不变过完备小波的构造方 器的构造方法，并给出了相应的构造表达式：然后利用 法，计算方法简单直观，并且构造的分数阶小波具有近 更为简捷的拓普利兹(Toeplitz)矩阵分解法，求解出对 似对称性、平移不变性及较高的消失矩 应的近似平移不变的高通滤波器，从而完成分数阶的 设山()(1∈R)是一小波函数，信号与这个小波函 具有近似对称性（平移不变性）的小波滤波器组的设 数的平移量b以及a倍尺度因子伸缩量√a地(al-b)的 计：最后，将构造出的小波变换应用于滚动轴承故障检 卷积就是连续小波变换，其中a,b∈R.将参数a,b离 测，通过实验比较（与dN系列小波比较）发现由于具 散化{2“，k}。:z即可得到严格采样的离散小波变换； 有对称、平移不变以及较高消失矩等良好特性，本文提 但是，这种小波变换不具有平移不变性.分数阶过完 出的分数阶对称平移不变过完备小波对故障的监测效 备小波√a山(at-b)的参数a和b离散化为{g/p",spk/ 果更为理想 q}。,:z,其中p,qs∈Z,经过这种离散化后的小波伸缩 滚动轴承是一般机械设备中非常重要的零部件， 由于在使用过程中经常会承受较高的冲击负荷，因此 因子是分数阶的，即为gp,并且具有上1 s q/p-1 倍冗余 容易产生磨损（监测信号中有噪声）或者轴承断裂等 度，从而使得该小波变换是分数阶的且是过完备的，这种 现象.当滚动轴承发生故障时，往往会引发高频冲击 过完备性质可以保证构造出的小波具有平移不变性 振动四，所以对故障振动信号高频部分的研究在轴承 本文中取q=3、p=2和s=1的情况为例进行讨 故障检测中显得尤为重要.小波分析技术由于具备良 论（其他情况类似），这是一个3/2阶的过完备小波变 好的时频局部分析能力，在信号处理方面已经被广泛 换.理想的3/2阶过完备小波变换的分解与重构过程 地应用：而基于小波技术的轴承故障检测-)也同样 可以用图1表示，其中H(z)是低通滤波器，G(z)是高 取得了良好的效果.但是，这些方法都没有考虑信号 通滤波器组. 自身带有延迟或者发生偏移的情况.本文利用构造的 具有对称性的分数阶过完备小波变换，分别对某大型风 2 Me) 12 机轴承故障信号以及电主轴轴承故障信号进行监测与诊 CE Ge) 断.实验结果表明，与其他小波变换处理的结果比较，本 图1理想的32阶过完备滤波器组 文构造的小波变换能更准确的监测出故障信号 Fig.1 Ideal 3/2 overcomplete filter bank 1分数阶过完备小波变换 从图中可以看出，信号通过低通滤波器分支是3/ 2倍采样，而通过高通滤波器没有进行采样.因此，整 分数阶离散小波变换作为对严格采样正交小波 个滤波器组也只能是一个近似平移不变的的滤波器 变换的一种改进，能更加有效地在时频域分析信号：它 组，因为只有高通部分是平移不变的.之所以说该滤 虽然具有更高的频域分辨率，但是时域分辨率却不够 波器组是理想的是因为除了一组特解（无实际意义） 好，并且也不具备对称性，从而影响了信号分析的精 外，图1中的滤波器组是不可能利用有限冲击响应 度.基于此，文献]首先提出了一种基于迭代法的分 (finite impulse response,FIR)滤波器来实现完全可重 数阶离散小波变换的构造方法，证明了随着分数阶小 构的(perfect reconstruction).然而，实际的信号处理中 波函数的正则性（光滑性）的提高，由平移导致的误差 需要满足完全可重构并且是有限冲击响应的滤波器 可以达到任意小.在文献0]中，作者使用Grobner基 组.因此，需要对图1中的滤波器组进行改进.基于 方法构造了分数阶正交滤波器组，但是该方法计算较 此，文献2]提出了如图2所示的滤波器组 耗时，且仅能求解一些简单实例，之后作者又利用多项 从图2中看出，图1中的高通滤波器G(z)被分解 式矩阵谱因式分解方法求出了分数阶小波紧框架.基 于之前的相关研究⑧四，文献2]详细介绍了分数阶 2 H(e) 过完备小波变换的特性，并与分数阶严格采样离散小 G日目 3→3 G) 波变换进行了详细对比：且在给定消失矩的情况下，通 过矩阵谱因式分解求出了具有最小长度的完全可重构 G同 03→3 G, 的滤波器。然而，文献2]在构造低通滤波器时没有考 3→③G) 虑对称性等滤波器特性，并且整个滤波器组远未达到 图2改进的3/2阶过完备滤波器组4 平移不变的要求：另外，在构造相应高通滤波器时方法 Fig.2 Improved 3/2 overcomplete filter bank

沈政伟等: 分数阶对称平移不变过完备小波的构造及其在轴承故障诊断中的应用 性． 在信号处理中，系统具有线性相位特性是必要的， 因为这样的系统非常适合处理已经发生平移或者变形 的信号． 本文在对分数阶过完备小波变换理论研究的 基础上，提出了一种分数阶对称过完备小波低通滤波 器的构造方法，并给出了相应的构造表达式; 然后利用 更为简捷的拓普利兹( Toeplitz) 矩阵分解法，求解出对 应的近似平移不变的高通滤波器，从而完成分数阶的 具有近似对称性( 平移不变性) 的小波滤波器组的设 计; 最后，将构造出的小波变换应用于滚动轴承故障检 测，通过实验比较( 与 dbN 系列小波比较) 发现由于具 有对称、平移不变以及较高消失矩等良好特性，本文提 出的分数阶对称平移不变过完备小波对故障的监测效 果更为理想． 滚动轴承是一般机械设备中非常重要的零部件， 由于在使用过程中经常会承受较高的冲击负荷，因此 容易产生磨损( 监测信号中有噪声) 或者轴承断裂等 现象． 当滚动轴承发生故障时，往往会引发高频冲击 振动［3］，所以对故障振动信号高频部分的研究在轴承 故障检测中显得尤为重要． 小波分析技术由于具备良 好的时频局部分析能力，在信号处理方面已经被广泛 地应用; 而基于小波技术的轴承故障检测［3--7］也同样 取得了良好的效果． 但是，这些方法都没有考虑信号 自身带有延迟或者发生偏移的情况． 本文利用构造的 具有对称性的分数阶过完备小波变换，分别对某大型风 机轴承故障信号以及电主轴轴承故障信号进行监测与诊 断． 实验结果表明，与其他小波变换处理的结果比较，本 文构造的小波变换能更准确的监测出故障信号． 1 分数阶过完备小波变换 分数阶离散小波变换［8］作为对严格采样正交小波 变换的一种改进，能更加有效地在时频域分析信号; 它 虽然具有更高的频域分辨率，但是时域分辨率却不够 好，并且也不具备对称性，从而影响了信号分析的精 度． 基于此，文献［9］首先提出了一种基于迭代法的分 数阶离散小波变换的构造方法，证明了随着分数阶小 波函数的正则性( 光滑性) 的提高，由平移导致的误差 可以达到任意小． 在文献［10］中，作者使用 Grobner 基 方法构造了分数阶正交滤波器组，但是该方法计算较 耗时，且仅能求解一些简单实例，之后作者又利用多项 式矩阵谱因式分解方法求出了分数阶小波紧框架． 基 于之前的相关研究［8--12］，文献［2］详细介绍了分数阶 过完备小波变换的特性，并与分数阶严格采样离散小 波变换进行了详细对比; 且在给定消失矩的情况下，通 过矩阵谱因式分解求出了具有最小长度的完全可重构 的滤波器． 然而，文献［2］在构造低通滤波器时没有考 虑对称性等滤波器特性，并且整个滤波器组远未达到 平移不变的要求; 另外，在构造相应高通滤波器时方法 非常繁琐，除了需要用计算复杂度高的矩阵谱因式分 解方法，还 需 要 在 此 基 础 上 对 仿 酉 矩 阵( paraunitary matrix) 进行频率响应修正． 本文在此基础上，提出了 一种改进的分数阶对称平移不变过完备小波的构造方 法，计算方法简单直观，并且构造的分数阶小波具有近 似对称性、平移不变性及较高的消失矩． 设 ψ( t) ( t∈Ｒ) 是一小波函数，信号与这个小波函 数的平移量 b 以及 a 倍尺度因子伸缩量槡aψ( at － b) 的 卷积就是连续小波变换，其中 a，b∈Ｒ． 将参数 a，b 离 散化{ 2n ，k} n，k∈Z即可得到严格采样的离散小波变换; 但是，这种小波变换不具有平移不变性． 分数阶过完 备小波槡aψ( at － b) 的参数 a 和 b 离散化为{ qn / pn ，spk / q} n，k∈Z，其中 p，q，s∈Z，经过这种离散化后的小波伸缩 因子是分数阶的，即为 q n / pn ，并且具有 1 s 1 q/ p － 1 倍冗余 度，从而使得该小波变换是分数阶的且是过完备的，这种 过完备性质可以保证构造出的小波具有平移不变性． 本文中取 q = 3、p = 2 和 s = 1 的情况为例进行讨 论( 其他情况类似) ，这是一个 3 /2 阶的过完备小波变 换． 理想的 3 /2 阶过完备小波变换的分解与重构过程 可以用图 1 表示，其中 H( z) 是低通滤波器，G( z) 是高 通滤波器组． 图 1 理想的 3 /2 阶过完备滤波器组 Fig． 1 Ideal 3 /2 overcomplete filter bank 图 2 改进的 3 /2 阶过完备滤波器组［2］ Fig． 2 Improved 3 /2 overcomplete filter bank［2］ 从图中可以看出，信号通过低通滤波器分支是 3 / 2 倍采样，而通过高通滤波器没有进行采样． 因此，整 个滤波器组也只能是一个近似平移不变的的滤波器 组，因为只有高通部分是平移不变的． 之所以说该滤 波器组是理想的是因为除了一组特解( 无实际意义) 外，图 1 中的滤波器组是不可能利用有限冲击响应 ( finite impulse response，FIＲ) 滤波器来实现完全可重 构的( perfect reconstruction) ． 然而，实际的信号处理中 需要满足完全可重构并且是有限冲击响应的滤波器 组． 因此，需要对图 1 中的滤波器组进行改进． 基于 此，文献［2］提出了如图 2 所示的滤波器组． 从图 2 中看出，图 1 中的高通滤波器 G( z) 被分解 · 973 ·

·380· 工程科学学报，第37卷，第3期 成了三个高通滤波器G。(z)、G(z)和G,(z),并随后进 点个数，且N≥K(i=0,1,2),并有(1+z)(1+z1+ 行3倍的采样.如果假设图2中的三个高通滤波器 z2)不能整除Q(z). G,(z)(i=0,1,2)满足 如果要得到冲击响应最短的低通滤波器，分数阶 go(n-2）=g,(n-1)=g2(n), (1) 过完备小波低通滤波器长度还应该满足下面的条件： 则图1和图2在本质上是相同的，其中g:(n)为高通滤 Lg≥3N+2K-1. (7) 波器G,()的冲激响应.因此，也同样没有满足完全平 式中，K=min(K,K,K2). 移不变的有限冲击响应解.但是，实际上图2对滤波 因此，由式(5)和(7)可知，为了得到阶次（滤波器 器组的限制和要求比图1少一些，是可以用有限冲击 最短)最低并且具有线性相位特性的低通滤波器 响应滤波器实现近似解的.下面给出图2中的滤波器 H(z),式(5)中的Q(z)因子应该是对称的，并且长度 组的完全可重构条件 应该小于2K-1;也就是低通滤波器的长度应该至少为 由文献2]知道，图2中的滤波器组满足的完全 3N+2K-1,从而使得其最短支集长度为3N+2K-1. 可重构条件为 下面给出具有线性相位特性以及最小长度的低通 H(z1)H(z)+H(-z)H(-z)+ 滤波器的构造方法.此处考虑N为奇数的情况，N为 2c,e)c,)=6, (2) 偶数时同理可得. 首先对于给定的N和K,假设低通滤波器H(z)是 H(2)H(Wz）+H(-z)H(-Wz）+ 对称的并且其冲激响应h(n)长度为L,其中L≥3N+ 2盒cG(=0. (3) 2K-1.令z=e“,那么H(z)可以表示成w的多项式. 为了求解H(z),首先进行下述转换： H(a)H(Wa)+H(-)H(-Wa)+ 2gce-y6,(=0 x=(e+2+》=im号=1-co).(8) 4 2 (4) 进而有 式中，W=e2m 构造图2中的滤波器组，并对其特性加以改进，使 1、4 =号+1+=宁1++).9) 得图2中的滤波器组具有半对称性并且满足近似平移 不变.首先，提出一种构造具有线性相位特性（即对称 1-(e+2+=宁1+)1+以. 性)的低通滤波器的方法.其次，在此基础上利用谱分 (10) 解求出对应的高通滤波器。本文对已有方法进行了改 因子1+z1可以用x表示： 进，使用拓普利兹矩阵分解方法，计算方法简单，因此 1+z1=2c0s%=2√1-x (11) 有效地提高了计算效率. 在进行一位平移后，1+z+z2也可以用x表示： 2滤波器的构造 +1+=4(层- (12) 2.1对称低通滤波器的构造 将H(z)转换成x的表达式T(x),其中A(x)对应 离散多项式信号由低通滤波器H(z)处理之后，通 Q(z),有 常情况下应该能够保留离散多项式信号中的阶次小于 或等于V-1的部分，而阶次大于N的部分可以由高 T=1-)(子-)广A. (13) 通滤波器G(z)(i=0,1,2)消除，即离散多项式信号 由z和x的零点对应关系，可以得到T(x)有下列等价 中阶次为K(i=0,1,2)的部分.因此，严格采样正交 关系成立： 小波系统应该满足N=K(i=0,1,2)的条件，而对于 T(0)=1: 分数阶过完备小波系统，为了获得更高阶的消失矩，通 T(0)=0,(i=1,2,,K-1): 常要求N≥K(i=0,1,2).基于此条件和图2中的滤 T(1）=0,(i=0,1/2,…,N/2-1): 波器组结构，本文中的32阶过完备小波低通滤波器 T0(3/4)=0,(i=0,1,…,N-1). 和高通滤波器组可以写成如下形式： 考虑H(z)在z=1处的特性，即T(x)在x=0处的特 H(z)=(1+z)(1+z1+z2)Q(z),(5) 性，有下式成立： G,(z)=(1-z)V(z),i=0,1,2. (6) 1-T(x)=B(x)x* (14) 式中，N是低通滤波器在点z=-l、z=e2=小和z=e4m 综合式(13)和(14)，得到 处的零点个数，因此可以保留多项式中阶次小于等于 W-1部分：K(i=0,1,2)为高通滤波器在点z=1处零 1-)(子-)广A)+B()=1,（15)

工程科学学报，第 37 卷，第 3 期 成了三个高通滤波器 G0 ( z) 、G1 ( z) 和 G2 ( z) ，并随后进 行 3 倍的采样． 如果假设图 2 中的三个高通滤波器 Gi ( z) ( i = 0，1，2) 满足 g0 ( n － 2) = g1 ( n － 1) = g2 ( n) ， ( 1) 则图1 和图2 在本质上是相同的，其中 gi ( n) 为高通滤 波器 Gi ( z) 的冲激响应． 因此，也同样没有满足完全平 移不变的有限冲击响应解． 但是，实际上图 2 对滤波 器组的限制和要求比图 1 少一些，是可以用有限冲击 响应滤波器实现近似解的． 下面给出图 2 中的滤波器 组的完全可重构条件． 由文献［2］知道，图 2 中的滤波器组满足的完全 可重构条件为 H( z － 1 ) H( z) + H( － z － 1 ) H( － z) + 2 ∑ 2 i = 0 Gi ( z － 2 ) Gi ( z 2 ) = 6， ( 2) H( z － 1 ) H( Wz) + H( － z － 1 ) H( － Wz) + 2 ∑ 2 i = 0 Gi ( z － 2 ) Gi ( W2 z 2 ) = 0， ( 3) H( z － 1 ) H( W2 z) + H( － z － 1 ) H( － W2 z) + 2 ∑ 2 i = 0 Gi ( z － 2 ) Gi ( Wz2 ) = 0． ( 4) 式中，W = e － 2π/3 ． 构造图 2 中的滤波器组，并对其特性加以改进，使 得图 2 中的滤波器组具有半对称性并且满足近似平移 不变． 首先，提出一种构造具有线性相位特性( 即对称 性) 的低通滤波器的方法． 其次，在此基础上利用谱分 解求出对应的高通滤波器． 本文对已有方法进行了改 进，使用拓普利兹矩阵分解方法，计算方法简单，因此 有效地提高了计算效率． 2 滤波器的构造 2. 1 对称低通滤波器的构造 离散多项式信号由低通滤波器 H( z) 处理之后，通 常情况下应该能够保留离散多项式信号中的阶次小于 或等于 N － 1 的部分，而阶次大于 N 的部分可以由高 通滤波器 Gi ( z) ( i = 0，1，2) 消除，即离散多项式信号 中阶次为 Ki ( i = 0，1，2) 的部分． 因此，严格采样正交 小波系统应该满足 N = Ki ( i = 0，1，2) 的条件，而对于 分数阶过完备小波系统，为了获得更高阶的消失矩，通 常要求 N≥Ki ( i = 0，1，2) ． 基于此条件和图 2 中的滤 波器组结构，本文中的 3 /2 阶过完备小波低通滤波器 和高通滤波器组可以写成如下形式: H( z) = ( 1 + z － 1 ) N ( 1 + z － 1 + z － 2 ) N Q( z) ， ( 5) Gi ( z) = ( 1 － z － 1 ) Ki Vi ( z) ，i = 0，1，2． ( 6) 式中，N 是低通滤波器在点 z = － 1、z = ei2π/3 和 z = ei4π/3 处的零点个数，因此可以保留多项式中阶次小于等于 N － 1 部分; Ki ( i = 0，1，2) 为高通滤波器在点 z = 1 处零 点个数，且 N≥Ki ( i = 0，1，2) ，并有( 1 + z － 1 ) ( 1 + z － 1 + z － 2 ) 不能整除 Q( z) ． 如果要得到冲击响应最短的低通滤波器，分数阶 过完备小波低通滤波器长度还应该满足下面的条件: LH≥3N + 2K － 1． ( 7) 式中，K = min( K0，K1，K2 ) ． 因此，由式( 5) 和( 7) 可知，为了得到阶次( 滤波器 最短) 最低并且具有线性相位特性的低通滤波器 H( z) ，式( 5) 中的 Q( z) 因子应该是对称的，并且长度 应该小于 2K － 1; 也就是低通滤波器的长度应该至少为 3N + 2K － 1，从而使得其最短支集长度为 3N + 2K － 1． 下面给出具有线性相位特性以及最小长度的低通 滤波器的构造方法． 此处考虑 N 为奇数的情况，N 为 偶数时同理可得． 首先对于给定的 N 和 K，假设低通滤波器 H( z) 是 对称的并且其冲激响应 h( n) 长度为 L，其中 L≥3N + 2K － 1． 令 z = eiω ，那么 H( z) 可以表示成 ω 的多项式． 为了求解 H( z) ，首先进行下述转换: x = 1 4 ( z － 1 + 2 + z) = sin2 ω 2 = 1 2 ( 1 － cosω) ． ( 8) 进而有 1 － 4 3 x = 1 3 ( z － 1 + 1 + z) = 1 3 z( 1 + z － 1 + z － 2 ) ，( 9) 1 － x = 1 4 ( z － 1 + 2 + z) = 1 4 ( 1 + z － 1 ) ( 1 + z) ． ( 10) 因子 1 + z － 1可以用 x 表示: 1 + z － 1 = 2cos ω 2 = 2 1 － 槡 x． ( 11) 在进行一位平移后，1 + z － 1 + z － 2也可以用 x 表示: z － 1 + 1 + z ( = 4 3 4 － ) x ． ( 12) 将 H( z) 转换 成 x 的 表 达 式 T ( x) ，其中 A ( x) 对 应 Q( z) ，有 T( x) = ( 1 － x) N/ ( 2 3 4 － ) x N A( x) ． ( 13) 由 z 和 x 的零点对应关系，可以得到 T( x) 有下列等价 关系成立: T( 0) = 1; T( i) ( 0) = 0，( i = 1，2，…，K － 1) ; T( i) ( 1) = 0，( i = 0，1 /2，…，N /2 － 1) ; T( i) ( 3 /4) = 0，( i = 0，1，…，N － 1) ． 考虑 H( z) 在 z = 1 处的特性，即 T( x) 在 x = 0 处的特 性，有下式成立: 1 － T( x) = B( x) xK ． ( 14) 综合式( 13) 和( 14) ，得到 ( 1 － x) N/ ( 2 3 4 － ) x N A( x) + B( x) xK = 1， ( 15) · 083 ·

沈政伟等：分数阶对称平移不变过完备小波的构造及其在轴承故障诊断中的应用 381· 则A(x)可以表示为 构条件的等价多相位矩阵表达式： A(x）=- 1-B(x)x* (16) [H(z) Gg81-4 (20) 式中 A(x)可以通过泰勒级数展开式得到： rHo(z） Ho(2)Ho.(2)1 =A(x)+0(1x).(17) H(）= 1-)(-) [Ho(2) H,(e)H2(a)] TGm(z） Go(=)Go (= 展开式子 G(z)= Gio(z) Gu(2) G2(z） LGm(a)G(z)G2(a)」 通过式(20)得到 和 G"(a-)G(z)=I,-H(2)H(a). (21) 下面将通过矩阵谱因式分解的方式求解上式中的 G(z).对于矩阵谱因式分解，文献3]给出了多种方 此处的 法的详细介绍，但这些方法仅限于简单例子的手工计 N-1 ]为组合运算，取泰勒展开式中的前 算.文献2]中使用的是对称因式分解方法，但计算结 K项，即可得到A(x).再通过标准化(H(1)=√6求 果还需要最后通过仿酉矩阵调整才能使高通滤波器组 出T(x),最后经过参数转换即得到如下结果： 满足近似平移不变特性，因此计算复杂度非常高.本 文将采用拓普利兹矩阵分解法对式(21)进行直接分 解，该方法是一种迭代计算方法，计算方法非常简单. q*2 -g+N-1 通过程序计算对比发现，本文采用拓普利兹矩阵 -1 N-1 分解法求高通滤波器组具有以下三点优势：(1)对称 因式分解法回和使用拓普利兹矩阵分解法迭代5000 (告)”(2 (18) 次计算时间开销虽然接近，但是本文方法在经过3000 次拓普利兹矩阵分解迭代运算（时间为5000次迭代 式中，z-”因子是用来调整平移量.同理当V为偶 的1/2)之后，其计算结果误差范围已经在0.0055%~ 数时，也可以得到相应的表达式： 0.0097%.(2)利用拓普利兹矩阵分解法得到的高通 滤波器组已经满足近似平移不变特性，无需再如文献 ]一样做额外的计算，因此省去繁杂的调整过程，降 +22 低了计算复杂度.(3)根据不同的数据精度需求和滤 p-q+N-1Y N-1 波器的设计复杂程度，采用拓普利兹矩阵分解法也具 2-2 有更好的灵活度，同时也更方便计算机实现. 表1是文献2]N=3和K=1高通滤波器组.对 (19) 比表2中由本文中提出的分解法经过3000次迭代运 因此给定N和K,通过式(18)和式(19)可以很容 算之后的结果，可以看出两种方法求出的高通滤波器 易求出相应的低通滤波器。需要强调说明的是，利用 组是近似一样的.这验证了本文方法的正确性. 式(18)和式(19)构造的分数阶低通滤波器是具有线 表1文献2]中的N=3和K=1的高通滤波器系数 性相位特性的（对称性）. Table 1 High-pass filter coefficients of N=3 and K=1 in Ref. 2.2高通滤波器的构造 50(n) &(n) &2(n) 构造了具有线性相位特性的低通滤波器之后，将 0.64917778505741 0 求解相应的高通滤波器，从而完成图2中滤波器组的 0 -0.482626543662260.63770868747435 0 构造.首先，将低通滤波器H(z)用多相位方法表示为 H(z)=H。()+z3H,(),则低通滤波器部分被分解 -0.15059130119969-0.466878032128120.64520631583316 成两个低通通道H。(z)和H,(z).为方便使用完全可 -0.01477135217528-0.14815175304968-0.49098922627425 重构条件(2)~(4)，再将H。(z)展开成多相位形式， -0.00118858802016-0.02267890229656-0.13149490989732 即H(z)=Ho(z3)+zHo(z3)+z2H2(z),H,(z）和 0 0 -0.02152627546889 G,(z)(i=0,1,2)按照同样的方式分解，得到完全可重 0 0 -0.00119590419272

沈政伟等: 分数阶对称平移不变过完备小波的构造及其在轴承故障诊断中的应用 则 A( x) 可以表示为 A( x) = 1 － B( x) xK ( 1 － x) N/ ( 2 3 4 － ) x N ． ( 16) A( x) 可以通过泰勒级数展开式得到: 1 ( 1 － x) N/ ( 2 3 4 － ) x N = A( x) + O( | x | K ) ． ( 17) 展开式子 1 ( 1 － x) N = ∑ ∞ n = 0 n + N － 1 ( N ) － 1 xn 和 1 ( 3 /4 － x) N = ∑ ∞ n = 0 ( 3 /4) N n + N － 1 ( N ) － 1 ( 3 /4) n xn ， 此处的 n + N － 1 ( N ) － 1 为组合运算，取泰勒展开式中的前 K 项，即可得到 A( x) ． 再通过标准化( H( 1) = 6槡) 求 出 T( x) ，最后经过参数转换即得到如下结果: H( z) = z － ( K － 1) 槡 ( 6 1 + z － 1 ) 2 ( N 1 + z － 1 + z － 2 ) 3 N · ∑ K －1 p = 0 ∑ p q = 0 q + N 2 － 1 N 2           － 1 p － q + N － 1 ( N ) － 1 ( · ) 4 3 p － ( q － z － 1 + 2 － z ) 4 p ． ( 18) 式中，z － ( K － 1) 因子是用来调整平移量． 同理当 N 为偶 数时，也可以得到相应的表达式: H( z) = z － ( K － 1) 槡 ( 6 1 + z － 1 ) 2 N ( + 1 1 + z － 1 + z － 2 ) 3 N · ∑ K －1 p = 0 ∑ p q = 0 q + N 2 － 1 2 N 2 －           1 2 p － q + N － 1 ( N ) － 1 ( · ) 4 3 p － ( q － z － 1 + 2 － z ) 4 p ． ( 19) 因此给定 N 和 K，通过式( 18) 和式( 19) 可以很容 易求出相应的低通滤波器． 需要强调说明的是，利用 式( 18) 和式( 19) 构造的分数阶低通滤波器是具有线 性相位特性的( 对称性) ． 2. 2 高通滤波器的构造 构造了具有线性相位特性的低通滤波器之后，将 求解相应的高通滤波器，从而完成图 2 中滤波器组的 构造． 首先，将低通滤波器 H( z) 用多相位方法表示为 H( z) = H0 ( z 2 ) + z － 3H1 ( z 2 ) ，则低通滤波器部分被分解 成两个低通通道 H0 ( z) 和 H1 ( z) ． 为方便使用完全可 重构条件( 2) ～ ( 4) ，再将 H0 ( z) 展开成多相位形式， 即 H0 ( z) = H00 ( z 3 ) + z － 1H01 ( z 3 ) + z － 2H02 ( z 3 ) ，H1 ( z) 和 Gi ( z) ( i = 0，1，2) 按照同样的方式分解，得到完全可重 构条件的等价多相位矩阵表达式: HT ( z － 1 ) GT ( z － 1 [ ) ] H( z) G( z [ ] ) = I3 ． ( 20) 式中 H( z) = H00 ( z) H01 ( z) H02 ( z) H10 ( z) H11 ( z) H12 [ ( z) ]， G( z) = G00 ( z) G01 ( z) G02 ( z) G10 ( z) G11 ( z) G12 ( z) G20 ( z) G21 ( z) G22 ( z        )  ． 通过式( 20) 得到 GT ( z － 1 ) G( z) = I3 － HT ( z － 1 ) H( z) ． ( 21) 下面将通过矩阵谱因式分解的方式求解上式中的 G( z) ． 对于矩阵谱因式分解，文献［13］给出了多种方 法的详细介绍，但这些方法仅限于简单例子的手工计 算． 文献［2］中使用的是对称因式分解方法，但计算结 果还需要最后通过仿酉矩阵调整才能使高通滤波器组 满足近似平移不变特性，因此计算复杂度非常高． 本 文将采用拓普利兹矩阵分解法对式( 21) 进行直接分 解，该方法是一种迭代计算方法，计算方法非常简单． 通过程序计算对比发现，本文采用拓普利兹矩阵 分解法求高通滤波器组具有以下三点优势: ( 1) 对称 因式分解法［2］和使用拓普利兹矩阵分解法迭代 5 000 次计算时间开销虽然接近，但是本文方法在经过 3 000 次拓普利兹矩阵分解迭代运算( 时间为 5 000 次迭代 的 1 /2) 之后，其计算结果误差范围已经在 0. 005 5% ～ 0. 009 7% ． ( 2) 利用拓普利兹矩阵分解法得到的高通 滤波器组已经满足近似平移不变特性，无需再如文献 ［2］一样做额外的计算，因此省去繁杂的调整过程，降 低了计算复杂度． ( 3) 根据不同的数据精度需求和滤 波器的设计复杂程度，采用拓普利兹矩阵分解法也具 有更好的灵活度，同时也更方便计算机实现． 表 1 是文献［2］N = 3 和 K = 1 高通滤波器组． 对 比表 2 中由本文中提出的分解法经过 3 000 次迭代运 算之后的结果，可以看出两种方法求出的高通滤波器 组是近似一样的． 这验证了本文方法的正确性． 表 1 文献［2］中的 N = 3 和 K = 1 的高通滤波器系数 Table 1 High-pass filter coefficients of N = 3 and K = 1 in Ｒef．［2］ g0 ( n) g1 ( n) g2 ( n) 0. 649 177 785 057 41 0 0 － 0. 482 626 543 662 26 0. 637 708 687 474 35 0 － 0. 150 591 301 199 69 － 0. 466 878 032 128 12 0. 645 206 315 833 16 － 0. 014 771 352 175 28 － 0. 148 151 753 049 68 － 0. 490 989 226 274 25 － 0. 001 188 588 020 16 － 0. 022 678 902 296 56 － 0. 131 494 909 897 32 0 0 － 0. 021 526 275 468 89 0 0 － 0. 001 195 904 192 72 · 183 ·

·382· 工程科学学报，第37卷，第3期 表2通过拓普利兹矩阵分解得到的N=3和K=1的高通滤波器系数 求出 Table 2 High-pass filter coefficients of N=3 and K=1 obtained by To- 例1N=3和K=1的低通滤波器： eplitz matrix decomposition B0(n) &i(n) &(n) o=6广 0.649214180021336 0 0 例2N=4,K=2的低通滤波器： -0.4825822417505970.637744293143371 0 -0.150580473087161 -0.466832735837093 0.645241656686216 -' -0.014770391219193-0.148141693288885 -0.490944166453673 (-47z2+118z1-47). -0.001188521365662 -0.022677502664602 -0.131485064578830 例3N=5和K=3的低通滤波器： 0 0 -0.021525096049414 o源告字) 0 0 -0.001195838647559 (3435z4-16380z-3+27042z-2+16380z1+3435). 下面给出三个具体构造例子.其中低通滤波器用 具体的底通和高通滤波器组系数在此省略，本文 式(18)或者式(19)求出，对应的高通滤波器组则可以 中只画出了上述三个例子的对应高通滤波器对应的小 利用本节中的拓普利兹矩阵分解法，即用式(21) 波函数图像，如图3所示 0.10 0.2 (a) --平 0.05- 01 0 -0.05- -0.1 0106 -0. 50 100 20 70 120 时间 时间 0.2 一2 01 0 -0.1 -0. 40 90 140 时间 图3对称分数阶过完备小波.(a)N=3,K=1:(b)N=4,K=2:(c)N=5,K=3 Fig.3 Symmetrical fractional overcomplete wavelet:(a)N=3,K=1;(b)N=4,K=2:(c)N=5,K=3 3在轴承故障诊断中的应用 74rmin和582r"min.数据采样来自于风机左侧 的1号轴承，采样频率为25600Hz,样本数据长度为 3.1风机系统滚动轴承故障信号检测 16384.通过计算可以得出外圈的故障频率为 图4为某风机系统中滚动轴承发生故障时的信号 63.96Hz. (含有噪声)时域波形图.电机设备转速和性能分别为 为了对比构造的具有线性相位特性和近似平移不 740r·min和1200kW,风机两侧的滚动轴承型号为 变性的过完备小波在故障诊断中的优势，选择d5作 22344CA,轴承滚柱的数量为13，输入和输出转速为 为对比小波.原因有两点：第一，dh5小波不具有线性

工程科学学报，第 37 卷，第 3 期 表2 通过拓普利兹矩阵分解得到的 N = 3 和 K = 1 的高通滤波器系数 Table 2 High-pass filter coefficients of N = 3 and K = 1 obtained by To￾eplitz matrix decomposition g0 ( n) g1 ( n) g2 ( n) 0. 649 214 180 021 336 0 0 － 0. 482 582 241 750 597 0. 637 744 293 143 371 0 － 0. 150 580 473 087 161 － 0. 466 832 735 837 093 0. 645 241 656 686 216 － 0. 014 770 391 219 193 － 0. 148 141 693 288 885 － 0. 490 944 166 453 673 － 0. 001 188 521 365 662 － 0. 022 677 502 664 602 － 0. 131 485 064 578 830 0 0 － 0. 021 525 096 049 414 0 0 － 0. 001 195 838 647 559 下面给出三个具体构造例子． 其中低通滤波器用 式( 18) 或者式( 19) 求出，对应的高通滤波器组则可以 利用本 节 中 的 拓 普 利 兹 矩 阵 分 解 法，即 用 式 ( 21 ) 求出． 例 1 N = 3 和 K = 1 的低通滤波器: H( z) 槡 ( = 6 1 + z － 1 ) 2 ( 3 1 + z － 1 + z － 2 ) 3 3 ． 例 2 N = 4，K = 2 的低通滤波器: H( z) =槡6 ( 24 1 + z － 1 ) 2 ( 5 1 + z － 1 + z － 2 ) 3 4 · ( － 47z － 2 + 118z － 1 － 47) ． 例 3 N = 5 和 K = 3 的低通滤波器: H( z) = 槡6 ( 1 152 1 + z － 1 ) 2 ( 5 1 + z － 1 + z － 2 ) 3 5 · ( 3 435z － 4 － 16 380z － 3 + 27 042z － 2 + 16 380z － 1 + 3 435) ． 具体的底通和高通滤波器组系数在此省略，本文 中只画出了上述三个例子的对应高通滤波器对应的小 波函数图像，如图 3 所示． 图 3 对称分数阶过完备小波． ( a) N = 3，K = 1; ( b) N = 4，K = 2; ( c) N = 5，K = 3 Fig． 3 Symmetrical fractional overcomplete wavelet: ( a) N = 3，K = 1; ( b) N = 4，K = 2; ( c) N = 5，K = 3 3 在轴承故障诊断中的应用 3. 1 风机系统滚动轴承故障信号检测 图 4 为某风机系统中滚动轴承发生故障时的信号 ( 含有噪声) 时域波形图． 电机设备转速和性能分别为 740 r·min － 1和 1 200 kW，风机两侧的滚动轴承型号为 22344CA，轴承滚柱的数量为 13，输入和输出转速为 74 r·min － 1和 582 r·min － 1 ． 数据采样来自于风机左侧 的 1 号轴承，采样频率为 25 600 Hz，样本数据长度为 16 384． 通 过 计 算 可 以 得 出 外 圈 的 故 障 频 率 为 63. 96 Hz． 为了对比构造的具有线性相位特性和近似平移不 变性的过完备小波在故障诊断中的优势，选择 db5 作 为对比小波． 原因有两点: 第一，db5 小波不具有线性 · 283 ·