英文地址https://ww2.mathworks.cn/help/predmaint/examples/Rolling-Element-Bearing-Fault-Diagnosis.html
這個例子展示了如何基于加速度信號對滾動軸承進(jìn)行故障診斷，特別是在存在來自其他機(jī)器部件的強(qiáng)掩蔽信號的情況下。該實(shí)例將演示如何應(yīng)用包絡(luò)譜分析和頻譜峰度來診斷軸承故障，并能夠擴(kuò)展到大數(shù)據(jù)應(yīng)用。
問題總覽
滾動軸承的局部故障可能發(fā)生在外圈、內(nèi)圈、保持架或滾動部件上。當(dāng)滾動元件撞擊外圈或內(nèi)圈的局部故障，或滾動元件上的故障撞擊外圈或內(nèi)圈[1]時，軸承和響應(yīng)傳感器之間的高頻共振被激發(fā)。下圖顯示了一個滾動部件在內(nèi)圈撞擊一個局部故障。問題是如何檢測和識別各種類型的故障。
機(jī)械故障預(yù)防技術(shù)(MFPT)挑戰(zhàn)數(shù)據(jù)
MFPT質(zhì)詢數(shù)據(jù)[4]包含從不同故障條件下的機(jī)器上收集的23個數(shù)據(jù)集。前20組數(shù)據(jù)采集于軸承試驗(yàn)臺，其中3組工況良好，3組外套圈故障為恒載，7組外套圈故障為各種載荷，7組內(nèi)套圈故障為各種載荷。剩下的3組數(shù)據(jù)來自真實(shí)的機(jī)器:一個油泵軸承，一個中速軸承和一個行星軸承。故障位置未知。在本例中，只使用從具有已知條件的試驗(yàn)臺收集的數(shù)據(jù)。
每個數(shù)據(jù)集包含加速度信號“gs”、采樣率“sr”、軸轉(zhuǎn)速“rate”、負(fù)載重量“l(fā)oad”，以及代表不同故障位置的四個關(guān)鍵頻率:球通頻率外圈(BPFO)、球通頻率內(nèi)圈(BPFI)、基本列車頻率(FTF)和球自旋頻率(BSF)。這是臨界頻率[1]的公式。
如圖所示，d為球直徑，d為節(jié)距直徑。變量fr軸轉(zhuǎn)速,n是滾動原件的數(shù)量,?是軸承接觸角[1]。
包絡(luò)譜分析用于軸承診斷
在MFPT數(shù)據(jù)集中，軸向速度是恒定的，因此不需要將跟蹤作為預(yù)處理步驟來消除軸向速度變化的影響。當(dāng)滾動元件撞擊外圈或內(nèi)圈的局部故障時，或當(dāng)滾動元件上的故障撞擊外圈或內(nèi)圈時，沖擊將調(diào)制相應(yīng)的臨界頻率，如BPFO、BPFI、FTF、BSF。因此，振幅解調(diào)產(chǎn)生的包絡(luò)信號傳遞了原始信號頻譜分析無法提供的更多診斷信息。以MFPT數(shù)據(jù)集中的內(nèi)部race故障信號為例。
dataInner = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
'predmaintdemos', 'bearingFaultDiagnosis', ...
'train_data', 'InnerRaceFault_vload_1.mat'));
在時域內(nèi)可視化原始的內(nèi)圈故障數(shù)據(jù)
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])
在頻域中可視化原始數(shù)據(jù)
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')
現(xiàn)在放大原始信號在低頻范圍的功率譜，以更近地觀察BPFI及其前幾次諧波的頻率響應(yīng)
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])
在BPFI及其諧波中看不到清晰的模式。原始信號的頻率分析不能提供有用的診斷信息。從時域數(shù)據(jù)可以看出，原始信號的幅值被調(diào)制到一定的頻率，調(diào)制的主頻約為1/0.009 Hz≈111 Hz。已知滾動單元內(nèi)圈局部故障即BPFI的頻率為118.875 Hz。這表明軸承可能存在內(nèi)圈故障
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)])
為了提取調(diào)制后的振幅，計(jì)算原始信號的包絡(luò)線，并將其顯示在底圖上
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')
現(xiàn)在計(jì)算包絡(luò)信號的功率譜，看看BPFI的頻率響應(yīng)及其諧波
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')
envelope spectrum包絡(luò)譜
BPFI Harmonics 內(nèi)圈頻率諧波
結(jié)果表明，大部分能量集中在BPFI及其諧波上。這表示軸承的內(nèi)圈故障，與數(shù)據(jù)的故障類型匹配。
將包絡(luò)譜分析應(yīng)用于其它故障類型
現(xiàn)在對正常數(shù)據(jù)和外圈故障數(shù)據(jù)重復(fù)相同的包絡(luò)譜分析。
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')
正如所料，正常軸承信號的包絡(luò)譜在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諧波也沒有明顯的峰值。包絡(luò)譜分析不能區(qū)分外圈故障軸承和正常軸承嗎?讓我們回過頭來看看不同條件下時域內(nèi)的信號。首先，讓我們再次在時域中形象化這些信號并計(jì)算它們的峰度。峰度是隨機(jī)變量的第四個標(biāo)準(zhǔn)化矩。它描述了信號的沖動性或隨機(jī)變量尾部的沉重。
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])
結(jié)果表明，內(nèi)圈故障信號具有較大的沖動性，使得包絡(luò)譜分析能夠有效地捕獲故障信號在BPFI處的特征。對于外環(huán)故障信號，BPFO處的調(diào)幅值略顯明顯，但被強(qiáng)噪聲掩蓋。正常信號不顯示任何調(diào)幅。在BPFO處進(jìn)行幅度調(diào)制提取脈沖信號(或提高信噪比)是包絡(luò)譜分析前的一個關(guān)鍵預(yù)處理步驟。下一節(jié)將介紹峰度圖和光譜峰度，提取峰度最高的信號，并對濾波后的信號進(jìn)行包絡(luò)譜分析。
峰圖和光譜峰度用于波段選擇
峰度圖和頻譜峰度計(jì)算在頻帶內(nèi)的局部峰度。它們是定位具有最高峰度(或最高信噪比)[2]的頻帶的強(qiáng)大工具。在確定峰度最高的頻段后，可以對原始信號進(jìn)行帶通濾波，獲得更強(qiáng)的脈沖信號進(jìn)行包絡(luò)譜分析。
level = 9;
figure
kurtogram(xOuter, fsOuter, level)
峰度圖表明，以2.67 kHz為中心，帶寬為0.763 kHz的頻帶峰度最高，為2.71?，F(xiàn)在使用最優(yōu)窗口長度所建議的kurtogram來計(jì)算光譜峰度。
figure
wc = 128;
pkurtosis(xOuter, fsOuter, wc)
為了使光譜圖上的頻帶可視化，計(jì)算光譜圖并將光譜峰度放在一邊。另一種解釋譜峰度的方法是，高譜峰度值表示對應(yīng)頻率下功率的高方差，這使得譜峰度成為定位信號[3]非平穩(wěn)分量的有用工具。
helperSpectrogramAndSpectralKurtosis(xOuter, fsOuter, level)
通過對中心頻率和帶寬建議的信號進(jìn)行帶通濾波，可以增強(qiáng)峰度，恢復(fù)外環(huán)故障的調(diào)制幅值。
[~, ~, ~, 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')
可以看出，帶通濾波后峰度值增大?，F(xiàn)在在頻域中形象化包絡(luò)信號。
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')
結(jié)果表明，通過對原始信號進(jìn)行帶通濾波，利用峰度和頻譜峰度所建議的頻帶，包絡(luò)譜分析能夠較好地揭示BPFO的故障特征及其諧波特征。
Batch Process
現(xiàn)在，讓我們使用文件集成數(shù)據(jù)存儲將該算法應(yīng)用于一批訓(xùn)練數(shù)據(jù)。工具箱中提供了數(shù)據(jù)集的有限部分。將數(shù)據(jù)集復(fù)制到當(dāng)前文件夾，并啟用寫入權(quán)限:
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')
對于完整的數(shù)據(jù)集，請轉(zhuǎn)到這個鏈接https://github.com/mathworks/rollingelementbearingfaultdiagnostics - data，以zip文件的形式下載整個存儲庫，并將其保存在與活動腳本相同的目錄中。使用以下命令解壓縮文件:
if exist('RollingElementBearingFaultDiagnosis-Data-master.zip', 'file')
unzip('RollingElementBearingFaultDiagnosis-Data-master.zip')
end
本例中的結(jié)果是從完整的數(shù)據(jù)集生成的。整個數(shù)據(jù)集包含一個包含14個mat文件的訓(xùn)練數(shù)據(jù)集(2個normal, 4個inner race fault, 7個outer race fault)和一個包含6個mat文件的測試數(shù)據(jù)集(1個normal, 2個inner race fault, 3個outer race fault)。通過將函數(shù)句柄分配給ReadFcn和WriteToMemberFcn，文件集成數(shù)據(jù)存儲將能夠?qū)Ш降轿募校运璧母袷綑z索數(shù)據(jù)。例如，MFPT數(shù)據(jù)具有存儲振動信號gs、采樣率sr等的結(jié)構(gòu)軸承。不返回軸承結(jié)構(gòu)本身，而是編寫readmfpt支座函數(shù)，使文件集成數(shù)據(jù)存儲返回軸承數(shù)據(jù)結(jié)構(gòu)內(nèi)部的振動信號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"]
ensembleTrainTable = tall(ensembleTrain)
從分析的最后一部分可以看出，BPFO和BPFI處的帶通濾波包絡(luò)譜振幅是軸承故障診斷的兩個條件指標(biāo)。因此，下一步就是從所有的訓(xùn)練數(shù)據(jù)中提取這兩個狀態(tài)指標(biāo)。使算法更加健壯,設(shè)置一個窄帶(帶寬= 10Δf,Δf功率譜的頻率分辨率)周圍BPFO BPFI,然后找到最大振幅在這個狹窄的區(qū)間。該算法包含在下面的bearingFeatureExtraction函數(shù)中。注意，在示例的其余部分中，BPFI和BPFO周圍的包絡(luò)頻譜振幅稱為“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
一旦將新的條件指示器添加到文件集成數(shù)據(jù)存儲中，指定SelectedVariables從文件集成數(shù)據(jù)存儲中讀取相關(guān)數(shù)據(jù)，并創(chuàng)建一個包含提取的條件指示器的特征表。
ensembleTrain.SelectedVariables = ["BPFIAmplitude", "BPFOAmplitude", "Label"];
featureTableTrain = tall(ensembleTrain);
featureTableTrain = gather(featureTableTrain);
featureTableTrain
可視化已創(chuàng)建的特性表
figure
gscatter(featureTableTrain.BPFIAmplitude, featureTableTrain.BPFOAmplitude, featureTableTrain.Label, [], 'ox+')
xlabel('BPFI Amplitude')
ylabel('BPFO Amplitude')
BPFI幅值與BPFO幅值的相對值可以作為判斷不同故障類型的有效指標(biāo)。這里創(chuàng)建了一個新特性，即兩個現(xiàn)有特性的日志比率，并以按不同故障類型分組的直方圖顯示。
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')
使用測試數(shù)據(jù)集進(jìn)行驗(yàn)證
現(xiàn)在，讓我們將工作流應(yīng)用到一個測試數(shù)據(jù)集中，并驗(yàn)證上一節(jié)中獲得的分類器。測試數(shù)據(jù)包括1個正常數(shù)據(jù)集、2個內(nèi)圈故障數(shù)據(jù)集和3個外圈故障數(shù)據(jù)集。
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')
測試數(shù)據(jù)集的BPFI和BPFO振幅的對數(shù)比與訓(xùn)練數(shù)據(jù)集的對數(shù)比分布一致。最后一節(jié)得到的樸素分類器在測試數(shù)據(jù)集上達(dá)到了很好的精度，需要注意的是單特征通常不足以得到一個很好的泛化分類器。通過將數(shù)據(jù)分成多個部分(創(chuàng)建更多的數(shù)據(jù)點(diǎn))，提取多個診斷相關(guān)特征，根據(jù)特征的重要程度選擇子集，使用Classification Learner App在Statistics &中訓(xùn)練各種分類器，可以得到更加復(fù)雜的分類器機(jī)器學(xué)習(xí)工具。有關(guān)此工作流的詳細(xì)信息，請參考示例“使用Simulink生成故障數(shù)據(jù)”。
總結(jié)
這個例子展示了如何使用kurtogram，頻譜峰度和包絡(luò)譜來識別不同類型的滾動軸承故障。將該算法應(yīng)用于磁盤上的一批數(shù)據(jù)集，表明帶通濾波包絡(luò)譜在BPFI和BPFO處的振幅是軸承診斷的兩個重要條件指標(biāo)。
