MATLAB环境下基于包络谱和谱峭度的一维振动信号分析

依旧是信号处理相关的东西，本文再次讲解如何应用包络谱和谱峭度分析一维振动信号进而诊断轴承故障，运行环境为MATLAB R2021B。

面包多第三方代码如下：

https://mianbaoduo.com/o/bread/Y5mTlp5q

滚动轴承的局部故障可能发生在外圈、内圈、保持架或滚动体中。 当滚动体撞击外圈或内圈上的局部故障，或者滚动体上的故障经过外圈或内圈时，轴承和传感器之间的高频共振会被激发， 下图显示了内圈处发生局部故障。

故障频率计算如下

包络谱分析

以轴承内圈故障信号dataInner 为例，在时域中可视化原始内圈故障数据

xInner = Not Found;
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])

信号的功率谱

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

放大低频范围内原始信号的功率谱，仔细观察 BPFI 的频率响应及其前几个谐波

在 BPFI 及其谐波处看不到清晰的谱线图。 看时域数据，观察到原始信号的幅度在一定的频率上被调制，调制的主频率在1/0.009Hz~111Hz左右。 已知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')

计算包络信号的功率谱，并查看 BPFI 的频率响应及其谐波

结果表明，大部分能量集中在 BPFI 及其谐波上， 这表明轴承发生了内圈故障，与数据的故障类型相匹配。

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

对正常数据dataNormal进行包络谱分析

正常轴承振动信号的包络谱在 BPFO 或 BPFI 处没有显示任何显著的峰值。

对外圈故障数据dataOuter进行包络谱分析

对于外圈故障信号，BPFO及其谐波也没有明显的峰值。

在时域中可视化信号并计算相应的峭度， 峭度是随机变量的4阶标准矩，表征了信号的冲击性

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

结果表明，内圈故障信号具有较大的冲击性，包络谱分析有效地捕捉到了BPFI处的故障特征。 对于外圈故障信号，BPFO 处的幅度调制略微明显，但被强噪声掩盖。 在 BPFO 处通过幅度调制提取脉冲信号（或提高信噪比）是包络谱分析之前的关键预处理步骤。

用于频带选择的峭度图和谱峭度图

峭度图和谱峭度计算频带内的局部峭度，在确定具有最高峭度的频带后，可以对原始信号应用带通滤波器，以获得更具冲击性的信号，用于包络谱分析。

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

峭度图表明以 2.67 kHz 为中心、带宽为 0.763 kHz 的频带具有 2.71 的最高峭度。使用峭度图建议的最佳窗口长度来计算谱峭度。

计算时频谱图并将谱峭度放在一边

通过使用建议的中心频率和带宽对信号进行带通滤波，可以增强峭度，并且可以检索到外圈故障的调制幅度

[~, ~, ~, 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')

可以看出，经过带通滤波后峭度值有所增加，可视化包络谱

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

此外，BPFO 和 BPFI 处的带通滤波包络谱幅度是轴承故障诊断的两个条件指标。 因此，下一步是从所有训练数据中提取两个条件指标。 为了使算法更加鲁棒，在BPFO和BPFI周围设置一个窄带，然后在这个窄带内找到最大幅度。 该算法包含在BearingFeatureExtraction 函数中。注意，BPFI 和 BPFO 周围的包络谱振幅被称为“BPFIAmplitude”和“BPFOAmplitude”。BPFI Amplitude 和 BPFO Amplitude 的相对值可能是不同故障类型的有效指标。 创建一个新特征，它是两个现有特征的对数比，并在按不同故障类型分组的直方图中可视化。

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

直方图清楚地显示了3种不同轴承故障之间的区别。 BPFI 和 BPFO 幅度之间的对数比是轴承故障分类的有效特征。

使用测试数据集进行验证，这里的测试数据包含 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
ensembleTest.SelectedVariables = ["BPFIAmplitude", "BPFOAmplitude", "Label"];
featureTableTest = tall(ensembleTest);
featureTableTest = gather(featureTableTest);

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 幅度的对数比与训练数据集的对数比显示出一致的分布。 需要注意的是，单个特征通常不足以得到一个泛化好的分类器。 可以通过将数据分成多个部分（以创建更多数据点）来获得更复杂的分类器，提取多个与诊断相关的特征，按重要性等级选择特征子集，并使用 Statistics & Machine 中的分类学习器应用程序训练各种分类器。

最后，BPFI 和 BPFO 处的带通滤波包络谱的幅度是轴承诊断的两个重要条件指标。