通过桥梁振动信号自动识别车辆(MATLAB)
只是简单参数建模,还没有实际场景应用。
Generation of the bridge response to multiple vehicles Initialisation
clearvars;close all;clc clf;close all; Nyy = 446; % Number of nodes to discretize the bridge structure. We need a spatial resolution of 1 m in our case. Nmodes = 6; % number of modes for each degree of freedoms, i.e. a total of 3x6 =18 nodes [Bridge] = LysefjordBridge(Nyy); % bridge structural properties [wn,phi] = eigenBridge(Bridge,Nmodes); % get the eigen frequencies and mode shapes [2,3] fn=wn/(2*pi); % get the frequencies in Hz zetaStruct = 5e-3*ones(3,Nmodes); % modal damping ratios % Define the structure bridge Bridge.zetaStruct = zetaStruct; Bridge.wn = 2*pi*fn; Bridge.phi = phi;
Sampling parameters
fs = 15; % sampling frequency M = 15; % M must be a power of 2 [t,f] = getSamplingPara(M,fs); N=numel(t);
Wind turbulence
% No wind in this example tmax = t(end); Wind.u = zeros(Nyy,N); Wind.w = zeros(Nyy,N); Wind.t = t;
Computation of the vehicle-induced response at midspan (target data)
yTarget = Bridge.L/2; % midspan [~,indY]=min(abs(Bridge.x.*Bridge.L-yTarget)); % get indice where the midspan is located Nvehicle = 4; % Number of vehicle selected tic %Create the structure variable vehicle_Target that contains the aprameters %of each vehicle rng(1) % To ensure reproducibility of this example clear vehicleTarget for ii=1:Nvehicle vehicleTarget(ii).mass = randi([1e3 50e3],1); % mass of the vehicle (kg) vehicleTarget(ii).speed = randi([30 80],1)/3.6; % speed of the vehicle (m/s) vehicleTarget(ii).direction = 1; % All the vehicle have the same direction. If "-1", the vehicle is moving in the opposite direction vehicleTarget(ii).tStart = 55+ii.*400; % arrival time (sec). t = 0s is the initial time of the simulation. end % Compute the wind and vehicle-induced bridge response % Newmark-beta algorithm is used for the solver [Do1,Ao1] = dynaResp_vehicle_TD(Bridge,Wind,'vehicle',vehicleTarget); % Get the vertical bridge displacement response at midspan rz = squeeze(Do1(2,indY,:)); % Some white noise is added to challenge the identification algorithm rz = rz + 0.03*std(rz).*randn(size(rz)); toc
Visualization of the vehicle-induced response at midspan
gcf;close all; figure plot(t,rz,'r') xlabel('Time (s)') ylabel('r_z (m)') set(gcf,'color','w') axis tight
Outlier-detection algorithm and clustering algorithm
tSim = t; minNp = 30; % minimum number of cluster points CO = 30*fs; % cutoff value for the clusters, constructed using a hierarchical cluster tree (cf Matlab function "cluster") fMin = 0.04; % cut-off frequency (high pass filter) -> must be above the lower frequency correctly measured by the accelerometer fMax = 0.15; % cut-off frequency (low pass filter) -> must be lower than the first eigen frequency. here the lowest vertical eigenfrequency is 0.21 Hz [rz_filtered,t_filtered,tImpact_guess,tCluster,rzCluster] = findVehicleID(rz,t,... 'minNp',minNp,'CuttOff',CO,'deltaTmax',inf,'fcLow',fMin,'fcUp',fMax); % tImpact_guess is the Initial guess estimate based on the cluster analysis % rz_filtered,t_filtered, tCluster and rzCluster are used only for visualization purpose! Ncluster = numel(tCluster); clf;close all; figure COLOR = {'r','b','g','m','c','y'}; if Ncluster>numel(COLOR), COLOR = repmat(COLOR,1,Ncluster); end for ii=1:Ncluster subplot(Ncluster,1,ii) plot(t_filtered,rz_filtered,'k') hold on box on plot(tCluster{ii},rzCluster{ii},'o','color',COLOR{ii},'markersize',2) xlabel('time (s)') ylabel('r_z (m)') axis tight end set(gcf,'color','w')
Identification of the vehicles' arrival time and speed
gradS = 0.3; % termination tolerance for the speed estimate, i.e. the speed is identified at +/- gradS (in m/s) gradT = 0.5; % termination tolerance for the arrival time estimate, i.e. it is identified at +/- gradT (in s) newN = 150; % number of data point for the time-domain simulation: less points means faster algorithm but also less accurate one posAcc = yTarget./Bridge.L; % relative position of the accelerometer tic [vehicleSpeed,tImpact] = findSpeed(Bridge,t,rz,posAcc,tImpact_guess,... 'gradS',gradS,'gradS',gradT,'newN',newN,'fcLow',fMin,'fcUp',fMax); toc fprintf('Target speed (km/h): \n') fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[vehicleTarget(:).speed]*3.6) fprintf('identified speed (km/h): \n') fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[vehicleSpeed]*3.6) fprintf('Target arrival time (s): \n') fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[vehicleTarget(:).tStart]) fprintf('identified arrival time (s): \n') fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[tImpact])
Identification of the vehicles' mass
[vehicleMass] = findMass(Bridge,t,rz,posAcc,tImpact,vehicleSpeed,'fcLow',fMin,'fcUp',fMax); fprintf('Target mass (kg): \n') fprintf([repmat('%4.0f \t',1,Nvehicle),' \n'],[vehicleTarget(:).mass]) fprintf('identified mass (kg): \n') fprintf([repmat('%4.0f \t',1,Nvehicle),' \n'],[vehicleMass])
Comparison between the target and simulated bridge response histories
% Create a structure vehicle with the identified parameters of the vehicles. clear vehicle Nvehicle = numel(vehicleSpeed); for ii=1:Nvehicle vehicle(ii).mass = abs(vehicleMass(ii)); vehicle(ii).speed = vehicleSpeed(ii); vehicle(ii).direction = 1; vehicle(ii).tStart = tImpact(ii); vehicle(ii).wn = 0; end % COmpute the relative erros on the identified mas, speed and arrival time errSpeed = nanmean(100.*([vehicle(:).speed]./[vehicleTarget(:).speed]-1)); errMass = nanmean(100.*([vehicle(:).mass]./[vehicleTarget(:).mass]-1)); errTime = nanmean(100.*([vehicle(:).tStart]./[vehicleTarget(:).tStart]-1)); fprintf(['Relative average error on speed: ',num2str(errSpeed,2),' percent \n']) fprintf(['Relative average error on speed: ',num2str(errMass,2),' percent \n']) fprintf(['Relative average error on speed: ',num2str(errTime,2),' percent \n']) % Compute the bridge response with the identified vehicle parameters [Do2] = dynaResp_vehicle_TD(Bridge,Wind,'vehicle',vehicle); % both wind and traffic rzFit = squeeze(Do2(2,indY,:)); clf;close all; figure plot(t,rz,'k',t,rzFit,'r'); legend('target','fitted') axis tight xlabel('Time (s)') ylabel('r_z (m)') set(gcf,'color','w') axis tight % PSD calculation [S1,f1] = pmtm(rz,7,numel(rz),fs); [S2,f2] = pmtm(rzFit,7,numel(rz),fs); figure loglog(f1,S1,'k',f2,S2,'r'); legend('target','fitted') axis tight ylim([1e-15,1]) xlabel('Frequencies (Hz)') ylabel('PSD vertical displacement response (m^2/Hz)') set(gcf,'color','w')
知乎学术咨询: https://www.zhihu.com/consult/people/792359672131756032?isMe=1
工学博士,担任《Mechanical System and Signal Processing》《中国电机工程学报》等期刊审稿专家,擅长领域:现代信号处理,机器学习,深度学习,数字孪生,时间序列分析,设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。
文章版权声明:除非注明,否则均为主机测评原创文章,转载或复制请以超链接形式并注明出处。