MATLAB環境下基于包絡譜和譜峭度的一維振動信号分析

依舊是信号處理相關的東西，本文再次講解如何應用包絡譜和譜峭度分析一維振動信号進而診斷軸承故障，運作環境為MATLAB R2021B。

面包多第三方代碼：

🍞正在為您運送作品詳情

滾動軸承的局部故障可能發生在外圈、内圈、保持架或滾動體中。 當滾動體撞擊外圈或内圈上的局部故障，或者滾動體上的故障經過外圈或内圈時，軸承和傳感器之間的高頻共振會被激發， 下圖顯示了内圈處發生局部故障。

故障頻率計算如下

包絡譜分析

以軸承内圈故障信号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 處的帶通濾波包絡譜的幅度是軸承診斷的兩個重要條件名額。