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

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

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

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

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

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

Когда прокрутка элементов поразила локальные отказы на внешних или внутренних гонках, или когда отказы на прокручивающемся элементе поразят внешние или внутренние гонки, удар будет модулировать соответствующие критические частоты, e.g. 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')

Как ожидалось огибающая спектра нормального сигнала подшипника не показывает значительного peaks в 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')

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

Kurtogram и Spectral Kurtosis для выбора полосы

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

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

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

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

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

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

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

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

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

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-файл и сохранить его в той же директории как live скрипт. Разархивируйте файл с помощью этой команды:

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

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

Путем присвоения указателей на функцию ReadFcn и WriteToMemberFcn, datastore ансамбля файла сможет перейти в файлы, чтобы получить данные в нужном формате. Например, данные MFPT имеют структуру bearing это хранит сигнал вибрации gs, частота дискретизации sr, и так далее. Вместо того, чтобы возвратить саму структуру подшипника readMFPTBearing функция записана так, чтобы datastore ансамбля файла возвратил сигнал вибрации 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 являются двумя индикаторами состояния для подшипника диагностики отказа. Поэтому следующий шаг должен извлечь эти два индикатора состояния из всех обучающих данных. Чтобы сделать алгоритм более устойчивым, установите узкую полосу (полоса пропускания = $10\Delta \mathit{f}$, где $\Delta \mathit{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

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

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 является допустимой функцией, чтобы классифицировать отказы подшипника. Чтобы упростить пример, очень простой классификатор выведен: если $\log \left(\frac{\mathrm{BPFIAmp}\mathrm{litude}}{\mathrm{BPFOAmp}\mathrm{litude}}\right)\le -1.5$, подшипник имеет внешний отказ гонки; если $-1.5<\log \left(\frac{\mathrm{BPFIAmp}\mathrm{litude}}{\mathrm{BPFOAmp}\mathrm{litude}}\right)\le 0.5$, подшипник нормален; и если $\log \left(\frac{\mathrm{BPFIAmp}\mathrm{litude}}{\mathrm{BPFOAmp}\mathrm{litude}}\right)>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 в Statistics & Machine Learning Toolbox. Для получения дополнительной информации этого рабочего процесса, обратитесь к примеру "Используя Simulink, чтобы сгенерировать данные об отказе".

Сводные данные

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

Ссылки

[1] Рэндалл, Роберт Б. и Джером Антони. "Диагностика подшипника качения — пример". Механические Системы и Обработка сигналов. Издание 25, Номер 2, 2011, стр 485–520.

[2] Антони, Жером. "Быстрый расчет kurtogram для обнаружения переходных отказов". Механические Системы и Обработка сигналов. Издание 21, Номер 1, 2007, стр 108–124.

[3] Антони, Жером. "Спектральный эксцесс: полезный инструмент для характеристики неустановившихся сигналов". Механические Системы и Обработка сигналов. Издание 20, Номер 2, 2006, стр 282–307.

[4] Bechhoefer, Эрик. "Основанная на условии База данных Отказа Обслуживания для Тестирования Диагностики и Предвещающих Алгоритмов". 2013. https://www.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