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

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

Каждый набор данных содержит сигнал ускорения «gs», частоту дискретизации «sr», скорость «вал скорости», «нагрузку веса нагрузки» и четыре критические частоты, представляющие различные местоположения отказа: внешнее кольцо с частотами баллпаса (BPFO), внутреннее кольцо с частотами баллпаса (BPFI), основная частота train (FTF F) и частота шарового спина (bin) Вот формулы для этих критических частот [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 (или увеличение отношения сигнал/шум) является ключевым шагом предварительной обработки перед огибающей спектра анализом. В следующем разделе введут куртограмму и спектральный куртоз для извлечения сигнала с наивысшим куртозом и проведут анализ огибающей спектра по отфильтрованному сигналу.

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

Куртограмма и спектральный куртоз вычисляют куртоз локально в полосах. Они являются мощными инструментами для определения местоположения полосы частот, которая имеет самый высокий куртоз (или самое высокое отношение сигнал/шум) [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 и его гармоники.

Пакетная обработка

Теперь применим алгоритм к пакету обучающих данных с помощью файла ensemble 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 файлами mat (2 файла normal, 4 отказа внутреннего кольца, 7 отказ внешнего кольца) и тестовый набор данных с 6 файлами mat (1 файл normal, 2 отказ внутреннего кольца, 3 отказа внешнего кольца).

Путем назначения указателей на функцию ReadFcn и WriteToMemberFcn, file ensemble datastore сможет перемещаться в файлы, чтобы получить данные в желаемом формате. Для примера данные MFPT имеют структуру bearing который сохраняет сигнал вибрации gs, частота дискретизации srи так далее. Вместо возврата самой несущей конструкции readMFPTBearing функция записывается так, что file ensemble 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 являются двумя индикаторами состояния для диагностики отказа подшипника. Поэтому следующим шагом является извлечение двух индикаторов состояния из всех обучающих данных. Чтобы сделать алгоритм более устойчивым, установите узкую полосу (bandwidth = $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 App в Statistics & Машинное Обучение 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