新手案例-基于变分模态分解的滚动轴承经典故障诊断研究

最新推荐文章于 2024-04-29 15:30:11 发布

在职研究僧

最新推荐文章于 2024-04-29 15:30:11 发布

阅读量1.3k

点赞数 27

文章标签： matlab 算法人工智能

本文链接：https://blog.csdn.net/ynn4818172/article/details/137565291

版权

1、文章概括

本文利用凯斯西楚大学轴承中心数据，采用了变分模态分解方法（variational mode decomposition，VMD）进行故障分量提取，然后采用了包络谱分析方法，提取故障特征频率，技术路线如下图所示。

提取的故障特征频率结果下图所示：

上图为凯斯西楚大学轴承数据的105.mat,该数据为内圈故障，理论故障特征频率为5.4152（内圈故障系数）*29.95（转频，1797/60）=162.18524Hz，从图中能明显发现各个IMF分量存在161.133Hz的特征峰以及2倍频，这与理论故障特频率接近，因此认为此时的轴承可能出现了内圈故障。

代码采用了Matlab 2023a进行运行，欢迎大家测试和提出问题！

2、基于变分模态分解的滚动轴承经典故障诊断方法

第一步：分析原始信号的时域波形，并利用FFT（快速傅里叶变换，matlab函数fft），分析原始信号的频域波形（频谱），观察各个特征谱峰，为后续VMD参数确定提供一定的判断依据；

第二步：将原始的振动信号输入变分模态分解方法(VMD)中,设置VMD的超参数（模态分解层数K为4，惩罚因子为2000，其他参数默认）；

第三步：分析VMD的获得本征模态函数（IMF），分析不同IMF分量的时域波形和频域波形，观察VMD分解过程中的各个IMF分量的中心频率的迭代过程。

第四步：对各个IMF分量进行进一步的包络谱分析，寻找相关的故障特征频率，确定滚动轴承的运行状态。

2.1 变分模态分解的理论

变分模态分解方法的理论可以参考文献①（见附录），它不是这篇软文的重点。作为一个工具，只需要了解它输入和输出，如下：

输入：

f:一维的时域振动信号（1XN，N为信号长度）

K：模态分解层数（1-10，经验所得）

alpha:惩罚因子（100~10000，经验）

输出：

u:IMF分量（多个时域的分量信号，K x N，N为信号长度，K为模态分解层数）

u_hat:IMF分量再进行FFT的结果（IMF的副产品）

omega:IMF分量中心频率的迭代过程（n X K,n为VMD的迭代次数，K为模态分解层数）

综上，输入的K和alpha需要人为预先定义，可以反复实验多次进行调整。输出仅关注u和omega，u_hat为u的衍生品，可忽略。

为了方便小白理解VMD方法的功能，可以将原始的一维信号理解为一个菜篮子，里面有各种蔬菜和水果，经过VMD方法后，可以将菜篮子中的西红柿、土豆、黄瓜、苹果、梨和香蕉分在不同的新篮子里。其中，新篮子的个数为模态分解层数，惩罚因子决定了蔬菜、水果分类的标准。

2.2 包络谱分析

包络谱分析是一种开展轴承故障诊断要用到的经典方法，几乎凡是提取故障特征频率，最后一步必须用到包络谱分析。这主要是因为故障信号往往多和原始信号发生调制现象，利用包络谱是进行解调，让故障特征峰更加明显和显眼，从而进一步确定轴承是否发生了故障。

包络谱分析主要用到Matlab自带的hilbert（希尔伯特变换函数）和fft（傅里叶变换）函数。

3、实验分析

实验数据采用凯斯西楚大学（CWRU）的轴承数据，相关数据说明可参考如下软文：

CWRU(凯斯西储大学)轴承数据

本文选择了0HP、1797rpm的105.mat数据，它为内圈故障，缺陷尺寸为0.007英寸。该信号的采样频率为12000Hz，选择了驱动端数据进行分析，即X105_DE_time。

第一步：原始信号的时域和频域分析，代码和结果如下：

clear all;
close all;
clc;

% Time Domain 0 to T
Fs = 12000;            % Sampling frequency  凯斯西储大学
load('105.mat');  %内圈
f_start= X105_DE_time;
f(1:2^12)=f_start(10001:10000+2^12);

T = 1/Fs;             % Sampling period       
L = length(f);             % Length of signal
t = (0:L-1)*T;        % Time vector
f_pin_yu=fft(f);
P2 = abs(f_pin_yu/L);
P1 = P2(1:L/2);  %%频谱图纵坐标
P1(2:end-1) = 2*P1(2:end-1);
frequ = Fs*(0:(L/2-1))/L; %%频谱图的横坐标频率
figure;
subplot(2,1,1);plot(t,f,'b');xlabel('t/s');ylabel('Amplitude/g');
subplot(2,1,2);plot(frequ,P1,'b');xlabel('f/Hz');ylabel('Amplitude');
sgtitle('原始信号的时域波形和频域波形')

此时，从上图原始信号的频谱中已经能找到内圈的故障特征频率166.133Hz，但是从频谱中能发现，该谱峰相对周围的谱峰较小。因此，为了进一步的验证猜想，需要开展进一步的振动信号分析。

第二步：将原始的振动信号输入变分模态分解方法(VMD)中,设置VMD的超参数（模态分解层数K为4,惩罚因子为2000，其他参数默认）

% some sample parameters for VMD
alpha = 2000;        % moderate bandwidth constraint
tau = 0;            % noise-tolerance (no strict fidelity enforcement)
K = 4;              % 3 modes
DC = 0;             % no DC part imposed
init = 1;           % initialize omegas uniformly
tol = 1e-7;

%--------------- Run actual VMD code

 [u, u_hat, omega] = VMD(f, alpha, tau, K, DC, init, tol);

第三步：分析VMD的获得本征模态函数（IMF），分析不同IMF分量的时域波形和频域波形，观察VMD分解过程中的各个IMF分量的中心频率的迭代过程。

%--------------- Visualization

% For convenience here: Order omegas increasingly and reindex u/u_hat
[~, sortIndex] = sort(omega(end,:));
omega = omega(:,sortIndex);
u_hat = u_hat(:,sortIndex);
u = u(sortIndex,:);
linestyles = {'b', 'g', 'm', 'c', 'r', 'k'};

% 原始信号和IMF分量信号的时域分析
figure;
for k = 1:K
    subplot(K,1,k);
    plot(t,u(k,:),linestyles{k});
    xlabel('t/s');ylabel('Amplitude/g');
    legend(['IMF',num2str(k)])
end
sgtitle('信号的时域分析')


% 原始信号和IMF分量信号的频域分析
figure;
plot(frequ, P1, 'k:');
hold on;
for k = 1:K
    plot(frequ(1:end), abs(u_hat(L/2+1:end,k))/L, linestyles{k});
end
legend('原始信号','IMF1','IMF2','IMF3','IMF4')
xlabel('f/Hz');ylabel('Amplitude');
title('信号的频谱分析')

% VMD方法的IMF分量的中心频率迭代过程
figure;
for k=1:K
    plot(Fs/2/0.5*omega(:,k), 1:size(omega,1), linestyles{k});
    hold on;
end
legend('IMF1','IMF2','IMF3','IMF4')
xlabel('f/Hz');ylabel('迭代次数');
title('中心频率的迭代过程')

上图为4个IMF分量的时域波形和频域波形，从时域波形中能发现类似缺陷导致的冲击信号的影子，从频域波形中4个IMF分量基本完全覆盖了原始信号的主要频带。

上图为随着VMD的迭代次数增加，IMF的中心频率的变化情况。从图中能发现，从第10次开始后，IMF的中心频率基本不发生变化，即原始信号的分解过程逐渐收敛。

第四步：对各个IMF分量进行进一步的包络谱分析，寻找相关的故障特征频率，确定滚动轴承的运行状态：

figure;
u_hi1=hilbert(u');%希尔伯特变换
u_h_bao_luo=abs(u_hi1);
u_mean=mean(u_h_bao_luo);
for i=1:(size(u,1))
    u_h_bao_luo(:,i)=u_h_bao_luo(:,i)-u_mean(1,i);
end
f_VMD_bao_luo=fft(u_h_bao_luo);
P2_VMD_bao_luo = abs(f_VMD_bao_luo/L);
P1_VMD_bao_luo = P2_VMD_bao_luo(1:L/2,:);  %%频谱图纵坐标
P1_VMD_bao_luo(2:end-1) = 2*P1_VMD_bao_luo(2:end-1);
for k = 1:K
    subplot(K,1,k);
    title(['IMF',num2str(k)]);
    plot(frequ(1:end),P1_VMD_bao_luo(:,k), linestyles{k});
    axis([0 1000 min(P1_VMD_bao_luo(:,k)) max(P1_VMD_bao_luo(:,k))]);
    legend(['IMF',num2str(k)])
    xlabel('f/Hz');ylabel('Amplitude');
end
sgtitle('各个IMF分量的包络谱分析')

上图为4个IMF分量的包络谱分析结果，从中我们能够发现，每个IMF中都包含了内圈的故障特征频率161.133Hz(f_i),在IMF3中甚至能发现2倍频，在IMF1中能找到转频的二倍频58.5939Hz(2f_r)。这进一步验证了轴承可能出现了内圈缺陷的可能性。

总结

上述仅仅是一个利用变分模态分解方法和包络谱分析开展滚动轴承的经典故障诊断的一个例子。针对不同的工程实际情况，采用不同的信号处理方法，进行信号降噪、故障信号提取等工作，最终提取关键部件的故障特征频率。

值得注意的是，这里的关键部件就我目前了解的一般有齿轮和轴承两类部件，欢迎小伙伴们进一步补充。

代码和参考文献

1、上述源代码及其数据文件：

①105.mat （凯斯西楚大学轴承数据）

②VMD.m (VMD的函数文件，可从mathwork的file exchange获得，网址：https://ww2.mathworks.cn/matlabcentral/fileexchange/44765-variational-mode-decomposition?s_tid=srchtitle_support_results_12_VMD)

③VMD_experminent.m (基于变分模态分解的滚动轴承经典故障诊断研究的源代码)

2、相关参考

①Dragomiretskiy K, Zosso D. Variational mode decomposition[J]. IEEE transactions on signal processing, 2013, 62(3): 531-544.（VMD的原文）

②凯斯西楚大学数据参考请点击过去的参考资料：CWRU(凯斯西楚大学)轴承数据

关注公众号“故障诊断与寿命预测工具箱”，每天进步一点点。