基于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是有效的
完整代码：https://mianbaoduo.com/o/bread/YZmbkppw
转载：https://blog.csdn.net/qq_41043389/article/details/117191687
查看评论