英文地址https://ww2.mathworks.cn/help/predmaint/examples/Rolling-Element-Bearing-Fault-Diagnosis.html
这个例子展示了如何基于加速度信号对滚动轴承进行故障诊断，特别是在存在来自其他机器部件的强掩蔽信号的情况下。该实例将演示如何应用包络谱分析和频谱峰度来诊断轴承故障，并能够扩展到大数据应用。

问题总览

滚动轴承的局部故障可能发生在外圈、内圈、保持架或滚动部件上。当滚动元件撞击外圈或内圈的局部故障，或滚动元件上的故障撞击外圈或内圈[1]时，轴承和响应传感器之间的高频共振被激发。下图显示了一个滚动部件在内圈撞击一个局部故障。问题是如何检测和识别各种类型的故障。

image.png

机械故障预防技术(MFPT)挑战数据

MFPT质询数据[4]包含从不同故障条件下的机器上收集的23个数据集。前20组数据采集于轴承试验台，其中3组工况良好，3组外套圈故障为恒载，7组外套圈故障为各种载荷，7组内套圈故障为各种载荷。剩下的3组数据来自真实的机器:一个油泵轴承，一个中速轴承和一个行星轴承。故障位置未知。在本例中，只使用从具有已知条件的试验台收集的数据。

每个数据集包含加速度信号“gs”、采样率“sr”、轴转速“rate”、负载重量“load”，以及代表不同故障位置的四个关键频率:球通频率外圈(BPFO)、球通频率内圈(BPFI)、基本列车频率(FTF)和球自旋频率(BSF)。这是临界频率[1]的公式。

image.png

如图所示，d为球直径，d为节距直径。变量fr轴转速,n是滚动原件的数量,ϕ是轴承接触角[1]。

包络谱分析用于轴承诊断

在MFPT数据集中，轴向速度是恒定的，因此不需要将跟踪作为预处理步骤来消除轴向速度变化的影响。当滚动元件撞击外圈或内圈的局部故障时，或当滚动元件上的故障撞击外圈或内圈时，冲击将调制相应的临界频率，如BPFO、BPFI、FTF、BSF。因此，振幅解调产生的包络信号传递了原始信号频谱分析无法提供的更多诊断信息。以MFPT数据集中的内部race故障信号为例。

dataInner = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingFaultDiagnosis', ...
    'train_data', 'InnerRaceFault_vload_1.mat'));

在时域内可视化原始的内圈故障数据

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])

image.png

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')

image.png

figure
plot(fpInner, pInner)
ncomb = 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])

image.png

在BPFI及其谐波中看不到清晰的模式。原始信号的频率分析不能提供有用的诊断信息。从时域数据可以看出，原始信号的幅值被调制到一定的频率，调制的主频约为1/0.009 Hz≈111 Hz。已知滚动单元内圈局部故障即BPFI的频率为118.875 Hz。这表明轴承可能存在内圈故障

figure
subplot(2, 1, 1)
plot(tInner, xInner)
xlim([0.04 0.06])
title('Raw Signal: Inner Race Fault')
ylabel('Acceleration (g)')
annotation('doublearrow', [0.37 0.71], [0.8 0.8])
text(0.047, 20, ['0.009 sec \approx 1/BPFI, BPFI = ' num2str(dataInner.BPFI)])

为了提取调制后的振幅，计算原始信号的包络线，并将其显示在底图上

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')

image.png

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')

image.png

将包络谱分析应用于其它故障类型

现在对正常数据和外圈故障数据重复相同的包络谱分析。

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')

image.png

正如所料，正常轴承信号的包络谱在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')

image.png

对于外圈故障信号，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])

image.png

结果表明，内圈故障信号具有较大的冲动性，使得包络谱分析能够有效地捕获故障信号在BPFI处的特征。对于外环故障信号，BPFO处的调幅值略显明显，但被强噪声掩盖。正常信号不显示任何调幅。在BPFO处进行幅度调制提取脉冲信号(或提高信噪比)是包络谱分析前的一个关键预处理步骤。下一节将介绍峰度图和光谱峰度，提取峰度最高的信号，并对滤波后的信号进行包络谱分析。

峰图和光谱峰度用于波段选择

峰度图和频谱峰度计算在频带内的局部峰度。它们是定位具有最高峰度(或最高信噪比)[2]的频带的强大工具。在确定峰度最高的频段后，可以对原始信号进行带通滤波，获得更强的脉冲信号进行包络谱分析。

level = 9;
figure
kurtogram(xOuter, fsOuter, level)

image.png

峰度图表明，以2.67 kHz为中心，带宽为0.763 kHz的频带峰度最高，为2.71。现在使用最优窗口长度所建议的kurtogram来计算光谱峰度。

figure
wc = 128;
pkurtosis(xOuter, fsOuter, wc)

image.png

为了使光谱图上的频带可视化，计算光谱图并将光谱峰度放在一边。另一种解释谱峰度的方法是，高谱峰度值表示对应频率下功率的高方差，这使得谱峰度成为定位信号[3]非平稳分量的有用工具。

helperSpectrogramAndSpectralKurtosis(xOuter, fsOuter, level)

image.png

通过对中心频率和带宽建议的信号进行带通滤波，可以增强峰度，恢复外环故障的调制幅值。

[~, ~, ~, fc, ~, BW] = kurtogram(xOuter, fsOuter, level);

bpf = designfilt('bandpassfir', 'FilterOrder', 200, 'CutoffFrequency1', fc-BW/2, ...
    'CutoffFrequency2', fc+BW/2, 'SampleRate', fsOuter);
xOuterBpf = filter(bpf, xOuter);
[pEnvOuterBpf, fEnvOuterBpf, xEnvOuterBpf, tEnvBpfOuter] = envspectrum(xOuter, fsOuter, ...
    'FilterOrder', 200, 'Band', [fc-BW/2 fc+BW/2]);

figure
subplot(2, 1, 1)
plot(tOuter, xOuter, tEnvOuter, xEnvOuter)
ylabel('Acceleration (g)')
title(['Raw Signal: Outer Race Fault, kurtosis = ', num2str(kurtOuter)])
xlim([0 0.1])
legend('Signal', 'Envelope')

subplot(2, 1, 2)
kurtOuterBpf = kurtosis(xOuterBpf);
plot(tOuter, xOuterBpf, tEnvBpfOuter, xEnvOuterBpf)
ylabel('Acceleration (g)')
xlim([0 0.1])
xlabel('Time (s)')
title(['Bandpass Filtered Signal: Outer Race Fault, kurtosis = ', num2str(kurtOuterBpf)])
legend('Signal', 'Envelope')

image.png

可以看出，带通滤波后峰度值增大。现在在频域中形象化包络信号。

figure
plot(fEnvOuterBpf, pEnvOuterBpf);
ncomb = 10;
helperPlotCombs(ncomb, dataOuter.BPFO)
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum of Bandpass Filtered Signal: Outer Race Fault ')
legend('Envelope Spectrum', 'BPFO Harmonics')

image.png

结果表明，通过对原始信号进行带通滤波，利用峰度和频谱峰度所建议的频带，包络谱分析能够较好地揭示BPFO的故障特征及其谐波特征。

Batch Process

现在，让我们使用文件集成数据存储将该算法应用于一批训练数据。工具箱中提供了数据集的有限部分。将数据集复制到当前文件夹，并启用写入权限:

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/rollingelementbearingfaultdiagnostics - data，以zip文件的形式下载整个存储库，并将其保存在与活动脚本相同的目录中。使用以下命令解压缩文件:

if exist('RollingElementBearingFaultDiagnosis-Data-master.zip', 'file')
    unzip('RollingElementBearingFaultDiagnosis-Data-master.zip')
end

本例中的结果是从完整的数据集生成的。整个数据集包含一个包含14个mat文件的训练数据集(2个normal, 4个inner race fault, 7个outer race fault)和一个包含6个mat文件的测试数据集(1个normal, 2个inner race fault, 3个outer race fault)。通过将函数句柄分配给ReadFcn和WriteToMemberFcn，文件集成数据存储将能够导航到文件中，以所需的格式检索数据。例如，MFPT数据具有存储振动信号gs、采样率sr等的结构轴承。不返回轴承结构本身，而是编写readmfpt支座函数，使文件集成数据存储返回轴承数据结构内部的振动信号gs。

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"]

image.png

ensembleTrainTable = tall(ensembleTrain)

image.png

从分析的最后一部分可以看出，BPFO和BPFI处的带通滤波包络谱振幅是轴承故障诊断的两个条件指标。因此，下一步就是从所有的训练数据中提取这两个状态指标。使算法更加健壮,设置一个窄带(带宽= 10Δf,Δf功率谱的频率分辨率)周围BPFO BPFI,然后找到最大振幅在这个狭窄的区间。该算法包含在下面的bearingFeatureExtraction函数中。注意，在示例的其余部分中，BPFI和BPFO周围的包络频谱振幅称为“BPFIAmplitude”和“BPFOAmplitude”。

% To process the data in parallel, use the following code
% 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

一旦将新的条件指示器添加到文件集成数据存储中，指定SelectedVariables从文件集成数据存储中读取相关数据，并创建一个包含提取的条件指示器的特征表。

ensembleTrain.SelectedVariables = ["BPFIAmplitude", "BPFOAmplitude", "Label"];
featureTableTrain = tall(ensembleTrain);
featureTableTrain = gather(featureTableTrain);

image.png

image.png

image.png

featureTableTrain

image.png

可视化已创建的特性表

figure
gscatter(featureTableTrain.BPFIAmplitude, featureTableTrain.BPFOAmplitude, featureTableTrain.Label, [], 'ox+')
xlabel('BPFI Amplitude')
ylabel('BPFO Amplitude')

image.png

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')

image.png

image.png

使用测试数据集进行验证

现在，让我们将工作流应用到一个测试数据集中，并验证上一节中获得的分类器。测试数据包括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"]

image.png

ensembleTest.DataVariables = [ensembleTest.DataVariables; "BPFIAmplitude"; "BPFOAmplitude"];
reset(ensembleTest)
while hasdata(ensembleTest)
    bearingFeatureExtraction(ensembleTest)
end

ensembleTest.SelectedVariables = ["BPFIAmplitude", "BPFOAmplitude", "Label"];
featureTableTest = tall(ensembleTest);
featureTableTest = gather(featureTableTest);

image.png

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')

image.png

测试数据集的BPFI和BPFO振幅的对数比与训练数据集的对数比分布一致。最后一节得到的朴素分类器在测试数据集上达到了很好的精度，需要注意的是单特征通常不足以得到一个很好的泛化分类器。通过将数据分成多个部分(创建更多的数据点)，提取多个诊断相关特征，根据特征的重要程度选择子集，使用Classification Learner App在Statistics &中训练各种分类器，可以得到更加复杂的分类器机器学习工具。有关此工作流的详细信息，请参考示例“使用Simulink生成故障数据”。

总结

这个例子展示了如何使用kurtogram，频谱峰度和包络谱来识别不同类型的滚动轴承故障。将该算法应用于磁盘上的一批数据集，表明带通滤波包络谱在BPFI和BPFO处的振幅是轴承诊断的两个重要条件指标。