CN102706555B - 一种复解析最优小波解调法 - Google Patents

Abstract

一种复解析最优小波解调法，其特征在于，包括如下步骤：获取振动信号频谱图，截取有效频带集中区域；设定Morlet小波的外形因子σ和初始频率ω₀的初始值，得到小波系数C_s(b，a)；得出连续小波幅值谱信息测度SH；小波峭度改进；计算平滑指数；计算最优小波包络谱；解调谱细化，对检测目标的故障定位。本发明不仅对分析小波参数进行了优化选择，而且对常规的小波峭度计算方法进行了改进，提高了小波系数筛选的准确性。同时，为了解决早期故障信号分解出的小波系数特征不明显的问题，增设了平滑指数，对故障频带所在的小波系数起到了更好的定位筛选作用。该方法不仅能更加准确有效地在多分量信号中提出复杂信号中微弱故障特征，具有很高的实用价值。

Description

一种复解析最优小波解调法

技术领域

本发明涉及机械工程信号处理领域，尤其是旋转机械振动信号的故障特征提取方法，主要是基于小波参数优化的一种复解析最优小波解调方法，可应用于在诸如齿轮箱、轴承等旋转机械的故障定位和故障程度的分析中。

背景技术

机械传动机构在运转过程中，若传动部件存在局部损伤，如齿轮箱，若某个齿轮存在局部损伤，会产生间断性冲击激振，使得在整个机械设备上测得的振动信号包含有迅速衰减的脉冲分量。

这种脉冲信号通常在时域上呈一定周期间隔出现，包含这种冲击的信号频谱中往往会出现以啮合频率或固有频率为中心的两侧间隔均匀的调制边频。凭借调制频率和幅值强度可以判断零件损伤的程度和部位，因此解调分析是齿轮箱等机械传动机构的状态监测与故障诊断的一种重要方法。

从齿轮箱表面采集的振动信号包含了多级传动齿轮及转轴与轴承等多个零部件振动信息，因此具有多载波多调制信息的特征，信号整体包络解调无法得出有效的故障信息，通常必须先滤波再进行解调分析。使用传统单一带通滤波器由于无法准确定位具体的故障特征信号所在频段，很难准确地提取出隐含故障信息所在的调制频带。国内外许多学者采用小波变换方法实现齿轮箱信号频带的任意划分，然后根据齿轮箱的理论啮合频率及倍频大致确定故障频带，这种方法虽然在单级齿轮传动系统或结构较简单的齿轮箱系统，且信噪比高的情况下是有效的。但对于工况复杂的多级耦合的齿轮箱系统来，这种方法往往会因为其他频率干扰而使齿轮啮合成分无法分辨，无法有效地提取出有用的故障调制信息。

发明内容

本发明的目的就是提供一种复解析最优小波解调方法，能更优化的调制信息提取方法，从而克服传统技术的缺陷，有效的解决多载波多调制，频段混叠的复杂信号的调制信息提取问题，以方便对齿轮箱等故障位置和故障程度做出准确的判断，为相关的机械设备状态辨识与故障诊断研究提供可靠的保障。

本发明的思想在于：针对齿轮故障信号特点，采用解析小波变换对信号分量的包络幅值进行提取，应用最小化时间尺度熵的方法对分析小波参数进行了优化选择，将修正峭度和平滑指数作为故障的度量指标对故障频带进行准确定位，继而又对解调谱进行了细化分析以提高复解析小波解调谱的频谱精度，因此，将各种优化方法相融合，共同应用于复解析小波变换中，可以使齿轮箱的故障调制信息更加准确地提取出来。

改进的小波参数的优化选择方式，就是采用了时频信息熵、小波峭度及平滑指数同时对小波中心频率、衰减参数及变换尺度进行了优化。

本发明的技术方案是：一种复解析最优小波解调法，其关键在于，包括如下步骤：

步骤一、通过传感器获取检测目标的振动信号时域波形，并对其进行傅里叶变换可得到振动信号频谱图，确定信号分析频率的上下界限f_xmin，f_xmax；

步骤二、截取所述振动信号频谱图中信号的有效频带集中区域；

步骤三、尺度范围内的解析小波变换：设定Morlet小波的外形因子σ和初始频率ω₀的初始值，Morlet小波函数的表达式为：

结合分析频率的上下界限f_xmin，f_xmax，根据小波变换的尺度a与基小波采样率、中心频率和信号分析频率关系确定出a的取值范围E：

$E=(\frac{{\mathrm{\ω}}_0{f}_s}{{2\mathrm{\πf}}_{x\max }{f}_{\mathrm{\ω}}},\frac{{\mathrm{\ω}}_0{f}_s}{{2\mathrm{\πf}}_{x\min }{f}_{\mathrm{\ω}}})$

其中：f_s为信号采样率；f_ω为小波采样率；

采用以上设置的具体参数值，在尺度范围E内，对被分析信号s(t)作关于小波函数的连续小波变换，得到小波系数C_s(b，a)，其中b是位置参数；

步骤四、小波参数迭代确定最佳小波基：由小波变换系数C_s(b，a)可得出在尺度a下被分析信号s(t)包络幅值A_S(b，a)＝|C_s(b，a)|，进一步可以得出连续小波幅值谱信息测度SH：

$\mathrm{SH}=-\underset{a\∈E}{\mathrm{\Σ}}\underset{b=i}{\mathrm{\Σ}}\left(\frac{{\left|A_s(b,a)\right|}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left|A_s(b,a)\right|}^2}\right)\lg \left(\frac{{\left|A_s(b,a)\right|}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left|A_s(b,a)\right|}^2}\right)$

其中N为时间的离散序列长度；i为小波系数离散序列号；对ω₀与σ进行迭代搜索：先选用ω₀的初始值，将σ从初始值开始增加一个步幅，求得小波系数与小波幅值谱信息测度SH，直至完成σ的搜索范围；然后将ω₀加上一个步幅，σ又从初始值开始增加一个步幅，计算小波系数与小波幅值谱信息测度SH，直至完成σ的搜索范围，重复前述两个步骤直至完成ω₀的搜索范围；通过前面的迭代搜索可确定出SH的最小值，其所对应的就是最优ω₀和σ，将其带入Morlet小波函数的表达式即可得到最佳小波函数；

步骤五、小波峭度改进：得到最佳小波函数后，在设定的尺度范围E内完成连续小波变换，为加大噪声和冲击之间的峭度差距，将小波峭度改进为小波系数的能量包络幅值的峭度K_sopt(a)：

$K_{\mathrm{sopt}}\left(a\right)=\frac{C_{4W}^{\′}\left(a\right)}{{(C_{2W}^{\′}\left(a\right)+\mathrm{\η}\max [C_{2W}^{\′}\left]\right)}^2}-3$

其中：C′_4W(a)表示小波包络幅值A_s(b，a)的四阶矩；C′_2W(a)表示A_s(b，a)的二阶矩；max[C′_2W]表述所有A_s(b，a)二阶矩的最大值；η为限制二阶矩幅度的常数，取决于最小二阶矩到最大二阶矩的动态波动范围；

搜索出最大小波峭度值，继而计算出每个尺度上的改进小波峭度K_sopt(a)；

步骤六、计算平滑指数：为解决早期故障信号峭度本身较小的问题，引入平滑指数对小波系数作进一步的筛选；对尺度范围E内的小波系数包络幅值计算平滑指数SI(a)：

$\mathrm{SI}\left(a\right)=\frac{{\left(\underset{i=1}{\overset{N}{\mathrm{\Π}}}{A}_s(b_i,a)\right)}^{1/N}}{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_s(b_i,a)}$

其中b_i为小波系数的时间离散序列；N为时间的离散序列长度；i为小波系数离散序列号；A_s(b_i，a)为不同时间离散序列下对应的包络幅值；

步骤七、计算最优小波包络谱：对尺度范围E内的与SI(a)均做归一化处理，得到小波系数筛选指标：

KSI(a)＝Normalized(K_sopt(a))-Normalized(SI(a))，

其中Normalized(K_sopt(a))是不同尺度a下的峭度；Normalized(SI(a))是不同尺度a下的平滑指数；按照KSI(a)最大原则对对小波系数进行筛选，即，在尺度范围E内，最大KSI(a)所对应的小波系数为最优小波系数，继而计算出最优小波包络幅值谱；

步骤八、解调谱细化：对最优小波包络幅值谱进行傅里叶变换得到最优小波系数的解调谱，为了提高包络谱线的分辨率，精确地诊断出缺陷类型，对解调谱进行细化分析，得到最优小波系数解调细化谱，然后根据调制频率，即边频带的间隔与检测目标故障所在的转频的对应关系，对检测目标的故障精确定位。

步骤八中所述，检测目标里每个部件的转频不同，所述细化分析就是对解调谱进行低频段细化，对包络信号进行重新抽样，即每取10个点抽取一个点重新组成新的数据序列，在对新的数据序列进行频谱分析，得到细化后的解调谱。

本发明的显著效果是：本发明结合基于复解析小波变换的信号分量的包络幅值提取原理，提出了一种基于峭度和平滑指数的复解析最优小波解调方法。这种方法不仅对分析小波参数进行了优化选择，而且对常规的小波峭度计算方法进行了改进，提高了小波系数筛选的准确性。

同时，为了解决齿轮箱早期故障信号分解出的小波系数特征不明显的问题，增设了平滑指数，对故障频带所在的小波系数起到了更好的定位筛选作用。为解决调制频带密集繁杂无法分辨问题，经频谱细化方法引入到了解调谱的细化中，能更准确地捕捉到信号的故障信息。该方法不仅能更加准确有效地提出齿轮箱复杂信号中微弱故障特征，在其他多分量信号的调制信息提取中也具有很高的实用价值。

附图说明

图1：本发明方法的基本流程图

图2：使用本方法完成的齿轮故障诊断的实施例波形组合

2a齿轮箱振动信号时域波形

2b齿轮箱振动信号频谱图

2c齿轮箱振动信号连续小波幅值谱

2d最优小波连续变换谱

2e最优小波系数的包络幅值谱

2f最优小波系数的解调谱

2g最优小波系数解调细化谱

具体实施方式：

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

用于实施的硬件环境是：AMD Athlon(tm)2.60G计算机、2.00G内存、128M显卡，运行的软件环境是：Labview 8.6和Windows XP。我们用Labview软件实现了本发明提出的方法。被分析信号来自于试验台上存在故障的齿轮变速箱。我们通过安装在箱体表面上的加速度传感器，用30kHz的采样频率为测得该变速箱在II档运行时的振动信号，为去除高次谐波的影响，模拟滤波器上限截止频率为3000Hz，然后利用本发明提出的方法对振动信号进行分析以诊断其故障所在。运行啮合齿轮的齿数比为41∶37，大齿轮所在轴的转速为600rpm，计算出齿轮的啮合频率为410Hz，大齿轮所在轴的转频为10.0Hz，小齿轮所在轴的转频为11.081Hz。

如图1所示：一种复解析最优小波解调法，包括如下步骤：

步骤一、通过传感器获取检测目标的振动信号时域波形，如图2a所示，并对其进行傅里叶变换可得到振动信号频谱图，如图2b所示，从图2b中可以看出，信号的有效频带大致集中在[750，2250]之间，确定信号分析频率的上下界限f_xmin＝750Hz，f_xmax＝2250Hz；

步骤二、截取所述振动信号频谱图中信号的有效频带集中区域[750，2250]；

步骤三、设定Morlet小波的外形因子σ和初始频率ω₀的初始值，Morlet小波函数的初始参数初始值选取越小越好，但是计算量也会增加，为了计算方便，将Morlet小波函数的初始参数选为σ＝10，ω₀＝20；

Morlet小波函数的表达式为：

结合分析频率的上下界限f_xmin，f_xmax，根据小波变换的尺度a与基小波采样率、中心频率和信号分析频率关系确定出a的取值范围E＝[2，16]：

$E=(\frac{{\mathrm{\ω}}_0{f}_s}{{2\mathrm{\πf}}_{x\max }{f}_{\mathrm{\ω}}},\frac{{\mathrm{\ω}}_0{f}_s}{{2\mathrm{\πf}}_{x\min }{f}_{\mathrm{\ω}}})$

其中：f_s为信号采样率；f_ω为小波采样率；

采用以上设置的具体参数值，在尺度范围E内，对被分析信号s(t)作关于小波函数的连续小波变换，得到小波系数C_s(b，a)，其中b是位置参数；

步骤四、小波参数迭代确定最佳小波基：由小波变换系数C_s(b，a)可得出在尺度a下被分析信号s(t)包络幅值A_S(b，a)＝|C_s(b，a)|，进一步可以得出连续小波幅值谱信息测度SH：

$\mathrm{SH}=-\underset{a\∈E}{\mathrm{\Σ}}\underset{b=i}{\mathrm{\Σ}}\left(\frac{{\left|A_s(b,a)\right|}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left|A_s(b,a)\right|}^2}\right)\lg \left(\frac{{\left|A_s(b,a)\right|}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left|A_s(b,a)\right|}^2}\right)$

其中N为时间的离散序列长度；i为小波系数离散序列号；对ω₀与σ进行迭代搜索：先选用ω₀的初始值，将σ从初始值开始增加一个步幅，求得小波系数与小波幅值谱信息测度SH，直至完成σ的搜索范围；然后将ω₀加上一个步幅，σ又从初始值开始增加一个步幅，计算小波系数与小波幅值谱信息测度SH，直至完成σ的搜索范围，重复前述两个步骤直至完成ω₀的搜索范围；通过前面的迭代搜索可确定出SH的最小值，其所对应的就是最优ω₀和σ，将其带入Morl et小波函数的表达式即可得到最佳小波函数；

如图2c所示，结合初始的连续小波复值谱，虽然搜索步长越小越好，但是搜索步长越小，计算量就会增大，所以在计算速度容许的情况下，将步长舍得适中即可，因而此处设定小波参数σ和ω₀的搜索步长分别为10，20；调节小波中心频率参数ω₀和衰减参数σ，经过迭代搜索，当SH达到最小值时，可确定最优化小波参数：σ＝470，ω₀＝7380，将其带入Morlet小波函数的表达式即可得到最佳小波函数，采用该最优小波函数得到的连续小波变换谱如图2d所示，比较图2d与图2c可知，经过优化后的小波连续谱的时间尺度分辨率得到了较大的提高，冲击信号的特征刻画得更加清楚；

步骤五、小波峭度改进：得到最佳小波函数后，在设定的尺度范围E内完成连续小波变换，为加大噪声和冲击之间的峭度差距，将小波峭度改进为小波系数的能量包络幅值的峭度K_sopt(a)：

$K_{\mathrm{sopt}}\left(a\right)=\frac{C_{4W}^{\′}\left(a\right)}{{(C_{2W}^{\′}\left(a\right)+\mathrm{\η}\max [C_{2W}^{\′}\left]\right)}^2}-3$

其中：C′_4W(a)表示小波包络幅值A_s(b，a)的四阶矩；C′_2W(a)表示A_s(b，a)的二阶矩；max[C′_2W]表述所有A_s(b，a)二阶矩的最大值；η为限制二阶矩幅度的常数，取决于最小二阶矩到最大二阶矩的动态波动范围；

搜索出最大小波峭度值，继而计算出每个尺度上的改进小波峭度K_sopt(a)；

步骤六、计算平滑指数：为解决早期故障信号峭度本身较小的问题，引入平滑指数对小波系数作进一步的筛选；对尺度范围E内的小波系数包络幅值计算平滑指数SI(a)：

$\mathrm{SI}\left(a\right)=\frac{{\left(\underset{i=1}{\overset{N}{\mathrm{\Π}}}{A}_s(b_i,a)\right)}^{1/N}}{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_s(b_i,a)}$

其中b_i为小波系数的时间离散序列；N为时间的离散序列长度；i为小波系数离散序列号；A_s(b_i，a)为不同时间离散序列下对应的包络幅值；

步骤七、计算最优小波包络谱：对尺度范围E内的与SI(a)均做归一化处理，得到小波系数筛选指标：

在图2d的基础上，设定尺度a的搜索步长为0.1，按照上述步骤五、六计算每个搜索尺度下的小波系数包络幅值的峭度和平滑指数，并计算小波系数筛选指数KSI(a)，在尺度范围[2，16]内，最大KSI(a)对应的尺度为a＝6.1，对应的最优小波系数的包络幅值谱如图2e所示；

KSI(a)＝Normalized(K_sopt(a))-Normalized(SI(a))，

其中Normalized(K_sopt(a))是不同尺度a下的峭度；Normalized(SI(a))是不同尺度a下的平滑指数；按照KSI(a)最大原则对对小波系数进行筛选，即，在尺度范围E内，最大KSI(a)所对应的小波系数为最优小波系数，继而计算出最优小波包络幅值谱；

步骤八、解调谱细化：对最优小波包络幅值谱进行傅里叶变换得到最优小波系数的解调谱，为了提高包络谱线的分辨率，精确地诊断出缺陷类型，对解调谱进行细化分析，得到最优小波系数解调细化谱，然后根据调制频率，即边频带的间隔与检测目标故障所在的转频的对应关系，对检测目标的故障精确定位。

对最优小波系数包络幅值进行解调，图2f为的最优小波系数的解调谱，从图2f中可看出解调出频率都密集在低频段，两个齿轮的故障频率分别为10Hz和11.08Hz相差不大，故障频率更加无法分辨，为了提高包络谱线的分辨率，精确地诊断出缺陷类型，有必要对解调谱进行低频段细化。

步骤八中所述，检测目标里每个部件的转频不同，所述细化分析就是对解调谱进行低频段细化，对包络信号进行重新抽样，即每取10个点抽取一个点重新组成新的数据序列，在对新的数据序列进行频谱分析，得到细化后的解调谱。

这里对图2f做10倍的细化，对图2f中的包络信号进行重新抽样，即每取10个点抽取一个点重新组成新的数据序列，在对新的数据序列进行频谱分析，得到细化后的解调谱如图2g所示；从图2g中可以看出，小齿轮所在轴的转频11.093Hz，与理论值11.081Hz非常接近，其倍频22.174Hz在解调谱上出现，由此可断定是小齿轮出现了故障，且从图g与图f的比较可以看出，细化后的包络谱所提取的特征频率，这里为转频11.093Hz，比未细化的包络谱所提取的特征频率的精度高，更接近于理论值，降低了误诊断的概率，大大提高了故障诊断的准确性。

Claims (2)

1.一种复解析最优小波解调法，其特征在于，包括如下步骤：

步骤一、通过传感器获取检测目标的振动信号时域波形，并对其进行傅里叶变换可得到振动信号频谱图，确定信号分析频率的上下界限f_xmin，f_xmax；

步骤二、截取所述振动信号频谱图中信号的有效频带集中区域；

步骤三、设定Morlet小波的外形因子σ和初始频率ω₀的初始值，Morlet小波函数的表达式为：

结合分析频率的上下界限f_xmin，f_xmax，根据小波变换的尺度a与基小波采样率、中心频率和信号分析频率关系确定出a的取值范围E：

$E=(\frac{{\mathrm{\ω}}_0{f}_s}{2\mathrm{\π}{f}_{x\max }{f}_{\mathrm{\ω}}},\frac{{\mathrm{\ω}}_0{f}_s}{2\mathrm{\π}{f}_{x\min }{f}_{\mathrm{\ω}}})$

其中：f_s为信号采样率；f_ω为小波采样率；

采用以上设置的具体参数值，在尺度范围E内，对被分析信号s(t)作关于小波函数的连续小波变换，得到小波系数C_s(b，a)，其中b是位置参数；

步骤四、小波参数迭代确定最佳小波基：由小波变换系数C_s(b，a)可得出在尺度a下被分析信号s(t)包络幅值A_S(b，a)＝|C_s(b，a)|，进一步可以得出连续小波幅值谱信息测度SH：

$\mathrm{SH}=-\underset{a\∈E}{\mathrm{\Σ}}\underset{b=i}{\mathrm{\Σ}}\left(\frac{{\left|A_s(b,a)\right|}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left|A_s(b,a)\right|}^2}\right)1g\left(\frac{{\left|A_s(b,a)\right|}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left|A_s(b,a)\right|}^2}\right)$

其中N为时间的离散序列长度；i为小波系数离散序列号；对ω₀与σ进行迭代搜索：先选用ω₀的初始值，将σ从初始值开始增加一个步幅，求得小波系数与小波幅值谱信息测度SH，直至完成σ的搜索范围；然后将ω₀加上一个步幅，σ又从初始值开始增加一个步幅，计算小波系数与小波幅值谱信息测度SH，直至完成σ的搜索范围，重复对ω₀与σ进行迭代搜索步骤直至完成ω₀的搜索范围；

通过前面的迭代搜索可确定出SH的最小值，其所对应的就是最优ω₀和σ，将其带入Morlet小波函数的表达式即可得到最佳小波函数；

步骤五、小波峭度改进：得到最佳小波函数后，在设定的尺度范围E内完成连续小波变换，为加大噪声和冲击之间的峭度差距，将小波峭度改进为小波系数的能量包络幅值的峭度K_sopt(a)：

$K_{\mathrm{sopt}}\left(a\right)=\frac{C_{4W}^{\′}\left(a\right)}{{(C_{2W}^{\′}\left(a\right)+\mathrm{\η}\max [C_{2W}^{\′}\left]\right)}^2}-3$

其中：C′_4W(a)表示小波包络幅值A_s(b，a)的四阶矩；C′_2W(a)表示A_s(b，a)的二阶矩；max[C′_2W]表述所有A_s(b，a)二阶矩的最大值；η为限制二阶矩幅度的常数，取决于最小二阶矩到最大二阶矩的动态波动范围；

搜索出最大小波峭度值，继而计算出每个尺度上的改进小波峭度K_sopt(a)；

步骤六、计算平滑指数：引入平滑指数对小波系数作进一步的筛选；对尺度范围E内的小波系数包络幅值计算平滑指数SI(a)：

$\mathrm{SI}\left(a\right)=\frac{{\left(\underset{i=1}{\overset{N}{\mathrm{\Π}}}{A}_s(b_i,a)\right)}^{1/N}}{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_s(b_i,a)}$

其中b_i为小波系数的时间离散序列；N为时间的离散序列长度；i为小波系数离散序列号；A_s(b_i，a)为不同时间离散序列下对应的包络幅值；

步骤七、计算最优小波包络谱：对尺度范围E内的与SI(a)均做归一化处理，得到小波系数筛选指标：

KSI(a)＝Normalized(K_sopt(a))-Normalized(SI(a))，

其中Normalized(K_sopt(a))是不同尺度a下的峭度；Normalized(SI(a))是不同尺度a下的平滑指数；按照KSI(a)最大原则对对小波系数进行筛选，即，在尺度范围E内，最大KSI(a)所对应的小波系数为最优小波系数，继而计算出最优小波包络幅值谱；

步骤八、解调谱细化：对最优小波包络幅值谱进行傅里叶变换得到最优小波系数的解调谱，为了提高包络谱线的分辨率，精确地诊断出缺陷类型，对解调谱进行细化分析，得到最优小波系数解调细化谱；

然后根据调制频率，即边频带的间隔与检测目标故障所在的转频的对应关系，对检测目标的故障精确定位。

2.根据权利要求1所述的一种复解析最优小波解调法，其特征在于：步骤八中所述，检测目标里每个部件的转频不同，所述细化分析就是对解调谱进行低频段细化，对包络信号进行重新抽样，即每取10个点抽取一个点重新组成新的数据序列，在对新的数据序列进行频谱分析，得到细化后的解调谱。