利用加速度计数据进行裂纹识别

这个示例使用:

打开实时脚本

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

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

数据-描述和必需的属性

本例中使用的数据是从Mendeley data开放数据存储库中检索的[2]．数据是在知识共享协议(CC) 4.0许可下分发的。在MATLAB®指定的临时目录中下载本例中使用的数据和模型tempdir命令。如果您选择在不同的文件夹中下载数据tempdir，请在后续说明中修改目录名称。将数据解压缩到指定的文件夹中TransverseCrackData．

数据URL=“https://ssd.mathworks.com/金宝appsupportfiles/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 ~ 13 mm，裂缝间距为7 ~ 35 mm。这些路段以30公里/小时、40公里/小时和50公里/小时三种不同的速度行驶。在1.28 kHz的采样频率下，从前乘客悬架转向节获得垂直加速度测量值。30、40和50公里/小时的速度对应于传感器在30公里/小时每6.5毫米、40公里/小时8.68毫米和50公里/小时10.85毫米的测量值。参见[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 =独特(tslen);Ng = histcounts (tslen);Ng = Ng(Ng > 0);栏(uLen, Ng, 0.5)网格在…上甘氨胆酸AX =;斧子。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)))包含(“样本长度”)伊拉贝尔(“系列数量”)标题(的时间序列长度）

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

countlabels (allroadLabel)

ans =2×3表标签计数百分比：CR109 33.333 UNCR 218 66.667

此数据集明显不平衡。没有裂缝的时间序列是正常时间序列的两倍(UNCR)作为包含裂纹(CR)这意味着，在没有任何学习的情况下，预测每个记录上“未开裂”的分类器将达到67%的准确率。

时间序列也有不同的长度。要使用卷积网络，需要一个通用的输入形状。在递归网络中，可以使用不等长的时间序列作为输入，但小批量中的所有时间序列都要根据训练选项进行填充或截断。这需要在为训练和测试创建小批量时小心谨慎，以确保填充序列的正确分布。此外，它要求您在训练期间不要打乱数据。对于这个小数据集，需要对每个epoch的训练数据进行变换。因此，使用普通的时间序列长度。

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

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

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

allroadData=equalents（allroadData）；all（cellfun（@numel，allroadData）==461）

ans =逻辑1

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

splt8020=拆分标签（allroadLabel，0.80）；计数标签（allroadLabel（splt8020{1}））

ans =2×3表标签计数百分比uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu

countlabels (allroadLabel (splt8020 {2}))

ans =2×3表标签计数百分比\uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu

创建培训和测试集。

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

训练前，将数据和标签洗一次。

idx = randperm(长度(TrainData));TrainData = TrainData (idx);TrainLabels = TrainLabels (idx);idx = randperm(长度(TrainData));TrainData = TrainData (idx);TrainLabels = TrainLabels (idx);

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

XTrainSCAT =细胞(大小(TrainData));为kk = 1:numel(TrainData) XTrainSCAT{kk} = helperscat(TrainData{kk});结束npaths = cellfun (@ (x)大小(x, 1), XTrainSCAT);inputSize = npaths (1);

培训- GRU网络

构建GRU网络层。使用两个带有30个隐藏单位的GRU层，以及两个dropout层。输入的大小是散射路径的数目。要在原始时间序列上训练这个网络，改变inputSize将每个时间序列转置为行向量(1 × 461)。如果您希望跳过网络培训，您可以直接到分类和分析部分。在那里，您可以加载经过训练的GRU网络以及预处理的训练和测试集。

numHiddenUnits1 = 30;numHiddenUnits2 = 30;numClasses = 2;GRUlayers = [...sequenceInputLayer (inputSize“名字”，“InputLayer”，...“归一化”，“零中心”) gruLayer (numHiddenUnits1“名字”，“GRU1”，“OutputMode”，“序列”)dropoutLayer（0.35，“名字”，“Dropout1”) gruLayer (numHiddenUnits2“名字”，“GRU2”，“OutputMode”，“最后一次”) dropoutLayer (0.2,“名字”，“Dropout2”) fullyConnectedLayer (numClasses“名字”，“FullyConnected”) softmaxLayer (“名字”，“smax”）;classificationLayer (“名字”，“ClassificationLayer”));

培训GRU网络。使用15个小批量，150个纪元。

maxEpochs = 150;miniBatchSize = 15;选择= trainingOptions (“亚当”，...“ExecutionEnvironment”，“图形”，...“L2Regularization”，1e-3，...“MaxEpochs”，maxEpochs，...“MiniBatchSize”miniBatchSize,...“洗牌”，“每个时代”，...“详细”0,...“阴谋”，“没有”)；iterPerEpoch=楼层（长度（XTrainSCAT）/小批量大小）；[scatGRUnet，infoGRU]=列车网络（XTrainSCAT，列车标签，GRUlayers，选项）；

绘制每个迭代的平滑训练准确性和损失。

figure subplot(2,1,1) smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1)，1，...infoGRU.TrainingAccuracy);smoothedLoss =过滤器(1 / iterPerEpoch *的(iterPerEpoch, 1), 1,...infoGRU.TrainingLoss);情节(smoothedACC)标题(“训练准确性(平滑)”...num2str（伊特佩雷波奇）“每个纪元的迭代”]) ylabel (的精度(%)ylim([0 100.1])网格在…上xlim([1 length(smoothedACC)]) subplot(2,1,2) plot(smoothedLoss) ylim([-0.01 1]) grid在…上xlim（[1长度（平滑度）]ylabel(“损失”)xlabel(“迭代”）

对保留测试数据进行小波散射变换，用于分类。

XTestSCAT =细胞(大小(TestData));为kk = 1:numel(TestData) XTestSCAT{kk} = helperscat(TestData{kk});结束

训练-一维卷积网络

用小波散射序列训练一维卷积网络。散射序列是28 × 58，其中58是时间步长数，28是散射路径数。为了在一维卷积网络中使用这个方法，将散射序列转置并重塑为58 × 1 × 28 × Nsig，其中Nsig是训练示例的数量。如果您希望跳过网络培训，您可以直接到分类和分析在这里，您可以加载经过训练的卷积网络以及预处理的训练集和测试集。

scatCONVTrain = 0(1, 28日,58长度(XTrainSCAT),“单一”）;为kk = 1:length(XTrainSCAT) scatCONVTrain(:，:，:，kk) =重塑(XTrainSCAT{kk}'，[58,1,28]);结束

构建并训练一维卷积网络。

conver1dlayers = [imageInputLayer([58 1 28])); / /转换图层64年convolution2dLayer (1 [3],“步”1);batchNormalizationLayer;reluLayer;2 maxPooling2dLayer ([1]);convolution2dLayer (1 [3], 32);reluLayer;2 maxPooling2dLayer ([1]);fullyConnectedLayer (500);fullyConnectedLayer (2);softmaxLayer; classificationLayer; ]; convoptions = trainingOptions(“个”，...“InitialLearnRate”,0.01,...“LearnRateSchedule”，“分段”，...“LearnRateDropFactor”,0.5,...“LearnRateDropPeriod”,5，...“阴谋”，“没有”，...“MaxEpochs”, 50岁,...“详细”0,...“阴谋”，“没有”，...“MiniBatchSize”15);[scatCONV1Dnet, infoCONV] =...trainNetwork (scatCONVTrain TrainLabels、conv1dLayers convoptions);

绘制每个迭代的平滑训练准确性和损失。

iterPerEpoch =地板(大小(scatCONVTrain 4) / 15);figure subplot(2,1,1) smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1)，1，...信息转换培训（准确性）；SmoothdLoss=过滤器（1/iterPerEpoch*个（iterPerEpoch，1），1，...infoCONV.TrainingLoss);情节(smoothedACC)标题(“训练准确性(平滑)”...num2str（伊特佩雷波奇）“每个纪元的迭代”]) ylabel (的精度(%)ylim([0 100.1])网格在…上xlim([1 length(smoothedACC)]) subplot(2,1,2) plot(smoothedLoss) ylim([-0.01 1]) grid在…上xlim（[1长度（平滑度）]ylabel(“损失”)xlabel(“迭代”）

对散射序列测试集进行重塑，使其与网络兼容，便于分类。

scatCONVTest=零（58,1,28，长度（XTestSCAT））；为kk = 1:length(XTestSCAT) scatCONVTest(:，:，:，kk) = shape(XTestSCAT{kk}'，58,1,28);结束

分类和分析

将训练过的门控循环单元(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”，“scatCONVTest.mat”）)

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

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

检查测试集中每个类的时间序列数。请注意，测试集明显不平衡，如中所述数据预处理部分。

countlabels (TestLabels)

ans =2×3表标签计数百分比\uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu

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

在模型培训中未使用的测试数据上显示GRU模型性能。

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

尽管类的不平衡性大，数据集小，但精度和召回值表明该网络在测试数据上表现良好。具体来说，“破解”数据的精度值接近98%，召回率接近96%。这些值特别好，因为训练集中有整整67%的记录是“未破解”的。尽管不平衡，网络还没有过度学习将时间序列分类为“未破解”。

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

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

对于这个网络，性能不好。具体来说，“破解”数据的召回率很低。“破解”数据的漏报率相当高。这正是当模型学习不好时，不平衡数据集的情况。

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

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

带有散射序列的卷积网络性能优异，超过了GRU网络的性能。在少数族裔的课堂上，精确和记忆证明了稳健的学习。

在原始序列集上训练一维卷积网络inputSize = [461 1 1]在图像输入层．您可以加载和测试该网络。

负载(fullfile (saveFolder“transverse_crack_latest”，“tsCONV1Dnet.mat”))；XTestData=cell2mat（TestData'）；XTestData=restrape（XTestData[461 1 66]）；miniBatchSize=15；YPredRAW=classify（tsCONV1Dnet，XTestData，...“MiniBatchSize”, miniBatchSize);confusionchart (TestLabels YPredRAW,“RowSummary”，“row-normalized”，...“ColumnSummary”，“column-normalized”）;标题({“1-D卷积网络——原始序列”；...“精确和回忆的混乱图表”});

用原始序列训练的一维卷积网络的性能优于GRU网络，但不如用小波散射序列训练的卷积网络。此外，随着散射变换数据的时间维数的降低，产生了一个网络，它大约比在原始数据上训练的卷积网络小8.5倍。

小波推理与分析

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

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

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

科幻小说= waveletScattering (“信号长度”, 461,...“OversamplingFactor”，1，“qualityfactors”8 [1],...“不变性刻度”, 0.05,“边界”，“反映”，“SamplingFrequency”, 1280);idx = 1;data = TestData {idx};[smat x] = featureVectors(数据,科幻小说);PredictedClass =分类(scatCONV1Dnet smat);如果isequal (PredictedClass“CR”)流(裂纹检测。计算小波变换模极大值。\n′)wtmm（数据、，“缩放”，“本地”）;结束

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

小波变换模极大值（WTMM）技术显示，样本205处的最大值线收敛到最细尺度。收敛到精细尺度的最大值线是时间序列中奇异点位置的良好估计。这使样本205成为裂纹位置的良好估计。

情节(x,数据)轴紧持有在…上情节([x (205) (205)], [min(数据)max(数据)],“k”网格)在…上标题([“裂缝位于”num2str（x（205））“米”])包含(“距离(m)”)伊拉贝尔(“振幅”）

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

[mra，chngpts]=helperMRA（数据，x）；

基于MRA的变点分析有助于确认WTMM分析，将1.8米左右的区域确定为裂纹的可能位置。

总结

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

参考文献

杨群，周诗诗。“基于车辆振动信号分析的沥青路面横向裂缝识别”，《道路材料与路面设计》，2020,1-19。https://doi.org/10.1080/14680629.2020.1714699．

[2] 周，石狮。“车辆振动数据。”https://data.mendeley.com/datasets/3dvpjy4m22/1．数据使用于CC 4.0.数据从原始Excel数据格式重新打包到MAT文件。删除速度标签，只保留“裂纹”或“nocrack”标签。

附录

本例中使用的辅助函数。

函数smat=helperscat（数据输入）%这个辅助函数只支持小波工具箱的例子。金宝app%它可能会在未来的版本中改变或被删除。datain =单(datain);sn = waveletScattering (“信号长度”、长度(datain),...“OversamplingFactor”，1，“qualityfactors”8 [1],...“不变性刻度”, 0.05,“边界”，“反映”，“SamplingFrequency”, 1280);smat = sn.featureMatrix (datain);结束%-----------------------------------------------------------------------函数dataUL = equalLenTS(数据)这个函数在只支持小波工具箱的例子。金宝app%它可能会在未来的版本中改变或被删除。N =长度(数据);dataUL =细胞(N, 1);为kk = 1:N L = length(data{kk});转换l情况下461 dataUL{kk} = data{kk};情况下369 Ndiff = 461-369;垫= Ndiff / 2;dataUL {kk} =[翻转(数据{kk}(1:垫);数据{kk};...翻转(数据{kk} (L-pad + 1: L)));否则Ndiff = l - 461;zr = Ndiff / 2;dataUL {kk} = {kk}数据(zr: end-zrs-1);结束结束结束%--------------------------------------------------------------------------函数(fmat x) = featureVectors(数据、科幻小说)%这个函数只支持小波工具箱的例子。金宝app%它可能会在未来的版本中改变或被删除。数据=单个（数据）；N=长度（数据）；dt=1/1280；如果N < 461 Ndiff = 461-N;垫= Ndiff / 2;dataUL =[翻转(数据(1:垫);数据;...翻转(数据(N-pad + 1: N)));率= 5 e4/3600;dx =速度* dt;x = 0: dx: dx (N * dx);elseifN > 461 Ndiff = N-461;zr = Ndiff / 2;dataUL =数据(zr: end-zrs-1);率= 3 e4/3600;dx =速度* dt;x = 0: dx: dx (N * dx);其他的dataUL =数据;率= 4 e4/3600;dx =速度* dt;x = 0: dx: dx (N * dx);结束fmat = sf.featureMatrix (dataUL);Fmat =重塑(Fmat '，[58 1 28]);结束%------------------------------------------------------------------------------函数[mra, chngpts] = helperMRA(数据,x)%这个函数只支持小波工具箱的例子。金宝app%它可能会在未来的版本中改变或被删除。mra = modwtmra (modwt(数据,“sym3”)，“sym3”)；mraLev=mra（4:6，：）；Ns=size（mraLev，1）；thresh=[2,4,8]；chngpts=false（size（mraLev））；%确定changepoints。我们需要不同的阈值%的分辨率水平。为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，：）如果Kk == 1 title(“MRA组件”）结束yyaxis正确的hs=阀杆（x，idx）；hs=“关闭”；海关。标志='^'；海关。markfacecolor = [1 0 0];结束网格在…上轴紧包含(“距离(m)”）结束