Crack Identification From Accelerometer Data（二）——实现GRU循环神经网络

本文基于matlab2020版官方网页DocumentationCrack Identification From Accelerometer Data及个人理解。

该示例显示了如何使用小波wavelet和深度学习技术来检测横向路面裂缝并确定其位置。该示例演示了将小波散射序列用作门控循环单元（GRU）和一维卷积网络的输入，以便根据是否存在裂缝对时间序列进行分类。数据是从安装在前排乘客座椅车轮的转向节上的传感器获得的垂直加速度测量值。及早发现发展中的横向裂缝对于评估和维护路面性能非常重要。可靠的自动检测方法可实现更频繁，更广泛的监视。

本文主要介绍其中的第二部分，GRU循环神经网络的设置、训练和测试。GRU（Gate Recurrent Unit）是循环神经网络（Recurrent Neural Network, RNN）的一种，是为了解决长期记忆和反向传播中的梯度等问题而提出来的深度学习算法。

GRU训练集

计算每个训练序列的散射序列。散射序列存储在一个单元阵列中，以与GRU网络兼容。这些序列随后被重塑成一个4-D数组输入，用于1-D卷积网络。

XTrainSCAT = cell(size(TrainData));% GRU训练集
for kk = 1:numel(TrainData)
    XTrainSCAT{kk} = helperscat(TrainData{kk});
end
npaths = cellfun(@(x)size(x,1),XTrainSCAT);
inputSize = npaths(1);

GRU网络结构

使用两个有30个隐藏单位的GRU图层和两个dropout图层。输入的大小是散射路径的数目。要在原始时间序列上训练这个网络，将inputSize更改为1，并将每个时间序列转置为一个行向量(1×461)。

numHiddenUnits1 = 30;% 隐含单元数
numHiddenUnits2 = 30;
numClasses = 2;% 类别

GRUlayers = [ ...
    sequenceInputLayer(inputSize,'Name','InputLayer',...
        'Normalization','zerocenter')
    gruLayer(numHiddenUnits1,'Name','GRU1','OutputMode','sequence')
    dropoutLayer(0.35,'Name','Dropout1')
    gruLayer(numHiddenUnits2,'Name','GRU2','OutputMode','last')
    dropoutLayer(0.2,'Name','Dropout2')
    fullyConnectedLayer(numClasses,'Name','FullyConnected')
    softmaxLayer('Name','smax');
    classificationLayer('Name','ClassificationLayer')];

进行GRU网络的训练

通过options设置训练相关参数

maxEpochs = 150;%最大轮次
miniBatchSize = 15;%小批量

options = trainingOptions('adam', ...
    'ExecutionEnvironment','cpu', ...
    'L2Regularization',1e-3,...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'Shuffle','every-epoch',...
    'Verbose',0,...
    'Plots','none');
iterPerEpoch = floor(length(XTrainSCAT)/miniBatchSize);%一个轮次的迭代次数
[scatGRUnet,infoGRU] = trainNetwork(XTrainSCAT,TrainLabels,GRUlayers,options);%训练，返回网络和训练信息

绘制训练过程中每个迭代的准确度和损失

figure
subplot(2,1,1)
smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoGRU.TrainingAccuracy);
smoothedLoss = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoGRU.TrainingLoss);
plot(smoothedACC)
title(['Training Accuracy (Smoothed) '...
    num2str(iterPerEpoch) ' iterations per epoch'])
ylabel('Accuracy (%)')
ylim([0 100.1])
grid on
xlim([1 length(smoothedACC)])
subplot(2,1,2)
plot(smoothedLoss)
ylim([-0.01 1])
grid on
xlim([1 length(smoothedLoss)])
ylabel('Loss')
xlabel('Iteration')

运行结果

GRU测试

计算测试集的小波散射序列，作为GRU的测试输入

XTestSCAT = cell(size(TestData));
for kk = 1:numel(TestData)
    XTestSCAT{kk} = helperscat(TestData{kk});
end

进行GRU测试，绘制混淆矩阵

miniBatchSize = 15;
ypredSCAT = classify(scatGRUnet,XTestSCAT, ...
    'MiniBatchSize',miniBatchSize);
figure
confusionchart(TestLabels,ypredSCAT,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'GRU network -- Wavelet Scattering Sequences'; ...
    'Confusion chart with precision and recall'});

尽管类的不平衡较大，数据集又小，但精度和召回值表明网络对测试数据的识别效果良好。具体来说，“CR”数据的精度值约为100%，召回率约为95%。考虑到训练集中有整整67%的记录是“UNCR”，这样的结果已经很好了。尽管存在不平衡，但网络并没有因为过度学习而将时间序列分类为“UNCR”。

原始数据的GRU测试结果

如果设置inputSize = 1并在训练数据中转置时间序列，就可以在原始时间序列数据上对GRU网络进行再训练。这是在训练集中的相同数据上完成的，可以加载该网络并检查测试集中的性能。

load(fullfile(saveFolder,"transverse_crack_latest","tsGRUnet.mat"));
rawTest = cellfun(@transpose,TestData,'UniformOutput',false);
miniBatchSize = 15;
YPredraw = classify(tsGRUnet,rawTest, ...
    'MiniBatchSize',miniBatchSize);
confusionchart(TestLabels,YPredraw,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'GRU network -- Raw Time Series'; ...
    'Confusion chart with precision and recall'});

对于这个网络，性能不是很好。具体来说，对“CR”数据的召回很糟糕。“UNCR”数据的假阴性数量相当大。当模型不能很好的学习不平衡数据集时，就会产生这样的结果。

References

[1] Yang, Qun and Shishi Zhou. “Identification of asphalt pavement transverse cracking based on vehicle vibration signal analysis.”, Road Materials and Pavement Design, 2020, 1-19. https://doi.org/10.1080/14680629.2020.1714699.

[2] Zhou,Shishi. “Vehicle vibration data.” https://data.mendeley.com/datasets/3dvpjy4m22/1. Data is used under CC BY 4.0. Data is repackaged from original Excel data format to MAT-files. Speed label removed and only “crack” or “nocrack” label retained.

附录

本示例中使用的辅助函数如下

function smat = helperscat(datain)
% This helper function is only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
datain = single(datain);

sf = waveletScattering('SignalLength',length(datain),...
    'OversamplingFactor',1,'qualityfactors',[8 1],...
    'InvarianceScale',0.05,'Boundary','reflect','SamplingFrequency',1280);
smat = sf.featureMatrix(datain);

end

%-----------------------------------------------------------------------
function dataUL = equalLenTS(data)
% This function in only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
N = length(data);
dataUL = cell(N,1);
for kk = 1:N
    L = length(data{kk});
    switch L
        case 461
            dataUL{kk} = data{kk};
        case 369
            Ndiff = 461-369;
            pad = Ndiff/2;
            dataUL{kk} = [flip(data{kk}(1:pad)); data{kk} ; ...
                flip(data{kk}(L-pad+1:L))];     
        otherwise
            Ndiff = L-461;
            zrs = Ndiff/2;
            dataUL{kk} = data{kk}(zrs:end-zrs-1);
    end
end       

end

%--------------------------------------------------------------------------
function [fmat,x] = featureVectors(data,sf)
% This function is only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
data = single(data);
N = length(data);
dt = 1/1280;
if N < 461
    Ndiff = 461-N;
    pad = Ndiff/2;
    dataUL = [flip(data(1:pad)); data ; ...
               flip(data(N-pad+1:N))];   
     rate = 5e4/3600;
     dx = rate*dt;
     x = 0:dx:(N*dx)-dx;     
    
elseif N > 461
    Ndiff = N-461;
    zrs = Ndiff/2;
    dataUL = data(zrs:end-zrs-1);
    rate = 3e4/3600;
    dx = rate*dt;
    x = 0:dx:(N*dx)-dx;
else
    dataUL = data;
    rate = 4e4/3600;
    dx = rate*dt;
    x = 0:dx:(N*dx)-dx;
end
fmat = sf.featureMatrix(dataUL);
fmat = reshape(fmat',[58 1 38]);
end

%------------------------------------------------------------------------------
function [mra,chngpts] = helperMRA(data,x)
% This function is only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
mra = modwtmra(modwt(data,'sym3'),'sym3');
mraLev = mra(4:6,:);
Ns = size(mraLev,1);
thresh = [2, 4, 8];
chngpts = false(size(mraLev));
% Determine changepoints. We want different thresholds for different
% resolution levels.
for ii = 1:Ns
    chngpts(ii,:) = ischange(mraLev(ii,:),"linear",2,"Threshold",thresh(ii));
end

for kk = 1:Ns
    idx = double(chngpts(kk,:));
    idx(idx == 0) = NaN;    
    subplot(Ns,1,kk)
    plot(x,mraLev(kk,:))
    if kk == 1
        title('MRA Components')
    end
    yyaxis right
    hs = stem(x,idx);
    hs.ShowBaseLine = 'off';
    hs.Marker = '^';
    hs.MarkerFaceColor = [1 0 0];
end
grid on
axis tight
xlabel('Distance (m)')
end