该案例展示了如何基于加速信号,尤其是有其他机械部件强掩蔽信号存在的情况下,进行故障分析诊断。该案例将演示如何使用包络谱分析和谱峰度来诊断轴承故障,同时也可被拓展到大数据应用方面
滚动轴承的局部故障可能发生在外圈、内圈、保持架和滚动元件中。当健康滚动元件撞击内圈或外圈的局部故障,或故障滚动元件撞击健康内圈或外圈时,会激起高频谐振。下图中展示了健康滚动元件撞击轴承内圈局部故障。待研究的问题是如何检测并识别不同类型的故障。
本案例数据来源于MFPT(Machinery Failure Prevention Technology)挑战数据,该数据包含了来自遭受不同故障类型的不同机械,它包括23组数据集。前20组数据情况如下(前20组数据集来自轴承试验台,后3组为实际采集的真实数据,分别来自油泵轴承、中速轴承和行星轴承,故障位置未知):
本案例分析过程中仅使用前20组带标签的数据
每组数据包含了一条时序加速信号“gs”,采样频率“sr”,轴转速“rate”,载荷“load”,以及代表不同故障位置的四种临界频率:外圈缺陷频率(BPFO),内圈缺陷频率(BPFI),轴承保持器损坏频率(FTF),轴承滚动件损坏频率(BSF)10。以下四个公式描述了这四种频率的计算:
公式中d为球直径,D为节径,变量fr为轴转速,n为滚动元件(滚珠)的个数,ϕ为轴承接触角
在MFPT数据中,轴转速都是恒定的,因此没有必要进行阶次跟踪作为预处理步骤以消除转速变化的影响
当健康滚动元件撞击故障内圈或外圈,或故障滚动元件撞击健康内圈或外圈时,所带来的冲击会调制出对应的临界频率如BFPO、BFPI、FTF、BSF等。因此幅度解调产生的包络信号,会传达出更多通过使用原始信号进行频谱分析无法得到的诊断信息。接下来以MFPT数据中的内圈故障数据为例:
可视化时域内原始内圈故障数据
xInner = dataInner.bearing.gs; % 时序数据
fsInner = dataInner.bearing.sr; % 采样频率
tInner = (0:length(xInner)-1)/fsInner; % 生成时间轴
figure
plot(tInner, xInner)
xlabel('Time, (s)')
ylabel('Acceleration (g)')
title('Raw Signal: Inner Race Fault')
xlim([0 0.1])
使用功率谱生成函数pspectrum
查看频域内的原始数据
figure
[pInner, fpInner] = pspectrum(xInner, fsInner); % 生成功率谱,返回功率谱和对应的频率
pInner = 10*log10(pInner);
plot(fpInner, pInner)
xlabel('Frequency (Hz)')
ylabel('Power Spectrum (dB)')
title('Raw Signal: Inner Race Fault')
legend('Power Spectrum')
下面将频率放大到0-1000 Hz范围,进一步查看内圈缺陷频率的响应以及它出现第一次谐振时的情形
figure
plot(fpInner, pInner)
ncomb = 10; % 定义10个谐波游标
helperPlotCombs(ncomb, dataInner.BPFI) % 画出这些谐波游标
xlabel('Frequency (Hz)')
ylabel('Power Spectrum (dB)')
title('Raw Signal: Inner Race Fault')
legend('Power Spectrum', 'BPFI Harmonics')
xlim([0 1000])
帮助函数helperPlotCombs
如下:
function helperPlotCombs(ncomb, f)
%HELPERPLOTCOMBS Plot harmonic cursors on a power spectrum plot
% Copyright 2017 The MathWorks, Inc.
ylimit = get(gca, 'YLim');
ylim(ylimit);
ycomb = repmat([ylimit nan], 1, ncomb);
hold(gca, 'on')
for i = 1:length(f)
xcomb = f(i)*(1:ncomb);
xcombs = [xcomb; xcomb; nan(1, ncomb)];
xcombs = xcombs(:)';
plot(xcombs, ycomb, '--')
end
hold(gca, 'off')
end
从图中可以看到BPFI和它的谐波中观察不到明显的模式。此时对原始信号的频域分析也没法提供有用的诊断信息
再回过头来看时序数据,原始信号的振幅按一定频率调制,并且主频大约是1/0.009Hz ≈ 111Hz。已知健康滚动元件撞击内圈局部故障时的频率,也就是BPFI是118.875 Hz。这就说明该轴承可能有潜在的内圈故障
下面提取一段连续的振幅变化的时序数据进行可视化:
同时在下方使用函数envspectrum
函数计算包络谱信息以提取调制振幅,并画出原始信号的包络
subplot(2, 1, 2)
[pEnvInner, fEnvInner, xEnvInner, tEnvInner] = envspectrum(xInner, fsInner);
plot(tEnvInner, xEnvInner)
xlim([0.04 0.06])
xlabel('Time (s)')
ylabel('Acceleration (g)')
title('Envelope signal')
之后再使用包络信号计算功率谱,并查看对应的BPFI和它谐振的响应
figure
plot(fEnvInner, pEnvInner)
xlim([0 1000])
ncomb = 10;
helperPlotCombs(ncomb, dataInner.BPFI)
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum: Inner Race Fault')
legend('Envelope Spectrum', 'BPFI Harmonics')
可以看出,大部分能量集中于BPFI和它的谐振上。这一现象指示了轴承内圈的故障,这也佐证了数据的故障标签
现在将同样的包络谱分析方法应用于正常数据和轴承外圈故障数据中
% 载入正常数据
dataNormal = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
'predmaintdemos', 'bearingFaultDiagnosis', ...
'train_data', 'baseline_1.mat'));
xNormal = dataNormal.bearing.gs;
fsNormal = dataNormal.bearing.sr;
tNormal = (0:length(xNormal)-1)/fsNormal;
[pEnvNormal, fEnvNormal] = envspectrum(xNormal, fsNormal); % 分析正常数据的包络谱
figure
plot(fEnvNormal, pEnvNormal)
ncomb = 10;
helperPlotCombs(ncomb, [dataNormal.BPFO dataNormal.BPFI])
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum: Normal')
legend('Envelope Spectrum', 'BPFO Harmonics', 'BPFI Harmonics')
可以看出,对正常数据而言,包络谱并未表现出在BPFO和BPFI谐振上有明显的峰值出现
% 载入外圈故障数据
dataOuter = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
'predmaintdemos', 'bearingFaultDiagnosis', ...
'train_data', 'OuterRaceFault_2.mat'));
xOuter = dataOuter.bearing.gs;
fsOuter = dataOuter.bearing.sr;
tOuter = (0:length(xOuter)-1)/fsOuter;
[pEnvOuter, fEnvOuter, xEnvOuter, tEnvOuter] = envspectrum(xOuter, fsOuter);
figure
plot(fEnvOuter, pEnvOuter)
ncomb = 10;
helperPlotCombs(ncomb, dataOuter.BPFO)
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum: Outer Race Fault')
legend('Envelope Spectrum', 'BPFO Harmonics')
从外圈故障数据的包络谱来看,也没有明显的BPFO谐波出现。那么就可以断定包络谱无法区分外圈故障和健康状态吗?对此问题,可以退一步回过头来看不同条件下的原始时域数据
首先再次可视化这些时域数据并计算它们的峰度。峰度是随机变量的第四个标准化矩,它表征了信号的脉冲或随机变量尾部的肥瘦程度
figure
subplot(3, 1, 1)
kurtInner = kurtosis(xInner); % 计算内圈故障数据的峰度
plot(tInner, xInner)
ylabel('Acceleration (g)')
title(['Inner Race Fault, kurtosis = ' num2str(kurtInner)])
xlim([0 0.1])
subplot(3, 1, 2)
kurtNormal = kurtosis(xNormal); % 计算健康数据的峰度
plot(tNormal, xNormal)
ylabel('Acceleration (g)')
title(['Normal, kurtosis = ' num2str(kurtNormal)])
xlim([0 0.1])
subplot(3, 1, 3)
kurtOuter = kurtosis(xOuter); % 计算外圈故障数据的峰度
plot(tOuter, xOuter)
xlabel('Time (s)')
ylabel('Acceleration (g)')
title(['Outer Race Fault, kurtosis = ' num2str(kurtOuter)])
xlim([0 0.1])
从上图可以看出,内圈故障数据有更大的峰度(即信号脉冲变化),这使得包络谱分析对它更有效。但是对于外圈故障而言,仅仅能看出轻微的BPFO谐波幅值重合,因为噪声较严重。对正常数据来说,甚至没有丝毫的幅值调制现象可以被观察到。所以综上,提取带有在BPFO处振幅调制的脉冲信号(或提高信噪比),在包络谱分析前属于关键的预处理步骤。下一节将介绍使用峰度图和谱峭度提取有着最高峰度值的信号,同时将包络谱分析方法应用在滤出的信号中
峰度图和谱峭度计算一定频率范围内的局部峰度。它们是找到拥有最高峰度(或最高信噪比)频率带信号的强力工具。在找到拥有最高峰度的频率带后,可以将带通滤波方法应用于原始信号来获得对包络谱分析更有帮助的强冲击信号
接下来,将该算法应用于批量训练数据中,该数据以文件集成数据的形式存储
工具箱中提供了有限的数据集,可将该数据集复制到当前文件夹下并开启写权限:
copyfile(...
fullfile(matlabroot, 'toolbox', 'predmaint', 'predmaintdemos', ...
'bearingFaultDiagnosis'), ...
'RollingElementBearingFaultDiagnosis-Data-master')
fileattrib(fullfile('RollingElementBearingFaultDiagnosis-Data-master', 'train_data', '*.mat'), '+w') % 训练数据
fileattrib(fullfile('RollingElementBearingFaultDiagnosis-Data-master', 'test_data', '*.mat'), '+w') % 测试数据
要想获得完整的数据集,可以使用本链接https://github.com/mathworks/RollingElementBearingFaultDiagnosis-Data来下载整个数据存储的压缩文件,可以用以下命令解压缩文件(压缩数据文件应位于当前程序位置):
if exist('RollingElementBearingFaultDiagnosis-Data-master.zip', 'file')
unzip('RollingElementBearingFaultDiagnosis-Data-master.zip')
end
本案例中的结果从完整的数据集中得到。完整数据集包括14个mat文件组成的训练数据集(2个正常数据集,4个内圈故障数据集,7个外圈故障数据集)以及6个mat文件组成的测试数据集(1个正常数据集,2个内圈故障数据集,3个外圈故障数据集)
通过将函数句柄分配给ReadFcn
以及WriteToMemberFcn
,文件集成数据能导航到文件中,并以所需的格式检索数据。例如,MFPT数据是一个struct
数据,存有振动信号gs,采样率sr等等。通过写入到readMFPTBearing
函数中,而不是返回一个struct
数据本身,可以使得文件集成数据可以直接返回振动信号gs,而不是struct
结构数据
fileLocation = fullfile('.', 'RollingElementBearingFaultDiagnosis-Data-master', 'train_data');
fileExtension = '.mat';
ensembleTrain = fileEnsembleDatastore(fileLocation, fileExtension); % 预测性维护工具箱中新增函数,用来存储数据文件,之后直接读取文件夹中符合文件扩展形式的所有数据文件
ensembleTrain.ReadFcn = @readMFPTBearing;
ensembleTrain.DataVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF"];
ensembleTrain.ConditionVariables = ["Label", "FileName"];
ensembleTrain.WriteToMemberFcn = @writeMFPTBearing;
ensembleTrain.SelectedVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF", "Label", "FileName"]
之后将其组合成tall
数据类型:
从上一节的分析看出,带通滤波后的包络谱在BPFO和BPFI幅值位置时,是用来指示轴承故障的指示器。因此,下一步可以提取从所有训练数据中提取这两个健康指示量。为了使算法结果具有更好的鲁棒性,可以在BPFO和BPFI附近设置窄频带(带宽=10Δf,Δf表示功率谱的频率分辨率),然后找到窄频带中的最大幅值。这一算法被写进了以下bearingFeatureExtraction
函数中。在之后的讨论中,包络谱在BPFO和BPFI附近的幅值分别被叫做“BPFOAmplitude”和“BPFIAmplitude”
function bearingFeatureExtraction(ensemble)
% 从轴承数据中提取指示量
data = read(ensemble);
x = data.gs{1};
fs = data.sr;
% 临界频率
BPFO = data.BPFO;
BPFI = data.BPFI;
level = 9;
[~, ~, ~, fc, ~, BW] = kurtogram(x, fs, level);
% 带通滤波后的包络谱
[pEnvpBpf, fEnvBpf] = envspectrum(x, fs, 'FilterOrder', 200, 'Band', [max([fc-BW/2 0]) min([fc+BW/2 0.999*fs/2])]);
deltaf = fEnvBpf(2) - fEnvBpf(1);
BPFIAmplitude = max(pEnvpBpf((fEnvBpf > (BPFI-5*deltaf)) & (fEnvBpf < (BPFI+5*deltaf))));
BPFOAmplitude = max(pEnvpBpf((fEnvBpf > (BPFO-5*deltaf)) & (fEnvBpf < (BPFO+5*deltaf))));
writeToLastMemberRead(ensemble, table(BPFIAmplitude, BPFOAmplitude, 'VariableNames', {'BPFIAmplitude', 'BPFOAmplitude'}));
end
% 以下被注释的代码用来进行并行处理
% ppool = gcp;
% n = numpartitions(ensembleTrain, ppool);
% parfor ct = 1:n
% subEnsembleTrain = partition(ensembleTrain, n, ct);
% reset(subEnsembleTrain);
% while hasdata(subEnsembleTrain)
% bearingFeatureExtraction(subEnsembleTrain);
% end
% end
ensembleTrain.DataVariables = [ensembleTrain.DataVariables; "BPFIAmplitude"; "BPFOAmplitude"];
reset(ensembleTrain)
while hasdata(ensembleTrain)
bearingFeatureExtraction(ensembleTrain)
end
一旦新的状态指示量被添加到文件集成数据中,可以指定SelectVariables
来从文件集成数据存储中读取相关数据,同时创建一个包含所提取状态指示量的特征表
ensembleTrain.SelectedVariables = ["BPFIAmplitude", "BPFOAmplitude", "Label"];
featureTableTrain = tall(ensembleTrain);
featureTableTrain = gather(featureTableTrain);
查看特征表
之后画散点图可视化创建的特征表
figure
gscatter(featureTableTrain.BPFIAmplitude, featureTableTrain.BPFOAmplitude, featureTableTrain.Label, [], 'ox+')
xlabel('BPFI Amplitude')
ylabel('BPFO Amplitude')
图中可以看到,BPFI和BPFO幅值的相对值可以作为判断不同故障类型的有效指标。因此,可以创造一个新的特征,即两个特征的对数比率,通过柱状图可画出不同故障类型时的该量的分布
% 创造新的特征,使用BPFI和BPFO两者的对数比率
featureTableTrain.IOLogRatio = log(featureTableTrain.BPFIAmplitude./featureTableTrain.BPFOAmplitude);
figure
hold on
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Inner Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Outer Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Normal"), 'BinWidth', 0.5)
plot([-1.5 -1.5 NaN 0.5 0.5], [0 3 NaN 0 3], 'k--')
hold off
ylabel('Count')
xlabel('log(BPFIAmplitude/BPFOAmplitude)')
legend('Inner Race Fault', 'Outer Race Fault', 'Normal', 'Classification Boundary')
柱状图显示了明显的划分三种轴承状态类型能力,因此BPFI和BPFO的对数比率是一个有效分类故障类型的特征。为了简化例子,简单的分类器可以通过图中虚线的绝对值确定,即对数值 ≤ − 1.5时判断为外圈故障, − 1.5<对数值 ≤ 0.5时判断为正常状态,对数值 ≥ 0.5时判断为内圈故障
现在,可对测试数据使用上述完整的工作流并评估上一部分中得到的分类器。以下测试数据包含1条正常数据集,2条内圈故障数据集,3条外圈故障数据集
fileLocation = fullfile('.', 'RollingElementBearingFaultDiagnosis-Data-master', 'test_data');
fileExtension = '.mat';
ensembleTest = fileEnsembleDatastore(fileLocation, fileExtension);
ensembleTest.ReadFcn = @readMFPTBearing;
ensembleTest.DataVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF"];
ensembleTest.ConditionVariables = ["Label", "FileName"];
ensembleTest.WriteToMemberFcn = @writeMFPTBearing;
ensembleTest.SelectedVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF", "Label", "FileName"]
设置变量名并提取特征
ensembleTest.DataVariables = [ensembleTest.DataVariables; "BPFIAmplitude"; "BPFOAmplitude"];
reset(ensembleTest)
while hasdata(ensembleTest)
bearingFeatureExtraction(ensembleTest)
end
将所有测试数据生成tall
数据并得到测试数据特征表
之后生成对数比率特征值并画出其柱状分布图
featureTableTest.IOLogRatio = log(featureTableTest.BPFIAmplitude./featureTableTest.BPFOAmplitude);
figure
hold on
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Inner Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Outer Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Normal"), 'BinWidth', 0.5)
histogram(featureTableTest.IOLogRatio(featureTableTest.Label=="Inner Race Fault"), 'BinWidth', 0.1)
histogram(featureTableTest.IOLogRatio(featureTableTest.Label=="Outer Race Fault"), 'BinWidth', 0.1)
histogram(featureTableTest.IOLogRatio(featureTableTest.Label=="Normal"), 'BinWidth', 0.1)
plot([-1.5 -1.5 NaN 0.5 0.5], [0 3 NaN 0 3], 'k--')
hold off
ylabel('Count')
xlabel('log(BPFIAmplitude/BPFOAmplitude)')
legend('Inner Race Fault - Train', 'Outer Race Fault - Train', 'Normal - Train', ...
'Inner Race Fault - Test', 'Outer Race Fault - Test', 'Normal - Test', ...
'Classification Boundary')
测试数据中BPFI和BPFO幅值计算得到的对数比率特征分布几乎同训练数据得到的该分布一致。说明上一节得到的这种简单的分类器在测试数据中可以实现完美分类。
当然需要注意的是,通常情况下单个特征不足以得到一个泛化能力好的分类器。更复杂的分类器可以通过将数据分割成多个块(创造更多的数据点),提取多个诊断相关的特征,然后根据它们的重要性排名选取一部分特征子集,最后使用机器学习工具箱中的Classification Learner APP训练模型来获得。案例6展示了这一工作流程的更多细节
本案例展示了如何使用峰度图、谱峭度和包络谱识别滚动轴承中不同的故障类型。然后将算法应用于批量数据集,这一过程显示了在BPFI和BPFO处带通滤波后的幅值,是两个应用于轴承诊断时重要的状态指示量
参考:https://ww2.mathworks.cn/help/predmaint/ug/Rolling-Element-Bearing-Fault-Diagnosis.html