Воздушное обнаружение отказа компрессора Используя рассеивание вейвлета

В этом примере показано, как классифицировать отказы на акустические записи воздушных компрессоров с помощью сети рассеивания вейвлета и машины опорных векторов (SVM). Пример обеспечивает возможность использовать графический процессор, чтобы ускориться, расчет рассеивания вейвлета преобразовывают. Если вы хотите использовать графический процессор, у вас должны быть Parallel Computing Toolbox™ и поддерживаемый графический процессор. Смотрите Поддержку графического процессора Релизом (Parallel Computing Toolbox) для деталей.

Набор данных

Набор данных состоит из акустических записей, собранных на одноступенчатом воздушном компрессоре типа оплаты [1]. Данные производятся на уровне 16 кГц. Технические требования воздушного компрессора следующие:

  • Область значений Давления воздуха: 0-500 lb/m2, 0-35 Kg/cm2

  • Асинхронный двигатель: 5 л. с., 415 В, 5:00, 50 Гц, 1440 об/мин

  • Датчик давления: введите PR-15, область значений 100-213 фунтов на квадратный дюйм

Каждая запись представляет одно из 8 состояний, которое включает здоровое состояние и 7 дефектных состояний. 7 дефектных состояний:

  1. Отказ утечки вставила клапан (LIV)

  2. Отказ клапана выхода утечки (LOV)

  3. Отказ обратного клапана (NRV)

  4. Поршневой кольцевой отказ

  5. Отказ маховика

  6. Отказ пояса наездника

  7. Подшипник отказа

Загрузите набор данных и разархивируйте файл данных в папке, где у вас есть разрешение записи. Этот пример принимает, что вы загружаете данные во временной директории, определяемой как tempdir в MATLAB®. Если вы принимаете решение использовать другую папку, замените той папкой tempdir в следующем. Записи хранятся как .wav файлы в папках, названных по имени их соответствующего состояния.

url = 'https://www.mathworks.com/supportfiles/audio/AirCompressorDataset/AirCompressorDataset.zip';
downloadFolder = fullfile(tempdir,'AirCompressorDataSet');
if ~exist(fullfile(tempdir,'AirCompressorDataSet'),'dir')
    loc = websave(downloadFolder,url);
    unzip(loc,fullfile(tempdir,'AirCompressorDataSet'))
end

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

datasetLocation = fullfile(tempdir,'AirCompressorDataSet','AirCompressorDataset');
ads = audioDatastore(datasetLocation,'IncludeSubfolders',true,...
    'LabelSource','foldernames');

Исследуйте количество примеров в каждом классе. Существует 225 записей в каждом классе для в общей сложности 1 800 записей.

countcats(ads.Labels)
ans = 8×1

   225
   225
   225
   225
   225
   225
   225
   225

Разделите данные в наборы обучающих данных и наборы тестов. Используйте 80% данных для обучения и протяните остающиеся 20% для тестирования. Переставьте данные, однажды разделяющие.

rng default
ads = shuffle(ads);
[adsTrain,adsTest] = splitEachLabel(ads,0.8,0.2);

Проверьте, что количество примеров в каждом классе является ожидаемым номером.

uniqueLabels = unique(adsTrain.Labels);
tblTrain = countEachLabel(adsTrain);
tblTest = countEachLabel(adsTest);
H = bar(uniqueLabels,[tblTrain.Count, tblTest.Count],'stacked');
legend(H,["Training Set","Test Set"],'Location','NorthEastOutside')

Выберите некоторые случайные примеры из набора обучающих данных для графического вывода. Каждая запись имеет 50 000 выборок, произведенных на уровне 16 кГц.

idx = randperm(numel(adsTrain.Files),8);
Fs = 16e3;
for n = 1:numel(idx)
    x = audioread(adsTrain.Files{idx(n)});
    t = (0:size(x,1)-1)/Fs;
    subplot(4,2,n);
    plot(t,x);
    if n == 7 || n == 8
        xlabel('Seconds');
    end
    title(string(adsTrain.Labels(idx(n))));
end

Сеть рассеивания вейвлета

Создайте сеть рассеивания вейвлета на основе характеристик данных. Установите шкалу инвариантности составлять 0,5 секунды.

N = 5e4;
Fs = 16e3;
IS = 0.5;
sn = waveletScattering('SignalLength',N,'SamplingFrequency',Fs,...
    'InvarianceScale',0.5);

С этими сетевыми настройками существует 330 рассеивающихся путей и 25 раз окна на пример. Вы видите это со следующим кодом.

[~,numpaths] = paths(sn);
Ncfs = numCoefficients(sn);
sum(numpaths)
ans = 330
Ncfs
Ncfs = 25

Обратите внимание, что это уже представляет 6-кратное сокращение размера данных для каждой записи. Мы уменьшали размер данных с 50 000 выборок до 8 250. В этом случае мы получаем дальнейшее сокращение подвыборкой окна времени на коэффициент 6, уменьшая размер каждого примера на коэффициент 25.

Получите функции рассеивания вейвлета наборов обучающих данных и наборов тестов. Если у вас есть подходящий графический процессор и Parallel Computing Toolbox, можно установить useGPU к true ускорять рассеивающееся преобразование. Функциональный helperBatchScatFeatures получает рассеивающееся преобразование и подпроизводит окна времени на коэффициент 6.

batchsize = 64;
useGPU = false; 
scTrain = [];
while hasdata(adsTrain)
    sc = helperBatchScatFeatures(adsTrain,sn,N,batchsize,useGPU);
    scTrain = cat(3,scTrain,sc);
end

Повторите процесс для протянутого набора тестов.

scTest = [];
while hasdata(adsTest)
   sc = helperBatchScatFeatures(adsTest,sn,N,batchsize,useGPU);
   scTest = cat(3,scTest,sc); 
end
scTest = gather(scTest);

Определите количество рассеивающихся путей и количество окон времени. Существует 330 рассеивающихся путей и пять раз окна.

[Npaths,numTimeWindows,~] = size(scTrain);

Измените рассеивающиеся преобразования для наборов обучающих данных и наборов тестов так, чтобы каждая строка была окном времени через все 330 рассеивающихся путей. Поскольку существует пять раз, когда окна в этом подпроизведенном рассеивании вейвлета преобразовывают, размер TrainFeatures 5*1440-by-330. Точно так же размер TestFeatures 5*360-by-330.

TrainFeatures = permute(scTrain,[2 3 1]);
TrainFeatures = reshape(TrainFeatures,[],Npaths,1);
TestFeatures = permute(scTest,[2 3 1]);
TestFeatures = reshape(TestFeatures,[],Npaths,1);

Для того, чтобы подбирать нашу модель SVM, реплицируйте метки так, чтобы была метка для каждой строки TrainFeatures.

trainLabels = adsTrain.Labels;
numTrainSignals = numel(trainLabels);
trainLabels = repmat(trainLabels,1,numTimeWindows);
trainLabels = reshape(trainLabels',numTrainSignals*numTimeWindows,1);

Классификация SVM

Используйте классификатор машины опорных векторов (SVM) мультикласса с ядром кубического полинома. Соответствуйте SVM к обучающим данным.

template = templateSVM(...
    'KernelFunction', 'polynomial', ...
    'PolynomialOrder', 3, ...
    'KernelScale', 'auto', ...
    'BoxConstraint', 1, ...
    'Standardize', true);

classificationSVM = fitcecoc(...
    TrainFeatures, ...
    trainLabels, ...
    'Learners', template, ...
    'Coding', 'onevsall','ClassNames',uniqueLabels);

Определите точность перекрестной проверки модели SVM с помощью 5-кратной перекрестной проверки.

partitionedModel = crossval(classificationSVM, 'KFold', 5);
validationAccuracy = (1 - kfoldLoss(partitionedModel))*100
validationAccuracy = single
    99.5278

Точность перекрестной проверки составляет более чем 99%, когда каждый раз окно отдельно классифицируется. Эта эффективность превосходна, но фактическая точность перекрестной проверки на самом деле выше. Поскольку существует пять окон на пример, мы должны использовать все пять раз окна, чтобы присвоить метку класса. Существуют многочисленные способы сделать это, но здесь мы используем голосование простого большинства по этим пяти окнам. Если нет никакого большинства, мы присваиваем класс "NoUniqueMode" и полагаем что ошибка классификации.

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

validationPredictions = kfoldPredict(partitionedModel);
TrainVotes = helperMajorityVote(validationPredictions,adsTrain.Labels,uniqueLabels);
crossvalAccuracy = sum(TrainVotes == adsTrain.Labels)/numel(adsTrain.Labels)*100;
crossvalAccuracy
crossvalAccuracy = 100

Точность перекрестной проверки составляет близкие 100%. Примените модель к протянутому набору тестов. Используйте голосование простого большинства, чтобы присвоить метку класса.

predLabels = predict(classificationSVM,TestFeatures);
[TestVotes,TestCounts] = helperMajorityVote(predLabels,adsTest.Labels,uniqueLabels);
testAccuracy = sum(TestVotes == adsTest.Labels)/numel(adsTest.Labels)*100;
testAccuracy
testAccuracy = 100

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

figure
confusionchart(adsTest.Labels,TestVotes)

Сводные данные

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

Ссылки

[1] Verma, Нищел К., Рахул Кумар Севэкула, Сонэл Диксит и Аль Сэлур. “Интеллектуальный основанный на условии Контроль Используя Акустические Сигналы для Воздушных Компрессоров”. Транзакции IEEE на Надежности 65, № 1 (март 2016): 291–309. https://doi.org/10.1109/TR.2015.2459684.

helperBatchScatFeatures - Эта функция возвращает время вейвлета, рассеивая матрицу функции для данного входного сигнала. Рассеивающиеся функции подпроизводятся на коэффициент 6. Если useGPU установлен в true, рассеивающееся преобразование вычисляется на графическом процессоре.

function sc = helperBatchScatFeatures(ds,sn,N,batchsize,useGPU)
% This function is only intended to support examples in the Wavelet
% Toolbox. It may be changed or removed in a future release.

% Read batch of data from audio datastore
batch = helperReadBatch(ds,N,batchsize);
if useGPU
    batch = gpuArray(batch);
end

% Obtain scattering features
S = sn.featureMatrix(batch,'transform','log');
gather(batch);
S = gather(S);
    
% Subsample the features
sc = S(:,1:6:end,:);
end

helperReadBatch - Эта функция читает пакеты заданного размера от datastore и возвращает выходной параметр в одинарной точности. Каждый столбец выхода является отдельным сигналом от datastore. Выход может иметь меньше столбцов, чем batchsize, если datastore не имеет достаточного количества записей.

function batchout = helperReadBatch(ds,N,batchsize)
% This function is only in support of Wavelet Toolbox examples. It may
% change or be removed in a future release.
%
% batchout = readReadBatch(ds,N,batchsize) where ds is the Datastore and
%   ds is the Datastore
%   batchsize is the batchsize

kk = 1;

while(hasdata(ds)) && kk <= batchsize
    tmpRead = read(ds);
    batchout(:,kk) = cast(tmpRead(1:N),'single'); %#ok<AGROW>
    kk = kk+1;
end

end

helperMajorityVote - Эта функция получает решение большинством голосов для меток класса. Если нет никакого большинства, "NoUniqueMode" возвращен и обработан как ошибка.

function [ClassVotes,ClassCounts] = helperMajorityVote(predLabels,origLabels,classes)
% This function is in support of wavelet scattering examples only. It may
% change or be removed in a future release.

% Make categorical arrays if the labels are not already categorical
predLabels = categorical(predLabels);
origLabels = categorical(origLabels);
% Expects both predLabels and origLabels to be categorical vectors
Npred = numel(predLabels);
Norig = numel(origLabels);
Nwin = Npred/Norig;
predLabels = reshape(predLabels,Nwin,Norig);
ClassCounts = countcats(predLabels);
[mxcount,idx] = max(ClassCounts);
ClassVotes = classes(idx);
tmpsum = sum(ClassCounts == mxcount);
ClassVotes(tmpsum > 1) = categorical({'NoUniqueMode'});
end

Смотрите также

Связанные примеры

Больше о