Данные вызова технологии предотвращения отказов оборудования (MFPT)

Данные MFPT Challenge [4] содержат 23 набора данных, собранных с машин при различных условиях отказа. Первые 20 наборов данных собраны с установки для испытания подшипников, 3 в хороших условиях, 3 с разломами наружного кольца при постоянной нагрузке, 7 с разломами наружного кольца при различных нагрузках и 7 с разломами внутреннего кольца при различных нагрузках. Остальные 3 набора данных получены от реальных машин: подшипник масляного насоса, подшипник промежуточной скорости и планетарный подшипник. Местоположения отказов неизвестны. В этом примере используются только данные, собранные с испытательной установки с известными условиями.

Каждый набор данных содержит сигнал ускорения «gs», частоту дискретизации «sr», скорость «частоты вращения вала», «вес нагрузки» и четыре критические частоты, представляющие различные места возникновения отказов: BPFO, BPFI, FTF и BSF. Вот формулы для этих критических частот [1].

Как показано на чертеже, $\mathit{d}$ - диаметр шарика, $\mathit{D}$ - диаметр шага. Переменный ${\mathit{}}_{\mathit{франк}}$ - скорость шахты, $\mathit{n}$ - число катящихся элементов, $\varphi$ - угол контакта отношения [1].

Анализ огибающей спектра для диагностики подшипников

В наборе данных MFPT скорость вала является постоянной, поэтому нет необходимости выполнять отслеживание заказа в качестве этапа предварительной обработки для устранения влияния изменений скорости вала.

Когда элементы качения ударяют по локальным разломам на внешнем или внутреннем кольце или когда разломы на элементе качения ударяют по внешнему или внутреннему кольцу, удар будет модулировать соответствующие критические частоты, например BPFO, BPFI, FTF, BSF. Следовательно, сигнал огибающей, полученный амплитудной демодуляцией, передает больше диагностической информации, которая недоступна из спектрального анализа исходного сигнала. Возьмем в качестве примера сигнал отказа внутренней расы в наборе данных MFPT.

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

Визуализация необработанных данных в частотной области.

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 и его первые несколько гармоник.

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 и его гармониках. Частотный анализ необработанного сигнала не дает полезной диагностической информации.

При рассмотрении данных временной области наблюдается, что амплитуда исходного сигнала модулируется на определенной частоте, а основная частота модуляции составляет около 1/0,009 Гц $\approx$ 111 Гц. Известно, что частота наезда элемента качения на локальный отказ во внутреннем кольце, то есть BPFI, составляет 118,875 Гц. Это указывает на то, что подшипник потенциально имеет неисправность внутреннего кольца.

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 и его гармоники.

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

Показано, что большая часть энергии сосредоточена на BPFI и его гармониках. Это указывает на отказ внутренней цепи подшипника, который соответствует типу отказа данных.

Применение анализа огибающей спектра к другим типам отказов

Теперь повторите один и тот же анализ огибающей спектра для нормальных данных и данных об отказах внешнего контура.

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

Как и ожидалось, спектр огибающей нормального несущего сигнала не показывает каких-либо значительных пиков в 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. Не позволяет ли анализ огибающей спектра отличить подшипник с неисправностью внешнего кольца от здоровых подшипников? Давайте сделаем шаг назад и снова посмотрим на сигналы во временной области при различных условиях.

Сначала давайте снова визуализируем сигналы во временной области и рассчитаем их куртоз. Куртоз - четвёртый стандартизированный момент случайной величины. Он характеризует импульсивность сигнала или тяжесть хвоста случайной величины.

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 (или усиление отношения сигнал/шум) является ключевым этапом предварительной обработки перед анализом огибающей спектра. В следующем разделе будут введены куртограмма и спектральный куртоз для выделения сигнала с наивысшим куртозом и для отфильтрованного сигнала будет выполнен анализ огибающей спектра.

Куртограмма и спектральный куртоз для выбора диапазона

Куртограмма и спектральный куртоз вычисляют куртоз локально в частотных диапазонах. Они представляют собой мощные инструменты для определения диапазона частот с наивысшим куртозом (или наивысшим отношением сигнал/шум) [2]. После определения полосы частот с наивысшим куртозом к исходному сигналу может быть применен полосовой фильтр для получения более импульсного сигнала для анализа спектра огибающей.

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

Куртограмма показывает, что полоса частот, центрированная на 2,67 кГц с шириной полосы 0,763 кГц, имеет наивысший куртоз 2,71.

Теперь используйте оптимальную длину окна, предложенную куртограммой, для вычисления спектрального куртоза.

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

Чтобы визуализировать полосу частот на спектрограмме, вычислите спектрограмму и поместите спектральный куртоз сбоку. Для интерпретации спектрального куртоза другим способом высокие значения спектрального куртоза указывают на высокую дисперсию мощности на соответствующей частоте, что делает спектральный куртоз полезным инструментом для нахождения нестационарных составляющих сигнала [3].

helperSpectrogramAndSpectralKurtosis(xOuter, fsOuter, level)

Посредством полосовой фильтрации сигнала с предлагаемыми центральной частотой и полосой пропускания может быть усилен куртоз и может быть извлечена модулированная амплитуда отказа внешнего кольца.

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

Можно видеть, что значение куртоза увеличивается после полосовой фильтрации. Теперь визуализируйте сигнал огибающей в частотной области.

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

Показано, что посредством полосовой фильтрации необработанного сигнала с диапазоном частот, предложенным куртограммой и спектральным куртозом, анализ спектра огибающей способен выявить сигнатуру неисправности в BPFO и ее гармоники.

Пакетный процесс

Теперь применим алгоритм к пакету обучающих данных, используя хранилище данных файлового ансамбля.

На панели инструментов доступна ограниченная часть набора данных. Скопируйте набор данных в текущую папку и включите разрешение на запись:

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/RollingElementBearingFaultDiagnosis-Data, чтобы загрузить весь репозиторий в виде zip-файла и сохранить его в том же каталоге, что и сценарий. Распакуйте файл с помощью следующей команды:

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

Результаты в этом примере генерируются из полного набора данных. Полный набор данных содержит обучающий набор данных с 14 файлами мата (2 нормальные, 4 ошибки внутренней гонки, 7 ошибок внешней гонки) и тестовый набор данных с 6 файлами мата (1 нормальная, 2 ошибки внутренней гонки, 3 ошибки внешней гонки).

Путем назначения дескрипторов функций ReadFcn и WriteToMemberFcnхранилище данных ансамбля файлов сможет перемещаться в файлы для извлечения данных в требуемом формате. Например, данные MFPT имеют структуру bearing который сохраняет вибросигнал gs, частота отбора проб srи так далее. Вместо возврата самой несущей конструкции readMFPTBearing функция записывается таким образом, что хранилище данных ансамбля файлов возвращает сигнал вибрации gs внутри bearing структура данных.

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

Из последнего раздела анализа следует, что амплитуды спектра огибающей с полосовой фильтрацией в BPFO и BPFI являются двумя индикаторами состояния для диагностики неисправности подшипника. Поэтому следующим шагом является извлечение двух индикаторов условий из всех данных обучения. Чтобы сделать алгоритм более надежным, установите узкую полосу (полоса пропускания = $\mathrm{}\mathit{10\Delta f}$, где $\mathit{\Delta 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);

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"          

Визуализация созданной таблицы элементов.

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

Относительные значения амплитуды 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')

Гистограмма показывает четкое разделение между тремя различными условиями подшипника. Логарифмическое отношение между амплитудами BPFI и BPFO является допустимой функцией для классификации неисправностей подшипников. Чтобы упростить пример, очень простой классификатор получен: если $\mathrm{регистрация}\frac{\mathrm{}\mathrm{}}{\mathrm{}\mathrm{}}(BPFIAmplitudeBPFOAmplitude)\mathrm{\le }\mathrm{}$-1.5, у отношения есть внешняя ошибка гонки; $если\mathrm{}\mathrm{}\mathrm{-1.5\; <}\frac{\mathrm{регистрация}\mathrm{}}{\mathrm{}\mathrm{}}\mathrm{}(BPFIAmplitudeBPFOAmplitude)\mathrm{}$≤0.5, отношение нормально; и $\mathrm{если}\frac{\mathrm{регистрация}\mathrm{}}{\mathrm{}\mathrm{}}\mathrm{}\mathrm{}$(BPFIAmplitudeBPFOAmplitude)> 0.5, у отношения есть внутренняя ошибка гонки.

Проверка с использованием наборов тестовых данных

Теперь применим рабочий процесс к набору тестовых данных и проверим классификатор, полученный в последнем разделе. Здесь тестовые данные содержат 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 = 
  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')

Логарифмическое отношение амплитуд BPFI и BPFO из наборов тестовых данных показывает согласованное распределение с логарифмическим отношением из наборов обучающих данных. Наивный классификатор, полученный в последнем разделе, достиг совершенной точности в наборе тестовых данных.

Следует отметить, что одной особенности обычно недостаточно для получения классификатора, который хорошо обобщается. Более сложные классификаторы могут быть получены путем разделения данных на несколько частей (для создания большего количества точек данных), извлечения нескольких диагностических функций, выбора подмножества функций по рангам важности и обучения различных классификаторов с помощью приложения Classification Learner App in Statistics & Machine Learning Toolbox. Дополнительные сведения об этом рабочем процессе см. в примере «Использование Simulink для генерации данных об отказах».

Резюме

В этом примере показано, как использовать куртограмму, спектральный куртоз и спектр огибающей для выявления различных типов неисправностей в подшипниках элементов качения. Затем алгоритм применяется к пакету наборов данных на диске, что помогает показать, что амплитуды полосового отфильтрованного огибающего спектра в BPFI и BPFO являются двумя важными индикаторами состояния для диагностики подшипников.

Ссылки

[1] Рэндалл, Роберт Б. и Джером Антони. «Диагностика подшипников качения - учебное пособие». Механические системы и обработка сигналов. Том 25, номер 2, 2011, стр. 485-520.

[2] Антони, Жером. «Быстрое вычисление куртограммы для обнаружения переходных неисправностей». Механические системы и обработка сигналов. Том 21, номер 1, 2007, стр. 108-124.

[3] Антони, Жером. «Спектральный куртоз: полезный инструмент для характеристики нестационарных сигналов». Механические системы и обработка сигналов. Том 20, номер 2, 2006, стр. 282-307.

[4] Беххёфер, Эрик. «База данных отказов технического обслуживания на основе состояния для тестирования диагностических и прогностических алгоритмов». 2013. https://mfpt.org/fault-data-sets/

Вспомогательные функции

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