MATLAB基于格拉姆角場與2DCNN-BiGRU的軸承故障診斷模型

本博客來源于CSDN機器魚，未同意任何人轉載。

更多內容，歡迎點擊本專欄目錄，查看更多內容。

0 引言

1 格拉姆角場原理

2 2DCNN-BiGRU網絡結構

3 應用實例

3.1 數據準備

3.2 格拉姆角場數據提取

3.3 網絡模型搭建-重中之重

3.4 使用訓練好的網絡

4 結語

0 引言

為確保對滾動軸承故障診斷的有效性,結合卷積神經網絡(CNN)在圖像特征提取與分類識別的優勢,利用格拉姆角場(CAF)將滾動軸承一維振動信號轉換為二維圖像數據,既保留了信號完整的信息,也保持著信號對于時間的依賴性。并由此提出基于卷積神經網絡與雙向門控循環單元(BiGRU)的診斷模型。首先，將轉換得到的二維圖像作為模型的輸人數據,通過卷積神經網絡提取圖像的空間特征,再由雙向門控循環單元篩選其時間特征,最終由分類器完成模式識別。

1 格拉姆角場原理

具體原理可以看下面這篇文章，另外【博客】提供了Python代碼。

2 2DCNN-BiGRU網絡結構

結構很簡單，就是先將振動信號轉換為GADF或者GASF圖片，然后resize成100*100*3輸入進2DCNN-BiGRU進行特征提取，最后接一個全連接層做分類層，訓練的時候損失函數為交叉熵函數。

在MATLAB中建立這種網絡比較復雜，自帶的深度學習工具箱沒有BiGRU這種函數，我是參考【這里】手動建立的BiGRUlayer，順便說一下MATLAB深度學習工具箱里面很多經典網絡，并且帶有數據，可直接跑。我代碼中具體實現是，分別建立一個CNN與一個BiGRU，先把數據輸入進CNN，再提取出來轉成BiGRU的格式輸入進BiGRU。

%% 模型建立
cnn = [imageInputLayer([100 100 3],'Normalization','none','Name','cnn_channel')convolution2dLayer(3,32,'Padding','same','Name','conv1')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool1')dropoutLayer(0.5,'Name','dropout1')reluLayer('Name','relu1')convolution2dLayer(3,64,'Padding','same','Name','conv2')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool2')dropoutLayer(0.5,'Name','dropout2')reluLayer('Name','relu2')convolution2dLayer(3,64,'Padding','same','Name','conv3')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool3')dropoutLayer(0.5,'Name','dropout3')reluLayer('Name','relu3')convolution2dLayer(3,64,'Padding','same','Name','conv4')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool4')dropoutLayer(0.5,'Name','dropout4')reluLayer('Name','relu4')];%     fullyConnectedLayer(384,'Name','fc1')
lgraph1 = layerGraph(cnn);
dlnet1 = dlnetwork(lgraph1);% bigru層與輸出層 bigru參考
% https://ww2.mathworks.cn/support/search.html/answers/1567853-is-a-bi-gru-available-bidirectional-gated-recurrent-unit-gru-or-a-way-to-implement-a-bi-gru.html?fq[]=asset_type_name:answer&fq[]=category:deeplearning/deep-learning-with-time-series-sequences-and-text&page=1lg=layerGraph();
lg = addLayers(lg, [sequenceInputLayer(448, "Name", "bigruinput")gruLayer( 512, 'OutputMode', 'sequence', ..."Name", "gru1")concatenationLayer(1, 2, "Name", "cat1")gruLayer( 512, 'OutputMode', 'last', ..."Name", "gru3")concatenationLayer(1, 2, "Name", "cat2")fullyConnectedLayer(10, 'Name', 'output')
] );
lg = addLayers( lg, [FlipLayer("flip1")gruLayer( 512, 'OutputMode', 'sequence', ..."Name", "gru2" )FlipLayer("flip2")] );
lg = addLayers(lg, [FlipLayer("flip3")gruLayer( 512, 'OutputMode', 'last', ..."Name", "gru4" )] );lg = connectLayers(lg, "bigruinput", "flip1");
lg = connectLayers(lg, "flip2", "cat1/in2");
lg = connectLayers(lg, "cat1", "flip3");
lg = connectLayers(lg, "gru4", "cat2/in2");
dlnet2 = dlnetwork(lg);

3 應用實例

3.1 數據準備

采用西儲大學軸承故障診斷數據集，48K/0HP數據，共10類故障（正常作為一類特殊的故障類型），劃分后每個樣本的采樣點為864，每類故障各200個樣本，一共得到2000個樣本。博客中用到的數據可在【這里】下載。最后按照7：2：1的比例劃分訓練集、驗證集、測試集。

%% 數據預處理（訓練集 驗證集 測試集劃分）
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%% 
N=200;
L=864;% 每種狀態取N個樣本  每個樣本長度為L
data=[];label=[];
for i=1:10if i==1;ori_data=a1;endif i==2;ori_data=a2;endif i==3;ori_data=a3;endif i==4;ori_data=a4;endif i==5;ori_data=a5;endif i==6;ori_data=a6;endif i==7;ori_data=a7;endif i==8;ori_data=a8;endif i==9;ori_data=a9;endif i==10;ori_data=a10;endfor j=1:Nstart_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
%% 標簽轉換 onehot編碼
N=size(data,1);
output=zeros(N,10);
for i = 1:Noutput(i,label(i))=1;
end
%% 劃分訓練集 驗證集與測試集 7:2:1比例
n=randperm(N);
m1=round(0.7*N);
m2=round(0.9*N);
train_X=data(n(1:m1),:);
train_Y=output(n(1:m1),:);valid_X=data(n(m1+1:m2),:);
valid_Y=output(n(m1+1:m2),:);test_X=data(n(m2+1:end),:);
test_Y=output(n(m2+1:end),:);save data_process train_X train_Y valid_X valid_Y test_X test_Y

3.2 格拉姆角場數據提取

%% 使用格拉姆角場(GAF)以將時間序列數據轉換為圖像
clear;clc;close all
%% 加載信號
load data_process.mat
s=train_X(1,:);
N=864;
fs=48000;
t=0:1/fs:(N-1)/fs;
%%
figure
plot(t,s)
xlabel('t/s')
ylabel('幅值/A')
grid on%% 
window_size = 2;%定義窗口大小為2
transformer = PPA(s,window_size); %分段聚合
transformer = MinMaxScaler(transformer);%歸一化
arccos_X = acos(transformer);%轉換成極坐標
X=GADF(arccos_X);
figure
imagesc(t,t,abs(X)/max(max(abs(X))));
set(gca,'YDir','normal')
axis ij;    %從上到下坐標遞增
xlabel('時間 t/s');
ylabel('時間 t/s');
title('GSDF');X=GASF(arccos_X);
figure
imagesc(t,t,abs(X)/max(max(abs(X))));
set(gca,'YDir','normal')
axis ij;    %從上到下坐標遞增
xlabel('時間 t/s');
ylabel('時間 t/s');
title('GASF');

下面是求角場需要的一些函數，建立4個函數分別保存，不然有時會報錯找不到函數。?

function y=PPA(x,window)
if size(x,1)==1x=x';
end
while mod(size(x,1),window)~=0x=[x; 0];
end
y=reshape(x,[window,length(x)/window])';
y=sum(y,2)/window;
endfunction y=MinMaxScaler(data)
max_data = max(data);
min_data = min(data);
y = ((data-max_data)+(data-min_data))/(max_data-min_data);  %得到了一個歸一化后的矩陣
end
function X=GADF(arccos_X)X = ones(length(arccos_X),length(arccos_X));
for i = 1:length(arccos_X)for j = 1:length(arccos_X)X(i,j) = arccos_X(i)-arccos_X(j); end
end
X=sin(X);
end
function X=GASF(arccos_X)X = ones(length(arccos_X),length(arccos_X));
for i = 1:length(arccos_X)for j = 1:length(arccos_X)X(i,j) = arccos_X(i)+arccos_X(j); end
end
X=cos(X);
end

?舉例得到的GADF、GASF如圖所示。

?然后再寫一個循環，分別提取訓練集、驗證集、測試集的所有數據格拉姆角場圖片并保存下來，這里我們提取的是GADF數據，并采用多線程的方式求每個圖片，速度快一點，代碼如下：

%% 批處理GSDF圖
clear;clc;close all;format compact
if exist('GSDF')rmdir('GSDF','s')
end
mkdir('GSDF');
mkdir('GSDF/train_img/');
mkdir('GSDF/valid_img/');
mkdir('GSDF/test_img/');
%% 加載信號
load data_process.mat
N=864;
fs=48000;
t=0:1/fs:(N-1)/fs;
window_size = 2;%定義窗口大小為2
%% 訓練集
[m,n]=size(train_X);
[~,label]=max(train_Y,[],2);
p=parpool(maxNumCompThreads);%圖片太多 我們將采用多線程生成圖片 大大節約時間
%如果電腦不支持多線程報錯 就把p=parpool(maxNumCompThreads) 與delete(p) delete(gcp('nocreate'))刪掉 把3個地方的parfor改成for
parfor  i=1:mitransformer = PPA(train_X(i,:),window_size); %分段聚合transformer = MinMaxScaler(transformer);%歸一化arccos_X = acos(transformer);%轉換成極坐標coefs=GADF(arccos_X);imagesc(t,t,abs(coefs)/max(max(abs(coefs))));
%     set(gcf,'PaperPositionMode','auto');set(gca,'position',[0 0 1 1])set(gcf,'position',[0,0,250,250]);fname=['GSDF/train_img/',num2str(i),'-',num2str(label(i)),'.jpg'];saveas(gcf,fname);
end%% 驗證集
[m,n]=size(valid_X);
[~,label]=max(valid_Y,[],2);
parfor i=1:mitransformer = PPA(valid_X(i,:),window_size); %分段聚合transformer = MinMaxScaler(transformer);%歸一化arccos_X = acos(transformer);%轉換成極坐標coefs=GADF(arccos_X);imagesc(t,t,abs(coefs)/max(max(abs(coefs))));set(gca,'position',[0 0 1 1])set(gcf,'position',[0,0,250,250]);fname=['GSDF/valid_img/',num2str(i),'-',num2str(label(i)),'.jpg'];saveas(gcf,fname);
end
%% 測試集
[m,n]=size(test_X);
[~,label]=max(test_Y,[],2);
parfor i=1:mitransformer = PPA(test_X(i,:),window_size); %分段聚合transformer = MinMaxScaler(transformer);%歸一化arccos_X = acos(transformer);%轉換成極坐標coefs=GADF(arccos_X);imagesc(t,t,abs(coefs)/max(max(abs(coefs))));set(gca,'position',[0 0 1 1])set(gcf,'position',[0,0,250,250]);fname=['GSDF/test_img/',num2str(i),'-',num2str(label(i)),'.jpg'];saveas(gcf,fname);
end
delete(p);
delete(gcp('nocreate'))

3.3 網絡模型搭建-重中之重

MATLAB搭建復合網絡比較復雜，主函數如下所示：

%% 1，清空變量
clc;clear;close all;format compact
%% 2.加載數據
disp('讀取訓練集圖片用于訓練模型')
[traindata,traintarget]=readimgdata('GSDF/train_img',[100 100],1);
disp('讀取驗證集圖片用于訓練過程中進行驗證')
[XValidation,YValidation]=readimgdata('GSDF/valid_img',[100 100],1);trainT = onehotencode(traintarget',2)';
validT = onehotencode(YValidation',2)';
%% 模型建立
cnn = [imageInputLayer([100 100 3],'Normalization','none','Name','cnn_channel')convolution2dLayer(3,32,'Padding','same','Name','conv1')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool1')dropoutLayer(0.5,'Name','dropout1')reluLayer('Name','relu1')convolution2dLayer(3,64,'Padding','same','Name','conv2')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool2')dropoutLayer(0.5,'Name','dropout2')reluLayer('Name','relu2')convolution2dLayer(3,64,'Padding','same','Name','conv3')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool3')dropoutLayer(0.5,'Name','dropout3')reluLayer('Name','relu3')convolution2dLayer(3,64,'Padding','same','Name','conv4')maxPooling2dLayer(2,'Stride',2,'Padding','same','Name','pool4')dropoutLayer(0.5,'Name','dropout4')reluLayer('Name','relu4')];%     fullyConnectedLayer(384,'Name','fc1')
lgraph1 = layerGraph(cnn);
dlnet1 = dlnetwork(lgraph1);% bigru層與輸出層 bigru參考
% https://ww2.mathworks.cn/support/search.html/answers/1567853-is-a-bi-gru-available-bidirectional-gated-recurrent-unit-gru-or-a-way-to-implement-a-bi-gru.html?fq[]=asset_type_name:answer&fq[]=category:deeplearning/deep-learning-with-time-series-sequences-and-text&page=1lg=layerGraph();
lg = addLayers(lg, [sequenceInputLayer(448, "Name", "bigruinput")gruLayer( 512, 'OutputMode', 'sequence', ..."Name", "gru1")concatenationLayer(1, 2, "Name", "cat1")gruLayer( 512, 'OutputMode', 'last', ..."Name", "gru3")concatenationLayer(1, 2, "Name", "cat2")fullyConnectedLayer(10, 'Name', 'output')
] );
lg = addLayers( lg, [FlipLayer("flip1")gruLayer( 512, 'OutputMode', 'sequence', ..."Name", "gru2" )FlipLayer("flip2")] );
lg = addLayers(lg, [FlipLayer("flip3")gruLayer( 512, 'OutputMode', 'last', ..."Name", "gru4" )] );lg = connectLayers(lg, "bigruinput", "flip1");
lg = connectLayers(lg, "flip2", "cat1/in2");
lg = connectLayers(lg, "cat1", "flip3");
lg = connectLayers(lg, "gru4", "cat2/in2");
dlnet2 = dlnetwork(lg);
%% 網絡訓練參數
numEpochs = 150;% 訓練代數
miniBatchSize = 128;%batchsize
learnRate = 0.01;%初始學習率
decay = 0.9;%學習率衰減
momentum = 0.9;%動量
%建一個圖 來展示實時loss曲線
figure
lineLossTrain = animatedline('Color',[0.85 0.325 0.098]);
ylim([0 inf])
xlabel("Iteration")
ylabel("Loss")
grid on
%% 開始訓練
train_again=1;% 為1就代碼重新訓練模型，為0就是調用訓練好的網絡
if train_again==1velocity1 = [];velocity2 = [];iteration = 0;start = tic;trainloss=zeros(numEpochs,1);for epoch = 1:numEpochs% 每一代 都先打亂數據，然后分批次遍歷所有數據[x,y]=shuffle_data(traindata,trainT);batchs=ceil(size(y,2)/miniBatchSize);% 計算總批次all_loss=0;for i=1:batchsiteration = iteration + 1;% 讀取第i個批次的數據[dX,T] = next_batch(x,y,miniBatchSize,i);% 計算梯度 lossdX=gpuArray(dX);T=gpuArray(T);[gradients1, gradients2,loss] = dlfeval(@modelGradients,dlnet1,dlnet2,dX,T,miniBatchSize);% 利用 SGDM 優化器更新網絡參數.[dlnet1.Learnables, velocity1] = sgdmupdate(dlnet1.Learnables, gradients1, velocity1, learnRate, momentum);[dlnet2.Learnables, velocity2] = sgdmupdate(dlnet2.Learnables, gradients2, velocity2, learnRate, momentum);all_loss = all_loss+loss;endif mod(epoch,100)==0% 學習率衰減learnRate = learnRate*decay;endtrainloss(epoch)=all_loss;% 實時顯示lossD = duration(0,0,toc(start),'Format','hh:mm:ss');addpoints(lineLossTrain,epoch,all_loss/size(y,2))title("Epoch: " + epoch + ", Elapsed: " + string(D))drawnowendsave result/net dlnet1 dlnet2 
elseload result/net
end%% 測試集正確率
[XValidation,YValidation]=readimgdata('GSDF/test_img',[100 100],1);
pred=modelPredict(dlnet1,dlnet2,XValidation,miniBatchSize);
label=onehotdecode(pred,[1,2,3,4,5,6,7,8,9,10],1);
acc=sum(label==YValidation)/length(label)
figure
cm = confusionchart(YValidation,label,'Normalization','row-normalized');
xlabel('真實標簽')
ylabel('預測標簽')
save result/cnnigruresult acc%% 可視化
xb = dlarray(XValidation,'SSCB');
xb=gpuArray(xb);
dlYPred1 = predict(dlnet1,xb);%predict
dlYPred1=reshape(dlYPred1,[7*64,7,size(xb,4)]);
dlYPred2=dlarray(dlYPred1,'CTB');
testFeatures= predict(dlnet2,dlYPred2);
testFeatures=gather(extractdata(testFeatures));
Y = tsne(testFeatures');figure
gscatter(Y(:,1),Y(:,2),YValidation)
legend('1','2','3','4','5','6','7','8','9','10')
title('測試集特征可視化')

然后定義需要用到的子函數：

1.讀取圖片數據的子函數

function [data,label]=readimgdata(path,size,regular)
tic
list=dir([path,'/','*.jpg']);
N=length(list);
for i=1:N% 提取圖片，并轉換為 64*64*3（size）filepath=[path,'/',list(i).name];img=imread(filepath);img=double(imresize(img,size));
%     img1(1,:,:,:)=img;data(:,:,:,i)=img;%提取標簽，10-1.ipg 10代表第10張，1代表其標簽為1類S1 = regexp(list(i).name, '\-', 'split');S=regexp(S1{2}, '\.', 'split');label(:,i)=str2num(S{1});
endif regulardata=data/255;
end
label=categorical(label);
disp('讀取圖片耗時')
toc

2.打亂圖像的子函數

function [x,y]=shuffle_data(x,y)
n=length(y);
a=randperm(n);x=x(:,:,:,a);
y=y(:,a);
end

3.將數據轉換為CNN輸入格式的函數

function [xb,yb] = next_batch(x,y,bs,i)
start=(i-1)*bs+1;
ennd=start+bs-1;
if ennd>length(y)start=length(y)-bs;ennd=start+bs-1;
end
xb=x(:,:,:,start:ennd);
yb=y(:,start:ennd);% 轉換為matlab深度學習工具箱的輸入格式，這個是必須的
xb = dlarray(xb,'SSCB');
yb = dlarray(yb,'CB');
end

4.FlipLayer這個函數是用來構建BiGRU的子函數，詳細原理看【這里】，代碼如下：

classdef FlipLayer < nnet.layer.Layermethodsfunction layer = FlipLayer(name)layer.Name = name;endfunction Y = predict(~, X)Y = flip(X,3);endend
end

5.計算梯度的函數，用于反向更新網絡參數，就是將數據輸入到CNN中，然后提取CNN的輸出轉換成BiGRU的輸入，再輸入進BiGRU中，最后提取BiGRU的輸出與真實值計算交叉熵損失，代碼如下：

function [gradients1, gradients2,loss] = modelGradients(dlnet1,dlnet2,X, T,miniBatchSize)
% CNN通道前向計算
dlYPred1 = forward(dlnet1,X);
% size(dlYPred1)
dlYPred1=reshape(dlYPred1,[7*64,7,miniBatchSize]);%time distributed 論文里面算出來是6*6*64  matlab出來是7*7*64
dlYPred1=dlarray(dlYPred1,'CTB');
% BiGRU通道前向計算
output= forward(dlnet2,dlYPred1);
% size(output)
output=dlarray(output,'CB');
output = softmax(output);
loss = crossentropy(output,T);% 計算梯度
[gradients1,gradients2] = dlgradient(loss,dlnet1.Learnables,dlnet2.Learnables);
loss = double(gather(extractdata(loss)));
end

6.模型預測函數

function pred=modelPredict(dlnet1,dlnet2,data,batchsize)
N_samples=size(data,4);
batchs=ceil(N_samples/batchsize);% 計算總批次
pred=[];
for i =1:batchsxb = dlarray(data(:,:,:, (i-1)*batchsize+1: min(i*batchsize,N_samples)),'SSCB');xb=gpuArray(xb);dlYPred1 = predict(dlnet1,xb);%predictdlYPred1=reshape(dlYPred1,[7*64,7,size(xb,4)]);dlYPred1=dlarray(dlYPred1,'CTB');output= predict(dlnet2,dlYPred1);output=extractdata(output);%提取數據 將dlarray轉為一般的數組pred =[pred softmax(output)];
end
end

3.4 使用訓練好的網絡

clc;clear;close all;format compact
%% 加載訓練好的網絡
load result/net
%% 讀取一個信號樣本作為待測樣本 轉換為圖片  并輸入訓練好的分類網絡中
load data_process
data=test_X(1,:);
N=864;
fs=48000;
t=0:1/fs:(N-1)/fs;%% 
window_size = 2;%定義窗口大小為2
transformer = PPA(data,window_size); %分段聚合
transformer = MinMaxScaler(transformer);%歸一化
arccos_X = acos(transformer);%轉換成極坐標
X=GADF(arccos_X);
imagesc(t,t,abs(X)/max(max(abs(X))));
set(gca,'position',[0 0 1 1])
set(gcf,'position',[0,0,250,250]);
saveas(gcf,'test.jpg');
%% 轉為file4中cnn的輸入形式
testD=double(imresize(imread('test.jpg'),[100,100]))/255;
pred=modelPredict(dlnet1,dlnet2,testD,1);
label=onehotdecode(pred,[1,2,3,4,5,6,7,8,9,10],1);
disp(['該待測樣本標簽為：',num2str(double(label))])