使用自编码器小波散射检测异常

这个例子使用了:

打开实时脚本

这个例子展示了如何在LSTM和卷积自编码器中使用小波散射变换来开发预测性维护的警报系统。该实例比较了小波散射变换+自编码器和原始数据+自编码器两种方法。

数据

数据集收集自一台由20齿小齿轮驱动的2MW风力涡轮机高速轴[1］．从2013-03-07到2013-04-25连续50天，每天以97,656 Hz的采样率采集6秒的振动信号。在2013-03-17上有两个测量值，在本例中视为两天。每个时间序列由585,936个样本组成。在50天的记录期内，出现了一个内圈故障，导致轴承失效。

数据下载

获取数据https://github.com/mathworks/WindTurbineHighSpeedBearingPrognosis-Data．您可以将整个存储库下载为zip文件，并将其保存在具有写权限的文件夹中。本例中的命令假设您已经下载了MATLAB指定为的文件夹中的数据tempdir．如果选择使用不同的文件夹，请更改的值parentDir在下面。下载zip文件后，使用此命令解压缩数据。

parentDir = tempdir;如果存在(fullfile (parentDir“WindTurbineHighSpeedBearingPrognosis-Data-main.zip”)，“文件”)解压缩(fullfile (parentDir,“WindTurbineHighSpeedBearingPrognosis-Data-main.zip”), parentDir)其他的错误("文件未找到。"+.．.从GitHub手动下载。zip文件的存储库。+.．."确认。zip文件在:\n%s"parentDir)结束

如果您喜欢使用Git命令克隆存储库，则该文件夹不包括“小镇命名和数据被放置在WindTurbineHighSpeedBearingPrognosis-Data文件夹中。

信号数据存储

解开WindTurbineHighSpeedBearingPrognosis-Data-main.zipFile创建一个包含50个.mat文件的文件夹。每个.mat文件包含两个变量:性心动过速而且振动．本例仅利用振动数据。因此，创建一个信号数据存储，只从.mat文件中读取振动数据。

filepath = fullfile(parentDir，“WindTurbineHighSpeedBearingPrognosis-Data-main”）;sds = signalDatastore(filepath,SignalVariableNames =“振动”）;

因为这个数据集中只有50条记录，所以要一次性将所有数据读入内存。

allSignals = readall(sds);

数据记录相对较长，有585,936个样本。作为数据缩减的第一次尝试，检查6个信号的随机样本的光谱特征。使用默认小波和级别的最大重叠小波包变换。这导致了信号能量的正交分解到宽度通带 $F 年代 / 2^{5}$ 其中Fs为采样率。画出每个信号的相对能量作为通带的函数。打印小于采样频率1/4的信号能量百分比。

rng默认的Fs = 97656;rdIdx = randperm(length(allSignals)，6);图tiledlayout(2、3)为nexttile [~，~，f，~，re] = modwpt(allSignals{ii});栏((f。* fs) / 1 e3,再保险。* 100)sprintf (' % 2.2 f '，sum(re(f <= 1/4))*100) record = rdIdx(ii);标题({的通频带能量:；［“记录”num2str(记录)]});包含(“赫兹”) ylabel (“能量”比例）结束

Ans = '96.73'

Ans = '95.84'

Ans = '96.00'

Ans = '95.61'

Ans = '95.77'

Ans = '95.95'

在这些样本中，大约96%的信号能量低于截止频率24.414 kHz，需要将信号长度降采样2倍。对数据集中所有信号的检查显示，所有记录的94.7%到98.2%的能量低于24.414 kHz。以原采样率的1/2重新采样所有信号，并相应地调整采样率。

allSignals = cellfun(@(x)resample(x,1,2)，allSignals,UniformOutput = false);Fs = Fs /2;

本例中的数据没有标记。要了解数据在50天内的演变情况，可以创建一个累积图。

Tstart = 0;数字为ii = 1:length(allSignals) t = (1:length(allSignals{ii}))。/fs + tstart;持有在plot(t,allSignals{ii}) tstart = t(end);结束持有从标题(《风力机振动累积监测》)包含(“时间(秒)——每天6秒，持续50天”) ylabel (“电压”）

从前面的图中可以明显看出，数据中的振动量似乎随着平均时间逐渐增加。在记录200秒后，这一点显得尤为明显，这相当于超过33天的时间。虽然这种平均趋势是明显的，但请注意，对于个别录音来说，情况并非总是如此。接近50天周期结束时的一些录音表现出与接近开始时的录音类似的行为。

数据准备与特征提取

将每段6秒的录音分成11段，每段1秒重叠0.5秒。将数据转换为单精度。

frameSize = 1*fs;frameRate = 0.5*fs;nframe = (length(allSignals{1})-frameSize)/frameRate + 1;nday =长度(allSignals);myXAll = 0 (frameSize,nframe*nday);XAll = 0 (frameSize,nframe*nday);colIdx = 1;为ii = 1:length(allSignals) XAll(:，colIdx:colIdx+nframe-1) = buffer(allSignals{ii}，frameSize，.．.帧速率,“nodelay”）;colIdx = colIdx+nframe;结束XAll =单(XAll);

由此产生的XAll有550列代表每50天11个记录。每一栏XAll有48,828个样品。

建立小波散射网络。将不变性尺度设置为0.2秒，每个八度的小波滤波器数量在第一层为4，在第二层为1。设置过采样因子为1，优化网络的路径或通道计算。

N = size(XAll,1);sn =小波散射(SignalLength = N, SamplingFrequency = fs，.．.InvarianceScale = 0.2, QualityFactors = [4 1]，.．.OptimizePath = true, overamplingfactor = 1, Precision =“单一”）;

显示每个路径(通道)的散射系数的数量和通道的总数。

[~，pathbyLev] =路径(sn);Ncfs = numCoefficients(sn)

Ncfs = 48

总和(pathbyLev)

Ans = 180

散射变换的输出有180个通道。然而，只有48个时间样本。这意味着每个信号的数据量减少了近6倍。

得到数据的所有小波散射特征。在这个例子中，我们省略了零阶散射系数。因此，通道的数量减少了1到179。如果你有一个支持的GPU，你可以通过设置加速散射变换useGPU来真正的．否则,设置useGPU = false．

useGPU = true;如果useGPU XAll = gpuArray(XAll);结束SAll = sn.featureMatrix(XAll);SAll = SAll(2:end，:，:);npaths = size(SAll,1);scatfeatures = squeeze(num2cell(SAll，[1,2]));

scatfeatures包含550个元素的单元格数组。每个元素包含一秒记录的散射变换。由于这些数据没有标记，我们无法确定哪些记录表明运行正常，哪些记录表明存在种族内部故障。我们从数据描述中所知道的是，在这50天的时间里，一个种族内部的故障会出现。

因此，我们从最早的录音中构建训练集和验证集，以最大限度地提高这些集只包含录音而不存在种族内部故障的概率。前66条记录(6天)用于形成训练集，后44条记录(4天)用于形成验证集。剩下的440个记录(40天)作为测试集。

Ntrain = 6;trainscatFeatures = scatfeatures(1:ntrain*nframe);validationscatFeatures = scatfeatures(ntrain*nframe+1:ntrain*nframe+4*nframe);testscatFeatures = scatfeatures((ntrain*nframe+4*nframe+1):end);

在这个例子中，我们将散射变换与拟合原始时间序列绝对值的网络进行比较。使用绝对值的依据是观察到，在50天期间，振动的振幅似乎平均增加了。

rawfeatures = num2cell(XAll,1)';rawfeatures = cellfun(@转置，rawfeatures,UniformOutput = false);rawABSfeatures = cellfun(@abs,rawfeatures,UniformOutput = false);Ntrain = 6;trainrFeatures = rawABSfeatures(1:ntrain*nframe);validationrFeatures = rawABSfeatures(ntrain*nframe+1:ntrain*nframe+4*nframe);testrFeatures = rawABSfeatures((ntrain*nframe+4*nframe+1):end);

深层网络

在这个例子中，使用了两个深度学习自编码器:一个LSTM和一个卷积自编码器。

LSTM自动编码器在输入端使用z评分归一化，然后在编码阶段使用两个LSTM层，包括179个通道和floor (179/2)渠道分别。编码器中的最后一个LSTM层只从最后一个时间步单元格输出，OutputMode = "last"．随后，一个自定义层，repeatVectorLayer，用于为下一个LSTM层复制这个示例。

Ntimesteps = nfs;lstmAutoEncoder = [sequenceInputLayer(npaths, normalized = .“zscore”，.．.Name =“输入”， MinLength = Ntimesteps) lstmLayer(npaths, Name =“lstm1a”) reluLayer(Name = .“relu1”) lstmLayer(floor(npaths/2)， Name =“lstm2a”， OutputMode =“最后一次”dropoutLayer(0.2, Name = .“drop1”) reluLayer(Name = .“relu2”) repeatVectorLayer(Ntimesteps) lstmLayer(floor(npaths/2)， Name =“lstm2b”dropoutLayer(0.2,Name = .“drop2”) reluLayer(Name = .“relu3”lstmLayer(npaths, Name = .“lstm1b”) reluLayer(Name = .“relu4”) regressionLayer(Name =“回归”));

训练自动编码器200秒。使用16个大小的小批处理，每个历元重新排列训练数据。输出具有最佳验证损耗的网络。

options = trainingOptions(“亚当”，.．.“MaxEpochs”, 200,.．.“MiniBatchSize”, 16岁,.．.“洗牌”，“every-epoch”，.．.“ValidationData”{validationscatFeatures, validationscatFeatures},.．.“阴谋”，“训练进步”，.．.“详细”假的,.．.“OutputNetwork”，“best-validation-loss”）;scatLSTMAutoencoder = trainNetwork(trainscatFeatures).．.lstmAutoEncoder选项);

阈值确定

在异常检测中，没有确定阈值的明确规则。在这个例子中，由于L1范数相对于异常值具有鲁棒性，所以使用了平均绝对误差(MAE)。首先，计算训练、验证和测试数据的MAE误差。

ypredTrain = cellfun(@(x)predict(scatLSTMAutoencoder,x)，trainscatFeatures“UniformOutput”、假);maeTrain = cellfun(@(x,y)maeLoss(x,y)，ypredTrain,trainscatFeatures);ypredValidation = cellfun(@(x)predict(scatLSTMAutoencoder,x)，validationscatFeatures“UniformOutput”、假);maeValid = cellfun(@(x,y)maeLoss(x,y)，ypredValidation,validationscatFeatures);= cellfun(@(x)predict(scatLSTMAutoencoder,x)，testscatFeatures，“UniformOutput”、假);maeTest = cellfun(@(x,y)maeLoss(x,y)，ypredTest,testscatFeatures);如果useGPU [maeTrain,maeValid,maeTest] = gather(maeTrain,maeValid,maeTest);结束

仅使用验证数据来确定异常检测的阈值。本例使用基于验证误差上四分位数加上1.5倍四分位数范围的非参数方法。注意，根据应用程序中假阳性和假阴性的代价，您可以选择一个或多或少保守的阈值。上四分位数加上四分位数范围的1.5倍是一个相当保守的估计，通常会将误报的风险降到最低。

thresh = quantile(maeValid,0.75)+1.5*iqr(maeValid);

确定阈值后，将训练、验证和测试误差绘制为一天的函数。黑色水平线标记疑似异常行为的阈值。

图绘制(.．.(1:长度(maeTrain)) / 11日maeTrain,“b”，.．.(长度(maeTrain) +(1:长度(maeValid))) / 11日maeValid,‘g’，.．.(长度(maeTrain) + (maeValid) +(1:长度(maeTest))) / 11日maeTest,“r”，.．.“线宽”, 1.5)在情节((1:550)/ 11,打*的(550 1),“k”)举行从包含(“天”) ylabel (“美”) xlim([1 50]) legend(“培训”，“确认”，“测试”，“位置”，“西北”)标题(小波散射序列LSTM自编码器网格)在

在第11天和第12天之间出现了异常行为的初步迹象。随后，风力涡轮机似乎显示异常行为几乎连续后30天。这与数据的累积图是一致的。

具有小波散射特性的LSTM自编码器的一个优点是小波散射在时间维度上提供了显著的数据减少，从48,828个样本下降到48个。这使我们能够在不到2分钟的时间内使用GPU训练自动编码器。因为训练健壮的深度学习模型通常涉及调优超参数，拥有一个快速训练的模型是一个显著的优势。

相反，使用GPU对原始数据进行自动编码器训练需要超过1.5小时。因此，在本例中不重复使用原始数据训练自动编码器。这里只介绍结果。下面的LSTM自动编码器在原始数据上进行训练。

Nt =长度(rawfeatures{1});lstmAutoEncoder = [sequenceInputLayer(1, normalized = .“zscore”，.．.Name =“输入”， MinLength = Nt) lstmLayer(32, Name =“lstm1a”) reluLayer(Name = .“relu1”lstmLayer(16, Name = .“lstm2a”， OutputMode =“最后一次”dropoutLayer(0.2, Name = .“drop1”) reluLayer(Name = .“relu2”) repeatVectorLayer(Nt) lstmLayer(16, Name =“lstm2b”dropoutLayer(0.2,Name = .“drop2”) reluLayer(Name = .“relu3”lstmLayer(32, Name = .“lstm1b”) reluLayer(Name = .“relu4”fulllyconnectedlayer (1) regressionLayer(Name = .“回归”));

由于计算上的考虑，LSTM层中隐藏单元的数量用原始数据减少了。此外，使用原始数据和小波散射序列的网络是相同的。

阈值的确定方法与前面描述的完全相同。下图显示了结果。

小波散射序列得到的结果与原始数据有一些相似之处。两者都显示第11天到第12天之间获得的数据是一个异常值。两者也表明异常行为在大约30天后增加规律性。然而，原始数据上的自动编码器似乎对训练数据拟合不足，少数训练记录出现异常。这与我们所知道的风力涡轮机在记录期开始时正常运转的情况不一致。在原始数据上使用自动编码器也会显示出更多的异常。鉴于训练数据的假阳性，有理由怀疑在这些检测中存在一些假阳性。鉴于我们使用的是保守估计，其他阈值方法可能会产生更多的检测。

卷积Autoencoder

虽然循环网络(RNNS)是异常检测的强大架构，但当数据的时间维变大时，RNNS的计算成本很高。需要重申的是，上面使用的LSTM自编码器计算效率高，因为使用了小波散射变换，将数据的时间维数从48,828个样本减少到48个样本。因此，使用GPU训练自动编码器需要不到两分钟。另一方面，训练一个LSTM的原始数据需要超过1.5小时。使用卷积网络可以减小两种方法之间的训练差异。卷积网络模仿LSTM自编码器，在编码阶段使用下采样的卷积层，然后在解码阶段使用转置的上采样的卷积层。

创建一个卷积网络。使用z分数规范化在输入层规范化数据。要对输入进行下采样，请指定重复的一维卷积、ReLU和dropout层块。为了对编码的输入进行上采样，需要包含相同数量的一维转置卷积、ReLU和dropout层块。

对于卷积层，指定大小为11的过滤器数量递减。要确保输出按2的因子均匀采样，请指定步长为2，并设置填充选项“相同”．对于转置卷积层，指定大小为11的过滤器的数量递增。为确保输出被均匀地上采样为2的因子，指定的步长为2，并设置裁剪选项“相同”．

对于退出层，指定退出概率为0.2。要输出与输入信道数相同的序列，指定一个一维转置卷积层，其中滤波器的数量与输入信道数匹配。为确保输出序列与层输入长度相同，设置裁剪选项“相同”．*在输出处使用回归层。

numDownsamples = 2;minLength = 48;filterSize = 11;numFilters = 16;dropoutProb = 0.2;numChannels = npaths;convlayers = [sequenceInputLayer(numChannels, normalized = .“zscore”，.．.最小长度=最小长度,Name =“InputLayer”）%编码器阶段convolution1dLayer (2 * numFilters filterSize,填充=“相同”，Stride=2) reluLayer dropoutLayer(dropoutProb,Name =“dropout1”) convolution1dLayer (filterSize numFilters填充=“相同”，Stride=2) reluLayer dropoutLayer(dropoutProb,Name =“dropout2”）%解码器阶段transposedConv1dLayer (filterSize numFilters裁剪=“相同”，Stride=2) reluLayer dropoutLayer(dropoutProb,Name =“transpdropout1”) transposedConv1dLayer (2 * numFilters filterSize,裁剪=“相同”，Stride=2) reluLayer dropoutLayer(dropoutProb,Name =“transpdropout2”）使通道同意输出transposedConv1dLayer (filterSize numChannels裁剪=“相同”，.．.Name =“FinalConvLayer”) regressionLayer (“名字”，“回归”));

训练卷积网络300个epoch。对每个历元的训练数据进行洗牌。输出具有最佳验证损耗的网络。

options = trainingOptions(“亚当”，.．.MaxEpochs = 300,.．.洗牌=“every-epoch”，.．.ValidationData = {validationscatFeatures, validationscatFeatures},.．.Verbose = 0,.．.OutputNetwork =“best-validation-loss”，.．.情节=“训练进步”）;convNetSCAT = trainNetwork(trainscatFeatures,trainscatFeatures,convlayers,options);

计算与LSTM自动编码器所做的MAE损失。

ypredCTrain = cellfun(@(x)predict(convNetSCAT,x)，trainscatFeatures“UniformOutput”、假);maeCTrain = cellfun(@(x,y)maeLoss(x,y)，ypredCTrain,trainscatFeatures);ypredCValidation = cellfun(@(x)predict(convNetSCAT,x)，validationscatFeatures“UniformOutput”、假);maeCValid = cellfun(@(x,y)maeLoss(x,y)，ypredCValidation,validationscatFeatures);= cellfun(@(x)predict(convNetSCAT,x)，testscatFeatures，“UniformOutput”、假);maeCTest = cellfun(@(x,y)maeLoss(x,y)，ypredCTest,testscatFeatures);如果useGPU [maeCTrain,maeCValid,maeCTest] = gather(maeCTrain,maeCValid,maeCTest);结束

仅使用验证数据来确定异常检测的阈值。使用与LSTM自动编码器相同的阈值确定方法。

threshCV = quantile(maeCValid,0.75)+1.5*iqr(maeCValid);

画出结果。

图绘制(.．.(1:长度(maeCTrain)) / 11日maeCTrain,“b”，.．.(长度(maeCTrain) +(1:长度(maeCValid))) / 11日maeValid,‘g’，.．.(长度(maeCTrain) + (maeCValid) +(1:长度(maeCTest))) / 11日maeTest,“r”，.．.“线宽”, 1.5)在情节((1:550)/ 11,打*的(550 1),“k”)举行从包含(“天”) ylabel (“美”) xlim([1 50]) legend(“培训”，“确认”，“测试”，“位置”，“西北”）;标题(带小波散射序列的卷积网络网格)在

使用相同的卷积网络处理原始数据。将通道数更改为1，以匹配原始数据维度。退出概率为0.2的训练导致训练数据严重不拟合。因此，将退出概率降低到0.1。另外，小波散射序列所使用的网络与原始数据相同。

lgraph = layerGraph(convlayers);inputlayer = sequenceInputLayer(1, normalized =“zscore”，.．.最小长度= 48名=“InputLayerRaw”）;doLayer = dropoutLayer(0.1);tconvLayer = transpsedconv1dlayer (filterSize,1，裁剪=“相同”，.．.Name =“FinalConvLayer”）;rawConvLayers = replaceLayer(lgraph，“InputLayer”, inputlayer);rawConvLayers = replaceLayer(rawConvLayers，“dropout1”, doLayer);rawConvLayers = replaceLayer(rawConvLayers，“dropout2”, doLayer);rawConvLayers = replaceLayer(rawConvLayers，“transpdropout1”, doLayer);rawConvLayers = replaceLayer(rawConvLayers，“transpdropout2”, doLayer);rawConvLayers = replaceLayer(rawConvLayers，“FinalConvLayer”, tconvLayer);

用原始数据训练网络。使用与小波散射序列相同的选项。

options = trainingOptions(“亚当”，.．.MaxEpochs = 300,.．.洗牌=“every-epoch”，.．.ValidationData = {validationrFeatures, validationrFeatures},.．.Verbose = 0,.．.OutputNetwork =“best-validation-loss”，.．.情节=“训练进步”）;convNetRAW = trainNetwork(trainrFeatures,trainrFeatures,rawConvLayers,options);

计算MAE损失并确定阈值。

ypredRTrain = cellfun(@(x)predict(convNetRAW,x)，trainrFeatures，“UniformOutput”、假);maeRTrain = cellfun(@(x,y)maeLoss(x,y)，ypredRTrain,trainrFeatures);ypredRValidation = cellfun(@(x)predict(convNetRAW,x)，validationrFeatures，“UniformOutput”、假);maeRValid = cellfun(@(x,y)maeLoss(x,y)，ypredRValidation,validationrFeatures);= cellfun(@(x)predict(convNetRAW,x)，testrFeatures，“UniformOutput”、假);maeRTest = cellfun(@(x,y)maeLoss(x,y)，ypredRTest,testrFeatures);如果useGPU [maeRTrain,maeRValid,maeRTest] = gather(maeRTrain,maeRValid,maeRTest);结束threshCVraw = quantile(maeRValid,0.75)+1.5*iqr(maeRValid);

画出结果。

图绘制(.．.(1:长度(maeRTrain)) / 11日maeRTrain,“b”，.．.(长度(maeRTrain) +(1:长度(maeRValid))) / 11日maeRValid,‘g’，.．.(长度(maeRTrain) + (maeRValid) +(1:长度(maeRTest))) / 11日maeRTest,“r”，.．.“线宽”, 1.5)在情节(threshCVraw(1:550) / 11日* 1 (550 1),“k”)举行从包含(“天”) ylabel (“美”) xlim([1 50]) legend(“培训”，“确认”，“测试”，“位置”，“西北”）;标题(原始数据卷积网络网格)在

注意，使用卷积自编码器将小波散射序列的训练时间从略低于1.5分减少到使用GPU的大约45秒。最重要的训练时间输入发生在原始数据上，训练LSTM自动编码器大约需要1.5小时，而卷积网络完成训练大约需要3分钟。

对于原始数据和小波散射序列，卷积网络得到的结果与LSTM自编码器的结果相似。与LSTM自编码器的结果一致，基于小波散射的卷积自编码器在第11天到第12天之间表现出异常行为。风力涡轮机的行为再次回到正常范围，直到第30天左右，异常行为的检测开始出现，频率越来越高。

原始数据的结果也与使用LSTM自动编码器获得的结果相似。在训练数据中有一些检测，这些检测很可能表明存在错误检测。这引起了对原始数据网络在接近第15天的记录的早期部分所发生的检测的怀疑。

值得注意的是，使用传统的信号分析技术，如短时傅里叶变换或连续小波变换，对两种方法显示为异常的记录的初步分析并没有显示出任何明显的差异。

讨论

在本例中，我们使用小波散射序列和两种深度自编码器的原始数据来检测风力发电机的内部故障。使用小波散射序列代替原始数据有一些优点。首先，它在时间维度上大大降低了问题的维度。这允许更快速的模型原型化，包括超参数的优化。因为这是深度学习成功应用的关键部分，所以很难夸大LSTM自动编码器的优势。其次，基于小波散射序列训练的深度网络对假检测具有较强的鲁棒性。

卷积自编码器对小波散射序列和原始数据都具有训练时间优势，但对原始数据的相对优势要显著得多。使用具有这两组特性的卷积自编码器在合理的时间内提供了充分的机会来优化超参数。

最后，使用这些不同的网络和不同的功能也是潜在的互补。具体地说，它们提供了通过更详细的信号分析技术更仔细地研究那些所有方法都指定为正常和异常的记录的可能性。在无监督学习问题中，这允许我们增加对数据的理解，并开发更强大的机器监控模型。

参考文献

[1]埃里克·贝赫弗，布兰登·范·赫克，大卫·何。2013。“改进光谱分析的处理”。预测与健康管理学会年会5(1)。

支持功能

函数mae = mean(abs(ypred-target)，“所有”）;结束

另请参阅

waveletScattering|signalDatastore(信号处理工具箱)