从加速度计数据中识别裂纹

这个例子使用了:

打开实时脚本

该实例说明了如何利用小波和深度学习技术检测路面横向裂缝并对其位置进行定位。该实例演示了使用小波散射序列作为门控循环单元(GRU)和一维卷积网络的输入，根据裂缝的存在或不存在对时间序列进行分类。这些数据是由安装在前乘客座椅车轮悬挂关节上的传感器获得的垂直加速度测量数据。横向裂缝的早期识别对路面性能评价和养护具有重要意义。可靠的自动检测方法能够进行更频繁和更广泛的监测。

请细阅数据-描述和必要的属性在运行此示例之前。所有的数据导入和预处理都在数据-导入和数据-预处理部分。如果您想跳过训练-测试集创建、预处理、特征提取和模型训练，您可以直接转到分类与分析部分。在那里，您可以加载预处理的测试数据以及提取的特征和训练的模型。

数据-描述和必要的属性

本例中使用的数据是从Mendeley data开放数据存储库中检索的[2］．数据是在知识共享(CC) 4.0许可下分发的。在MATLAB®指定的临时目录中下载本例中使用的数据和模型tempdir命令。如果选择在不同的文件夹中下载数据tempdir，在后续提示中修改目录名称。将数据解压缩到指定为的文件夹中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);解压缩(zipFile saveFolder)

解压缩数据后，将TransverseCrackData文件夹包含子文件夹称为transverse_crack_latest．所有后续命令都必须在此文件夹中运行，或者可以将此文件夹放在matlabpath．

文本文件，vehiclevibrationdata.rights，包含在zip文件中包含CC BY 4.0许可证的文本。数据已经从原始的Excel格式重新打包到mat文件中。

数据采集的详细描述见[1］．12个4米长的沥青段包含一个中心位置的横向裂缝和12个等长未裂缝段被使用。这些数据来自三条不同的道路。横向裂缝宽度为2- 13mm，裂缝间距为7- 35mm。这些路段以三种不同的速度行驶:30公里/小时、40公里/小时和50公里/小时。在1.28 kHz的采样频率下，获得了来自前乘客悬架关节的垂直加速度测量。30、40和50公里/小时的速度对应的传感器测量每6.5毫米在30公里/小时，8.68毫米在40公里/小时，10.85毫米在50公里/小时。参见[1]，以便对这些数据进行详细的小波分析，这与本例中的分析不同。

为了与CC BY 4.0许可证的规定保持一致，我们注意到在本例中使用的数据中没有保留原始数据的速度信息。速度和道路信息保存在road1.mat，road2.mat,road3.mat为了完整性，数据文件包含在数据文件夹中。

数据-下载和导入

加载加速度计数据及其对应的标签。有327个加速度计记录。

负载(fullfile (saveFolder“transverse_crack_latest”，“allroadData.mat”));负载(fullfile (saveFolder“transverse_crack_latest”，“allroadLabel.mat”));

数据-预处理

得到所有时间序列的长度。显示每个长度的时间序列数量的柱状图。

tslen = cellfun(@length,allroadData);uLen = unique(tslen);Ng = histcounts(tslen);Ng = Ng(Ng > 0);栏(uLen, Ng, 0.5)网格在AX = gca;斧子。YLim = [0 max(Ng)+15];文本(uLen (1), Ng (1) + 10, num2str (Ng(1)))文本(uLen (2), Ng (2) + 10, num2str (Ng(2)))文本(uLen (3), Ng (3) + 10, num2str (Ng(3)))包含(“样本长度”) ylabel (“系列数”)标题(“时间序列长度”）

图中包含一个axes对象。标题为Time Series Length的axes对象包含4个类型为bar、text的对象。

数据集中有三个惟一长度:369、461和615样本。车辆以三种不同的速度行驶，但经过的距离和采样率是恒定的，导致不同的数据记录长度。确定有多少记录在“已破解”(CR)及“未破解”(UNCR)类。

countlabels (allroadLabel)

ans =2×3表标签计数百分比_____ _____ _______ CR 109 33.333 UNCR 218 66.667

这个数据集明显不平衡。有两倍的时间序列没有裂缝(UNCR)作为包含裂纹(CR)。这意味着在每个记录上预测“未破解”的分类器不需要任何学习就能达到67%的准确率。

时间序列也有不同的长度。要使用小波散射变换，需要一个公共的输入长度。在循环网络中，可以使用不等长度的时间序列作为输入，但是小批中的所有时间序列都根据训练选项进行填充或截断。这需要在为训练和测试创建小批次时谨慎，以确保填充序列的适当分布。此外，它要求您在训练期间不要打乱数据。对于这个小数据集，需要对每个时代的训练数据进行洗牌。因此，使用一个公共的时间序列长度。

最常见的长度是461个样本。此外，裂纹，如果存在，是在中央位置的记录。因此，通过反映初始和最终的46个样本，我们可以将369个样本的序列对称扩展到461的长度。在615个样本的记录中，去掉最初的77个和最后的77个样本。

训练-特征提取和网络训练

下面的部分生成训练集和测试集，创建小波散射序列，训练门控循环单元(GRU)和一维卷积网络。首先，扩展或截断两个集合的时间序列，得到461个样本的共同长度。

allroadData = equalLenTS(allroadData);(cellfun (@numel allroadData) = = 461)

ans =逻辑1

现在每个时间序列在破解和未破解的数据集中都有461个样本。将每个类中80%的时间序列组成的训练集中的数据拆分，并保留每个类中剩下的20%进行测试。验证不平衡的比例保留在每一组。

splt8020 = splitlabels(allroadLabel,0.80);countlabels (allroadLabel (splt8020 {1}))

ans =2×3表标签计数百分比_____ _____ _______ CR 87 33.333 UNCR 174 66.667

countlabels (allroadLabel (splt8020 {2}))

ans =2×3表标签计数百分比_____ _____ _______ CR 22 33.333 UNCR 44 66.667

创建培训和测试集。

TrainData = allroadData(splt8020{1});TrainLabels = allroadLabel(splt8020{1});TestData = allroadData(splt8020{2});TestLabels = allroadLabel(splt8020{2});

训练前对数据和标签进行一次洗牌。

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

计算每个训练序列的散射序列。散射序列存储在单元阵列中，以兼容GRU和一维卷积网络。

XTrainSCAT = cell(size(TrainData));为kk = 1: numl (TrainData) XTrainSCAT{kk} = helperscat(TrainData{kk});结束= cellfun(@(x)size(x,1)，XTrainSCAT);inputSize = npaths(1);

培训- GRU网络

构建GRU网络层。使用两个GRU层，每个有30个隐藏单位，以及两个dropout层。因为类是明显不平衡的，使用加权分类层，类权重与逆类频率成正比。输入的大小是散射路径的数量。若要在原始时间序列上训练此网络，请更改inputSize到1，并将每个时间序列转置为行向量(1 × 461)。如果您想跳过网络培训，您可以直接访问分类与分析部分。在那里，您可以加载经过训练的GRU网络以及预处理的训练和测试集。

numhiddenunit1 = 30;numhiddenunit2 = 30;numClasses = 2;classfrequency = countcats(allroadLabel);Nsamp = sum(classfrequency);weightCR = 1/ classfrequency (1)*Nsamp/2;weightUNCR = 1/ classfrequency (2)*Nsamp/2;GRUlayers = [.．.sequenceInputLayer (inputSize“名字”，“InputLayer”，.．.“归一化”，“zerocenter”) gruLayer (numHiddenUnits1“名字”，“GRU1”，“OutputMode”，“序列”) dropoutLayer (0.35,“名字”，“Dropout1”) gruLayer (numHiddenUnits2“名字”，“GRU2”，“OutputMode”，“最后一次”) dropoutLayer (0.2,“名字”，“Dropout2”) fullyConnectedLayer (numClasses“名字”，“FullyConnected”) softmaxLayer (“名字”，“smax”）;classificationLayer (“名字”，“ClassificationLayer”，“类”，[“CR”“UNCR”),.．.“ClassWeights”, (weightCR weightUNCR]);];

训练GRU网络。使用15个小批次，150个周期。

maxEpochs = 150;miniBatchSize = 15;options = trainingOptions(“亚当”，.．.“ExecutionEnvironment”，“图形”，.．.“L2Regularization”1 e - 3,.．.“MaxEpochs”maxEpochs,.．.“MiniBatchSize”miniBatchSize,.．.“洗牌”，“every-epoch”，.．.“详细”0,.．.“阴谋”，“没有”）;iterPerEpoch = floor(length(XTrainSCAT)/miniBatchSize);[scatGRUnet,infoGRU] = trainNetwork(XTrainSCAT,TrainLabels,GRUlayers,options);

绘制平滑训练精度和每次迭代的损失。

图subplot(2,1,1) smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1)，1，.．.infoGRU.TrainingAccuracy);smooththeloss = filter(1/iterPerEpoch*ones(iterPerEpoch,1)，1，.．.infoGRU.TrainingLoss);情节(smoothedACC)标题(“训练准确度(平滑)”.．.num2str (iterPerEpoch)“每个纪元的迭代”]) ylabel (的精度(%)ylim([0 100.1])网格在xlim([1 length(smoothedACC)]) subplot(2,1,2) plot(smoothedLoss) ylim([-0.01 1])网格在xlim([1 length(smoothedLoss)])“损失”)包含(“迭代”）

图中包含2个轴对象。axis对象1具有标题Training Accuracy (smooened)，每个纪元17次迭代，包含一个类型为line的对象。Axes对象2包含一个类型为line的对象。

得到搁置试验数据的小波散射变换，用于分类。

XTestSCAT = cell(size(TestData));为kk = 1:numel(TestData) XTestSCAT{kk} = helperscat(TestData{kk});结束

训练-一维卷积网络

训练小波散射序列的一维卷积网络。如果您想跳过网络培训，您可以直接访问分类与分析部分。在那里可以加载经过训练的卷积网络以及预处理的训练集和测试集。

构造和训练一维卷积网络。散射网络中有28条路径。

conv1dLayers = [sequenceInputLayer(28，的最小长度58岁的“归一化”，“zerocenter”）;24岁的convolution1dLayer (3“步”2);batchNormalizationLayer;reluLayer;maxPooling1dLayer (4);convolution1dLayer(3, 16岁,“填充”，“相同”）;batchNormalizationLayer;reluLayer;maxPooling1dLayer (2);fullyConnectedLayer (150);fullyConnectedLayer (2);globalAveragePooling1dLayer;softmaxLayer;classificationLayer (“名字”，“ClassificationLayer”，“类”，[“CR”，“UNCR”),.．.“ClassWeights”, (weightCR weightUNCR]);];convoptions = trainingOptions(“亚当”，.．.“InitialLearnRate”, 0.01,.．.“ExecutionEnvironment”，“cpu”，.．.“LearnRateSchedule”，“分段”，.．.“LearnRateDropFactor”, 0.5,.．.“LearnRateDropPeriod”5,.．.“阴谋”，“没有”，.．.“MaxEpochs”, 50岁,.．.“详细”0,.．.“阴谋”，“没有”，.．.“MiniBatchSize”, miniBatchSize);[scatCONV1Dnet, infoCONV] =.．.trainNetwork (XTrainSCAT TrainLabels、conv1dLayers convoptions);

绘制平滑训练精度和每次迭代的损失。

iterPerEpoch = floor(length(XTrainSCAT)/miniBatchSize);图subplot(2,1,1) smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1)，1，.．.infoCONV.TrainingAccuracy);smooththeloss = filter(1/iterPerEpoch*ones(iterPerEpoch,1)，1，.．.infoCONV.TrainingLoss);情节(smoothedACC)标题(“训练准确度(平滑)”.．.num2str (iterPerEpoch)“每个纪元的迭代”]) ylabel (的精度(%)ylim([0 100.1])网格在xlim([1 length(smoothedACC)]) subplot(2,1,2) plot(smoothedLoss) ylim([-0.01 1])网格在xlim([1 length(smoothedLoss)])“损失”)包含(“迭代”）

图中包含2个轴对象。axis对象1具有标题Training Accuracy (smooened)，每个纪元17次迭代，包含一个类型为line的对象。Axes对象2包含一个类型为line的对象。

分类与分析

加载训练过的门控循环单元(GRU)和一维卷积网络以及测试数据和散射序列。所有的数据、特性和网络都是在训练-特征提取和网络训练部分。

负载(fullfile (saveFolder“transverse_crack_latest”，“TestData.mat”)加载(fullfile (saveFolder“transverse_crack_latest”，“TestLabels.mat”)加载(fullfile (saveFolder“transverse_crack_latest”，“XTestSCAT.mat”)加载(fullfile (saveFolder“transverse_crack_latest”，“scatGRUnet”)加载(fullfile (saveFolder“transverse_crack_latest”，“scatCONV1Dnet.mat”）)

如果您还需要训练数据的预处理数据、标签和小波散射序列，您可以使用以下命令加载它们。如果您希望跳过以下内容，则在本示例的其余部分中不会使用这些数据和标签负载命令。

负载(fullfile (saveFolder“transverse_crack_latest”，“TrainData.mat”)加载(fullfile (saveFolder“transverse_crack_latest”，“TrainLabels.mat”)加载(fullfile (saveFolder“transverse_crack_latest”，“XTrainSCAT.mat”）)

检查测试集中每个类的时间序列的数量。注意，测试集是显著不平衡的数据-预处理部分。

countlabels (TestLabels)

ans =2×3表标签计数百分比_____ _____ _______ CR 22 33.333 UNCR 44 66.667

XTestSCAT包含为原始时间序列计算的小波散射序列TestData．

在未用于模型训练的测试数据上显示GRU模型的性能。

miniBatchSize = 15;ypredSCAT = classification (scatGRUnet,XTestSCAT，.．.“MiniBatchSize”, miniBatchSize);图confusionchart (TestLabels ypredSCAT,“RowSummary”，“row-normalized”，.．.“ColumnSummary”，“column-normalized”）;标题({“GRU网络—小波散射序列”；.．.“精确和回忆的混淆图表”})；

MATLAB图

尽管类之间存在较大的不平衡，数据集较小，但精度和召回值表明该网络在测试数据上表现良好。具体来说，“破解”数据的精度和召回值非常好。尽管训练集中67%的记录是“未破解的”，但仍然实现了这一目标。尽管存在不平衡，但电视台并没有过度学习将时间序列归类为“未破解”。

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

负载(fullfile (saveFolder“transverse_crack_latest”，“tsGRUnet.mat”));rawTest = cellfun(@转置，TestData，“UniformOutput”、假);miniBatchSize = 15;YPredraw = category (tsGRUnet,rawTest，.．.“MiniBatchSize”, miniBatchSize);confusionchart (TestLabels YPredraw,“RowSummary”，“row-normalized”，.．.“ColumnSummary”，“column-normalized”）;标题({“GRU网络—原始时间序列”；.．.“精确和回忆的混淆图表”})；

MATLAB图

对于这种网络，性能不是很好。具体来说，“破解”数据的召回率很低。“破解”数据的假阴性数量相当大。当模型没有很好地学习时，这正是不平衡数据集所期望的结果。

测试用小波散射序列训练的一维卷积网络。

miniBatchSize = 15;YPredSCAT = classification (scatCONV1Dnet,XTestSCAT，.．.“MiniBatchSize”, miniBatchSize);图confusionchart (TestLabels YPredSCAT,“RowSummary”，“row-normalized”，.．.“ColumnSummary”，“column-normalized”）;标题({“一维卷积网络——小波散射序列”；.．.“精确和回忆的混淆图表”})；

MATLAB图

具有散射序列的卷积网络性能优异，与GRU网络性能一致。在少数民族的课堂上，精确和回忆证明了稳健的学习。

训练原始序列集上的一维卷积网络inputSize到1sequenceInpuLayer．设置的最小长度到461年。您可以使用相同的数据和相同的网络架构加载和测试训练过的网络。

负载(fullfile (saveFolder“transverse_crack_latest”，“tsCONV1Dnet.mat”));miniBatchSize = 15;TestData = cellfun(@转置，TestData，“UniformOutput”、假);YPredRAW = class (tsCONV1Dnet,TestDataT，.．.“MiniBatchSize”, miniBatchSize);confusionchart (TestLabels YPredRAW,“RowSummary”，“row-normalized”，.．.“ColumnSummary”，“column-normalized”）;标题({“一维卷积网络——原始序列”；.．.“精确和回忆的混淆图表”})；

MATLAB图

用原始序列训练的卷积网络性能较好，但不如用小波散射序列训练的卷积网络。

小波推理与分析

本节演示如何使用小波分析和预训练的模型对单个时间序列进行分类。所使用的模型是小波散射序列训练的一维卷积网络。加载经过训练的网络和一些测试数据(如果您还没有在上一节中加载这些数据)。

负载(fullfile (saveFolder“transverse_crack_latest”，“scatCONV1Dnet.mat”));负载(fullfile (saveFolder“transverse_crack_latest”，“TestData.mat”));

构造小波散射网络对数据进行变换。从测试数据中选择一个时间序列并对数据进行分类。如果模型将时间序列分类为“裂缝”，研究该序列在波形中的裂缝位置。

sf =小波散射(“SignalLength”, 461,.．.“OversamplingFactor”, 1“qualityfactors”8 [1],.．.“InvarianceScale”, 0.05,“边界”，“反映”，“SamplingFrequency”, 1280);Idx = 22;data = TestData{idx};[smat,x] = featureVectors(data,sf);PredictedClass = classid (scatCONV1Dnet,smat);如果isequal (PredictedClass“CR”)流(裂纹检测。计算小波变换模极大值。\n') wtmm(数据,“缩放”，“本地”）;结束

裂纹检测。计算小波变换模极大值。

小波变换模极大值(WTMM)技术显示极大值线在样本225处收敛到最精细的尺度。收敛到细尺度的极大值线可以很好地估计时间序列中奇点的位置。这使样品225可以很好地估计裂缝的位置。

情节(x,数据)轴紧持有在图([x(225) x(225)]，[min(data) max(data)]，“k”网格)在标题([“裂缝位于”num2str (x (225))“米”])包含(“距离(m)”) ylabel (“振幅”）

图中包含一个axes对象。标题为Crack的轴对象位于1.9444米，包含2个类型为线的对象。

通过使用多分辨率分析(MRA)技术和在长尺度小波MRA序列中识别坡度变化，可以增加对该位置的信心。看到多分辨率分析实用导论参阅MRA技术的介绍。在[1“破裂”系列和“未破裂”系列之间的能量差异出现在低频波段，特别是在[10,20]Hz区间。因此，以下MRA主要关注[10,80]Hz频段内的信号成分。在这些波段中，识别数据的线性变化。将变更点与MRA组件一起绘制出来。

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

图中包含3个轴对象。标题为MRA Components的Axes对象1包含一个类型为stem的对象。Axes对象2包含一个stem类型的对象。Axes对象3包含一个stem类型的对象。

基于磁共振成像的变化点分析有助于确认WTMM分析，确定1.94米左右的区域为可能的裂缝位置。

总结

这个例子展示了如何使用小波散射序列与循环网络和卷积网络来分类时间序列。该示例进一步说明了小波技术如何在与原始数据相同的空间(时间)尺度上帮助定位特征。

参考文献

[1]杨群和周诗诗。“基于车辆振动信号分析的沥青路面横向裂缝识别”，《路面工程学报》，2020,1-19。https://doi.org/10.1080/14680629.2020.1714699．

[2]周,石狮。“车辆振动数据。”https://data.mendeley.com/datasets/3dvpjy4m22/1．资料使用于Cc / 4.0．数据从原始的Excel数据格式重新打包到mat文件中。速度标签被移除，只保留“裂缝”或“不裂缝”标签。

附录

本例中使用的Helper函数。

函数Smat = helperscat(datain)此辅助函数仅支持小波工具箱示例。它可能会在未来的版本中更改或删除。。Datain = single(Datain);sn =小波散射(“SignalLength”、长度(datain),.．.“OversamplingFactor”, 1“qualityfactors”8 [1],.．.“InvarianceScale”, 0.05,“边界”，“反映”，“SamplingFrequency”, 1280);smat = sn.featureMatrix(datain);结束%-----------------------------------------------------------------------函数dataUL = equalLenTS(数据)此函数仅在支持小波工具箱的示例中。它可能会在未来的版本中更改或删除。。N =长度(数据);dataUL = cell(N,1);为kk = 1:N L =长度(data{kk});开关l情况下461 dataUL{kk} = data{kk};情况下369 Ndiff = 461-369;pad = Ndiff/2;dataUL{kk} = [flip(data{kk}(1:pad));数据{kk};.．.翻转(数据{kk} (L-pad + 1: L)));否则Ndiff = L-461;zrs = Ndiff/2;dataUL{kk} = data{kk}(zrs:end-zrs-1);结束结束结束%--------------------------------------------------------------------------函数[fmat,x] = featureVectors(data,sf)此函数仅支持小波工具箱示例。它可能会在未来的版本中更改或删除。。数据=单(数据);N =长度(数据);Dt = 1/1280;如果N < 461 Ndiff = 461-N;pad = Ndiff/2;dataUL = [flip(data(1:pad));数据;.．.翻转(数据(N-pad + 1: N)));速率= 5e4/3600;Dx =速率*dt;x = 0:dx:(N*dx)-dx;elseifN > 461 Ndiff = N-461;zrs = Ndiff/2;dataUL = data(zrs:end-zrs-1);速率= 3e4/3600;Dx =速率*dt;x = 0:dx:(N*dx)-dx;其他的dataUL =数据;速率= 4e4/3600;Dx =速率*dt;x = 0:dx:(N*dx)-dx;结束fmat = sf.featureMatrix(dataUL);结束%------------------------------------------------------------------------------函数[mra,chngpts] = helperMRA(data,x)此函数仅支持小波工具箱示例。它可能会在未来的版本中更改或删除。。Mra = modwtmra(modwt(data，“sym3”)，“sym3”）;mraLev = mra(4:6，:);Ns = size(mraLev,1);Thresh = [2,4,8];chngpts = false(size(mraLev));确定变更点。我们想要不同的阈值%分辨率级别。为ii = 1:Ns chngpts(ii，:) = ischange(mraLev(ii，:)，“线性”2,“阈值”打(ii));结束为kk = 1:Ns idx = double(chngpts(kk，:));idx(idx == 0) = NaN;次要情节(Ns, 1 kk)情节(x, mraLev (kk,:))如果k == 1 title(“MRA组件”）结束yyaxis正确的Hs =茎(x,idx);海关。ShowBaseLine =“关闭”；海关。标志=“^”；海关。MarkerFaceColor = [1 0 0];结束网格在轴紧包含(“距离(m)”）结束