本文基于matlab2020版官方网页DocumentationCrack Identification From Accelerometer Data及个人理解。
该示例显示了如何使用小波wavelet和深度学习技术来检测横向路面裂缝并确定其位置。该示例演示了将小波散射序列用作门控循环单元(GRU)和一维卷积网络的输入,以便根据是否存在裂缝对时间序列进行分类。数据是从安装在前排乘客座椅车轮的转向节上的传感器获得的垂直加速度测量值。及早发现发展中的横向裂缝对于评估和维护路面性能非常重要。可靠的自动检测方法可实现更频繁,更广泛的监视。
本文主要介绍其中的第四部分Wavelet Inference and Analysis,即小波处理的相关分析。
小波散射数据分析waveletScattering
本节演示如何使用预先训练好的小波分析对单一时间序列进行分类。所使用的模型是基于小波散射序列训练的一维卷积网络。加载训练过的网络和一些测试数据(如果在前一节中没有加载这些数据)。
使用经过训练的一维CNN模型对一个时间序列进行分类,首先加载数据。
load(fullfile(saveFolder,"transverse_crack_latest","scatCONV1Dnet.mat"));
load(fullfile(saveFolder,"transverse_crack_latest","TestData.mat"));
构造小波散射网络对数据进行变换。从测试数据中选择一个时间序列并对数据进行分类。如果模型将该序列归类为“CR”,则研究波形中裂纹的位置。
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
运行结果
小波变换模极大值技术在样本205处有一条收敛到最细尺度的极大值线。收敛到细尺度的极大值线是对时间序列中奇点位置的一个很好的估计。这使得样本205很好地估计了裂纹的位置。
绘制样本的波形图,并标定205处的点。
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')
运行结果
可以通过使用多分辨率分析(MRA)技术并确定大尺度小波MRA系列中斜率的变化来增强对该位置的信心。有关MRA技术的介绍,请参见Practical Introduction to Multiresolution Analysis。 在[1]中,“CR”和“UNCR”系列之间的能量差发生在低频频段,特别是在[10,20] Hz的间隔内。 因此,以下的MRA集中在[10,80] Hz频带内的信号分量上。在这些波段中,确定数据中的线性变化。绘制变更点以及MRA组件。
[mra,chngpts] = helperMRA(data,x);
基于MRA的变化点分析帮助确认了WTMM的分析,确定了1.8米左右的区域可能是裂纹的位置。
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