【时频分析】轴承故障诊断的集中式时频分析工具【附MATLAB代码】

文章来源：微信公众号：EW Frontier 关注可了解更多的雷达、通信、人工智能相关代码。

摘要

在工业旋转机械中，瞬态信号通常对应于主要元件（例如轴承或齿轮）的故障。然而，面对实际工程的复杂性和多样性，提取瞬态信号是一项极具挑战性的任务。本文提出了一种新的时频分析方法--暂态提取变换法，该方法能有效地表征和提取故障信号中的暂态分量。该方法基于短时傅里叶变换，不需要扩展参数或先验信息。量化指标，如Rényi熵和峰度，所提出的方法与其他经典和先进的方法的性能进行比较。实验结果表明，该方法可以提供能量更集中的时频表示，并能以更大的峰度提取瞬态分量。数值和实验信号的使用表明我们的方法的有效性。

引言

在旋转机械的故障诊断领域中，信号处理方法被广泛应用于发现与机械故障密切相关的特征[1]，[2]。在记录的振动和声音信号中，故障通常表现出在短时间内出现的瞬态特征[3]。考虑到不同的故障信号占据不同的频带，联合时频（TF）分析（TFA）是表征具有非平稳TF特征的瞬态故障的有效工具[4]。虽然在许多研究中已经报道了TFA方法在故障诊断中的直接应用，但经典TFA方法的固有缺点从未得到有效解决。线性TFA方法，例如，利用短时傅里叶变换（STFT）和小波变换（WT）计算信号与具有TF特征定位能力的基函数之间的内积。然而，由于没有TF基函数可以被压缩支持，同时，线性TFA方法表现出很差的能力来表征精确的TF特征。双线性TFA方法，如Wigner-Ville分布和Cohen类分布，用于计算局部信号相关的傅里叶变换。然而，双线性TFA方法中的交叉项极大地限制了其应用。传统的故障树分析方法存在的这些缺陷，降低了诊断系统对早期弱故障和强噪声故障的敏感性。为了增强TFA方法在复杂环境中检测故障的能力，在过去十年中已经提出并引入了一些先进的方法，例如，经验模态分解（EMD）[5]、谱峰度（SK）方法[6]、[7]和同步压缩变换（SST）[8]-[10]。

经验模态分解是一种数据驱动的方法，可以将一维信号分解为一系列的本征模态函数（IMF）。由于不同的IMF占据不同的频带，包含故障频带的IMF的瞬态特征与原始信号相比可以被高度增强。由于这种优越性，许多基于EMD的故障诊断方法已经开发出来，并且可以在[11]中找到全面的综述。虽然我们不能很好地理解这种方法的数学基础，但一些研究表明，EMD在处理高斯噪声时表现为一个并矢滤波器组。它表示在执行时间序列信号处理时，EMD使用固定的二进滤波器组分解信号的因素。由于实际信号中的故障分量的频带不能预先知道，因此分解结果使得一些IMF可能包含期望的故障分量，或者一个故障分量可能被分解为多个IMF，这通常被称为模式混合。由于经验模态分解的过程行为难以控制，使得基于经验模态分解的故障诊断方法有时具有不可预测性和不稳定性。最近，建立了更先进的方法来改善EMD的性能，例如，局部均值分解[12]，集合EMD [13]和极值点加权模式分解[14]。

SK方法是一种基于峰度指标提取最瞬态分量的技术。峰度是一个用来度量时间序列信号时间离散度的统计量，它也可以用来检测故障信号中包含的暂态。SK方法首先需要基于STFT或带通滤波器将一维信号扩展到二维TF平面，可以重建或选择与具有最大峰度的故障最相关的分量。得益于峰度指标对瞬态故障的敏感性，SK方法在诊断机械故障方面显示出其有效性[15]，[16]。

SST方法被引入作为线性TFA方法的后处理工具，并已应用于旋转机械的故障诊断[10]。SST的目的是获得更清晰的TF表示，它可以在高TF分辨率下表征故障。同时，瞬态分量可以从更尖锐的TF结果中提取。为了从SST结果中提取信号，必须首先估计与瞬态分量相对应的IF轨迹。然而，它是具有挑战性的，以准确地估计的瞬时分量的IF，因为故障信号通常不满足弱时变的SST框架的要求。此外，不期望的背景噪声会给SST结果带来严重的干扰，这可能导致IF不能被准确地表征。为了进一步提高SST的性能，提出了一些先进的方法，例如，解调制SST [17]、匹配SST [18]、高阶SST [19]和同步提取变换（SET）[4]。

从上面的介绍中，我们可以看到，已经引入了许多先进的技术来提取原始信号中的瞬态分量，这是提高诊断系统故障检测能力的一个重要问题。在本文中，我们提出了一种新的TFA方法，可以在TF平面上精确地描述瞬态特征，并在时域中提取它。将该方法与SK、EMD、SST等先进的故障诊断方法及其改进方法进行了比较。本文的其余部分组织如下。第二节详细介绍了我们提出的方法的理论。在第三节中，使用Rényi熵和峰度指标来说明不同TFA方法产生的TF结果的量化比较。实验验证见第IV节和第V节。结论见第VI节。

MATLAB代码

function [tfr Te GD TEO] = TET_Y(x,hlength);
%   Transient-extracting transform
%   Input
%  x       : Signal.
%  hlength : Window length.
​
%   Output
%  tfr   :  STFT Representation.
%  Te    :  TET Representation.
%  GD    :  2D  group delay.
%  TEO   :  Transient-extracting operator.
​
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
%
%   Written at 2017.3.17.
​
[xrow,xcol] = size(x);
​
N=xrow;
​
if (xcol~=1),
 error('X must be column vector');
end;
​
if (nargin < 2),
 hlength=round(xrow/8);
end;
​
t=1:N;
%ft = 1:round(N/2);
​
[~,tcol] = size(t);
​
hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';
​
% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
​
% Window tg(t)
th=h.*ht;
[hrow,hcol]=size(h); Lh=(hrow-1)/2;
​
tfr1= zeros (N,tcol);
tfr2= zeros (N,tcol);
​
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr2(indices,icol)=rSig.*conj(th(Lh+1+tau));
end;
​
tfr1=fft(tfr1);
tfr2=fft(tfr2);
​
tfr1=tfr1(1:round(N/2),:);
tfr2=tfr2(1:round(N/2),:);
​
E=mean(abs(x));
omega= zeros(round(N/2),tcol);
​
for a=1:round(N/2)
omega(a,:) = t+(hlength-1)*real(tfr2(a,t)./tfr1(a,t));
end
​
GD=omega;
TEO=zeros(round(N/2),tcol);
​
for i=1:round(N/2)%frequency
for j=1:N%time
     if abs(tfr1(i,j))>12*E%
      if abs(omega(i,j)-j)<0.5%
        TEO(i,j)=1;
      end
    end
end
end
tfr=tfr1/(xrow/2);%the amplitude of tfr result has been pre-rectified.
%tfr=tfr1/(sum(h)/2);%
Te=TEO.*tfr;
​

MATLAB仿真结果