使用小波散射与Autoencoders检测异常

这个示例使用:

这个例子展示了如何使用小波散射变换LSTM和卷积autoencoders发展预见性维护的警报系统。散射的例子比较了小波变换+ + autoencoder autoencdoer和原始数据的方法。

数据

收集的数据集是一个2兆瓦风力涡轮机高速轴由20-tooth小齿轮(1]。6秒收购的振动信号采样率为97656 Hz每天从2013-03-07到2013-03-07连续50天。有两个测量2013-03-17,这两天被视为在这个例子。每个timeseries由585936个样本。内部种族断层发达,造成轴承在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)其他的错误(“文件未找到。”+…“\ nManually下载从GitHub库作为. zip文件。”+…“确认. zip文件:\ n % s”parentDir)结束

如果你喜欢使用Git命令克隆存储库,不包括”的文件夹小镇“名称和数据是放置在一个WindTurbineHighSpeedBearingPrognosis-Data文件夹中。

信号数据存储

解开WindTurbineHighSpeedBearingPrognosis-Data-main.zip文件创建一个文件夹和50个.mat文件。每个.mat文件包含两个变量:性心动过速和振动。这个例子只利用振动数据。因此,创建一个信号数据存储只读.mat振动数据文件。

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

因为只有50记录在这个数据集,所有的数据读入内存。

allSignals = readall (sds);

在585936个样本数据记录相对较长。第一次尝试数据简化,检查6信号的随机样本的光谱特征。使用最大重叠与默认的小波和小波包变换的水平。这导致正交分解的信号能量的通频带宽度 $F 年代 / 2^{5}$ Fs是采样率。阴谋每个信号的相对能量的函数的通频带。打印信号能量的百分比低于1/4的采样频率。

rng默认的fs = 97656;rdIdx = randperm(长度(allSignals), 6);图tiledlayout (2、3)为2 = 1:6 nexttile [~ ~ ~ f,再保险]= modwpt (allSignals {2});栏((f。* fs) / 1 e3,再保险。* 100)sprintf (' % 2.2 f '总和(re (f < = 1/4)) * 100)记录= 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)重新取样(x, 1, 2), allSignals, UniformOutput = false);fs = f / 2;

在这个例子中没有标记的数据。获得一个想法的发展50日期间的数据,创建一个累积的阴谋。

tstart = 0;图为2 = 1:长度(allSignals) t =(1:长度(allSignals {2}))。/ fs + tstart;持有在情节(t, allSignals {2}) tstart = t(结束);结束持有从标题(“累积风力涡轮机振动监测”)包含(“时间(sec)——每天6秒50天”)ylabel (“电压”)

从之前的情节,很明显,似乎是振动的数据量逐渐增加的平均时间。经过200秒的录音,这显得尤其明显,对应于时间超出了33天。虽然这平均趋势是明显的,注意个人录音这并非总是如此。附近的一些录音结束的50日期间表现出的行为类似于录音接近开始。

数据准备和特征提取

将每6秒记录分为11个1秒每个重叠0.5秒的录音。把数据单精度。

frameSize = 1 * fs;帧速率= 0.5 * fs;nframe =(长度(allSignals {1}) -frameSize) /帧速率+ 1;nday =长度(allSignals);myXAll = 0 (frameSize, nframe * nday);XAll = 0 (frameSize, nframe * nday);colIdx = 1;为2 = 1:长度(allSignals) XAll (:, colIdx: colIdx + nframe-1) =缓冲区(allSignals {2}, frameSize,…帧速率,“nodelay”);colIdx = colIdx + nframe;结束XAll =单(XAll);

由此产生的XAll有550列代表11为每个50天的录音。每一列的XAll有48828个样本。

建立了小波散射的网络。设置不变性规模为0.2秒,小波滤波器的数量每倍频程4在第一层和第二层1例。设置的过采样因素1和优化路径,或通道计算网络。

N =大小(XAll, 1);sn = waveletScattering (SignalLength = N, SamplingFrequency = fs,…InvarianceScale = 0.2, QualityFactors = (4 - 1),…OptimizePath = true, OversamplingFactor = 1,精度=“单一”);

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

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

Ncfs = 48

总和(pathbyLev)

ans = 180

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

获取所有数据的小波散射特性。在这个例子中,我们忽略第零散射系数。因此,由1到179频道的数量减少。如果你有一个支持GPU,可以万博1manbetx加速通过设置的散射变换useGPU来真正的。否则,设置useGPU = false。

useGPU = true;如果useGPU XAll = gpuArray (XAll);结束萨尔= sn.featureMatrix (XAll);萨尔=萨尔(2:,:,);npaths =大小(萨尔,1);scatfeatures =挤压(num2cell(萨尔,[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):结束);

在这个例子中,我们比较了散射变换网络适合原始时间序列的绝对值。使用的绝对值是基于这样的观察:振动的振幅增加平均50日期间出现。

rawfeatures = num2cell (XAll, 1)”;rawfeatures = cellfun (@transpose 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):结束);

深层网络

在这个例子中,两个深度学习autoencoders使用:一个LSTM和卷积autoencoder。

的LSTM autoencoder使用z分数归一化的输入之后,在编码阶段两个LSTM层组成的179通道和floor (179/2)渠道分别。最后LSTM层只编码器输出从过去的步伐细胞,OutputMode = "最后"。随后,一个自定义的层,repeatVectorLayer,复制这个示例用于下一个LSTM层。

Ntimesteps = Ncfs;lstmAutoEncoder = [sequenceInputLayer (npaths正常化=“zscore”,…Name =“输入”,最小长度= Ntimesteps) lstmLayer (npaths Name =“lstm1a”)reluLayer (Name =“relu1”)lstmLayer(地板(npaths / 2)、Name =“lstm2a”OutputMode =“最后一次”)dropoutLayer (0.2, Name =“drop1”)reluLayer (Name =“relu2”)repeatVectorLayer (Ntimesteps) lstmLayer(地板(npaths / 2),名称=“lstm2b”)dropoutLayer (0.2, Name =“drop2”)reluLayer (Name =“relu3”)lstmLayer (npaths Name =“lstm1b”)reluLayer (Name =“relu4”)regressionLayer (Name =“回归”));

火车200年的autoencoder时代。使用minibatch大小为16,洗牌训练数据每一时代。最好的网络输出验证的损失。

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

阈值确定

没有明确的规则来确定阈值用于异常检测。在这个例子中,平均绝对误差(MAE)由于强劲的行为的L1范数对离群值。首先,计算出美错误的培训、验证和测试数据。

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

只使用验证数据来确定异常检测的阈值。这个例子使用一个基于非参数法上四分位数的验证错误+四分位范围的1.5倍。注意,根据假阳性和假阴性的成本在你的应用程序中,您可以选择一个或多或少的保守的阈值。上四分位数和四分位范围是一个相当保守的估计的1.5倍和意志,在一般情况下,减少误报的风险失踪的真实事件。

打=分位数(maeValid, 0.75) + 1.5 *差(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 (“美”50)xlim([1])传说(“培训”,“确认”,“测试”,“位置”,“西北”)标题(“LSTM Autoencoder与小波散射序列”网格)在

有一个初步迹象天11和12之间的反常行为。随后,风力涡轮机似乎显示反常行为几乎不断后30天。这是与数据的累积图一致。

的一个优点LSTM autoencoder与小波散射特性是重要的数据减少小波散射提供了时间维度,从48828年到48个样品。这使我们能够培养autoencoder使用GPU在不到2分钟。因为训练通常需要调优hyperparameters健壮的深度学习模型,在模型火车很快就有一个显著的优势。

相反,培训一个autoencoder所需原始数据使用GPU超过1.5小时。因此,autoencoder的训练与原始数据不重复这个例子。这里只给出结果。以下LSTM autoencoder训练的原始数据。

元=长度(rawfeatures {1});lstmAutoEncoder = [sequenceInputLayer(1,标准化=“zscore”,…Name =“输入”,最小长度= Nt) lstmLayer(32名=“lstm1a”)reluLayer (Name =“relu1”Name =) lstmLayer(16日“lstm2a”OutputMode =“最后一次”)dropoutLayer (0.2, Name =“drop1”)reluLayer (Name =“relu2”)repeatVectorLayer (Nt) lstmLayer (Name =“lstm2b”)dropoutLayer (0.2, Name =“drop2”)reluLayer (Name =“relu3”)lstmLayer(32名=“lstm1b”)reluLayer (Name =“relu4”)fullyConnectedLayer (1) regressionLayer (Name =“回归”));

由于计算考虑,隐藏单位的数量LSTM层减少了原始数据。否则,网络使用与原始数据和小波散射序列是相同的。

相同的阈值确定方法如前所述。下面的图显示了结果。

有一些相似之处与小波散射的结果序列和原始数据。天11和12之间都显示获得的数据作为一个局外人。也显示了增加的反常行为规律之后大约30天。然而,原始数据上的autoencoder似乎underfit的训练数据和一些培训记录显示为异常。这并不同意我们的知识,风力涡轮机运转正常的记录。还有更多的异常显示通过使用autoencoder原始数据。给定的训练数据上的假阳性,有理由怀疑其中有很多假阳性检测。鉴于我们已经使用一个保守的估计,很可能,其他的阈值方法,将产生更多的检测。

卷积Autoencoder

而周期性网络异常检测(RNNS)强大的架构,RNNS计算昂贵的时间维度数据变得很大。重申,LSTM autoencoder上面使用的计算是有效的,因为小波散射变换的使用减少了时间维度的数据从48828年样本48个样品。因此,培训所需的autoencoder使用GPU不到两分钟。另一方面,培训一个LSTM所需原始数据超过1.5小时。培训两种方法之间的差异可以减少通过使用卷积网络。卷积网络模拟LSTM autoencoder利用卷积层将采样在编码阶段之后,转置卷积与upsampling层在解码阶段。

创建一个卷积网络。正常使用z分数归一化数据在输入层。downsample输入,指定重复的一维卷积,ReLU,辍学层。upsample编码输入,包括一维卷积转置的相同数量的块,ReLU,辍学层。

对于卷积层,指定大小11个数量减少的过滤器。确保输出downsampled均匀的2倍,指定的步伐2,并设置填充选项“相同”。转置的卷积层,指定越来越多的过滤器与大小11。确保输出upsampled均匀的因素2,指定的步伐2,并设置裁剪选项“相同”。

辍学的层,指定一个辍学概率为0.2。与相同数量的渠道输出序列作为输入,指定一个一维卷积转置层的过滤器匹配的输入通道的数量。确保输出序列是相同的长度作为输入层,设置裁剪选项“相同”。*在输出使用回归层。

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

火车300年的卷积网络时代。每个时代洗牌的训练数据。最好的网络输出验证的损失。

选择= trainingOptions (“亚当”,…MaxEpochs = 300,…洗牌=“every-epoch”,…ValidationData = {validationscatFeatures, validationscatFeatures},…Verbose = 0,…OutputNetwork =“best-validation-loss”,…情节=“训练进步”);convNetSCAT = trainNetwork (trainscatFeatures trainscatFeatures、convlayers选项);

美损失计算完成了LSTM autoencoder。

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

只使用验证数据来确定异常检测的阈值。利用相同的阈值确定方法用于LSTM autoencoder。

threshCV =分位数(maeCValid, 0.75) + 1.5 *差(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 (“美”50)xlim([1])传说(“培训”,“确认”,“测试”,“位置”,“西北”);标题(“与小波散射序列卷积网络”网格)在

使用相同的卷积网络原始数据。改变通道1的数量匹配原始数据维度。培训辍学概率为0.2导致了严重的underfitting训练数据。因此,减少辍学概率为0.1。否则,网络使用小波散射序列和原始数据是相同的。

lgraph = layerGraph (convlayers);inputlayer = sequenceInputLayer(1,标准化=“zscore”,…最小长度= 48名=“InputLayerRaw”);doLayer = dropoutLayer (0.1);tconvLayer = transposedConv1dLayer (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);

火车上的网络原始数据。使用相同的选项是使用小波散射序列。

选择= trainingOptions (“亚当”,…MaxEpochs = 300,…洗牌=“every-epoch”,…ValidationData = {validationrFeatures, validationrFeatures},…Verbose = 0,…OutputNetwork =“best-validation-loss”,…情节=“训练进步”);convNetRAW = trainNetwork (trainrFeatures trainrFeatures、rawConvLayers选项);

计算美损失和确定阈值。

ypredRTrain = cellfun (@ (x)预测(convNetRAW x), trainrFeatures,“UniformOutput”、假);maeRTrain = cellfun (@ (x, y) maeLoss (x, y), ypredRTrain, trainrFeatures);ypredRValidation = cellfun (@ (x)预测(convNetRAW x), validationrFeatures,“UniformOutput”、假);maeRValid = cellfun (@ (x, y) maeLoss (x, y), ypredRValidation, validationrFeatures);ypredRTest = cellfun (@ (x)预测(convNetRAW x), testrFeatures,“UniformOutput”、假);maeRTest = cellfun (@ (x, y) maeLoss (x, y), ypredRTest, testrFeatures);如果useGPU [maeRTrain、maeRValid maeRTest] =收集(maeRTrain, maeRValid maeRTest);结束threshCVraw =分位数(maeRValid, 0.75) + 1.5 *差(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 (“美”50)xlim([1])传说(“培训”,“确认”,“测试”,“位置”,“西北”);标题(“卷积网络原始数据”网格)在

注意,使用卷积autoencoder减少了小波散射的训练时间序列从不到半分钟使用GPU大约45秒。发生最重要的训练时间输入的原始数据,在训练LSTM autoencoder需要大约1.5小时,而卷积网络完成了培训大约3分钟。

结果,同卷积网络类似于LSTM autoencoder原始数据和小波散射序列。在协议LSTM autoencoder结果,基于小波散射的卷积autoencoder展览一天11和12之间的反常行为。风力涡轮机的行为再次返回到正常范围,直到大约30天当检测异常行为开始出现越来越频繁。

原始数据的结果也类似于同LSTM autoencoder。有检测训练数据,这可能表明错误的检测。这将导致一些怀疑的检测过程中发生的原始数据网络的早期部分记录附近15天。

重要的是要注意,记录显示为异常的初步分析这两种方法都没有透露任何明确的差异使用传统的信号分析技术,如短时傅里叶变换或连续小波变换。

讨论

在这个例子中,我们使用小波散射序列和原始数据检测内座圈故障的两种类型的深autoencoders风力涡轮机。使用小波散射序列代替原始数据提供了一些优势。首先,它大大降低了问题的维数在时间维度。这允许更快速原型模型包括hyperparameters的优化。因为这是这样一个深入学习的成功应用的关键部分,可以毫不夸张地说,这种优势对LSTM autoencoders。其次,深层网络训练小波散射序列似乎更健壮的反对错误的检测。

提供的卷积autoencoder训练时间优势为小波散射序列和原始数据,但相对优势与原始数据更重要。使用两组的卷积autoencoder特性为优化提供了充足的机会hyperparameters在合理的时间内。

最后,这些不同网络、不同的使用功能也可能互补。具体来说,他们提供的可能性更紧密地调查这些记录所有方法指定为正常与反常的更详细的信号分析技术。在一个无监督学习问题,这让我们增加我们对数据的理解和发展更强大的模型机监控。

引用

[1]Bechhoefer,埃里克,布兰登·范·Hecke他和大卫。2013。“处理改善光谱分析”。预测和健康管理协会的年会5 (1)。

万博1manbetx支持功能

函数美= maeLoss (ypred、目标)美=意味着(abs (ypred-target),“所有”);结束

另请参阅

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