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

Загрузка данных о вибрации, сохраненных в pdmBearingConditionMonitoringData.mat (это большой набор данных ~ 88 МБ, который загружается с сайта файлов поддержки MathWorks). Данные хранятся в массиве ячеек, которая была сгенерирована с помощью симулятора мяча сигнала с одной точкой дефектом на внешнем кольце подшипника. Он содержит несколько сегментов сигналов вибрации для подшипников, моделируемых в различных условиях здоровья (глубина дефекта увеличивается с 3 до 3 мм). Каждый сегмент хранит сигналы, собранные в течение 1 секунды, со частотой дискретизации 20 кГц. В pdmBearingConditionMonitoring Data.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,\omega )$$. Пиковая частота в каждом образце времени определяется как:

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

Вычислите среднюю пиковую частоту для исправных сигналов мяча.

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

Для прогностики и мониторинга условия необходимо задать порог, чтобы решить, когда остановить машину. В этом примере используйте статистические данные исправного и неисправного подшипника, сгенерированного при симуляции, чтобы определить порог. PDM- BearingConditionMonitoringStatistics.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), рассчитанной как

где $$x_{true}$$ является истинным значением, $$x_{pred}$$ - предсказанное значение.

Когда данные оценки являются комбинацией постоянного уровня и шума, $$x_{pred}\approx mean(x_{true})$$, придавая 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%, и тренд правильно фиксируется.

Заключения

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