Main Content

Crack Identification From Accelerometer Data

This example shows how to use wavelet and deep learning techniques to detect transverse pavement cracks and localize their position. The example demonstrates the use of wavelet scattering sequences as inputs to a gated recurrent unit (GRU) and 1-D convolutional network to classify time series based on the presence or absence of a crack. The data are vertical acceleration measurements obtained from a sensor mounted on the suspension knuckle of the front passenger seat wheel. Early identification of developing transverse cracks is important for pavement performance evaluation and maintenance. Reliable automatic detection methods enable more frequent and extensive monitoring.

Please read the Data - Description and Required Attributions before running this example. All data import and preprocessing is described in the Data — Import and the Data — Preprocessing sections. If you want to skip the training-test set creation, preprocessing, feature extraction, and model training, you can go directly to the Classification and Analysis section. There you can load the preprocessed test data as well as the extracted features and trained models.

Data — Description and Required Attributions

The data used in this example was retrieved from the Mendeley Data open data repository [2]. The data is distributed under a Creative Commons (CC) BY 4.0 license. Download the data and models used in this example in your temporary directory specified by MATLAB® tempdir command. If you choose to download the data in a folder different from tempdir, change the directory name in the subsequent instructions. Unzip the data into a folder specified as TransverseCrackData.

dataURL = 'https://ssd.mathworks.com/supportfiles/wavelet/crackDetection/transverse_crack.zip';
saveFolder = fullfile(tempdir,'TransverseCrackData'); 
zipFile = fullfile(tempdir,'transverse_crack.zip');
websave(zipFile,dataURL);
unzip(zipFile,saveFolder)

After you unzip the data, the TransverseCrackData folder contains a subfolder called transverse_crack_latest. All subsequent commands must be run in this folder, or you can place this folder on the matlabpath.

The text file, vehiclevibrationdata.rights, included in the zip file contains the text of the CC BY 4.0 license. The data has been repackaged from the original Excel format into MAT-files.

Data acquisition is described in detail in [1]. Twelve four meter long sections of asphalt containing a centrally-located transverse crack and twelve equal-length uncracked sections were used. The data is obtained from three separate roads. The transverse cracks ranged in width from 2-13 mm with crack spacings from 7-35 mm. The sections were driven at three different speeds: 30 km/hr, 40 km/hr, and 50 km/hr. Vertical acceleration measurements from the front passenger suspension knuckle are acquired at a sampling frequency of 1.28 kHz. The speeds of 30, 40, and 50 km/hr correspond to sensor measurements every 6.5 mm at 30 km/hr, 8.68 mm at 40 km/hr, and 10.85 mm at 50 km/hour. See [1] for a detailed wavelet analysis of these data distinct from the analyses in this example.

In keeping with the stipulations of the CC BY 4.0 license, we note that the speed information of the original data is not retained in the data used in this example. The speed and road information are retained in the road1.mat, road2.mat, road3.mat data files included in the data folder for completeness.

Data — Download and Import

Load the accelerometer data and their corresponding labels. There are 327 accelerometer recordings.

load(fullfile(saveFolder,"transverse_crack_latest","allroadData.mat"));
load(fullfile(saveFolder,"transverse_crack_latest","allroadLabel.mat"));

Data — Preprocessing

Obtain the length of all the time series. Display a bar graph of the number of time series per length.

tslen = cellfun(@length,allroadData);
uLen = unique(tslen);
Ng = histcounts(tslen);
Ng = Ng(Ng > 0);
bar(uLen,Ng,0.5)
grid on
AX = gca;
AX.YLim = [0 max(Ng)+15];
text(uLen(1),Ng(1)+10,num2str(Ng(1)))
text(uLen(2),Ng(2)+10,num2str(Ng(2)))
text(uLen(3),Ng(3)+10,num2str(Ng(3)))
xlabel('Length in Samples')
ylabel('Number of Series')
title('Time Series Length')

There are three unique lengths in the dataset: 369, 461, and 615 samples. The vehicle is traveling at three different speeds but the distance traversed and sample rate is constant resulting in different data record lengths. Determine how many records are in the "Cracked" (CR) and "Uncracked" (UNCR) classes.

tabulate(allroadLabel)
  Value    Count   Percent
     CR      109     33.33%
   UNCR      218     66.67%

This dataset is significantly imbalanced. There are twice as many time series without a crack (UNCR) as series containing a crack (CR). This means that a classifier which predicts "Uncracked" on each record would achieve an accuracy of 67% without any learning.

The time series are also of different lengths. To use a convolution network, a common input shape is needed. In recurrent networks it is possible to use unequal length time series as inputs, but all time series in a mini-batch are padded or truncated based on the training options. This requires care in creating mini-batches for both training and testing to ensure the proper distribution of padded sequences. Further, it requires that you do not shuffle the data during training. With this small dataset, shuffling the training data for each epoch is desirable. Accordingly, a common time series length is used.

The most common length is 461 samples. Further, the crack, if present, is centrally located in the recording. Accordingly, we can symmetrically extend the series with 369 samples to length 461 by reflecting the initial and final 46 samples. In the recordings with 615 samples, remove the initial 77 and final 77 samples.

Training — Feature Extraction and Network Training

The following sections generate the training and test sets, create the wavelet scattering sequences, and train both gated recurrent unit (GRU) and 1-D convolutional networks. First, split the data and corresponding labels into "Cracked" vs "Uncracked" sets.

[~,~,pop] = unique(allroadLabel);
CRidx = find(pop == 1);
UCRidx = find(pop == 2);
CRdata = allroadData(CRidx);
CRlabel = allroadLabel(CRidx);
UCRdata = allroadData(UCRidx);
UCRlabel = allroadLabel(UCRidx);

Extend or truncate the series in both sets to obtain a common length of 461 samples.

CRdataUL = equalLenTS(CRdata);
UCRdataUL = equalLenTS(UCRdata);

Now each time series in both the cracked and uncracked datasets has 461 samples. Split the data in a training set consisting of 80% of the time series in each class and hold out the remaining 20% of each class for testing.

NtrainCR = floor(0.8*length(CRlabel));
NtrainUCR = floor(0.8*length(UCRlabel));
idxTrainCR = randperm(length(CRlabel),NtrainCR);
idxTrainUCR = randperm(length(UCRlabel),NtrainUCR);
idxTestCR = setdiff(1:length(CRlabel),idxTrainCR);
idxTestUCR = setdiff(1:length(UCRlabel),idxTrainUCR);

Create the training and test sets and print a summary of each class.

TrainData = [CRdataUL(idxTrainCR) ; UCRdataUL(idxTrainUCR)];
TrainLabels = [CRlabel(idxTrainCR) ;  UCRlabel(idxTrainUCR)];
TestData = [CRdataUL(idxTestCR) ; UCRdataUL(idxTestUCR)];
TestLabels = [CRlabel(idxTestCR) ; UCRlabel(idxTestUCR)];
tabulate(TrainLabels)
  Value    Count   Percent
     CR       87     33.33%
   UNCR      174     66.67%
tabulate(TestLabels)
  Value    Count   Percent
     CR       22     33.33%
   UNCR       44     66.67%

Shuffle the data and labels once before training.

idxS = randperm(length(TrainData));
TrainData = TrainData(idxS);
TrainLabels = TrainLabels(idxS);
idxS = randperm(length(TrainData));
TrainData = TrainData(idxS);
TrainLabels = TrainLabels(idxS);

Compute the scattering sequences for each of the training series. The scattering sequences are stored in a cell array to be compatible with the GRU network. These sequences are subsequently reshaped into a 4-D array input for use with a 1-D convolutional network.

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

Training — GRU Network

Construct the GRU network layers. Use two GRU layers with 30 hidden units each as well as two dropout layers. The input size is the number of scattering paths. To train this network on the raw time series, change the inputSize to 1 and transpose each time series to a row vector (1-by-461). If you wish to skip the network training, you may go directly to the Classification and Analysis section. There you can load the trained GRU network as well as the preprocessed training and test sets.

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')];

Train the GRU network. Use a mini-batch size of 15 with 150 epochs.

maxEpochs = 150;
miniBatchSize = 15;

options = trainingOptions('adam', ...
    'ExecutionEnvironment','gpu', ...
    '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);

Plot the smoothed training accuracy and loss per iteration.

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')

Obtain the wavelet scattering transforms of the held-out test data for classification.

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

Training — 1-D Convolutional Network

Train a 1-D convolutional network with wavelet scattering sequences. The scattering sequences are 38-by-58 where 58 is the number of time steps and 38 is the number of scattering paths. To use this in a 1-D convolutional network, transpose and reshape the scattering sequences to be 58-by-1-by-38-by-Nsig where Nsig is the number of training examples. If you wish to skip the network training, you may go directly to the Classification and Analysis section. There you can load the trained convolutional network as well as the preprocessed training and test sets.

scatCONVTrain = zeros(58,1,38,length(XTrainSCAT),'single');
for kk = 1:length(XTrainSCAT)
    scatCONVTrain(:,:,:,kk) = reshape(XTrainSCAT{kk}',[58,1,38]);
end

Construct and train the 1-D convolutional network.

conv1dLayers = [
    imageInputLayer([58 1 38]);
    convolution2dLayer([3 1],64,'stride',1);
    batchNormalizationLayer;
    reluLayer;
    maxPooling2dLayer([2 1]);
    convolution2dLayer([3 1],32);
    reluLayer;
    maxPooling2dLayer([2 1]);
    fullyConnectedLayer(500);
    fullyConnectedLayer(2);
    softmaxLayer;
    classificationLayer;
    ];

convoptions = trainingOptions('sgdm', ...
    'InitialLearnRate',0.01, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropFactor',0.5, ...
    'LearnRateDropPeriod',5, ...
    'Plots','training-progress',...
    'MaxEpochs',50,...
    'Verbose',0,...
    'Plots','none',...
    'MiniBatchSize',15);
[scatCONV1dnet,infoCONV] = ...
    trainNetwork(scatCONVTrain,TrainLabels,conv1dLayers,convoptions);

Plot the smoothed training accuracy and loss per iteration.

iterPerEpoch = floor(size(scatCONVTrain,4)/15);
figure
subplot(2,1,1)
smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoCONV.TrainingAccuracy);
smoothedLoss = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoCONV.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')

Reshape the test set of scattering sequences to be compatible with the network for classification.

scatCONVTest = zeros(58,1,38,length(XTestSCAT));
for kk = 1:length(XTestSCAT)
    scatCONVTest(:,:,:,kk) = reshape(XTestSCAT{kk}',58,1,38);
end

Classification and Analysis

Load the trained gated recurrent unit (GRU) and 1-D convolutional networks along with the test data and scattering sequences. All data, features, and networks were created in the Training — Feature Extraction and Network Training section.

load(fullfile(saveFolder,"transverse_crack_latest","TestData.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","TestLabels.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","XTestSCAT.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","scatGRUnet"))
load(fullfile(saveFolder,"transverse_crack_latest","scatCONV1Dnet.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","scatCONVTest.mat"))

If you additionally want the preprocessed data, labels, and wavelet scattering sequences for the training data, you can load those with the following commands. These data and labels are not used in the remainder of this example if you wish to skip the following load commands.

load(fullfile(saveFolder,"transverse_crack_latest","TrainData.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","TrainLabels.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","XTrainSCAT.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","scatCONVTrain.mat"))

Examine the number of time series per class in the test set. Note the test set is significantly imbalanced as discussed in Data — Preprocessing section.

tabulate(TestLabels)
  Value    Count   Percent
     CR       22     33.33%
   UNCR       44     66.67%

XTestSCAT contains the wavelet scattering sequences computed for the raw time series in TestData.

Show the GRU model performance on the test data not used in model training.

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'});

In spite of the large imbalance in the classes and the small dataset, the precision and recall values indicate the network performs well on the test data. Specifically, the precision value for "Cracked" data is approximately 100% and the recall is approximately 95%. These values are especially good given that a full 67% of the records in the training set were "Uncracked". The network has not overlearned to classify the time series as "Uncracked" in spite of the imbalance.

If you set the inputSize = 1 and transpose the time series in the training data, you can retrain the GRU network on the raw time series data. This was done on the same data in the training set. You can load that network and check the performance on the test set.

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'});

For this network, the performance is not good. Specifically, the recall for the "Cracked" data is poor. The number of false negatives for the "Crack" data is quite large. This is exactly what you would expect with an imbalanced dataset when the model has not learned well.

Test the 1-D convolutional network trained with the wavelet scattering sequences.

miniBatchSize = 15;
YPredSCAT = classify(scatCONV1Dnet,scatCONVTest,...
    'MiniBatchSize',miniBatchSize);
figure
confusionchart(TestLabels,YPredSCAT,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'1-D Convolutional Network-- Wavelet Scattering Sequences'; ...
    'Confusion chart with precision and recall'});

The performance of the convolutional network with scattering sequences is excellent and exceeds the performance of the GRU network. Precision and recall on the minority class demonstrate robust learning.

To train the 1-D convolutional network on the raw sequences set inputSize = [461 1] in the imageInputLayer. You can load and test that network.

load(fullfile(saveFolder,"transverse_crack_latest","tsCONV1Dnet.mat"));
XTestData = cell2mat(TestData');
XTestData = reshape(XTestData,[461 1 1 66]);
miniBatchSize = 15;
YPredRAW = classify(tsCONV1Dnet,XTestData,...
    'MiniBatchSize',miniBatchSize);
confusionchart(TestLabels,YPredRAW,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'1-D Convolutional Network-- Raw Sequences'; ...
    'Confusion chart with precision and recall'});

The 1-D convolution network with the raw sequences performs better than the GRU network, but not as well as the convolutional network trained with the wavelet scattering sequences. Further, the reduction along the time dimension of the data with the scattering transform has resulted in a network which is approximately 8.5 times smaller than the convolutional network trained on the raw data.

Wavelet Inference and Analysis

This section demonstrates how to classify a single time series using wavelet analysis with a pretrained model. The model used is the 1-D convolutional network trained on wavelet scattering sequences. Load the trained network and some test data if you have not already loaded these in the previous section.

load(fullfile(saveFolder,"transverse_crack_latest","scatCONV1Dnet.mat"));
load(fullfile(saveFolder,"transverse_crack_latest","TestData.mat"));

Construct the wavelet scattering network to transform the data. Select a time series from the test data and classify the data. If the model classifies the time series as "Cracked", investigate the series for the position of the crack in the waveform.

sf = waveletScattering('SignalLength',461,...
    'OversamplingFactor',1,'qualityfactors',[8 1],...
    'InvarianceScale',0.05,'Boundary','reflect','SamplingFrequency',1280);
idx = 10;
data = TestData{idx};
[smat,x] = featureVectors(data,sf);
PredictedClass = classify(scatCONV1Dnet,smat);
if isequal(PredictedClass,'CR') 
    fprintf('Crack detected. Computing wavelet transform modulus maxima.\n')
    wtmm(data,'Scaling','local');
end
Crack detected. Computing wavelet transform modulus maxima.

The wavelet transform modulus maxima (WTMM) technique shows a maxima line converging to the finest scale at sample 205. Maxima lines that converge to fine scales are a good estimate of where singularities are in a time series. This makes sample 205 a good estimate of the location of the crack.

plot(x,data)
axis tight
hold on
plot([x(205) x(205)],[min(data) max(data)],'k')
grid on
title(['Crack located at ' num2str(x(205)) 'meters'])
xlabel('Distance (m)')
ylabel('Amplitude')

You can increase your confidence in this location by using multiresolution analysis (MRA) techniques and identifying changes in slope in long-scale wavelet MRA series. See Practical Introduction to Multiresolution Analysis for an introduction to MRA techniques. In [1] the difference in energy between "Cracked" and "Uncracked" series occurred in the low frequency bands, specifically in the interval of [10,20] Hz. Accordingly, the following MRA is focused on signal components in the frequency bands from [10,80] Hz. In these bands, identify linear changes in the data. Plot the change points along with the MRA components.

[mra,chngpts] = helperMRA(data,x);

The MRA-based changepoint analysis has helped to confirm the WTMM analysis in identifying the region around 1.8 meters as the probable location of the crack.

Summary

This example showed how to use wavelet scattering sequences with both recurrent and convolutional networks to classify time series. The example further demonstrated how wavelet techniques can help to localize features on the same spatial (time) scale as the original data.

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.

Appendix

Helper functions used in this example.

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

See Also

Related Topics