rolling element bearing fault diagnosis -凯发k8网页登录

machinery failure prevention technology (mfpt) challenge data

mfpt challenge data [4] contains 23 data sets collected from machines under various fault conditions. the first 20 data sets are collected from a bearing test rig, with 3 under good conditions, 3 with outer race faults under constant load, 7 with outer race faults under various loads, and 7 with inner race faults under various loads. the remaining 3 data sets are from real-world machines: an oil pump bearing, an intermediate speed bearing, and a planet bearing. the fault locations are unknown. in this example, only the data collected from the test rig with known conditions are used.

each data set contains an acceleration signal "gs", sampling rate "sr", shaft speed "rate", load weight "load", and four critical frequencies representing different fault locations: ballpass frequency outer race (bpfo), ballpass frequency inner race (bpfi), fundamental train frequency (ftf), and ball spin frequency (bsf). here are the formulae for those critical frequencies [1].

as shown in the figure, d is the ball diameter, d is the pitch diameter. the variable fr is the shaft speed, n is the number of rolling elements, ϕ is the bearing contact angle [1].

envelope spectrum analysis for bearing diagnosis

in the mfpt data set, the shaft speed is constant, hence there is no need to perform order tracking as a pre-processing step to remove the effect of shaft speed variations.

when rolling elements hit the local faults at outer or inner races, or when faults on the rolling element hit the outer or inner races, the impact will modulate the corresponding critical frequencies, e.g. bpfo, bpfi, ftf, bsf. therefore, the envelope signal produced by amplitude demodulation conveys more diagnostic information that is not available from spectrum analysis of the raw signal. take an inner race fault signal in the mfpt dataset as an example.

datainner = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingfaultdiagnosis', ...
    'train_data', 'innerracefault_vload_1.mat'));

visualize the raw inner race fault data in the time domain.

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

visualize the raw data in frequency domain.

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

now zoom in the power spectrum of the raw signal in low frequency range to take a closer look at the frequency response at bpfi and its first several harmonics.

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

no clear pattern is visible at bpfi and its harmonics. frequency analysis on the raw signal does not provide useful diagnosis information.

looking at the time-domain data, it is observed that the amplitude of the raw signal is modulated at a certain frequency, and the main frequency of the modulation is around 1/0.009 hz ≈ 111 hz. it is known that the frequency the rolling element hitting a local fault at the inner race, that is bpfi, is 118.875 hz. this indicates that the bearing potentially has an inner race fault.

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

to extract the modulated amplitude, compute the envelope of the raw signal, and visualize it on the bottom subplot.

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

now compute the power spectrum of the envelope signal and take a look at the frequency response at bpfi and its harmonics.

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

it is shown that most of the energy is focused at bpfi and its harmonics. that indicates an inner race fault of the bearing, which matches the fault type of the data.

applying envelope spectrum analysis to other fault types

now repeat the same envelope spectrum analysis on normal data and outer race fault data.

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

as expected, the envelope spectrum of a normal bearing signal does not show any significant peaks at bpfo or 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')

for an outer race fault signal, there are no clear peaks at bpfo harmonics either. does envelope spectrum analysis fail to differentiate bearing with outer race fault from healthy bearings? let's take a step back and look at the signals in time domain under different conditions again.

first let's visualize the signals in time domain again and calculate their kurtosis. kurtosis is the fourth standardized moment of a random variable. it characterizes the impulsiveness of the signal or the heaviness of the random variable's tail.

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

it is shown that inner race fault signal has significantly larger impulsiveness, making envelope spectrum analysis capture the fault signature at bpfi effectively. for an outer race fault signal, the amplitude modulation at bpfo is slightly noticeable, but it is masked by strong noise. the normal signal does not show any amplitude modulation. extracting the impulsive signal with amplitude modulation at bpfo (or enhancing the signal-to-noise ratio) is a key preprocessing step before envelope spectrum analysis. the next section will introduce kurtogram and spectral kurtosis to extract the signal with highest kurtosis, and perform envelope spectrum analysis on the filtered signal.

kurtogram and spectral kurtosis for band selection

kurtogram and spectral kurtosis compute kurtosis locally within frequency bands. they are powerful tools to locate the frequency band that has the highest kurtosis (or the highest signal-to-noise ratio) [2]. after pinpointing the frequency band with the highest kurtosis, a bandpass filter can be applied to the raw signal to obtain a more impulsive signal for envelope spectrum analysis.

level = 9;
figure
kurtogram(xouter, fsouter, level)

the kurtogram indicates that the frequency band centered at 2.67 khz with a 0.763 khz bandwidth has the highest kurtosis of 2.71.

now use the optimal window length suggested by the kurtogram to compute the spectral kurtosis.

figure
wc = 128;
pkurtosis(xouter, fsouter, wc)

to visualize the frequency band on a spectrogram, compute the spectrogram and place the spectral kurtosis on the side. to interpret the spectral kurtosis in another way, high spectral kurtosis values indicates high variance of power at the corresponding frequency, which makes spectral kurtosis a useful tool to locate nonstationary components of the signal [3].

helperspectrogramandspectralkurtosis(xouter, fsouter, level)

by bandpass filtering the signal with the suggested center frequency and bandwidth, the kurtosis can be enhanced and the modulated amplitude of the outer race fault can be retrieved.

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

it can be seen that the kurtosis value is increased after bandpass filtering. now visualize the envelope signal in frequency domain.

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

it is shown that by bandpass filtering the raw signal with the frequency band suggested by kurtogram and spectral kurtosis, the envelope spectrum analysis is able to reveal the fault signature at bpfo and its harmonics.

batch process

now let's apply the algorithm to a batch of training data using a file ensemble datastore.

a limited portion of the dataset is available in the toolbox. copy the dataset to the current folder and enable the write permission:

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

for the full dataset, go to this link to download the entire repository as a zip file and save it in the same directory as the live script. unzip the file using this command:

if exist('rollingelementbearingfaultdiagnosis-data-master.zip', 'file')
    unzip('rollingelementbearingfaultdiagnosis-data-master.zip')
end

the results in this example are generated from the full dataset. the full dataset contains a training dataset with 14 mat files (2 normal, 4 inner race fault, 7 outer race fault) and a testing dataset with 6 mat files (1 normal, 2 inner race fault, 3 outer race fault).

by assigning function handles to readfcn and writetomemberfcn, the file ensemble datastore will be able to navigate into the files to retrieve data in the desired format. for example, the mfpt data has a structure bearing that stores the vibration signal gs, sampling rate sr, and so on. instead of returning the bearing structure itself the readmfptbearing function is written so that file ensemble datastore returns the vibration signal gs inside of the bearing data structure.

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

ensembletrain = 
  fileensembledatastore with properties:
                 readfcn: @readmfptbearing
        writetomemberfcn: @writemfptbearing
           datavariables: [8×1 string]
    independentvariables: [0×0 string]
      conditionvariables: [2×1 string]
       selectedvariables: [10×1 string]
                readsize: 1
              nummembers: 14
          lastmemberread: [0×0 string]
                   files: [14×1 string]

ensembletraintable = tall(ensembletrain)

starting parallel pool (parpool) using the 'local' profile ...
connected to 6 workers.
ensembletraintable =
  m×10 tall table
           gs             sr      rate    load     bpfo      bpfi      ftf       bsf           label                   filename        
    _________________    _____    ____    ____    ______    ______    ______    _____    __________________    ________________________
    [146484×1 double]    48828     25       0     81.125    118.88    14.838    63.91    "inner race fault"    "innerracefault_vload_1"
    [146484×1 double]    48828     25      50     81.125    118.88    14.838    63.91    "inner race fault"    "innerracefault_vload_2"
    [146484×1 double]    48828     25     100     81.125    118.88    14.838    63.91    "inner race fault"    "innerracefault_vload_3"
    [146484×1 double]    48828     25     150     81.125    118.88    14.838    63.91    "inner race fault"    "innerracefault_vload_4"
    [146484×1 double]    48828     25     200     81.125    118.88    14.838    63.91    "inner race fault"    "innerracefault_vload_5"
    [585936×1 double]    97656     25     270     81.125    118.88    14.838    63.91    "outer race fault"    "outerracefault_1"      
    [585936×1 double]    97656     25     270     81.125    118.88    14.838    63.91    "outer race fault"    "outerracefault_2"      
    [146484×1 double]    48828     25      25     81.125    118.88    14.838    63.91    "outer race fault"    "outerracefault_vload_1"
            :              :       :       :        :         :         :         :              :                        :
            :              :       :       :        :         :         :         :              :                        :

from the last section of analysis, notice that the bandpass filtered envelope spectrum amplitudes at bpfo and bpfi are two condition indicators for bearing fault diagnosis. therefore, the next step is to extract the two condition indicators from all the training data. to make the algorithm more robust, set a narrow band (bandwidth = 10δf, where δf is the frequency resolution of the power spectrum) around bpfo and bpfi, and then find the maximum amplitude inside this narrow band. the algorithm is contained in the bearingfeatureextraction function listed below. note that the envelope spectrum amplitudes around bpfi and bpfo are referred to as "bpfiamplitude" and "bpfoamplitude" in the rest of the example.

% 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

once the new condition indicators are added into the file ensemble datastore, specify selectedvariables to read the relevant data from the file ensemble datastore, and create a feature table containing the extracted condition indicators.

ensembletrain.selectedvariables = ["bpfiamplitude", "bpfoamplitude", "label"];
featuretabletrain = tall(ensembletrain);
featuretabletrain = gather(featuretabletrain);

evaluating tall expression using the parallel pool 'local':
- pass 1 of 1: 0% complete
evaluation 0% complete

- pass 1 of 1: completed in 3 sec
evaluation completed in 3 sec

featuretabletrain

featuretabletrain=14×3 table
    bpfiamplitude    bpfoamplitude          label       
    _____________    _____________    __________________
        0.33918         0.082296      "inner race fault"
        0.31488         0.026599      "inner race fault"
        0.52356         0.036609      "inner race fault"
        0.52899         0.028381      "inner race fault"
        0.13515         0.012337      "inner race fault"
       0.004024          0.03574      "outer race fault"
      0.0044918           0.1835      "outer race fault"
      0.0074993          0.30166      "outer race fault"
       0.013662          0.12468      "outer race fault"
      0.0070963          0.28215      "outer race fault"
      0.0060772          0.35241      "outer race fault"
       0.011244          0.17975      "outer race fault"
      0.0036798        0.0050208      "normal"          
        0.00359        0.0069449      "normal"          

visualize the feature table that has been created.

figure
gscatter(featuretabletrain.bpfiamplitude, featuretabletrain.bpfoamplitude, featuretabletrain.label, [], 'ox ')
xlabel('bpfi amplitude')
ylabel('bpfo amplitude')

the relative values of bpfi amplitude and bpfo amplitude might be an effective indicator of different fault types. here a new feature is created, which is the log ratio of the two existing features, and is visualized in a histogram grouped by different fault types.

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

the histogram shows a clear separation among the three different bearing conditions. the log ratio between the bpfi and bpfo amplitudes is a valid feature to classify bearing faults. to simplify the example, a very simple classifier is derived: if log(bpfiamplitudebpfoamplitude)≤-1.5, the bearing has an outer race fault; if -1.5<log(bpfiamplitudebpfoamplitude)≤0.5, the bearing is normal; and if log(bpfiamplitudebpfoamplitude)>0.5, the bearing has an inner race fault.

validation using test data sets

now, let's apply the workflow to a test data set and validate the classifier obtained in the last section. here the test data contains 1 normal data set, 2 inner race fault data sets, and 3 outer race fault data sets.

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 = 
  fileensembledatastore with properties:
                 readfcn: @readmfptbearing
        writetomemberfcn: @writemfptbearing
           datavariables: [8×1 string]
    independentvariables: [0×0 string]
      conditionvariables: [2×1 string]
       selectedvariables: [10×1 string]
                readsize: 1
              nummembers: 6
          lastmemberread: [0×0 string]
                   files: [6×1 string]

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

evaluating tall expression using the parallel pool 'local':
- pass 1 of 1: completed in 1 sec
evaluation completed in 1 sec

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

the log ratio of bpfi and bpfo amplitudes from test data sets shows consistent distribution with the log ratio from training data sets. the naive classifier obtained in the last section achieved perfect accuracy on the test data set.

it should be noted that single feature is usually not enough to get a classifier that generalizes well. more sophisticated classifiers could be obtained by dividing the data into multiple pieces (to create more data points), extract multiple diagnosis related features, select a subset of features by their importance ranks, and train various classifiers using the classification learner app in statistics & machine learning toolbox. for more details of this workflow, please refer to the example "using simulink to generate fault data".

summary

this example shows how to use kurtogram, spectral kurtosis and envelope spectrum to identify different types of faults in rolling element bearings. the algorithm is then applied to a batch of data sets in disk, which helped show that the amplitudes of bandpass filtered envelope spectrum at bpfi and bpfo are two important condition indicators for bearing diagnostics.

references

[1] randall, robert b., and jerome antoni. "rolling element bearing diagnostics—a tutorial." mechanical systems and signal processing. vol. 25, number 2, 2011, pp. 485–520.

[2] antoni, jérôme. "fast computation of the kurtogram for the detection of transient faults." mechanical systems and signal processing. vol. 21, number 1, 2007, pp. 108–124.

[3] antoni, jérôme. "the spectral kurtosis: a useful tool for characterising non-stationary signals." mechanical systems and signal processing. vol. 20, number 2, 2006, pp. 282–307.

[4] bechhoefer, eric. "condition based maintenance fault database for testing diagnostics and prognostic algorithms." 2013.

helper functions

function bearingfeatureextraction(ensemble)
% extract condition indicators from bearing data
data = read(ensemble);
x = data.gs{1};
fs = data.sr;
% critical frequencies
bpfo = data.bpfo;
bpfi = data.bpfi;
level = 9;
[~, ~, ~, fc, ~, bw] = kurtogram(x, fs, level);
% bandpass filtered envelope spectrum
[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