Описание данных

Данные о вибрации нагрузки, хранящиеся в pdmBearingConditionMonitoringData.mat (это большой набор данных ~ 88 МБ, который загружается с сайта файлов поддержки MathWorks). Данные хранятся в матрице ячеек, которая была сформирована с помощью имитатора сигнала шарикоподшипника с одиночным точечным дефектом на внешнем кольце подшипника. Содержит множество сегментов вибрационных сигналов для подшипников, моделируемых в различных условиях здоровья (глубина дефекта увеличивается с 3 мкм до более 3 мм). Каждый сегмент хранит сигналы, собранные в течение 1 секунды при частоте дискретизации 20 кГц. В pdmBearingConditionMonitoringData.mat, defectDepthVec сохраняет, как изменяется глубина дефекта в зависимости от времени, в то время как expTime сохраняет соответствующее время в минутах.

url = 'https://www.mathworks.com/supportfiles/predmaint/condition-monitoring-and-prognostics-using-vibration-signals/pdmBearingConditionMonitoringData.mat';
websave('pdmBearingConditionMonitoringData.mat',url);
load pdmBearingConditionMonitoringData.mat

% Define the number of data points to be processed.
numSamples = length(data);

% Define sampling frequency.
fs = 20E3;          % unit: Hz

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

plot(expTime,defectDepthVec);
xlabel('Time (min)');
ylabel('Defect depth (m)');

Постройте график здоровых и неисправных данных.

time = linspace(0,1,fs)';

% healthy bearing signal
subplot(2,1,1);
plot(time,data{1});
xlabel('Time(s)');
ylabel('Acceleration(m/s^2)');
legend('Healthy bearing signal');

% faulty bearing signal
subplot(2,1,2);
plot(time,data{end});
xlabel('Time(s)');
ylabel('Acceleration(m/s^2)');
legend('Faulty bearing signal');

Извлечение элементов

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

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

Сначала вычислите спектрограмму данных здорового подшипника. Используйте размер окна 500 точек данных и коэффициент перекрытия 90% (эквивалентный 450 точкам данных). Установите количество точек для БПФ равным 512. fs представляет частоту выборки, определенную ранее.

[~,fvec,tvec,P0] = spectrogram(data{1},500,450,512,fs);

P0 - спектрограмма, fvec - частотный вектор и tvec - вектор времени.

Постройте график спектрограммы сигнала здорового подшипника.

clf;
imagesc(tvec,fvec,P0)
xlabel('Time(s)');
ylabel('Frequency(Hz)');
title('Healthy Bearing Signal Spectrogram');
axis xy

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

[~,fvec,tvec,Pfinal] = spectrogram(data{end},500,450,512,fs);
imagesc(tvec,fvec,Pfinal)
xlabel('Time(s)');
ylabel('Frequency(Hz)');
title('Faulty Bearing Signal Spectrogram');
axis xy

Поскольку спектрограммы для данных от здоровых и неисправных подшипников различны, репрезентативные признаки могут быть извлечены из спектрограмм и использованы для мониторинга состояния и прогностики. В этом примере извлекают средние пиковые частоты из спектрограмм в качестве показателей здоровья. Обозначим спектрограмму как $$P(t,$$λ). Пиковая частота в каждом моменте времени определяется как:

Средняя пиковая частота представляет собой среднее из пиковых частот, определенных выше.

Рассчитайте среднюю пиковую частоту для здоровых сигналов шарикоподшипника.

[~,I0] = max(P0);               % Find out where the peak frequencies are located.
meanPeakFreq0 = mean(fvec(I0))  % Calculate mean peak frequency.

meanPeakFreq0 = 666.4602

Здоровые сигналы вибрации подшипника имеют среднюю пиковую частоту около 650 Гц. Теперь рассчитайте среднюю пиковую частоту для неисправных сигналов подшипника. Средняя пиковая частота сдвигается выше 2500 Гц.

[~,Ifinal] = max(Pfinal);
meanPeakFreqFinal = mean(fvec(Ifinal))

meanPeakFreqFinal = 2.8068e+03

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

[~,fvec,tvec,Pmiddle] = spectrogram(data{end/2},500,450,512,fs);
imagesc(tvec,fvec,Pmiddle)
xlabel('Time(s)');
ylabel('Frequency(Hz)');
title('Bearing Signal Spectrogram');
axis xy

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

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

dataMiddleFilt = medfilt1(data{end/2},3);

Постройте график-спектрограмму после медианной фильтрации. Высокочастотные компоненты подавляются.

[~,fvec,tvec,Pmiddle] = spectrogram(dataMiddleFilt,500,450,512,fs);
imagesc(tvec,fvec,Pmiddle)
xlabel('Time(s)');
ylabel('Frequency(Hz)');
title('Filtered Bearing Signal Spectrogram');
axis xy

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

% Define a progress bar.
h = waitbar(0,'Start to extract features');
% Initialize a vector to store the extracted mean peak frequencies.
meanPeakFreq = zeros(numSamples,1);
for k = 1:numSamples
    % Get most up-to-date data.
    curData = data{k};
    % Apply median filter.
    curDataFilt = medfilt1(curData,3);
    % Calculate spectrogram.
    [~,fvec,tvec,P_k] = spectrogram(curDataFilt,500,450,512,fs);
    % Calculate peak frequency at each time instance.
    [~,I] = max(P_k);
    meanPeakFreq(k) = mean(fvec(I));
    % Show progress bar indicating how many samples have been processed.
    waitbar(k/numSamples,h,'Extracting features');
end
close(h);

Постройте график зависимости выделенных средних пиковых частот от времени.

plot(expTime,meanPeakFreq);
xlabel('Time(min)');
ylabel('Mean peak frequency (Hz)');
grid on;

Мониторинг и прогнозирование состояния

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

Прогнозирование помогает нам лучше подготовиться к потенциальной неисправности или даже остановить машину до отказа. Рассмотрим среднюю пиковую частоту как временной ряд. Мы можем оценить модель временных рядов для средней пиковой частоты и использовать модель для прогнозирования будущих значений. Используйте первые 200 средних пиковых значений частоты для создания начальной модели временного ряда, затем, как только будут доступны 10 новых значений, используйте последние 100 значений для обновления модели временного ряда. Этот пакетный режим обновления модели временных рядов фиксирует мгновенные тренды. Обновленная модель временных рядов используется для расчета прогноза на 10 шагов вперед.

tStart = 200;               % Start Time
timeSeg = 100;              % Length of data for building dynamic model
forecastLen = 10;           % Define forecast time horizon
batchSize = 10;             % Define batch size for updating the dynamic model

Для прогностики и мониторинга состояния необходимо установить порог, чтобы решить, когда остановить машину. В этом примере для определения порогового значения следует использовать статистические данные о работоспособном и неисправном подшипнике, полученные при моделировании. pdmBearingConditionMonitoringStatistics.mat сохраняет распределения вероятности средних пиковых частот для здоровых и неисправных подшипников. Распределения вероятности вычисляются путем возмущения глубины дефекта здоровых и неисправных подшипников.

url = 'https://www.mathworks.com/supportfiles/predmaint/condition-monitoring-and-prognostics-using-vibration-signals/pdmBearingConditionMonitoringStatistics.mat';
websave('pdmBearingConditionMonitoringStatistics.mat',url);
load pdmBearingConditionMonitoringStatistics.mat

Постройте график распределения вероятности средней пиковой частоты для здоровых и неисправных подшипников.

plot(pFreq,pNormal,'g--',pFreq,pFaulty,'r');
xlabel('Mean peak frequency(Hz)');
ylabel('Probability Distribution');
legend('Normal Bearing','Faulty Bearing');

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

threshold = 2000;

Рассчитайте время выборки и преобразуйте ее в секунды.

samplingTime = 60*(expTime(2)-expTime(1));                  % unit: seconds
tsFeature = iddata(meanPeakFreq(1:tStart),[],samplingTime);

Постройте график исходных данных 200 средних пиковых частот.

plot(tsFeature.y)

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

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

past_sys = ssest(tsFeature,2,'Ts',samplingTime,'Form','canonical')

past_sys =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + K e(t)
       y(t) = C x(t) + e(t)
 
  A = 
            x1       x2
   x1        0        1
   x2   0.9605  0.03942
 
  C = 
       x1  x2
   y1   1   0
 
  K = 
              y1
   x1  -0.003899
   x2   0.003656
 
Sample time: 600 seconds
  
Parameterization:
   CANONICAL form with indices: 2.
   Disturbance component: estimate
   Number of free coefficients: 4
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                               
Estimated using SSEST on time domain data "tsFeature".
Fit to estimation data: 0.2763% (prediction focus)    
FPE: 640, MSE: 602.7                                  

Первоначальная расчетная динамическая модель имеет низкое качество посадки. Достоверность метрики соответствия - нормированная среднеквадратическая ошибка (NRMSE), рассчитанная как

где $${}_{\mathrm{xtrue}}$$ - истинное значение, $${}_{\mathrm{xpred}}$$ - прогнозируемое значение.

Когда оценочные данные представляют собой комбинацию постоянного уровня и шума, $${}_{}\mathrm{xpred\approx mean}{(}_{}xtrue$$), давая NRMSE, близкий к 0. Чтобы проверить модель, постройте график автокорреляции остатков .

resid(tsFeature,past_sys)

Как видно, остатки не коррелируются, и сгенерированная модель является действительной.

Использовать идентифицированную модель past_sys прогнозировать средние пиковые значения частоты и вычислять стандартное отклонение прогнозируемых значений.

[yF,~,~,yFSD] = forecast(past_sys,tsFeature,forecastLen);

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

tHistory = expTime(1:tStart);
forecastTimeIdx = (tStart+1):(tStart+forecastLen);
tForecast = expTime(forecastTimeIdx);

% Plot historical data, forecast value and 95% confidence interval.
plot(tHistory,meanPeakFreq(1:tStart),'b',...
     tForecast,yF.OutputData,'kx',...
     [tHistory; tForecast], threshold*ones(1,length(tHistory)+forecastLen), 'r--',...
     tForecast,yF.OutputData+1.96*yFSD,'g--',...
     tForecast,yF.OutputData-1.96*yFSD,'g--');

ylim([400, 1.1*threshold]);
ylabel('Mean peak frequency (Hz)');
xlabel('Time (min)');
legend({'Past Data', 'Forecast', 'Failure Threshold', '95% C.I'},...
    'Location','northoutside','Orientation','horizontal');
grid on;

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

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

for tCur = tStart:batchSize:numSamples
    %  latest features into iddata object.
    tsFeature = iddata(meanPeakFreq((tCur-timeSeg+1):tCur),[],samplingTime);
    
    % Update system parameters when new data comes in. Use previous model
    % parameters as initial guesses.
    sys = ssest(tsFeature,past_sys);
    past_sys = sys;
    
    % Forecast the output of the updated state-space model. Also compute
    % the standard deviation of the forecasted output.
    [yF,~,~,yFSD]  = forecast(sys, tsFeature, forecastLen);
    
    % Find the time corresponding to historical data and forecasted values.
    tHistory = expTime(1:tCur);
    forecastTimeIdx = (tCur+1):(tCur+forecastLen);
    tForecast = expTime(forecastTimeIdx);
    
    % Plot historical data, forecasted mean peak frequency value and 95%
    % confidence interval.
    plot(tHistory,meanPeakFreq(1:tCur),'b',...
              tForecast,yF.OutputData,'kx',...
              [tHistory; tForecast], threshold*ones(1,length(tHistory)+forecastLen), 'r--',...
              tForecast,yF.OutputData+1.96*yFSD,'g--',...
              tForecast,yF.OutputData-1.96*yFSD,'g--');

    ylim([400, 1.1*threshold]);
    ylabel('Mean peak frequency (Hz)');
    xlabel('Time (min)');
    legend({'Past Data', 'Forecast', 'Failure Threshold', '95% C.I'},...
           'Location','northoutside','Orientation','horizontal');
    grid on;
    
    % Display an alarm when actual monitored variables or forecasted values exceed
    % failure threshold.
    if(any(meanPeakFreq(tCur-batchSize+1:tCur)>threshold))
        disp('Monitored variable exceeds failure threshold');
        break;
    elseif(any(yF.y>threshold))
        % Estimate the time when the system will reach failure threshold.
        tAlarm = tForecast(find(yF.y>threshold,1));
        disp(['Estimated to hit failure threshold in ' num2str(tAlarm-tHistory(end)) ' minutes from now.']);
        break;
    end
end

Estimated to hit failure threshold in 80 minutes from now.

Изучите последнюю модель временных рядов.

sys

sys =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + K e(t)
       y(t) = C x(t) + e(t)
 
  A = 
           x1      x2
   x1       0       1
   x2  0.2624   0.746
 
  C = 
       x1  x2
   y1   1   0
 
  K = 
           y1
   x1  0.3902
   x2  0.3002
 
Sample time: 600 seconds
  
Parameterization:
   CANONICAL form with indices: 2.
   Disturbance component: estimate
   Number of free coefficients: 4
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                               
Estimated using SSEST on time domain data "tsFeature".
Fit to estimation data: 92.53% (prediction focus)     
FPE: 499.3, MSE: 442.7                                

Благость посадки увеличивается до выше 90%, и тренд фиксируется правильно.

Заключения

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