一种改进解卷积算法在旋转机械故障诊断中的应用（MATLAB）

news2024/10/5 15:32:01

轴承振动是随机振动。在不同的时刻，轴承振动值是不尽相同的，不能用一个确定的时间函数来描述。这是由于滚动体除了有绕轴承公转运动以外，还有绕自身轴线的自旋运动，且在轴承运转时，滚动接触表面形貌是不断变化的，以及外来杂质进入运行区域的不确定性形成的。轴承振动也是宽带振动，其振动能量分布在一个宽带范围内，反映在幅值谱中就是轴承的振动频率范围为零到数万赫兹的连续谱。

当轴承正常时，其振动信号近似服从正态分布，此时信号中非噪声成分主要源于保持架转频和轴频，以及保持架对轴的旋转产生的微弱调制，信号能量在频域分布比较均匀，且各个频率分量的不确定性都较大且比较接近；而当轴承出现故障时，故障冲击将激励系统共振，产生脉冲冲击串，这将改变轴承振动信号的概率分布和能量分布，信号能量就会逐渐集中到故障特征频率（带）上。且随着轴承性能退化，信噪水平增加，信号的频率分布变得越来越规则，轴承振动信号的能量也会更加集中。

当轴承故障时，其原始脉冲冲击信号受随机噪声信号和周期信号干扰，并经过一定的传递路径衰减等因素影响，传感器采集的振动信号往往不能够反映真实的轴承故障信息，为设备故障诊断和状态监测增加了难度。解卷积类算法可以突出信号中的脉冲冲击成分，使得滤波后信号是原来脉冲冲击成分的一个近似。因此，可以认为经解卷积滤波后的信号是一个更加接近于轴承脉冲冲击信号的信号，对该信号做包络谱分析等即可找出轴承故障频率，进而反映轴承的故障类型。而当轴承正常时，轴承的脉冲冲击信号主要为轴的转频成分。经解卷积滤波处理后的信号，其在包络谱中对应的峰值频率为轴的转频，为判断轴承正常提供理论依据。

鉴于此，采用一种改进解卷积算法对旋转机械进行故障诊断，该算法不需要先验知识，效果较好。

function [CoSi, Bm] = CoSiEn(Sig, varargin)


narginchk(1,11)
Sig = squeeze(Sig);


p = inputParser;
Chk = @(x) isnumeric(x) && isscalar(x) && (x > 0) && (mod(x,1)==0);
Chk2 = @(x) isnumeric(x) && isscalar(x) && (x >= 0);
Chk3 = @(x) isnumeric(x) && isscalar(x) && (x > 1) && (mod(x,1)==0);
addRequired(p,'Sig',@(x) isnumeric(x) && isvector(x) && (length(x) > 10));
addParameter(p,'m',2,Chk3);
addParameter(p,'tau',1,Chk);
addParameter(p,'r',.1,@(x) isnumeric(x) && (x > 0) && (x < 1));
addParameter(p,'Logx',2,Chk2);
addParameter(p,'Norm',0,@(x) ismember(x,[0:4]));
parse(p,Sig,varargin{:})
m = p.Results.m; tau = p.Results.tau; r = p.Results.r;
Logx = p.Results.Logx; Norm = p.Results.Norm;


N = length(Sig);
if Logx == 0
    Logx = exp(1);
end
if Norm == 1
    Xi = Sig - median(Sig);
elseif Norm == 2
    Xi = Sig - mean(Sig);
elseif Norm == 3
    Xi = (Sig - mean(Sig))/std(Sig,1);
elseif Norm == 4
    Xi = (2*(Sig - min(Sig))/range(Sig)) - 1;
else
    Xi = Sig;
end
Nx = N-((m-1)*tau);
Zm = zeros(Nx,m);
for n = 1:m
    Zm(:,n) = Xi((n-1)*tau+1:Nx+(n-1)*tau);
end


Num = Zm*Zm'; 
Mag = sqrt(diag(Num));
Den = Mag*Mag';
AngDis = acos(Num./Den)/pi;
if max(imag(AngDis(:))) < (10^-5) %max(max(imag(AngDis))) < (10^-5)
    Bm = sum(sum(triu(round(AngDis,6) < r,1)))/(Nx*(Nx-1)/2);
else
    Bm = sum(sum(triu(real(AngDis) < r,1)))/(Nx*(Nx-1)/2);
    warning('Complex values ignored')
end
if Bm == 1 || Bm == 0
    CoSi = 0;
else
    CoSi = -(Bm*log(Bm)/log(Logx)) - ((1-Bm)*log(1-Bm)/log(Logx));
end
end
代码可通过知乎学术咨询获得：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1