一种基于变分模式分解滤波的包络分析方法与流程

本发明涉及旋转机械状态监测与故障诊断领域，具体涉及一种基于变分模式分解滤波的包络分析方法。

背景技术：

包络分析技术广泛应用于齿轮和滚动轴承的故障诊断中。现有的包络分析技术有下面三个缺陷：①现有的包络分析技术或者是直接对原始信号进行分析，或者是仅对原始信号进行简单的滤波后再进行分析，因此现有的方法容易受到噪声、趋势及其它成分的干扰，从而导致现有技术的分析精度较低；②现有的包络分析技术是以Hilbert变换为基础，而Hilbert变换要求被分析的信号必须是单分量的窄带信号，否则信号的频率调制部分将要污染信号的幅值包络分析结果，但是目前待分析的信号都不严格满足单分量且窄带的条件，这样就会导致现有技术因精度不高而容易出现误判问题；③ 由传统方法得到的包络谱存在着端点效应。

技术实现要素：

本发明要解决的问题是针对以上不足，提出一种基于变分模式分解滤波的包络分析方法，采用本发明的包络分析方法后，具有分析结果准确度和精确度高，并能准确地检测出旋转机械故障类型的优点。

为解决以上技术问题，本发明采取的技术方案如下：一种基于变分模式分解滤波的包络分析方法，其特征在于，包括以下步骤：

步骤1：利用加速度传感器以采样频率fs测取旋转机械的振动信号x（k）, （k=1, 2, …,N），N为采样信号的长度；

步骤2：采用变分模式分解(Variational Mode Decomposition)算法将信号x（k）分解成n个分量之和，即 ，其中，c_i（k）代表由变分模式分解算法得到的第i个分量；

步骤3：对c_i（k）执行重排操作和替代操作，经重排操作得到的数据用c_i^shuffle（k）表示，替代操作后得到数据用c_i^FTran（k）表示；

步骤4：对c_i（k）、c_i^shuffle（k）和c_i^FTran（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，c_i（k）的广义Hurst指数曲线用H_i（q）表示；c_i^shuffle（k）的广义Hurst指数曲线用H_i^shuffle（q）表示；c_i^FTran（k）的广义Hurst指数曲线用H_i^FTran（q）表示；

步骤5：如果H_i（q） 与H_i^shuffle（q）或H_i（q） 与H_i^FTran（q）之间的相对误差小于5%，或者H_i（q） 、H_i^shuffle（q） 和H_i^FTran（q）三者都不随q而变化，则抛弃对应的c_i（k）分量；

步骤6：对剩余的c_i（k）分量求和，将该和记为信号经重排和替代滤波后的结果x_f1（k）；

步骤7：对x_f1（k）执行谱峭度分析，求出信号峭度最大处所对应的中心频率f₀和带宽B；

步骤8： 根据中心频率f₀和带宽B对x_f1（k）进行带通滤波，得到x_f2（k）；

步骤9：对信号x_f2（k）执行有理样条迭代平滑包络分析，得到信号包络eov（k）；

步骤10：对得到的信号包络eov（k）执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。

一种优化方案，所述步骤3中数据重排操作包括以下步骤：

随机打乱分量c_i（k）的排列顺序。

进一步地，所述步骤3中数据替代操作包括以下步骤：

1） 对分量c_i（k）执行离散傅里叶变换，获得分量c_i（k）的相位；

2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量c_i（k）的原始相位；

3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据c_i^IFFT（k），求取数据c_i^IFFT（k）的实部。

进一步地，所述步骤4中MFDFA方法包括以下步骤：

1）构造x(k)(k=1,2,…,N)的轮廓Y(i)：

；

x (k )代表权利要求1所述步骤4中的c_i（k）或c_i^shuffle（k）或c_i^FTran（k）；

2）将信号轮廓Y (i )分成不重叠的N_s 段长度为s 的数据，由于数据长度N 通常不能整除s，所以会剩余一段数据不能利用；

为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N_s段数据；

3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：

；

y_v(i )为拟合的第v 段数据的趋势，若拟合的多项式趋势为m 阶，则记该去趋势过程为（MF-）DFAm ；本例中，m=1；

4） 计算第q 阶波动函数的平均值：

；

5）如果x (k )存在自相似特征，则第q 阶波动函数的平均值F_q(s )和时间尺度s 之间存在幂律关系：

F_q(s)~s^H(q)；

当q =0时，步骤4）中的公式发散，这时H (0)通过下式所定义的对数平均过程来确定：

；

6）对步骤5）中的公式两边取对数可得ln[F_q(s)]= H (q)ln(s )+c (c为常数)，由此可以获得直线的斜率H (q )。

进一步地，所述步骤7中的谱峭度方法包括以下步骤：

1）构造一个截止频率为f_c =0.125+ε 的低通滤波器h(n )；ε >0，本例中f_c=0.3；

2）基于h (n )构造通频带为[0;0.25]的准低通滤波器h₀(n )和通频带为[0.25;0.5]的准高通滤波器h₁(n )，

；

3）信号经h₀(n )、h₁(n)滤波并降采样后分解成低频部分和高频部分，降采样的因子为2，再经多次迭代滤波后形成滤波器树，第k层有2^k个频带，其中表示滤波器树中第k 层上的第i 个滤波器的输出信号，i =0,…,2^k-1，0≤k≤K-1，本例中K=8；c₀(n )代表权利要求1所述步骤7中x_f1（k）；

4）分解树中第k 层上的第i 个滤波器的中心频率f_ki和带宽B_k分别为

f_ki=(i +2^-1)2^-k-1

B_k=2^-k-1

5）计算每一个滤波器结果(i =0,…,2^k-1)的峭度；

6）将所有的谱峭度汇总，得到信号总的谱峭度。

进一步地，所述步骤9中的有理样条迭代平滑包络分析方法包括以下步骤：

1）计算信号z (k )的绝对值∣z(k )∣的局部极值；在第1次迭代中，z (k )代表权利要求1所述步骤9中x_f2（k）；

2）采用有理样条曲线拟合局部极值点得到包络线eov₁（k）；

3）对z(k )进行归一化处理得到；

4）第2次迭代：把z₁(k)重新作为新数据，重复执行上述步骤1）~3），得到；

5）第i次迭代：把z_i-1(k) 重新作为新数据，重复执行上述步骤1）~3)，得到；

6) 如果第n 次迭代得到的z_n(k)的幅值小于或等于1，则迭代过程停止，最后得到信号z(k )的包络为。

本发明采用以上技术方案，与现有技术相比，本发明具有以下优点：

1)利用变分模式分解对原始信号进行分解，然后利用数据的重排和替代操作排除其中的噪声和趋势分量，仅仅保留信号分量中的有用成分，从而避免了噪声和趋势分量对包络分析结果的影响，分析结果准确度和精确度高。

2)利用有理样条迭代平滑包络分析方法将信号包络与频率调制部分完全分离，能够避免频率调制部分对信号包络分析结果的影响，从而提高包络分析的精度。

3) 能够准确地检测出旋转机械的故障类型。

4) 由传统方法得到的包络谱存在端点效应，而由本发明得到的包络谱能够避免端点效应。

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

附图说明

附图1为本发明实施例中本发明方法的流程图；

附图2为本发明实施例中采用低通滤波器和高通滤波器对信号进行初步分解的示意图；

附图3为本发明实施例中采用树状滤波器结构快速计算谱峭度的示意图；

附图4为本发明实施例中具有内圈故障的滚动轴承振动信号；

附图5为本发明实施例中采用传统包络分析方法对内圈故障滚动轴承振动信号的分析结果；

附图6为本发明实施例中本发明对内圈故障滚动轴承振动信号的分析结果；

附图7为本发明实施例中具有外圈故障的滚动轴承振动信号；

附图8为本发明实施例中采用传统包络分析方法对外圈故障滚动轴承振动信号的分析结果；

附图9为本发明实施例中本发明对外圈故障滚动轴承振动信号的分析结果。

具体实施方式

实施例，如图1、图2、图3所示，一种基于变分模式分解滤波的包络分析方法，其特征在于，包括以下步骤：

步骤1：利用加速度传感器以采样频率fs测取旋转机械的振动信号x（k）, （k=1, 2, …,N），N为采样信号的长度；

步骤2：采用变分模式分解(Variational Mode Decomposition)算法将信号x（k）分解成n个分量之和，即 ，其中，c_i（k）代表由变分模式分解算法得到的第i个分量，变分模式分解已公知，见文献

Konstantin Dragomiretskiy, Dominique Zosso. Variational Mode Decomposition, IEEE TRANSACTIONS ON SIGNAL PROCESSING, 2014, 62(3) : 531-544；本例中，设置模态数目为10；

步骤3：对c_i（k）执行重排操作和替代操作，经重排操作得到的数据用c_i^shuffle（k）表示，替代操作后得到数据用c_i^FTran（k）表示；

步骤4：对c_i（k）、c_i^shuffle（k）和c_i^FTran（k）分别执行多重分形去趋势波动分析（Multifractal Detrended Fluctuation Analysis, MFDFA），得到广义Hurst指数曲线，c_i（k）的广义Hurst指数曲线用H_i（q）表示；c_i^shuffle（k）的广义Hurst指数曲线用H_i^shuffle（q）表示；c_i^FTran（k）的广义Hurst指数曲线用H_i^FTran（q）表示；

步骤5：如果H_i（q） 与H_i^shuffle（q）或H_i（q） 与H_i^FTran（q）之间的相对误差小于5%，或者H_i（q） 、H_i^shuffle（q） 和H_i^FTran（q）三者都不随q而变化，则抛弃对应的c_i（k）分量；

步骤6：对剩余的c_i（k）分量求和，将该和记为信号经重排和替代滤波后的结果x_f1（k）；

步骤7：对x_f1（k）执行谱峭度分析，求出信号峭度最大处所对应的中心频率f₀和带宽B；

步骤8： 根据中心频率f₀和带宽B对x_f1（k）进行带通滤波，得到x_f2（k）；

步骤9：对信号x_f2（k）执行有理样条迭代平滑包络分析，得到信号包络eov（k）；

步骤10：对得到的信号包络eov（k）执行离散傅里叶变换得到包络谱，根据包络谱特征频率判断机器的故障类型。

步骤3中数据重排操作包括以下步骤：

随机打乱分量c_i（k）的排列顺序。

步骤3中数据替代操作包括以下步骤：

1） 对分量c_i（k）执行离散傅里叶变换，获得分量c_i（k）的相位；

2） 用一组位于(-π,π)区间内的伪独立同分布数来代替分量c_i（k）的原始相位；

3） 对经过相位替代后的频域数据执行离散傅里叶逆变换得到数据c_i^IFFT（k），求取数据c_i^IFFT（k）的实部。

步骤4中MFDFA方法包括以下步骤：

1）构造x (k )(k =1,2,…,N )的轮廓Y (i )：

；

x (k )代表权利要求1所述步骤4中的c_i（k）或c_i^shuffle（k）或c_i^FTran（k）；

2）将信号轮廓Y (i )分成不重叠的N_s段长度为s 的数据，由于数据长度N 通常不能整除s，所以会剩余一段数据不能利用；

为了充分利用数据的长度，再从数据的反方向以相同的长度分段，这样一共得到2N_s段数据；

3）利用最小二乘法拟合每段数据的多项式趋势，然后计算每段数据的方差：

；

y_v(i )为拟合的第v 段数据的趋势，若拟合的多项式趋势为m 阶，则记该去趋势过程为（MF-）DFAm ；本例中，m=1；

4） 计算第q 阶波动函数的平均值：

；

5）如果x (k )存在自相似特征，则第q 阶波动函数的平均值F_q(s )和时间尺度s 之间存在幂律关系：

F_q(s)~s^H(q)；

当q=0时，步骤4）中的公式发散，这时H (0)通过下式所定义的对数平均过程来确定：

；

6）对步骤5）中的公式两边取对数可得ln[F_q(s )]=H (q )ln(s )+c (c 为常数)，由此可以获得直线的斜率H (q )。

步骤7中的谱峭度方法包括以下步骤：

1）构造一个截止频率为f_c=0.125+ε 的低通滤波器h (n )；ε >0，本例中f_c=0.3；

2）基于h (n )构造通频带为[0;0.25]的准低通滤波器h₀(n )和通频带为[0.25;0.5]的准高通滤波器h₁(n )，

；

3）信号经h₀(n )、h₁(n )滤波并降采样后分解成低频部分和高频部分，降采样的因子为2，再经多次迭代滤波后形成滤波器树，第k层有2^k个频带，其中表示滤波器树中第k 层上的第i 个滤波器的输出信号，i =0,…,2^k-1，0≤k≤K-1，本例中K=8；c₀ (n )代表权利要求1所述步骤7中x_f1（k）；

4）分解树中第k 层上的第i 个滤波器的中心频率f_ki和带宽B_k分别为

fki =(i +2^-1)2^-k-1

B_k=2^-k-1

5）计算每一个滤波器结果(i=0,…,2^k-1)的峭度；

6）将所有的谱峭度汇总，得到信号总的谱峭度。

步骤9中的有理样条迭代平滑包络分析方法包括以下步骤：

1）计算信号z(k )的绝对值∣z(k )∣的局部极值；在第1次迭代中，z(k )代表权利要求1所述步骤9中x_f2（k）；

2）采用有理样条曲线拟合局部极值点得到包络线eov₁（k）；

3）对z(k )进行归一化处理得到；

4）第2次迭代：把z₁(k)重新作为新数据，重复执行上述步骤1）~3），得到；

5）第i次迭代：把z_i-1(k) 重新作为新数据，重复执行上述步骤1）~3)，得到；

6) 如果第n 次迭代得到的z_n(k)的幅值小于或等于1，则迭代过程停止，最后得到信号z(k )的包络为。

试验1，利用具有内圈故障的滚动轴承振动数据对本发明所述算法的性能进行验证。

实验所用轴承为6205-2RS JEM SKF，利用电火花加工方法在轴承内圈上加工深度为0.2794mm、宽度为0.3556mm的凹槽来模拟轴承内圈故障，本实验负载约为0.7457kW，驱动电机转频约为29.5Hz，轴承内圈故障特征频率约为160Hz，采样频率为4.8KHz，信号采样时长为1s。

采集到的内圈故障信号如图4所示。

首先采用传统的包络分析方法对图4所示的信号进行分析，得到的分析结果如图5所示。从图5可以看出，轴承的故障特征完全被掩盖，因此传统的包络分析方法不能有效地提取轴承的故障特征；此外，从图5可以看出，包络谱的左端点存在着异常高值，这说明由传统方法得到的包络谱存在着端点效应。

采用本发明所提出的方法对图4所示的信号进行分析，得到的分析结果如图6所示。从图6可以看出，160Hz和320Hz所对应的谱线明显高于其它谱线，这两个频率分别对应轴承内圈故障特征频率的1倍频和2倍频，据此可以判断轴承具有内圈故障；从图6可以看出，由本发明得到包络谱没有端点效应。

经多次实验表明，在负载和故障尺寸深度不变的情况下，本发明能够可靠识别的最小内圈故障尺寸宽度约为0.22 mm，而传统方法能够可靠识别的最小内圈故障尺寸宽度约为0.53mm，精度提高58.5%。

试验2，利用具有外圈故障的滚动轴承振动数据对本发明所述算法的性能进行验证。

实验所用轴承为6205-2RS JEM SKF，利用电火花加工方法在轴承外圈上加工深度为0.2794mm、宽度为0.5334mm的凹槽来模拟轴承外圈故障，本实验负载约为2.237 kW，驱动电机转频约为28.7Hz，轴承外圈故障特征频率约为103Hz，采样频率为4.8KHz，信号采样时长为1s。

采集到的外圈故障信号如图7所示。

首先采用传统的包络分析方法对图7所示的信号进行分析，得到的分析结果如图8所示。从图8可以看出，轴承的故障特征完全被掩盖，因此传统的包络分析方法不能有效地提取轴承的故障特征；此外，从图8可以看出，包络谱的左端点存在着异常高值，这说明由传统方法得到的包络谱存在着端点效应。

采用本发明所提出的方法对图7所示的信号进行分析，得到的分析结果如图9所示。从图9可以看出，103Hz和206Hz所对应的谱线明显高于其它谱线，这两个频率分别对应轴承外圈故障特征频率的1倍频和2倍频，据此可以判断轴承具有外圈故障；从图9可以看出，由本发明得到包络谱没有端点效应。。

经多次实验表明，在负载和故障尺寸深度不变的情况下，本发明能够可靠识别的最小外圈故障尺寸宽度约为0.28mm，而传统方法能够可靠识别的最小外圈故障尺寸宽度约为0.68mm，精度提高58.8%。

根据试验结果，分析后认为：

1) 传统的包络分析方法直接对原始信号进行包络分析，或者对仅经过简单处理后的原始信号进行包络分析，与传统的包络分析方法不同，本发明首先利用变分模式分解对原始信号进行分解，然后利用数据的重排和替代操作排除其中的噪声和趋势分量，仅仅保留信号分量中的有用成分，从而避免了噪声和趋势分量对包络分析结果的影响，提高了准确度和精确度。

2) 传统的包络分析方法以Hilbert变换为基础，而Hilbert变换要求被分析的信号必须是单分量的窄带信号，否则信号的频率调制部分将要污染信号的包络分析结果，但是目前待分析的信号都不严格满足单分量且窄带的条件，这样就会导致现有技术因精度不高而容易出现误判问题，与传统包络分析方法不同，本发明利用有理样条迭代平滑包络分析方法将信号包络与频率调制部分完全分离，能够避免频率调制部分对信号包络分析结果的影响，从而提高包络分析的精度。

3)能够准确地检测出旋转机械的故障类型。

4) 由传统方法得到的包络谱存在端点效应，而由本发明得到的包络谱能够避免端点效应。

本领域技术人员应该认识到，上述的具体实施方式只是示例性的，是为了使本领域技术人员能够更好的理解本发明内容，不应理解为是对本发明保护范围的限制，只要是根据本发明技术方案所作的改进，均落入本发明的保护范围。