基于连续小波变换的旋转机械部件的自适应故障诊断方法

技术领域

本发明涉及机械故障诊断技术领域，尤其涉及对旋转机械部件(如轴承或齿轮等)的故 障诊断技术。

背景技术

旋转机械部件(如滚动轴承)发生局部故障(如剥落、点蚀、裂缝等)时，其振动信 号通常伴随着冲击的产生，这些低频段的低频冲击信号通过与机械系统和传感器的激励在 高频段发生共振产生高频振动信号。因此，通过对高频振动信号的分析能够有效对轴承进 行故障诊断，同时也消除了低频噪声的干扰。在过去很长一段时间内，高频共振解调技术 (又称为包络解调技术)成为旋转机械部件故障诊断的有效方法，它能够较准确的获取旋 转机械部件的故障特征频率。

在实际的故障诊断中，高频共振解调技术往往要求准确的获取高频共振频带，也即确 定带通滤波器的中心频率及带宽，然而中心频率与带宽的确定往往是比较困难的，这是因 为高频共振解调技术中高频共振频带的中心频率通常是事先未知的，这就导致了高频共振 频带的获取不完全或者带宽过大，从而影响到最终故障诊断的准确性。

发明内容

本发明的目的是针对现有的高频共振解调技术中带通滤波器中心频率难以确定的问 题，提出了基于连续小波变换的旋转机械部件的自适应故障诊断方法，该方法尤其适合带 有冲击性损伤的旋转机械部件的故障诊断。

本发明的内容是：基于连续小波变换的旋转机械部件的自适应故障诊断方法，如图1 所示，其步骤包括：

步骤1：对获取的离散初始振动信号s(n)进行零均值预处理，得到消除直流分量的预处 理信号s₁(n)；

步骤2：将步骤1中所得到的预处理信号s₁(n)进行连续小波变换，得到各尺度参数对应 的小波系数WT_a(n)；

步骤3：分别计算步骤2中各尺度参数所对应的小波系数WT_a(n)的峭度；

步骤4：利用一种自适应算法找出步骤3中小波尺度参数中较大峭度所对应的小波尺度 参数，即最优小波尺度，然后选择这些尺度参数所对应的小波系数构造最优分析信号s_a(n)；

步骤5：对步骤4中得到的最优分析信号s_a(t)进行包络解调，获取其包络信号s_e(t)；

步骤6：将步骤5中的包络信号s_e(t)做FFT变换，得到其包络频谱S_e(f)，通过对包络解 调谱的分析实现对旋转机械部件的故障诊断。

本发明的有益效果是：由于不需要预先估计滤波器的中心频率和带宽，避免了估计误 差对诊断结果的影响，因此与传统的高频共振解调技术相比，本发明提高了故障诊断的准 确性，尤其适合带有冲击性损伤的机械零部件的故障诊断。

附图说明

图1是本发明的主流程图。

图2是零均值处理后的预处理信号的时域图。

图3是自适应获取最优尺度系数的流程示意图。

图4是各尺度a＝[1，2，...，32]所对应的峭度值。

图5是自适应获取尺度系数结果a_opt＝[7，8，...，16]所对应的峭度值。

图6是重构最优尺度系数得到的分析信号的时域图。

图7是分析信号解调后的包络信号时域图。

图8是分析信号的包络解调谱。

图9是图2中零均值处理后的预处理信号的包络时域图。

图10是图9中包络信号的FFT谱

具体实施方式

下面结合具体实施例和附图对本发明做进一步的说明。

如图1所示，基于连续小波变换的旋转机械部件的自适应故障诊断方法，如图1所示， 其步骤包括：

步骤1：对获取的离散初始振动信号s(n)进行零均值预处理，得到消除直流分量的预处 理信号s₁(n)。

本实施例以齿轮箱作为故障诊断的对象，该齿轮箱的振动信号通过加速度传感器进行 采集，对采集到的振动信号利用信号调理器、模数转换后得到离散初始振动信号s(n)，并将 该信号送入计算机进行零均值预处理，零均值预处理是为了消除初始振动信号s(n)中直流分 量的影响，这里对初始振动信号进行零均值处理得到消除直流分量的预处理信号s1(n)，其 具体算法如下：

$s_1\left(n\right)=s\left(n\right)-\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}s\left(n\right)---\left(1\right)$

这里s(n)为离散初始振动信号，n为离散时间点，N为信号采样点数，s₁(n)为零均值处 理后的预处理信号。

如图2所示，为零均值处理后的预处理信号的时域图形，该信号消除了直流分量。

步骤2：将步骤1中所得到的预处理信号s₁(n)进行连续小波变换，得到各尺度参数对应 的小波系数WT_a(n)；

连续小波变换定义如下：

其中x(t)为信号函数，t为时间自变量；为小波母函数，本发明中使用的小波母函数 是Morlet小波函数；a为小波尺度参数；b为小波时间参数。这里的小波系数WT(a，b)为一 个二元函数，当小波尺度参数a给定时，其对应的小波系数就变为WT_a(b)，将b离散化可表 示为WT_a(n)，n为自然数；从等式(2)中可以看出，连续小波变换相当于将信号函数x(t)与小 波函数作卷积，上标^*表示对小波母函数取共轭，根据信号处理中技术中卷积定理 的性质，连续小波变换又可表示为：

其中X(f)、ψ(f)分别为x(t)、的傅里叶变换，F^-1[□]表示傅里叶变换的逆变换，b为 逆变换后的时间参数。由等式(3)可见，连续小波变换相当于利用一组带通滤波器对信号进 行带通滤波，每一尺度a对应于一个带通滤波器即单个滤波器的带宽和中心频率 由尺度参数a决定。以Morlet小波为例：

其中为小波母函数，σ为衰减参数常量，f₀为小波母函数的中心频率，其频域 表示：

$\psi \left(f\right)=e^{-({\pi }^2/{\sigma }^2){(f-f_0)}^2}---\left(5\right)$

故其半功率带宽可计算为：；因此，小波母函数ψ(f)对应的滤波器通带为： 对应小波尺度参数a的小波带通滤波器为：

$\psi \left(\mathrm{af}\right){=e}^{-({\pi }^2/{\sigma }^2){\left(a(f-\frac{f_0}{a})\right)}^2}---\left(6\right)$

因此将步骤2中得到的预处理信号s₁(n)作连续小波变换可以得到分布在不同频带内的 小波系数WT_a(n)，这些小波系数与小波变换的尺度参数a和各频带一一对应，发明中小波变 换的尺度参数表示为a＝[a₁，a₂，...，a_l]，l为小波变换中尺度参数的个数。

步骤3：分别计算步骤2中l个尺度参数所对应的小波系数WT_a(n)的峭度(kurtosis)；

自从1970年Stewart等提出峭度(kurtosis)的概念起，过去很长一段时间里其一直 被用来衡量机械故障的严重程度。峭度ku_a是一个对冲击信号非常敏感的统计量，它的大小 衡量了信号的冲击强度，其值越大信号的冲击就越强，反之则越弱。因此，本发明利用该 统计来衡量步骤2中各尺度小波系数的冲击强弱，计算如下：

${\mathrm{ku}}_a=\frac{E\left[{\mathrm{WT}}_a{\left(n\right)}^4\right]}{{\left\{E\right[{\mathrm{WT}}_a{\left(n\right)}^2\left]\right\}}^2}---\left(7\right)$

其中WT_a(n)为步骤2中得到的小波系数，E[□]表示求期望。

步骤4：利用一种自适应算法找出步骤3中小波尺度参数a＝[a₁，a₂，...，a_l]中较大峭度ku_a所 对应的小波尺度参数a_opt＝[a_n，a_n+1，...，a_m]，即最优小波尺度，然后选择这些尺度参数所对应的 小波系数构造最优分析信号s_a(n)；

步骤2中提到利用不同尺度小波对信号做小波变换相当于利用一滤波器组对信号滤波， 而包含有轴承故障特征的信号往往峭度ku较大，因此可以选择ku值较大的几个尺度所对应 的系数来构成分析信号s_a(n)进行故障诊断，s_a(n)为具有较大峭度的小波系数的叠加，其计 算方法如下：

$s_a=\underset{a=a_m}{\overset{a_n}{\Sigma }}\mathrm{WT}(a,b)=\underset{a=a_m}{\overset{a_n}{\Sigma }}\sqrt{a}{F}^{-1}\left[X\left(f\right){\psi }^{*}\left(\mathrm{af}\right)\right]=F^{-1}\left[X\left(f\right)\underset{a=a_m}{\overset{a_n}{\Sigma }}\sqrt{a}{\psi }^{*}\left(\mathrm{af}\right)\right]=F^{-1}\left[X\left(f\right)\phi \left(f\right)\right]---\left(8\right)$

式中，φ(f)相当于一个新的滤波器，最优尺度a_opt＝[a_m，a_m+1，...，a_n]的选择利用了一种自适 应的算法，其流程如图3所示，为了能够选择出具有最大峭度ku的分析信号s_a(n)，这里将 合并后的小波系数s_a作为分析信号，并对其峭度ku进行比较，它的具体算法为：选择第i尺 度系数作为上一个系数WT_pre(n)，和第i尺度系数与第i+1尺度系数的 合并系数WT_cur(n)进行比较(i为自然数，初始值为1)，即WT_pre(n)和WT_cur(n)进行比较，如果 合并系数WT_cur(n)的ku_cur值大于上一个尺度系数WT_pre(n)的ku_pre值，那么将当前尺度系数 WT_cur(n)作为上一个尺度系数WT_pre(n)，将当前尺度系数WT_cur(n)与第i+2个尺度系数的合并 系数作为当前尺度系数WT_cur(n)，并进行比较；反之，则重新选择第i+1尺度系数作为上一 个尺度系数WT_pre(n)，并和第i+1尺度系数与第i+2尺度系数的合并系数WT_cur(n)进行比较。 以此类推，直到所有尺度系数都被合并比较完，选择出最大ku值的合并系数作为分析信号， 所合并的各组系数的尺度作为最优尺度a_opt＝[a_m，a_m+1，...，a_n]。

通过本步骤的自适应算法从尺度a＝[a₁，a₂，...，a_l]中选取具有较大峭度的小波系数WT_a(n) 所对应的小波尺度参数a_opt＝[a_n，a_n+1，...，a_m]作为最优的尺度，该最优的尺度参数的合并小波系 数s_a(n)就是我们需要的高频共振带信号，这一过程相当于对若干个小波滤波器进行筛选， 留下最可能是本发明中真实的共振频带所对应的小波滤波器，这样，在对高频共振带信号 进行提取的过程中就不必要首先确定滤波器的中心频率和带宽，从而提高了诊断的准确性。

如图4和图5所示，为本步骤选取较大峭度ku_a的一个具体实施例，图4中，横坐标为 小波变换的尺度a＝[1，2，...，32]，纵坐标为各尺度所对应的小波系数WT_a(n)的峭度ku，可以 发现在尺度为12左右，小波系数的峭度较大。通过上述自适应方法处理后，图5中展示的 是合并后的小波系数的峭度与尺度的对应关系，可以发现尺度a_opt＝[7，8，...，16]的合并系数 所对应的峭度取得最大值。因此选择尺度a_opt＝[7，8，...，16]为最优尺度，并通过式(8)构造 最优分析信号s_a(n)，最优分析信号s_a(n)的时域图如图6所示。从图6中可以看出，该最优 信号具有明显的周期冲击性，与轴承局部损伤产生的冲击相符合，实现了对高频共振带的 提取。

步骤5：对步骤4中得到的最优分析信号s_a(t)进行包络解调，获取其包络信号s_e(t)；其 具体计算方法如下：

$s_e\left(t\right)=\sqrt{s_a^2\left(t\right)+H^2\left[s_a\left(t\right)\right]}---\left(9\right)$

$H\left[s_a\left(t\right)\right]=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{s_a\left(\tau \right)}{t-\tau }\mathrm{d\tau }---\left(10\right)$

其中H[□]表示希尔伯特变换，τ为变换参数。包络解调为现有技术，这里不再详述。

如图7所示，为解调后最优分析信号s_a(t)的包络信号，其周期与轴承的内圈故障特征频 率BPFI＝257Hz相一致。

步骤6：将步骤5中的包络信号s_e(t)做FFT变换，得到其包络频谱s_e(f)，通过对包络解 调谱的分析实现对旋转机械部件(如轴承、齿轮等机械零部件)的故障诊断。

由于FFT变换为现有技术，这里不再详述。通过步骤6得到了最优分析信号的包络解 调谱，对包络解调谱进行分析，如果其包含有轴承、齿轮的故障特征频率，则表明轴承、 齿轮发生了故障；反之则运行正常。

本实施例中，我们分别给出了通过直接解调预处理信号s₁(n)的包络时域图形和频域图 形，如图9和图10所示。从图9可以看出直接解调后的包络信号虽然也具有较明显的周期 性，但是可以发现，该周期T₁＝0.0075s与轴承内圈的转频的2倍频2f_r＝133.33Hz相一致， 从图10的频域图形中可以得到验证，因此并不能反映轴承的故障特征。然从步骤6中得到 的结果(如图8所示)可以看出，本发明所得到的结果明显包含了轴承的故障特征频率 BPFI＝257Hz及其谐波频率(2BPFI、3BPFI、…)。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的 原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通 技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体 变形和组合，这些变形和组合仍然在本发明的保护范围内。