# A Crack Identification Method for Bridge Type Structures under Vehicular Load Using Wavelet Transform and Particle Swarm Optimization

Authors: Hakan Gökdağ
Journal: Advances in Acoustics and Vibration
 (2013)

---

## Abstract

In this work a crack identification method is proposed for bridge type structures carrying moving vehicle. The bridge is modeled as an Euler-Bernoulli beam, and open cracks exist on several points of the beam. Half-car model is adopted for the vehicle. Coupled equations of the beam-vehicle system are solved using Newmark-Beta method, and the dynamic responses of the beam are obtained. Using these and the reference displacements, an objective function is derived. Crack locations and depths are determined by solving the optimization problem. To this end, a robust evolutionary algorithm, that is, the particle swarm optimization (PSO), is employed. To enhance the performance of the method, the measured displacements are denoised using multiresolution property of the discrete wavelet transform (DWT). It is observed that by the proposed method it is possible to determine small cracks with depth ratio 0.1 in spite of 5% noise interference.


---

## Body

## 1. Introduction

Structures under vehicular loads have significant applications such as railway tracks, bridges, and roadways. Moving load yields larger deflections and higher stresses than equivalent static loading; thus, dynamics of such structures has received significant attention [1, 2]. If the carrying structure includes crack-like defects, then the impact of moving load becomes more pronounced [3–6]. Various damage detection methods have been developed for such structures using the continuous wavelet transform (CWT) [7–11]. They are based on the fact that CWT coefficients of beam dynamic response demonstrate local peaks at crack locations, and magnitudes of these peaks are proportional to crack depths. The significant deficiency of these methods is that CWT coefficients become insensitive to crack at higher moving load speeds. The reason is that time-dependent structural response becomes smoother at higher vehicle speeds. In this case crack-induced singularities disappear, and CWT coefficients of the response signal fail to extract damage info.Some of the methods in the area of structural damage detection are based on model updating [12]. The basic idea is to update mathematical or finite element model of the structure to match the calculated response to that measured from damaged structure, so that damage parameters are estimated. This is achieved solving an optimization problem. To this end metaheuristic methods are generally preferred, as they do not require gradient info and have less possibility of being trapped by local minima in comparison with the gradient-based methods. The particle swarm optimization (PSO) is one of these algorithms. PSO is a stochastic optimization technique inspired by natural flocking and swarm behavior of birds and insects [13]. It is known to have less parameters and rapid convergence compared with the genetic algorithms (GA) [14] and has been successfully employed in model-updating-based damage detection applications [15–18].Bilello and Bergman [5] state that changes in the time response of the carrying beam due to damage are more perceptible in comparison to the changes in the natural frequencies. Motivated by this conclusion and the approach of Buezas et al. [19], Gökdağ [20] introduced a model-update-based damage detection approach for beam type structures subject to vehicular load. In this respect, time-dependent displacements from several points on a cracked beam are obtained, and an objective function is defined by subtracting these from the ones calculated by the mathematical model of the structure. Then, PSO algorithm is employed to minimize this objective function, so that crack parameters are determined. In the present work, the performance of this approach is improved employing the multiresolution property of the discrete wavelet transform (DWT). Also, it was demonstrated that the proposed method is efficient for small cracks and better than CWT coefficients method in [7] at higher vehicle speeds.


## 2. Material and Method

### 2.1. Dynamic Response of the Vehicle-Beam System

Figure1 illustrates the bridge-vehicle system. Euler-Bernoulli beam model is considered for the bridge, and half-car model is adopted for the vehicle moving with speed V. An open crack with depth h1 is located at z1 on the beam. Surface unevenness of the beam is disregarded, and tires are assumed to be always in contact with the beam [9]. Linear deflections of the beam are considered. Longitudinal deformations and bending-longitudinal coupling are neglected. Under these assumptions the equations of motion for the vehicle-beam system can be derived as follows:
(1)[Iv0000mv0000m10000m2]{d¨1d¨2d¨3d¨4}+[c1L12+c2L22L2c2-L1c1L1c1-L2c2L2c2-L1c1c2+c1-c1-c2L1c1-c1c3+c10-L2c2-c20c2+c4]{d˙1d˙2d˙3d˙4}+[k1L12+k2L22L2k2-L1k1L1k1-L2k2L2k2-L1k1k2+k1-k1-k2L1k1-k1k3+k10-L2k2-k20k2+k4]{d1d2d3d4}={00k3δIu+c3δIu˙k4δIIu+c4δIIu˙},(2)ρA∂2y(z,t)∂t2+C∂y(z,t)∂t+EI∂4y(z,t)∂x4=-PI-PII,
where Iv,mv,m1,m2,L1,L2,k1,k2,k3,c1,c2,c3,c4 are the vehicle parameters shown in Figure 1 and C is the damping of beam material. di, i=1:4, denote vehicle degrees of freedom, and δI and δII are the Dirac delta functions defined as follows:
(3)δII=δ(z-Vt),δI=δ(z-Vt-(L1+L2)),
where PI and PII are the interaction forces acting on the beam through the contact points I and II as follows:(4a)PI=(m1+L2L1+L2mv)9.81+k3(δIu-d3)+c3(δIu˙-d˙3),(4b)PII=(m2+L1L1+L2mv)9.81+k4(δIIu-d4)+c4(δIIu˙-d˙4),where u is the vertical displacement at the tire contact point, that is, u=y(z,t), and its derivative with respect to time is u˙=dy(z,t)/dt=V∂y(z,t)/∂z+∂y(z,t)/∂t and ′=d/dz. At the crack location vibration modes should satisfy the following compatibility equations [7, 11]:(5a)Y1(z1)=Y2(z1),Y1′(z1)+EIkY1′′(z1)=Y2′(z1),(5b)Y1′′(z1)=Y2′′(z1),Y1′′′(z1)=Y2′′′(z1),where
(6)Yi(z)=Ci1cos(λz)+Ci2sin(λz)+Ci3cosh(λz)+Ci4sinh(λz),i=1,2,λ=(ρAω2(EI))0.25,
where λ and ω in (6) are eigenvalue and natural frequency parameters, respectively. Cij are the constants to be determined by solving the relevant eigenvalue problem. The factor “1/k” in (5a) is called crack compliance. This is, for a rectangular beam with crack depth ratio of h-1, defined as follows [3, 6, 7, 11]:
(7)1k=2hEI(h-11-h-1)2×(5.93-19.69h-1+37.14h-12-35.84h-13+13.12h-14),h-1=h1h.Figure 1
Bridge-vehicle system.To solve the coupled equations (1) and (2) vibration modes of the cracked beam are obtained first by solving the relevant eigenvalue problem (for details see [11]). Then, the series solution y(z,t)=Y(z)Tq(t) is formed to determine dynamic response of the beam, where Y(z) is the vector of size Nm×1 containing vibration modes of the cracked beam, Nm is the number of modes used, T means transpose, and q(t) stands for the modal coordinates. Substituting y(z,t)=Y(z)Tq(t) into (2), premultiplying the equation by Y(z), integrating from 0 to L, and finally combining with (1) lead to the coupled beam-vehicle differential equations [20, 21]:
(8)[MNm×Nm0Nm×404×NmM-4×4]{q¨d¨}+[C+A-8A-10-A6C-]{q˙d˙}+[K+A-7A-9-A5K-]{qd}={F10},
where(9a)M=∫0LρAYYTdz,(9b)C=∫0LCYYTdz,(9c)K=∫0LEIYY′′′′Tdz,
where M-,C-, and K- are the mass, damping, and stiffness matrices in (1), that is,
(9d)M-=[Iv0000mv0000m10000m2],(9e)C-=[c1L12+c2L22L2c2-L1c1L1c1-L2c2L2c2-L1c1c2+c1-c1-c2L1c1-c1c3+c10-L2c2-c20c2+c4],(9f)K-=[k1L12+k2L22L2k2-L1k1L1k1-L2k2L2k2-L1k1k2+k1-k1-k2L1k1-k1k3+k10-L2k2-k20k2+k4],where d is the vector including the coordinates di in (1). Additionally,(10a)A1=[Y(z)Tk3+Vc3Y′(z)T]δI,(10b)A2=[Y(z)Tc3]δI,(10c)A3=[Y(z)Tk4+Vc4Y′(z)T]δII,(10d)A4=[Y(z)Tc4]δII,(10e)A5=[02×NmA1A2],A6=[02×NmA2A4],(10f)A7=A1+A3,(10g)A8=A2+A4,(10h)A9=[00-k3-k4],(10i)A10=[00-c3-c4],(10j)A-j=∫0LY(z)Ajdz,j=7,8,9,10,(10k)F1=∫0LY(z)(W1δI+W2δII)dz,(10l)WI=(m1+L2L1+L2mv)9.81,(10m)WII=(m2+L1L1+L2mv)9.81.Newmark-Beta method (with β=1/6 and γ=1/2 [22]) was employed in this work to solve this set of equations. For this purpose, computer codes were written in MATLAB environment. To check the reliability of the codes the midspan response of the beam subject to moving vehicle is obtained as shown in Figure 2, employing the same parameters in [21]. One can verify the agreement by comparing Figure 2 with the relevant figure in [21].Figure 2
Midspan deflection of beam (vehicle and beam parameters are from the work of Wu and Law [21]).


### 2.2. The Objective Function and Constrains

The aim is to correlate the response of the damaged beam to that calculated by the mathematical model of the structure. By this way damage configuration of the carrying structure will be determined. To achieve this, it is proposed to adjust crack sizes and locations by solving an optimization problem. The objective function of the problem is introduced as follows:(11)f1(x)=∑n=1Nmp∫0T|y(zn,t)-y-(zn,t)|max(|y-(zn,t)|)dt,
where Nmp is the number of measurement points on the beam, zn is the location of nth measurement point on the beam, y- denotes the reference displacements (RDs) measured from damaged beam, y stands for the corresponding displacements computed by the mathematical model of the structure, T is the total time for the vehicle to move across the beam, and x is the vector containing crack location and depth parameters, that is,
(12)x={h-1,h-2,…,h-Nc,z-1,z-2,…,z-Nc},h-i=hih,z-i=ziL,
where h is the height of beam cross-section and Nc denotes the number of cracks. As to the measurement points, it is better to choose them close to the midpoint, since maximum deflections occur around this point. Thus, five points on the beam were determined as {0.30.40.50.60.7}L, that is, Nmp=5 in (11). Since number of cracks may not be known in advance, it is required to make initial assumption about their number. To this end, the beam span is divided into five equal intervals each of which is accompanied by a crack (see Figure 3). Thus, ten design parameters, that is, five crack locations and five depth ratios, should be determined to estimate damage condition of the structure. In general, if the beam span is divided into N intervals, then the number of variables to be determined is 2N. Ideally, the number of beam spans should be as large as possible. But this leads to more computational effort. Thus, based on experience, the beam is divided into five spans each of which is associated with a crack. With these explanations in mind, the optimization problem can be formulated as follows:
(13)minf(x)subjectto:0≤h-j<1,0.2(j-1)≤z-j<0.2j,j=1,2,…,Nc.
Crack locations and depths are determined by solving (13). In this work the PSO algorithm is employed for this purpose, and its details are given in the next section.Figure 3
Five beam spans with cracks.


### 2.3. The Particle Swarm Optimization Algorithm

PSO algorithm is initialized with a “swarm” composed ofN particles. Particles refer to the candidate points in the search space of the optimization problem. To obtain best solution each particle adjusts its trajectory toward its own previous best position and toward the previous best position of swarm. By this way, each particle moves in the search space with an adaptive velocity and stores the best position of the search space. Location (x) and velocity (v) of a particle are updated with the following equations [13, 14, 23]:
(14)vijk+1=wvijk+c1R1(pijk-xijk)+c2R2(pgjk-xijk),xijk+1=xijk+vijk+1,i=1,2,…,N,j=1,2,…,m,k=1,2,…,Kmax,
where k is the iteration counter, Kmax denotes the maximum number of iterations, m is dimension of the problem, pijk and pgjk are, respectively, the best positions of the ith particle and the swarm found until the kth iteration, R1 and R2∈U(0,1), where U means the uniform random distribution, and c1 and c2 are positive weighting constants called cognitive and social coefficients, respectively. These two constants regulate the relative velocity toward global and local best points. w is the inertia coefficient employed to prevent swarm explosion. In this work the values of these parameters are selected as w=0.6, c1=c2=1.7 [23]. Initial values of particles and their velocities were obtained drawing random numbers within the range of each dimension, that is,
(15)xij1=xLB,j+R·(xUB,j-xLB,j),vij1=xLB,j+R·(xUB,j-xLB,j),
where xLB,j and xUB,j are, respectively, the lower and upper boundary values of the jth dimension and R∈U(0,1). Besides, a particle moving beyond ranges was bounced back to the search space in the following way:
(16)xijk=xUB,j-R·(xijk-xUB,j),ifxijk>xUB,j,xijk=xLB,j+R·(xLB,j-xijk),ifxijk<xLB,j.


### 2.4. The Proposed Objective Function

In practice, the RDs in (11) should be measured from damaged structure. However, in this work they are produced by the mathematical model due to the lack of such experimental data. To simulate the effect of experimental noise arising from instrumentation, environmental conditions, numerical errors, and so forth, a sequence of random numbers are added to RDs in the following way [7, 9]:
(17)y-(z,t)noisy=y-(z,t)+Np·G·σ,
where y-(z,t) is the calculated RD, which is free of noise, at point z of the damaged beam. Np is the noise percentage, G is Gaussian distribution with zero mean and unit standard deviation, and σ is the standard deviation of y-(z,t). If noisy data is employed in (11), then erroneous results may be obtained. Thus, the noise should be alleviated by a proper approach. To this end, multiresolution property of the DWT is exploited. By the DWT a time function f(t) can be decomposed to an approximation function (AM) at Mth level and sum of details (Dm) up to that level as follows [24, 25]:
(18)f(t)=AM+∑m=-∞MDm,AM=∑n=-∞∞SM,nθM,n(t),Dm=∑n=-∞∞Wm,nψm,n(t),Sm,n=∫-∞∞f(t)θm,n(t)dtWm,n=∫-∞∞f(t)ψm,n(t)dt,ψm,n(t)=2-m/2ψ(t-n2m2m),
where Sm,n is the scale coefficient at mth level, Wm,n is the wavelet coefficient, ψ is the wavelet function, and θ is the scale function associated with ψ. AM represents smooth part of the signal, whereas details include high-frequency components such as noise. If suitable wavelet and decomposition level (M) are determined, then the approximation function of y- in (11) can be used instead of y- itself. In this case, the objective function in (11) can be revised as follows:
(19)f2(x)=∑n=1Nmp∫0T|y(zn,t)-y-M(zn,t)|max(|y-M(zn,t)|)dt,
where y-M is the approximation function of y- at Mth decomposition level.


## 2.1. Dynamic Response of the Vehicle-Beam System

Figure1 illustrates the bridge-vehicle system. Euler-Bernoulli beam model is considered for the bridge, and half-car model is adopted for the vehicle moving with speed V. An open crack with depth h1 is located at z1 on the beam. Surface unevenness of the beam is disregarded, and tires are assumed to be always in contact with the beam [9]. Linear deflections of the beam are considered. Longitudinal deformations and bending-longitudinal coupling are neglected. Under these assumptions the equations of motion for the vehicle-beam system can be derived as follows:
(1)[Iv0000mv0000m10000m2]{d¨1d¨2d¨3d¨4}+[c1L12+c2L22L2c2-L1c1L1c1-L2c2L2c2-L1c1c2+c1-c1-c2L1c1-c1c3+c10-L2c2-c20c2+c4]{d˙1d˙2d˙3d˙4}+[k1L12+k2L22L2k2-L1k1L1k1-L2k2L2k2-L1k1k2+k1-k1-k2L1k1-k1k3+k10-L2k2-k20k2+k4]{d1d2d3d4}={00k3δIu+c3δIu˙k4δIIu+c4δIIu˙},(2)ρA∂2y(z,t)∂t2+C∂y(z,t)∂t+EI∂4y(z,t)∂x4=-PI-PII,
where Iv,mv,m1,m2,L1,L2,k1,k2,k3,c1,c2,c3,c4 are the vehicle parameters shown in Figure 1 and C is the damping of beam material. di, i=1:4, denote vehicle degrees of freedom, and δI and δII are the Dirac delta functions defined as follows:
(3)δII=δ(z-Vt),δI=δ(z-Vt-(L1+L2)),
where PI and PII are the interaction forces acting on the beam through the contact points I and II as follows:(4a)PI=(m1+L2L1+L2mv)9.81+k3(δIu-d3)+c3(δIu˙-d˙3),(4b)PII=(m2+L1L1+L2mv)9.81+k4(δIIu-d4)+c4(δIIu˙-d˙4),where u is the vertical displacement at the tire contact point, that is, u=y(z,t), and its derivative with respect to time is u˙=dy(z,t)/dt=V∂y(z,t)/∂z+∂y(z,t)/∂t and ′=d/dz. At the crack location vibration modes should satisfy the following compatibility equations [7, 11]:(5a)Y1(z1)=Y2(z1),Y1′(z1)+EIkY1′′(z1)=Y2′(z1),(5b)Y1′′(z1)=Y2′′(z1),Y1′′′(z1)=Y2′′′(z1),where
(6)Yi(z)=Ci1cos(λz)+Ci2sin(λz)+Ci3cosh(λz)+Ci4sinh(λz),i=1,2,λ=(ρAω2(EI))0.25,
where λ and ω in (6) are eigenvalue and natural frequency parameters, respectively. Cij are the constants to be determined by solving the relevant eigenvalue problem. The factor “1/k” in (5a) is called crack compliance. This is, for a rectangular beam with crack depth ratio of h-1, defined as follows [3, 6, 7, 11]:
(7)1k=2hEI(h-11-h-1)2×(5.93-19.69h-1+37.14h-12-35.84h-13+13.12h-14),h-1=h1h.Figure 1
Bridge-vehicle system.To solve the coupled equations (1) and (2) vibration modes of the cracked beam are obtained first by solving the relevant eigenvalue problem (for details see [11]). Then, the series solution y(z,t)=Y(z)Tq(t) is formed to determine dynamic response of the beam, where Y(z) is the vector of size Nm×1 containing vibration modes of the cracked beam, Nm is the number of modes used, T means transpose, and q(t) stands for the modal coordinates. Substituting y(z,t)=Y(z)Tq(t) into (2), premultiplying the equation by Y(z), integrating from 0 to L, and finally combining with (1) lead to the coupled beam-vehicle differential equations [20, 21]:
(8)[MNm×Nm0Nm×404×NmM-4×4]{q¨d¨}+[C+A-8A-10-A6C-]{q˙d˙}+[K+A-7A-9-A5K-]{qd}={F10},
where(9a)M=∫0LρAYYTdz,(9b)C=∫0LCYYTdz,(9c)K=∫0LEIYY′′′′Tdz,
where M-,C-, and K- are the mass, damping, and stiffness matrices in (1), that is,
(9d)M-=[Iv0000mv0000m10000m2],(9e)C-=[c1L12+c2L22L2c2-L1c1L1c1-L2c2L2c2-L1c1c2+c1-c1-c2L1c1-c1c3+c10-L2c2-c20c2+c4],(9f)K-=[k1L12+k2L22L2k2-L1k1L1k1-L2k2L2k2-L1k1k2+k1-k1-k2L1k1-k1k3+k10-L2k2-k20k2+k4],where d is the vector including the coordinates di in (1). Additionally,(10a)A1=[Y(z)Tk3+Vc3Y′(z)T]δI,(10b)A2=[Y(z)Tc3]δI,(10c)A3=[Y(z)Tk4+Vc4Y′(z)T]δII,(10d)A4=[Y(z)Tc4]δII,(10e)A5=[02×NmA1A2],A6=[02×NmA2A4],(10f)A7=A1+A3,(10g)A8=A2+A4,(10h)A9=[00-k3-k4],(10i)A10=[00-c3-c4],(10j)A-j=∫0LY(z)Ajdz,j=7,8,9,10,(10k)F1=∫0LY(z)(W1δI+W2δII)dz,(10l)WI=(m1+L2L1+L2mv)9.81,(10m)WII=(m2+L1L1+L2mv)9.81.Newmark-Beta method (with β=1/6 and γ=1/2 [22]) was employed in this work to solve this set of equations. For this purpose, computer codes were written in MATLAB environment. To check the reliability of the codes the midspan response of the beam subject to moving vehicle is obtained as shown in Figure 2, employing the same parameters in [21]. One can verify the agreement by comparing Figure 2 with the relevant figure in [21].Figure 2
Midspan deflection of beam (vehicle and beam parameters are from the work of Wu and Law [21]).


## 2.2. The Objective Function and Constrains

The aim is to correlate the response of the damaged beam to that calculated by the mathematical model of the structure. By this way damage configuration of the carrying structure will be determined. To achieve this, it is proposed to adjust crack sizes and locations by solving an optimization problem. The objective function of the problem is introduced as follows:(11)f1(x)=∑n=1Nmp∫0T|y(zn,t)-y-(zn,t)|max(|y-(zn,t)|)dt,
where Nmp is the number of measurement points on the beam, zn is the location of nth measurement point on the beam, y- denotes the reference displacements (RDs) measured from damaged beam, y stands for the corresponding displacements computed by the mathematical model of the structure, T is the total time for the vehicle to move across the beam, and x is the vector containing crack location and depth parameters, that is,
(12)x={h-1,h-2,…,h-Nc,z-1,z-2,…,z-Nc},h-i=hih,z-i=ziL,
where h is the height of beam cross-section and Nc denotes the number of cracks. As to the measurement points, it is better to choose them close to the midpoint, since maximum deflections occur around this point. Thus, five points on the beam were determined as {0.30.40.50.60.7}L, that is, Nmp=5 in (11). Since number of cracks may not be known in advance, it is required to make initial assumption about their number. To this end, the beam span is divided into five equal intervals each of which is accompanied by a crack (see Figure 3). Thus, ten design parameters, that is, five crack locations and five depth ratios, should be determined to estimate damage condition of the structure. In general, if the beam span is divided into N intervals, then the number of variables to be determined is 2N. Ideally, the number of beam spans should be as large as possible. But this leads to more computational effort. Thus, based on experience, the beam is divided into five spans each of which is associated with a crack. With these explanations in mind, the optimization problem can be formulated as follows:
(13)minf(x)subjectto:0≤h-j<1,0.2(j-1)≤z-j<0.2j,j=1,2,…,Nc.
Crack locations and depths are determined by solving (13). In this work the PSO algorithm is employed for this purpose, and its details are given in the next section.Figure 3
Five beam spans with cracks.


## 2.3. The Particle Swarm Optimization Algorithm

PSO algorithm is initialized with a “swarm” composed ofN particles. Particles refer to the candidate points in the search space of the optimization problem. To obtain best solution each particle adjusts its trajectory toward its own previous best position and toward the previous best position of swarm. By this way, each particle moves in the search space with an adaptive velocity and stores the best position of the search space. Location (x) and velocity (v) of a particle are updated with the following equations [13, 14, 23]:
(14)vijk+1=wvijk+c1R1(pijk-xijk)+c2R2(pgjk-xijk),xijk+1=xijk+vijk+1,i=1,2,…,N,j=1,2,…,m,k=1,2,…,Kmax,
where k is the iteration counter, Kmax denotes the maximum number of iterations, m is dimension of the problem, pijk and pgjk are, respectively, the best positions of the ith particle and the swarm found until the kth iteration, R1 and R2∈U(0,1), where U means the uniform random distribution, and c1 and c2 are positive weighting constants called cognitive and social coefficients, respectively. These two constants regulate the relative velocity toward global and local best points. w is the inertia coefficient employed to prevent swarm explosion. In this work the values of these parameters are selected as w=0.6, c1=c2=1.7 [23]. Initial values of particles and their velocities were obtained drawing random numbers within the range of each dimension, that is,
(15)xij1=xLB,j+R·(xUB,j-xLB,j),vij1=xLB,j+R·(xUB,j-xLB,j),
where xLB,j and xUB,j are, respectively, the lower and upper boundary values of the jth dimension and R∈U(0,1). Besides, a particle moving beyond ranges was bounced back to the search space in the following way:
(16)xijk=xUB,j-R·(xijk-xUB,j),ifxijk>xUB,j,xijk=xLB,j+R·(xLB,j-xijk),ifxijk<xLB,j.


## 2.4. The Proposed Objective Function

In practice, the RDs in (11) should be measured from damaged structure. However, in this work they are produced by the mathematical model due to the lack of such experimental data. To simulate the effect of experimental noise arising from instrumentation, environmental conditions, numerical errors, and so forth, a sequence of random numbers are added to RDs in the following way [7, 9]:
(17)y-(z,t)noisy=y-(z,t)+Np·G·σ,
where y-(z,t) is the calculated RD, which is free of noise, at point z of the damaged beam. Np is the noise percentage, G is Gaussian distribution with zero mean and unit standard deviation, and σ is the standard deviation of y-(z,t). If noisy data is employed in (11), then erroneous results may be obtained. Thus, the noise should be alleviated by a proper approach. To this end, multiresolution property of the DWT is exploited. By the DWT a time function f(t) can be decomposed to an approximation function (AM) at Mth level and sum of details (Dm) up to that level as follows [24, 25]:
(18)f(t)=AM+∑m=-∞MDm,AM=∑n=-∞∞SM,nθM,n(t),Dm=∑n=-∞∞Wm,nψm,n(t),Sm,n=∫-∞∞f(t)θm,n(t)dtWm,n=∫-∞∞f(t)ψm,n(t)dt,ψm,n(t)=2-m/2ψ(t-n2m2m),
where Sm,n is the scale coefficient at mth level, Wm,n is the wavelet coefficient, ψ is the wavelet function, and θ is the scale function associated with ψ. AM represents smooth part of the signal, whereas details include high-frequency components such as noise. If suitable wavelet and decomposition level (M) are determined, then the approximation function of y- in (11) can be used instead of y- itself. In this case, the objective function in (11) can be revised as follows:
(19)f2(x)=∑n=1Nmp∫0T|y(zn,t)-y-M(zn,t)|max(|y-M(zn,t)|)dt,
where y-M is the approximation function of y- at Mth decomposition level.


## 3. Results and Discussion

### 3.1. Application

The following data [9] is employed for the numerical simulation: E=210 GPa, ρ=7860 kgm−3, L=50 m, b=1 m, h=2 m, mv=12404 kg, m1=m2=725.4 kg, Iv=172160 kgm2, k1=1969034 Nm−1, k2=727812 Nm−1, k3=4735000 Nm−1, k4=1972900 Nm−1, c1=7182 Nsm−1, c2=2190 Nsm−1, c3=c4=0, and L1=L2=3 m. The first five vibration modes of the beam, for which the natural frequencies are 1.9, 7.5, 16.9, 30, and 46.9 Hz, are considered. The sampling frequency of the simulation is 300 Hz which can capture the response of the first five vibration modes of the beam. Two percent modal damping is considered for each mode [7]. The damage scenarios in Table 1 are considered. In case 1 single crack with depth ratio of 0.1 is available at 25% of the beam length. In case 2 and case 3 multiple cracks exist on the beam. Note that the considered crack depths are small according to the relevant works [7–9]. Also, one may verify by the relevant works that the considered vehicle speeds are relatively high for these small cracks. In other words, the wavelet methods in the cited works cannot extract such small cracks through processing dynamic response recorded at higher speeds such as those in Table 1. However, it will be demonstrated that crack identification is possible for the cases in Table 1 by the method proposed in this work. The beam midspan deflections are compared in Figure 4. It is clear that deflections are relatively small and not bigger than 3 mm.Table 1
Damage scenarios.
Case
CL (zi/L)
CDR (hi/h)
V (m/s)
1
0.25
0.1
30
2
0.33, 0.67
0.1, 0.15
20
3
0.3, 0.5, 0.7
0.15, 0.1, 0.1
20
CL: crack location, CDR: crack depth ratio,V: vehicle speed.Comparison of beam midspan deflections for the cases in Table1. U: undamaged, D: damaged.
(a)
(b)
(c)For DWT a widely used wavelet, that is, Daubechies wavelet with four vanishing moments (db4) [25] shown in Figure 5, is employed. It is experienced that the approximation function at the fourth level (M=4) is suitable since it is well correlated with noiseless data. This is demonstrated in Figure 6, considering Case 1 in Table 1. The midspan deflection of the beam is obtained for Case 1, and it is contaminated adding 5% Gaussian distribution. Then, the approximation (AM) of this noisy data at the first level (M=1) is obtained, and its comparison with the original data that is free of noise is illustrated in Figure 6(a). It is clear that AM at this level is not suitable since it is not denoised sufficiently. This is true for those AM at the levels M=2 and 3 (hence their graphs are not given). However, the AM at M=4 is in well agreement with the noiseless data (see Figure 6(b)); thus, it can be employed as reference data in (19), that is, y-M=4. At higher levels the AM deviates from the noiseless data, since not only high-frequency components are removed but also some low-frequency parts that affect shape of the signal are removed. This observation is valid for other RDs on other points of the beam. Generally the suitable level depends on the length of reference data. That is, it may be higher when data is recorded at higher sampling frequencies.Figure 5
The db4 wavelet (dashed) and its scale function (solid).Comparison ofAM at various levels with the original noiseless response of beam midspan.
(a)
(b)
(c)
(d)After some trials the swarm size and maximum iteration number are determined as 30 and 100, respectively. Then the optimization problem is solved using the objective functionsf1 and f2 in (11) and (19), and damage detection results are given in Figures 7 and 8, respectively. As the PSO algorithm is stochastic, the mean of 30 runs [16, 18] of the algorithm is displayed in the figures. Ideally, only the crack depth ratios at the true damage locations are expected to be nonzero. However, this is not sufficiently satisfied by the objective function f1 in (11), the reason of which is the noise in RDs. For example, for Case 1 it seems there is a crack with depth ratio about 0.16 around the 0.07L on the beam (see Figure 7) though the only crack in this case is at 0.25L of the beam length and with depth ratio of 0.1. On the other hand, the standard deviation (SD) of crack depth ratio at true crack location is low compared with the value of relevant parameter. However, SD at a misleading crack point is high compared to the relevant crack depth ratio, which indicates that this is not real damage. Generally, SD at points close to beam ends is higher. This is because the curvature of vibration modes around the beam ends is small for simply supported beam [26].Damage detection results (the objective function in (11) is employed).
(a)
(b)
(c)Damage detection results (the objective function in (19) is employed).
(a)
(b)
(c)Thus cracks close to beam ends do not have much effect on the dynamic response of the beam. According to Figure8, all real crack points are predicted more accurately, and SD values are generally lower when compared to Figure 7. This implies that the proposed method produces more reliable results.


### 3.2. Comparison with CWT Coefficients Approach

Zhu and Law [7] demonstrated that magnitudes of CWT coefficients of beam dynamic response become large when moving load passes over crack point. The reason is that crack induces localized singularities in the structural response, and the CWT coefficients of the response are sensitive to these singularities. However, the deficiency of this approach is that CWT coefficients become insensitive to damage with increasing moving load speed. Hester and González [8] report that it is possible to detect crack depth ratios of 0.2 and 0.3 at the vehicle speed 7 m/s by their approach. However, their method fails to detect depth ratio of 0.1 at the same speed. Nguyen and Tran [9] state that cracks as small as 10% of beam height can be detected when the vehicle is moving at low velocity such as 2 m/s. The authors conclude that these small cracks are very difficult to be revealed at high vehicle speeds. Thus, one may verify that in their work the authors consider crack depth ratio of 0.5 to detect crack at the vehicle speed of 30 m/s. From the relevant literature one can conclude that crack depth ratio of 0.1 cannot be detected at a vehicle speed such as 20 or 30 m/s by the wavelet transform methods. In order to prove this, Case 2 in Table 1 is considered. Midspan deflection and its CWT coefficients are demonstrated in Figure 9. The scale interval is determined as in [8] using the relation Fa=Fc/(aΔ), where Fc denotes the centre frequency of the Mexican Hat wavelet (0.25 Hz), a is the scale, Δ is the sampling period, and Fa is the frequency corresponding to the scale a. The maximum scale corresponds to the first natural frequency of the healthy beam, that is, 1.9 Hz. At a singularity point the CWT coefficients form a cone of influence [27]; that is, absolute values of the CWT coefficients become extreme and the region affected by the singularity increases along with scale. However, this variation is not observed in Figure 9; hence, one can conclude that the CWT coefficients fail to extract damage signature for the considered case.The midspan deflection for Case 2 in Table1 and its CWT coefficients.
(a)
(b)


## 3.1. Application

The following data [9] is employed for the numerical simulation: E=210 GPa, ρ=7860 kgm−3, L=50 m, b=1 m, h=2 m, mv=12404 kg, m1=m2=725.4 kg, Iv=172160 kgm2, k1=1969034 Nm−1, k2=727812 Nm−1, k3=4735000 Nm−1, k4=1972900 Nm−1, c1=7182 Nsm−1, c2=2190 Nsm−1, c3=c4=0, and L1=L2=3 m. The first five vibration modes of the beam, for which the natural frequencies are 1.9, 7.5, 16.9, 30, and 46.9 Hz, are considered. The sampling frequency of the simulation is 300 Hz which can capture the response of the first five vibration modes of the beam. Two percent modal damping is considered for each mode [7]. The damage scenarios in Table 1 are considered. In case 1 single crack with depth ratio of 0.1 is available at 25% of the beam length. In case 2 and case 3 multiple cracks exist on the beam. Note that the considered crack depths are small according to the relevant works [7–9]. Also, one may verify by the relevant works that the considered vehicle speeds are relatively high for these small cracks. In other words, the wavelet methods in the cited works cannot extract such small cracks through processing dynamic response recorded at higher speeds such as those in Table 1. However, it will be demonstrated that crack identification is possible for the cases in Table 1 by the method proposed in this work. The beam midspan deflections are compared in Figure 4. It is clear that deflections are relatively small and not bigger than 3 mm.Table 1
Damage scenarios.
Case
CL (zi/L)
CDR (hi/h)
V (m/s)
1
0.25
0.1
30
2
0.33, 0.67
0.1, 0.15
20
3
0.3, 0.5, 0.7
0.15, 0.1, 0.1
20
CL: crack location, CDR: crack depth ratio,V: vehicle speed.Comparison of beam midspan deflections for the cases in Table1. U: undamaged, D: damaged.
(a)
(b)
(c)For DWT a widely used wavelet, that is, Daubechies wavelet with four vanishing moments (db4) [25] shown in Figure 5, is employed. It is experienced that the approximation function at the fourth level (M=4) is suitable since it is well correlated with noiseless data. This is demonstrated in Figure 6, considering Case 1 in Table 1. The midspan deflection of the beam is obtained for Case 1, and it is contaminated adding 5% Gaussian distribution. Then, the approximation (AM) of this noisy data at the first level (M=1) is obtained, and its comparison with the original data that is free of noise is illustrated in Figure 6(a). It is clear that AM at this level is not suitable since it is not denoised sufficiently. This is true for those AM at the levels M=2 and 3 (hence their graphs are not given). However, the AM at M=4 is in well agreement with the noiseless data (see Figure 6(b)); thus, it can be employed as reference data in (19), that is, y-M=4. At higher levels the AM deviates from the noiseless data, since not only high-frequency components are removed but also some low-frequency parts that affect shape of the signal are removed. This observation is valid for other RDs on other points of the beam. Generally the suitable level depends on the length of reference data. That is, it may be higher when data is recorded at higher sampling frequencies.Figure 5
The db4 wavelet (dashed) and its scale function (solid).Comparison ofAM at various levels with the original noiseless response of beam midspan.
(a)
(b)
(c)
(d)After some trials the swarm size and maximum iteration number are determined as 30 and 100, respectively. Then the optimization problem is solved using the objective functionsf1 and f2 in (11) and (19), and damage detection results are given in Figures 7 and 8, respectively. As the PSO algorithm is stochastic, the mean of 30 runs [16, 18] of the algorithm is displayed in the figures. Ideally, only the crack depth ratios at the true damage locations are expected to be nonzero. However, this is not sufficiently satisfied by the objective function f1 in (11), the reason of which is the noise in RDs. For example, for Case 1 it seems there is a crack with depth ratio about 0.16 around the 0.07L on the beam (see Figure 7) though the only crack in this case is at 0.25L of the beam length and with depth ratio of 0.1. On the other hand, the standard deviation (SD) of crack depth ratio at true crack location is low compared with the value of relevant parameter. However, SD at a misleading crack point is high compared to the relevant crack depth ratio, which indicates that this is not real damage. Generally, SD at points close to beam ends is higher. This is because the curvature of vibration modes around the beam ends is small for simply supported beam [26].Damage detection results (the objective function in (11) is employed).
(a)
(b)
(c)Damage detection results (the objective function in (19) is employed).
(a)
(b)
(c)Thus cracks close to beam ends do not have much effect on the dynamic response of the beam. According to Figure8, all real crack points are predicted more accurately, and SD values are generally lower when compared to Figure 7. This implies that the proposed method produces more reliable results.


## 3.2. Comparison with CWT Coefficients Approach

Zhu and Law [7] demonstrated that magnitudes of CWT coefficients of beam dynamic response become large when moving load passes over crack point. The reason is that crack induces localized singularities in the structural response, and the CWT coefficients of the response are sensitive to these singularities. However, the deficiency of this approach is that CWT coefficients become insensitive to damage with increasing moving load speed. Hester and González [8] report that it is possible to detect crack depth ratios of 0.2 and 0.3 at the vehicle speed 7 m/s by their approach. However, their method fails to detect depth ratio of 0.1 at the same speed. Nguyen and Tran [9] state that cracks as small as 10% of beam height can be detected when the vehicle is moving at low velocity such as 2 m/s. The authors conclude that these small cracks are very difficult to be revealed at high vehicle speeds. Thus, one may verify that in their work the authors consider crack depth ratio of 0.5 to detect crack at the vehicle speed of 30 m/s. From the relevant literature one can conclude that crack depth ratio of 0.1 cannot be detected at a vehicle speed such as 20 or 30 m/s by the wavelet transform methods. In order to prove this, Case 2 in Table 1 is considered. Midspan deflection and its CWT coefficients are demonstrated in Figure 9. The scale interval is determined as in [8] using the relation Fa=Fc/(aΔ), where Fc denotes the centre frequency of the Mexican Hat wavelet (0.25 Hz), a is the scale, Δ is the sampling period, and Fa is the frequency corresponding to the scale a. The maximum scale corresponds to the first natural frequency of the healthy beam, that is, 1.9 Hz. At a singularity point the CWT coefficients form a cone of influence [27]; that is, absolute values of the CWT coefficients become extreme and the region affected by the singularity increases along with scale. However, this variation is not observed in Figure 9; hence, one can conclude that the CWT coefficients fail to extract damage signature for the considered case.The midspan deflection for Case 2 in Table1 and its CWT coefficients.
(a)
(b)


## 4. Conclusions

In this work a damage detection method is proposed for bridge type structures under vehicular load. An objective function is defined by the time-dependent deflections on several points of damaged beam and the mathematical model of the structure. The reference data is denoised using the multiresolution property of the DWT. Then, the optimization problem is solved to determine crack locations and depths. It was observed that the proposed approach achieves identification of small cracks with depth ratio of 0.1 in spite of 5% noise interference and is superior than the CWT coefficients method. The present study should be verified by experimental data. This issue is planned to be dealt with in later studies.



---

*Source: 634217-2013-06-15.xml*