基于MATLAB的深度自动编码器的无监督轴承异常检测
作者：互联网
基于MATLAB的DeepLearnToolbox工具箱（https://github.com/rasmusbergpalm/DeepLearnToolbox），本文在此基础上改成深度自动编码器用于无监督学习，即含有多个隐含层的自动编码器，其输入=输出，简称DAE。后续将其用于轴承故障的异常检测中。
1.相关原理
图1 DAE初始化示意图
对于一个含n个隐含层的DAE，训练时分为预训练与微调两个阶段。预训练阶段时，相邻两层可以通过采用一个单隐层的自动编码器(AE)进行预训练，下一个AE的输入来自上一个AE的隐含层输出，即DAE的上一个隐含层的输出，相关理论可以看这个论文(https://kns.cnki.net/kcms/detail/detail.aspx?FileName=1020832508.nh&DbName=CMFD2021)。上面称之为逐层预训练，训练完毕之后，得到n个AE，利用这些AE的输入权重与输出权重替换DAE中的初始权重，初始化示意图如图1所示，接着在这个基础上利用bp算法进行微调，得到最终的DAE。
2.DAE用于MNIST手写字体重构
%% DAE
clc;clear;close all;format compact
addpath(genpath('../utils/'))
rng(0)
%% 1.导入数据
load mnist_uint8
data=double(train_x);
% method=@mapstd;
method=@mapminmax;
[xs,mapping]=method(data');
train_x=xs';
%% 无监督训练多个自动编码器 假设输入维度为 sizes=[n n1 n2]，则训练两个ae 一个是n-n1-n 一个是n1-n2-n1
sizes=[size(train_x,2) 200 100 50];
learningRate1=0.01;%各ae的学习率
numepochs1=10;%各ae的训练次数
batchsize1=256;%%各ae的batchsize
activation_function='sigm';%各ae的激活函数
dae = daesetup(sizes,activation_function,learningRate1);
opts.numepochs = numepochs1;
opts.batchsize = batchsize1;
opts.show=1;% 0就不显示训练过程
dae = daetrain(dae, train_x, opts);
%% 采用上面无监督预训练好的几个AE来堆栈DAE，构建深度自编码器
learningRate2=0.01;%sae的学习率
numepochs2=100;%dae的训练次数
batchsize2=256;%dae的batchsize
nn = nnsetup([size(train_x,2) 200 100 50 100 200 size(train_x,2)]);%对称的结构，和ae一样，只不过dae是包含多个隐含层的ae
nn.activation_function= 'sigm';
nn.learningRate= learningRate2;
nn=daeunfoldnn(nn,dae);%利用AE的参数初始化DAE
opts.numepochs = numepochs2;
opts.batchsize = batchsize2;
nn.output = 'linear';
opts.show=1;% 0就不显示训练过程
opts.plot = 1;%0就不显示loss曲线
nn = nntrain(nn, train_x, train_x, opts);%输入=输出
disp('训练完毕')
%%
re = nnpredict(nn, train_x);
% 取几个重构的画图
xs=method('reverse',re',mapping);
pred_x=uint8(xs');
pred_x=pred_x(1:25,:);
n=1;
convaspred_x=uint8(zeros(28*5,28*5));
for i = 1:5
for j =1:5
convaspred_x((i-1)*28+1:i*28,(j-1)*28+1:j*28)=reshape(pred_x(n,:),[28,28])';
n=n+1;
end
end
% 取对应的几个真实的画图
xs=method('reverse',train_x',mapping);
real_x=uint8(xs');
real_x=real_x(1:25,:);
n=1;
convasreal_x=uint8(zeros(28*5,28*5));
for i = 1:5
for j =1:5
convasreal_x((i-1)*28+1:i*28,(j-1)*28+1:j*28)=reshape(real_x(n,:),[28,28])';
n=n+1;
end
end
figure
imshow(convaspred_x)
title('reconstruct')
figure
imshow(convasreal_x)
title('original')
训练次数设置的100次，目测loss曲线还在继续下降，可以继续训练，这里就不做展示了。
3 基于DAE的轴承故障无监督异常检测
背景：自然界中含有大量正常样本 而少量异常样本，如果直接构建2分类模型，则容易造成数据不均衡，为此有专家提出无监督异常检测。作为异常检测的经典方法，DAE仅用正常样本进行训练，其输入=输出，以输入输出之间的重构误差最小化作为损失函数，学习的是正常样本的数据分布，当输入异常样本时，由于DAE没有用异常数据进行训练，其重构数据与本身将会产生较大的重构误差，基于此思想，DAE本广泛用于异常检测中。
3.1 数据集构建
采用西储大学轴承故障诊断数据集，48K/0HP数据，共9类故障+一类正常。为模拟自然界中的【大量正常样本 而少量异常样本】这一情况。本文从正常数据中取1000个正常样本，标签记为1；从9类故障中每类取10个样本，共90个故障样本，记为2。
%% 数据预处理（训练集 测试集划分）
clc;clear;close all
%% 加载原始数据
load 0HP/48k_Drive_End_B007_0_122; a1=X122_DE_time'; %1
load 0HP/48k_Drive_End_B014_0_189; a2=X189_DE_time'; %2
load 0HP/48k_Drive_End_B021_0_226; a3=X226_DE_time'; %3
load 0HP/48k_Drive_End_IR007_0_109; a4=X109_DE_time'; %4
load 0HP/48k_Drive_End_IR014_0_174 ; a5=X173_DE_time';%5
load 0HP/48k_Drive_End_IR021_0_213 ; a6=X213_DE_time';%6
load 0HP/48k_Drive_End_OR007@6_0_135 ;a7=X135_DE_time';%7
load 0HP/48k_Drive_End_OR014@6_0_201 ;a8=X201_DE_time';%8
load 0HP/48k_Drive_End_OR021@6_0_238 ;a9=X238_DE_time';%9
load 0HP/normal_0_97 ;a10=X097_DE_time';%10
%% 自然界中含有大量正常样本 而少量异常样本，因此下面生成数据集的时候正常弄多点
N1=1000;
N2=10;
L=1024;% 异常状态（故障类数据）每种状态取N2个样本 每个样本长度为L，而正常状态取N1个样本 每个样本长度为L
data=[];label=[];
for i=1:10
if i==1;ori_data=a1;end
if i==2;ori_data=a2;end
if i==3;ori_data=a3;end
if i==4;ori_data=a4;end
if i==5;ori_data=a5;end
if i==6;ori_data=a6;end
if i==7;ori_data=a7;end
if i==8;ori_data=a8;end
if i==9;ori_data=a9;end
if i==10;ori_data=a10;end
for j=1:N1
if i<=9
if j==N2+1
break
end
end
start_point=randi(length(ori_data)-L);%随机取一个起点
end_point=start_point+L-1;
data=[data ;ori_data(start_point:end_point)];
label=[label;i];
end
end
%% 异常检测标签，除了正常之外，其他的都是异常,正常用1来表示，异常用标签2表示
% 无监督的异常检测方法，训练时仅用正常样本来进行训练
data_normal=data(label==10,:);%注意上面，加载的时候正常是a10
data_abnormal=data(label~=10,:);%不是10的都是异常
% 正常样本取200个以及所有的异常数据作为测试集，用来测试网络性能
index=randperm(size(data_normal,1));
train_x=data_normal(index(1:end-200),:);
test_x=[data_normal(index(end-200+1:end),:);
data_abnormal];
test_y=[ones(size(data_normal(index(end-200+1:end),:),1),1);2*ones(size(data_abnormal,1),1)];%测试集的标签，训练集全是正常的，这里就不做标签了
save result/data_process train_x test_x test_y
3.2 DAE无监督异常检测，DAE结构为1024-10-8-5-8-10-1024（参数是随笔设的，并没有最优）
%% DAE 深度自动编码器 无监督异常检测
clc;clear;close all;format compact
addpath(genpath('utils/'))
rng(0)
%% 1.导入数据
load result/data_process
data=train_x;
% method=@mapstd;
method=@mapminmax;
[xs,mapping]=method(data');
train_x=xs';
%% 无监督训练多个自动编码器 假设输入维度为 sizes=[n n1 n2]，则训练两个ae 一个是n-n1-n 一个是n1-n2-n1
sizes=[size(train_x,2) 10 8 5];
learningRate1=0.01;%各ae的学习率
numepochs1=10;%各ae的训练次数
batchsize1=20;%%各ae的batchsize
activation_function='sigm';%各ae的激活函数
dae = daesetup(sizes,activation_function,learningRate1);
opts.numepochs = numepochs1;
opts.batchsize = batchsize1;
opts.show=1;% 0就不显示训练过程
dae = daetrain(dae, train_x, opts);
%% 采用上面无监督预训练好的几个AE来堆栈DAE，构建深度自编码器
learningRate2=0.01;%sae的学习率
numepochs2=100;%dae的训练次数
batchsize2=20;%dae的batchsize
nn = nnsetup([size(train_x,2) 10 8 5 8 10 size(train_x,2)]);%对称的结构，和ae一样，只不过dae是包含多个隐含层的ae
nn.activation_function= 'sigm';
nn.learningRate= learningRate2;
nn=daeunfoldnn(nn,dae);%利用AE的参数初始化DAE
opts.numepochs = numepochs2;
opts.batchsize = batchsize2;
nn.output = 'linear';
opts.show=1;% 0就不显示训练过程
opts.plot = 1;%0就不显示loss曲线
disp('Train DAE ')
nn = nntrain(nn, train_x, train_x, opts);%输入=输出
disp('训练完毕')
%% 预测 并计算指标
labels0 = nnpredict(nn, train_x);
% 计算每一个重构误差
error0=sum((labels0-train_x).^2,2);
% 以均值+3方差为阈值
threshold=mean(error0)+3*std(error0);
%% 对测试集进行分类
data=test_x;
xs=method('apply',data',mapping);
test_x=xs';
labels = nnpredict(nn, test_x);
% 计算每一个重构误差
error=sum((labels-test_x).^2,2);
figure;plot(error);
%% 计算检测的正确率
pred=[];
pred(error<=threshold,1)=1;
pred(error>threshold,1)=2;
real=test_y;
MatrixClassTable=[real pred];
[EA,OA,AA,Kappa]=ClassifiEvaluate(MatrixClassTable);
disp('----------结果-------------')
disp('各类精度');EA
disp('总体精度');OA
disp('平均精度');AA
disp('kappa系数');Kappa
figure
plot(real,'ro')
hold on ;grid on
plot(pred,'b*')
legend('真实标签','预测标签')
xlabel('测试集样本')
ylabel('标签')
loss curve：
测试集中每个样本的重构误差（头200个为正常数据，后面90个为异常数据）：
以训练集样本的重构误差的均值+3倍方差为阈值，对测试集样本进行判别，大于该阈值的为异常，小于该预测的为正常，得出以下：
混淆矩阵
ConfuMatrix =
199 1
0 90
各类精度
EA =
0.9950
1.0000
总体精度
OA =
0.9966
平均精度
AA =
0.9975
kappa系数
Kappa =
0.9920
可见DAE是有效的
标签：轴承,编码器,nn,DAE,%%,train,MATLAB,data,opts 来源： https://blog.csdn.net/qq_41043389/article/details/117191687