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

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

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

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

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

meanPeakFreq0 = 666.4602

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

[~,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');

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

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%, и тренд правильно получен.

Заключения

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