Документация

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

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

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

Как показано в фигуре, $\mathit{d}$ диаметр шара, $\mathit{D}$ диаметр подачи. Переменная ${\mathit{f}}_{\mathit{r}}$ скорость вала, $\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')

Как ожидалось спектр конверта нормального сигнала переноса не показывает значительного 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 является допустимой функцией, чтобы классифицировать отказы переноса. Чтобы упростить пример, очень простой классификатор выведен: если $\mathrm{журнал}\left(\frac{\mathrm{BPFIAmp}\mathrm{litude}}{\mathrm{BPFOAmp}\mathrm{litude}}\right)\le -1.5$, перенос имеет внешний отказ гонки; если $-1.5<\mathrm{журнал}\left(\frac{\mathrm{BPFIAmp}\mathrm{litude}}{\mathrm{BPFOAmp}\mathrm{litude}}\right)\le 0.5$, перенос нормален; и если $\mathrm{журнал}\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://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