# Simultaneous Identification of Moving Vehicles and Bridge Damages Considering Road Rough Surface

Authors: Qingxia Zhang; Łukasz Jankowski; Zhongdong Duan
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A method for the simultaneous identification of moving vehicles and the damages of the supporting structure from measured responses is presented. A two-axle vehicle model with two degrees of freedom (DOF) is adopted. The extent of the damage and the vehicle parameters were chosen as the optimisation variables, which allow ill conditioning to be avoided and decrease the number of sensors required. The identification is performed by minimising the distance between the measured responses and the computed responses to given optimisation variables. The virtual distortion method (VDM) was used, such that the response of the damaged structure can be computed from comparison with the intact structure subjected to the same vehicle excitation and to the response-coupled virtual distortions. These are related to the optimisation variables by the system impulse response matrix and are expressed by a linear system, which allowed both types of optimisation variables to be treated in a unified way. The numerical cost is reduced by using a moving influence matrix. The adjoint variable method is used for fast sensitivity analysis. A three-span bridge numerical example is presented, where the identification was verified with 5% root mean square (RMS) measurement, and model, error whilst also considering the surface roughness of the road.


---

## Body

## 1. Introduction

In structural health monitoring, accurate load and damage identification are indispensable to the maintenance of structural integrity, as well as providing the evidence for forensic engineering. Particularly, moving loads or moving vehicle parameter identification is important not only for prediction and analysis of the dynamic responses of bridges, but also in traffic studies, in design code calibration, for traffic control, and so forth. In recent decades, many investigations have been performed on either load identification or damage identification. However, in practice, both unknown system damage and unknown (moving) excitations usually coexist and together influence the system response in a coupled way. Therefore, their simultaneous identification is worthy of further exploration.As opposed to local, high-frequency ultrasonic scanning, this paper considers only global or low-frequency damage identification methods which are based on the structural vibrations. These methods are mainly categorised into two groups: frequency domain methods or time domain methods. In the frequency domain, damage is detected, located, and identified through measured changes in the related modal parameters or dynamic signature [1] such as natural frequency and mode shape. However, the modal parameters are obtained from the responses to a certain kind of excitation and these methods are limited to time-invariant systems. Under the action of a moving vehicle, a coupled bridge-vehicle system is time variant, and the system parameters are changing as the vehicles move [2]. The vibration caused by traffic is a nonstationary process that strengthens with decreasing span [3]. In this case, the analysis is most often performed in the time domain by the direct comparison of the simulated and measured responses [4]. However, this usually requires a known moving load or other such vehicle parameters. Sieniawska et al. [5] identify parameters of a linear structure from its responses to a moving load using a static substitution of the equation of motion, where the moving load is a known constant magnitude. Chang et al. [3, 6] present, and experimentally verify, a methodology for damage identification in bridges using a pseudostatic formulation based on coupling vibration with a moving vehicle, where the moving force is provided as a vector or calculated using the measured acceleration responses of the axles. Ling and Haldar [7] propose a linear system identification approach with unknown input, but they require a known dynamic response of all the degrees of freedoms; this is hard to perform on real-life, complex structures.Moving load identification has been studied extensively in the past two decades [8]. Techniques of indirect identification from measured responses have been investigated and can be performed more easily and at lower cost compared to methods requiring direct measurements of moving loads. Chan et al. and Law et al. have proposed four methods for indirect identification: the time-domain method (TDM) [9], the frequency-time domain method (FTDM) [10], interpretive method I (IM-I) [11], and interpretive method II (IM-II) [12]. All of them require, a priori, known model parameters for the bridge. Each method has its merits and limitations, which are compared elsewhere [13]. Load identification is an inverse problem, and numerical ill-conditioning seems to be the main factor that influences the accuracy of the identification results. The accuracy can be improved by some regularisation methods, such as Tikhonov regularisation [14]. However, the determination of the optimal value of the regularisation parameter is numerically costly and requires lengthy computation. Pinkaew and Asnachinda [15, 16] point out that the regularisation parameter turns out to be sensitive to properties of both the vehicle and bridge and is difficult to be precisely assigned: they propose, and experimentally verify, an iterative regularisation method called the updated static component (USC) technique to decrease the sensitivity of the regularisation parameter. Moving forces are usually treated as unknown quantities to be identified. This requires a number of sensors equal to or exceeding the number of unknown moving forces to obtain a unique solution. Jiang et al. [17] present a moving vehicle parameter identification method based on genetic algorithms, where each moving vehicle is modelled as a two-DOF system that comprises mass, spring stiffness, and damping. Zhang et al. [18] provide an approach of identifying moving vehicle parameters based on VDM using a dynamic influence matrix. In general, the identification of moving forces or vehicle(s) requires a known, well-defined, bridge model to establish the relationship between load and response.For unknown coexistent load and structural damage, it is generally difficult to decouple the related identification problems and solve any one of them independently. Since the two factors have essentially different natures, Zhang et al. [19] present a method to identify load and damage simultaneously using Chebyshev polynomials to parameterise the unknown force, and thus all parameters related to the damage and excitation can be updated simultaneously in each iteration. Zhang et al. [20] present, and experimentally verify, a method for the simultaneous identification of nonmoving excitations and damage, which uses a VDM [21] to model structural damage and thus improve the identification efficiency thereof. In case of unknown moving loads, the vehicle-bridge system is a coupled time-varying system. Zhu and Law [22] present a two-step iteration procedure to simultaneously identify the moving loads and the damage caused to an Euler-Bernoulli simply supported beam; the number of the sensors is one less than the number of the beam elements.Zhang et al. [23] present a method for the simultaneous identification of structural damage and moving masses; moving masses and damage extents are used as its optimisation variables, which avoids the ill-conditioning problem present in traditional moving force identification and decreases the number of sensors required; damage is modelled using virtual distortions, and a dynamic influence matrix is introduced to reduce numerical analysis cost. This paper simulates the vehicle using a two-axle, two-DOF, vehicle model, and the related parameters including mass, rotary inertia, spring stiffness, and damping are treated as unknowns together with the structural damage induced. The damaged bridge is assumed to be modelled by the decay of its mechanical properties and is often represented as a decreasing stiffness [24]. Based on VDM, the response of the damaged structure is modelled by the intact structure subjected to the same vehicle excitation and to the response-coupled virtual distortions. Through the related system impulse response matrix (dynamic influence matrix in the terminology of VDM), the moving loads and virtual distortions are connected with the optimisation variables in a linear system. In this way, the optimisation variables related to the vehicle parameter, and to the damage, are treated in a unified way. Then, they can be optimised using any standard optimisation algorithm. The numerical cost is reduced by using the moving dynamic influence matrix. Moreover, rapid sensitivity analysis is performed using the adjoint variable method.Sections2 and 3 discuss the equation of motion of the system and simulation of the system dynamic response by moving vehicle and virtual distortions; Section 4 discusses the inverse identification problem. Section 5 verifies the proposed approach using a numerical example of a three-span frame structure considering road surface roughness. The approach and the results are discussed in Section 6.


## 2. Dynamic Analysis of a Coupled Vehicle-Beam System

### 2.1. Equation of Motion for a Moving Vehicle on a Bridge

A popular two-axle, two-DOF, vehicle model [3] was adopted; see Figure 1, whereuv andθv, respectively, denote the bounce and pitching motions of the vehicle. Exactly speaking,uvis the vertical displacement of the vehicle, which is measured vertically upwards from its vertical static equilibrium position before the vehicle reached the bridge andθivare its rotary displacement relative to the static position. Denoted by m and J, respectively, are the vehicle mass and the rotary inertia; denoted by k and c, respectively, are the spring constant and damping coefficient of the vehicle.Figure 1
A coupled vehicle-bridge system.For theith moving vehicle on the bridge,  the additional elongation of the left and right spring except the elongation caused by the weight is respectively denoted by Δil, and Δir, then
(1)Δil=uiv+θi1eil-xil-ril;Δir=uiv-θi1eir-xir-rir,
whereeil,  eirdenote, respectively, the distance from the left and right axle to the vehicle’s centre of mass;xil,  xir, respectively, denote the vertical displacements of the bridge at the current position of the vehicle’s left and right axle;ril,  rirdenote the roadway surface roughness at the current respective position, which is also measured vertically upwards with reference to the deck level.The forces of theith moving vehicle are shown in Figure 2, which includes the inertial forces. Based on the dynamic balance, the equation of motion of theith moving vehicle can be formulated as
(2)mivY¨iv(t)+cieikY˙iv(t)+kieikYiv(t)=kiSi(xi(t)+ri(t))+ciSi(x˙i(t)+r˙i(t)),
wheremiv=[mi00Ji],Yiv(t)=[uiv(t)θiv(t)],ri(t)=[ril(t)rir(t)],si=[11eil-eir],xi(t)=[xil(t)xir(t)],eik=[2eil-eireil-eireil2+eir2].Figure 2
Forces of theith moving vehicle.In this study, the random roughness of the road surfacer(x)is simulated by a zero-mean, real-value, stationary Gaussian process as described in the following [17, 25]:
(3)r(x)=∑k=1NTαkcos(2πnkx+φk),
whereαkis the amplitude of the cosine wave,nkis the spatial frequency (in cycle/m),φkis a random phase angle with uniform probability distribution in the interval [0,2π],xis the position measured from the left end of the bridge, andNTis the total number of the cosine waves which are used to construct the roughness. Parameterαkis computed by
(4)αk2=4Gd(nk)Δn,
whereGd(nk)is the power spectral density function of the road surface and related to the ground flatness coefficientGd(n0)in the form as [26]
(5)Gd(nk)=Gd(n0)(nkn0)-2,nk=nl+kΔn,Δn=nu-nlNT,
whereGd(n0)depends on the classification of road surface condition.nlandnu, are respectively, lower and upper cut-off spatial frequencies, and the power spectral density functionGd(nk)is defined within the interval of[nl,nu]. For the vehicle velocity between 70 km/h and 120 km/h on expressway, the power spectral density is most meaningful in the frequency interval of [0.0221 m−1, 1.4142 m−1].


### 2.2. Equation of Motion for a Bridge under Moving Vehicles

The bridge is idealized as a beam, which suffers the traffic excitations. Denoted bypil(t)andpir(t), respectively, are the interaction forces of theith moving vehicle at its left and right wheel, and through the dynamic balance analysis, they can be expressed using the additional elongation of the spring and the gravity of the vehicle,
(6)pil(t)=kiΔil(t)+ciΔ˙il(t)+migeireil+e1r,pir(t)=kiΔir(t)+ciΔ˙ir(t)+migeileil+eir.The equation of motion for the bridge can be formulated as follows:(7)Mu¨b(t)+Cu˙b(t)+Kub(t)=∑i=1nmbi(t)pi(t),
whereM,C, andKare, respectively, the mass, stiffness, and damping matrices of the bridge,ubrepresents the displacement vector, and pi(t)=[pil(t)pir(t)]T,bi(t)=[bil(t)bir(t)],bil(t),  bir(t), are respectively, the time-varying load distribution vector in global coordinates with regard to the left and right wheels. Thus,bil(t),  bir(t)can be formulated using the shape function of the beam element in finite element (FE) analysis, and their elements are zeroed when the corresponding wheel is off the bridge during the sampling time.


### 2.3. Dynamic Response of the Coupled Vehicle-Bridge System

In the equation of motion for the moving vehicle and bridge ((2) and (7)), it can be seen that the system is coupled: the vibration of the bridge influenced the vibration of the moving vehicle (2), and the vehicle in return interactively influenced the vibration of the bridge (7). The equation of motion for the coupled system can be expressed as follows:
(8)[Mm1v⋱mnmv][u¨b(t)Y¨1v(t)⋮Y¨nmv(t)]+[K+ΔK(t)-b1(t)S1Tk1⋯-bnm(t)SnmTkd-S1b1T(t)k1k1e1k⋮⋱-SnmbnmT(t)knmknmenmk]×[ub(t)Y1v(t)⋮Ynmv(t)]+[C+ΔC(t)-b1(t)S1Tc1⋯-bnm(t)SnmTcd-S1b1T(t)c1c1e1k⋮⋱-SnmbnmT(t)cnmcnmenmk]×[u˙b(t)Y˙1v(t)⋮Y˙nmv(t)]=[∑ibi(mieig-kiri+cir˙i)k1S1r1+c1S1r˙1⋮knmSnmrnm+cnmSnmr˙nm],
whereΔK(t)=∑ikibi(t)biT(t),ΔC(t)=∑icibi(t)biT(t), ei=[eil/(eil+eir)eir/(eil+eir)]T.Equation (8) shows that the system is time variant with the vehicle’s motion, and this indicated that the vibration of the bridge caused by the traffic was nonstationary. Therefore, the performance of damage identification in this case may have similar limitations to those pertaining to methods subjected to a stationary assumption, for example, modal based methods in the frequency domain. In contrast, time-domain analysis for identification is more accurate and appropriate. Then, the key procedure is the estimation of the system response to damage. However, if the response is computed in the traditional way by using the system motion equation, it will need to reassemble the varying system parameter matrix at each sampling time-step, which will be time-consuming, especially for bridges with complex structures. Aiming at solving this problem, the following section presents a fast reanalysis, VDM-based method.


## 2.1. Equation of Motion for a Moving Vehicle on a Bridge

A popular two-axle, two-DOF, vehicle model [3] was adopted; see Figure 1, whereuv andθv, respectively, denote the bounce and pitching motions of the vehicle. Exactly speaking,uvis the vertical displacement of the vehicle, which is measured vertically upwards from its vertical static equilibrium position before the vehicle reached the bridge andθivare its rotary displacement relative to the static position. Denoted by m and J, respectively, are the vehicle mass and the rotary inertia; denoted by k and c, respectively, are the spring constant and damping coefficient of the vehicle.Figure 1
A coupled vehicle-bridge system.For theith moving vehicle on the bridge,  the additional elongation of the left and right spring except the elongation caused by the weight is respectively denoted by Δil, and Δir, then
(1)Δil=uiv+θi1eil-xil-ril;Δir=uiv-θi1eir-xir-rir,
whereeil,  eirdenote, respectively, the distance from the left and right axle to the vehicle’s centre of mass;xil,  xir, respectively, denote the vertical displacements of the bridge at the current position of the vehicle’s left and right axle;ril,  rirdenote the roadway surface roughness at the current respective position, which is also measured vertically upwards with reference to the deck level.The forces of theith moving vehicle are shown in Figure 2, which includes the inertial forces. Based on the dynamic balance, the equation of motion of theith moving vehicle can be formulated as
(2)mivY¨iv(t)+cieikY˙iv(t)+kieikYiv(t)=kiSi(xi(t)+ri(t))+ciSi(x˙i(t)+r˙i(t)),
wheremiv=[mi00Ji],Yiv(t)=[uiv(t)θiv(t)],ri(t)=[ril(t)rir(t)],si=[11eil-eir],xi(t)=[xil(t)xir(t)],eik=[2eil-eireil-eireil2+eir2].Figure 2
Forces of theith moving vehicle.In this study, the random roughness of the road surfacer(x)is simulated by a zero-mean, real-value, stationary Gaussian process as described in the following [17, 25]:
(3)r(x)=∑k=1NTαkcos(2πnkx+φk),
whereαkis the amplitude of the cosine wave,nkis the spatial frequency (in cycle/m),φkis a random phase angle with uniform probability distribution in the interval [0,2π],xis the position measured from the left end of the bridge, andNTis the total number of the cosine waves which are used to construct the roughness. Parameterαkis computed by
(4)αk2=4Gd(nk)Δn,
whereGd(nk)is the power spectral density function of the road surface and related to the ground flatness coefficientGd(n0)in the form as [26]
(5)Gd(nk)=Gd(n0)(nkn0)-2,nk=nl+kΔn,Δn=nu-nlNT,
whereGd(n0)depends on the classification of road surface condition.nlandnu, are respectively, lower and upper cut-off spatial frequencies, and the power spectral density functionGd(nk)is defined within the interval of[nl,nu]. For the vehicle velocity between 70 km/h and 120 km/h on expressway, the power spectral density is most meaningful in the frequency interval of [0.0221 m−1, 1.4142 m−1].


## 2.2. Equation of Motion for a Bridge under Moving Vehicles

The bridge is idealized as a beam, which suffers the traffic excitations. Denoted bypil(t)andpir(t), respectively, are the interaction forces of theith moving vehicle at its left and right wheel, and through the dynamic balance analysis, they can be expressed using the additional elongation of the spring and the gravity of the vehicle,
(6)pil(t)=kiΔil(t)+ciΔ˙il(t)+migeireil+e1r,pir(t)=kiΔir(t)+ciΔ˙ir(t)+migeileil+eir.The equation of motion for the bridge can be formulated as follows:(7)Mu¨b(t)+Cu˙b(t)+Kub(t)=∑i=1nmbi(t)pi(t),
whereM,C, andKare, respectively, the mass, stiffness, and damping matrices of the bridge,ubrepresents the displacement vector, and pi(t)=[pil(t)pir(t)]T,bi(t)=[bil(t)bir(t)],bil(t),  bir(t), are respectively, the time-varying load distribution vector in global coordinates with regard to the left and right wheels. Thus,bil(t),  bir(t)can be formulated using the shape function of the beam element in finite element (FE) analysis, and their elements are zeroed when the corresponding wheel is off the bridge during the sampling time.


## 2.3. Dynamic Response of the Coupled Vehicle-Bridge System

In the equation of motion for the moving vehicle and bridge ((2) and (7)), it can be seen that the system is coupled: the vibration of the bridge influenced the vibration of the moving vehicle (2), and the vehicle in return interactively influenced the vibration of the bridge (7). The equation of motion for the coupled system can be expressed as follows:
(8)[Mm1v⋱mnmv][u¨b(t)Y¨1v(t)⋮Y¨nmv(t)]+[K+ΔK(t)-b1(t)S1Tk1⋯-bnm(t)SnmTkd-S1b1T(t)k1k1e1k⋮⋱-SnmbnmT(t)knmknmenmk]×[ub(t)Y1v(t)⋮Ynmv(t)]+[C+ΔC(t)-b1(t)S1Tc1⋯-bnm(t)SnmTcd-S1b1T(t)c1c1e1k⋮⋱-SnmbnmT(t)cnmcnmenmk]×[u˙b(t)Y˙1v(t)⋮Y˙nmv(t)]=[∑ibi(mieig-kiri+cir˙i)k1S1r1+c1S1r˙1⋮knmSnmrnm+cnmSnmr˙nm],
whereΔK(t)=∑ikibi(t)biT(t),ΔC(t)=∑icibi(t)biT(t), ei=[eil/(eil+eir)eir/(eil+eir)]T.Equation (8) shows that the system is time variant with the vehicle’s motion, and this indicated that the vibration of the bridge caused by the traffic was nonstationary. Therefore, the performance of damage identification in this case may have similar limitations to those pertaining to methods subjected to a stationary assumption, for example, modal based methods in the frequency domain. In contrast, time-domain analysis for identification is more accurate and appropriate. Then, the key procedure is the estimation of the system response to damage. However, if the response is computed in the traditional way by using the system motion equation, it will need to reassemble the varying system parameter matrix at each sampling time-step, which will be time-consuming, especially for bridges with complex structures. Aiming at solving this problem, the following section presents a fast reanalysis, VDM-based method.


## 3. Fast Reanalysis of the Damaged System

### 3.1. Dynamic Response of the Damaged Structure under Known Excitation

Denoted byμithe damage extent of theith finite element, which is quantified by the ratio between its original stiffness matrixKiand the damaged stiffness matrixK~i,
(9)K~i=μiKi.Then, the stiffness matrix of the damaged bridge can be expressed by(10)K~=∑iμiKi.Substitute (10) into the equation of motion of the damaged beam to obtain
(11)Mu¨b(t)+Cu˙b(t)+∑iμiKiub(t)=∑i=1nmbi(t)pi(t).Compared with the intact structure, remove the modified part of the stiffness matrix on the right-hand side of (11); then the equation is transformed into the equation of motion of the distorted structure(in the terminology of VDM); that is, the intact structure was subjected to the same external excitation and to a certain response coupled pseudoloadp0(t),
(12)Mu¨b(t)+Cu˙b(t)+Kub(t)=∑i=1nmbi(t)pi(t)+p0(t),
where the pseudo-loadp0(t)is related to the extent of the damage by
(13)p0(t)=∑i(1-μi)Kiu(t).In VDM, virtual distortion is more often adopted to simulate stiffness-degrading damage, which is related to the pseudo-load as follows:(14)pe,i0(t)=Ke,i∑jκij0(t)φij,
wherepe,i0is the local pseudo-load applied on theith element expressed in the local DOFs,Ke,iis the stiffness matrix of theith element in its local coordinate,κij0(t)φijis thejth virtual distortion of the  ith element,φijis thejth nonzero eigenvector to thejth nonzero positive eigenvalueλijof the stiffness matrixKe,i=∑jλijφijφijT, andφijrepresents thejth local unit distortion,κij0(t)=(1-μi)φijTue,i(t)is the time-varying combination coefficient of the correspondingjth local unit distortionφij.ue,i(t)is the nodal displacement of theith element in its local coordinate system, which shows thatκij0(t)is coupled with the system vibration.Denoted bype,iis the local nodal load applied on theith element; then it can be similarly expressed by
(15)pe,i(t)=Ke,i∑jκij(t)φij,
whereκij(t)φijis thejth actual distortion of theith damaged element and the combination coefficientκij(t)=φijTue,i(t)has the relation to the coefficientκij0(t)as follows [23]:
(16)κij0(t)=(1-μi)κij(t).Based on (12), by virtual distortions, the response of the damaged bridge under moving vehicles can be expressed by
(17)yα(t)=∑i∫0tDαim(t,τ)pi(τ)dτ+∑i,j∫0tDαijκ(t-τ)κij0(τ)dτ.Here, it assumed zero value initial conditions andyα(t)was the response of theαth measurement, which is modelled as the sum of the linear responses of the intact structure to a known excitationpi(t)and to certain virtual distortions.Dαim(t,τ)denoted the impulse response of the intact structure at time t to unit impulsive excitation applied at timeτat the respective locations of the excitationpi(t). It is worth noticing that the excitation location is changing with the vehicles.Dαijκ(t)denotes the impulse response (in the scope of the VDM, it is called the dynamic influence matrix) to an impulse unit distortionφijof theith element which is equivalent to a local impulsive loadKe,iφij.Similarly to (17), the distortion responseκij(t)of the damaged structure can be formulated as
(18)κij(t)=∑k∫0tDijkκm(t,τ)pk(τ)dτ+∑k,l∫0tDijklκκ(t-τ)κkl0(τ)dτ,
whereDijkκm(t,τ)is thejth distortion of theith damaged element at timetto unit impulsive excitation applied at timeτat the respective location of the excitationpk(t), andDijklκκ(t)is thejth distortion of theith damaged element to an impulse unit distortionφklof thekth damaged element.Substitute (16) into (18) to obtain
(19)(1-μi)∑k∫0tDijkκm(t,τ)pk(τ)dτ=κij0(t)-(1-μi)∑k,l∫0tDijklκκ(t-τ)κkl0(τ)dτ.If combined, all the damaged elementsiand the distortionsj form a system of Volterra integral equations of the second kind, which is always well posed and thus uniquely solvable (Kress 1989 [27]). Then, the responseyα(t)can be computed by (17) with the obtained virtual distortions.


### 3.2. Response of the Damaged Structure under Moving Vehicles

In practice, the interaction between the moving vehicle and the bridge, (see (6)), is time varying and coupled with the system vibration, which is usually hard to measure in advance for monitoring and identification of the extent of damage in a bridge. While the vehicle parameters, such as mass, stiffness, and damping, can be obtained more easily through traffic count statistics on the bridge, in this case, the gain of the moving excitation, besides the virtual distortion, is the key procedure, by which the response of the damaged structure in (17) is computed. Equation (6) shows that the interaction was related to the vibration of the moving vehicle; that is,uiv(t), θiv(t), and the vertical response of the bridgexi(t)at the contact location of the vehicle wheels, as well as their corresponding vibration velocity: the solutions are discussed in the following paragraphs.For the responsexi(t), it can be formulated in a manner similar to
(20)xi(t)=∑β∫0tDiβmm(t,τ)pβ(τ)dτ+∑j,k∫0tDijkmκ(t-τ)κjk0(τ)dτ,
whereβ=2l+lw-2is the number of the contact point of thelwth wheel(lw=1,2)of thelth moving vehicle.Diβmm(t,τ)is the impulse response at the location of contact pointβto unit moving impulsive excitation applied at timeτat the location ofpβ(τ)and to an impulse unit distortionφjkof thejth damaged element. It is worth noticing that position ofxi(t)changes with vehicle move. MatrixDiβmm consisting ofDiβmm(t,τ)is thus called moving dynamic influence matrix.Dijkmκ(t)is the impulse response corresponding to an impulse unit distortionφjkof thejth damaged element.Similarly to (20), the velocity response of the bridgex˙i(t)can be expressed as follows:
(21)x˙i(t)=∑β∫0tD˙iβmm(t,τ)pβ(τ)dτ+∑j,k∫0tD˙ijkmκ(t-τ)κjk0(τ)dτ,
whereD˙iβmm(t,τ)andD˙ijkmκ(t)are the respective velocity impulse responses.Using the equation of motion of the moving vehicle (2), the dynamic response of the moving vehicle can be expressed by the convolution integral between the impulse response and the excitation (the right-hand side of (2))
(22)Yiv(t)=∫0tHiv(t-τ)si(ki(xi(τ)+ri(τ))+ci(x˙i(τ)+r˙i(τ)))dτ,
whereHiv(t)=[huu,i(t)huθ,i(t)hθu,i(t)hθθ,i(t)]consists of the impulse response with respect to each DOF of the vehicle. For example,huθ,i(t)is the impulse response along the vertical direction to unit impulsive excitation along the rotational direction: the velocity of the vehicle vibrationYiv(t)is formulated as
(23)Y˙iv(t)=∫0tH˙iv(t-τ)si×(ki(xi(τ)+ri(τ))+ci(x˙i(τ)+r˙i(τ)))dτ,
where matrixH˙iv(t)consists of the corresponding velocity impulse response.Then, given damage extents and the moving vehicle parameters, the moving excitationpi(t), and virtual distortionsκij0(t), can be computed by solving the combined equations (6), (19), (20), (21), (22), and (23), where the relevant impulse matrices of the intact bridge and of the vehicles are required. Then, substitute excitationpi(t)and virtual distortionsκij0(t)into (17); the responses of the damaged structure can be rapidly computed.


### 3.3. Discretization

In applications, the responses are usually measured or obtained by numerical simulation and are thus discrete. Collect the dynamic responsesyα(t)for all the time steps and arrange them in a certain sequence, and the discrete matrix form of (17) becomes
(24)y=Dmp+Dκκ0,
where vectorsy,  p,  κ0, respectively, consist of the discrete responses (of all the considered sensors), discrete excitations, discrete virtual distortions (of all the damaged elements) for all time steps, and thus they are of respective lengthsnant, 2nmnt, and ndnt, wherentdenotes the measured time step, andnadenotes the number of sensors, andnd denotes the considered number of virtual distortions. MatricesDmandDκare block matrices and, respectively, consist of the corresponding discrete integral operators in (17) with lower-triangularnt×ntblocks. Thereinto, the submatrix ofDκis Toeplitz matrix.Similarly, the discrete responsex,x˙and the discrete distortionsκdepend on the discrete excitationpand discrete virtual distortionsκ0in the following way:
(25)[xx˙κ]=[DmmDmκD˙mmD˙mκDκmDκκ][pκ0],
which is an aggregated discrete version of (20), (21), and (18).The discrete version of vehicle vibration in (22) and (23) is as follows:
(26)Yv=HvS(k(x+r)+c(x˙+r˙)),Y˙v=H˙vS(k(x+r)+c(x˙+r˙)),
whereHv=∑i=1nmIi⊗Hiv, S=∑i=1nmIi⊗S-i, S-i=Si⊗Int×nt,k=∑i=1nmkiI^i,c=∑i=1nmciI^i,I^i=Ii⊗I2nt,2nt,I2nt,2ntis a unit matrix with dimension2nt.Iiis a square matrix with dimensionnm, of which only elementIi(i,i)has value 1, and all other elements have zero values.Equation (27) shows the expression of discrete excitationp:
(27)p=kSTYv+cSTY˙v-k(x+r)-c(x˙+r˙)+me-g,
wherem=∑i=1nmmiI^i, andgis the vector of Earth’s gravitiesg.e-=[e1T,…enmT]T⊗1nt×1, ei=[eir/(eir+eil)eil/(eir+eil)]T.Further, the discrete version of (16) is as follows:
(28)κ0=(I-μ)κ,
whereμis a block diagonal matrix of respective dimensions with diagonal blocksμiI2nt,2nt.Finally, the previous discrete equations (25), (26), and (27) are combined and generate the following large linear system:
(29)[I-[0HcvSΘ10DallmmΘ2-DallmmΘ1Dallmκ(I-μ)DκmΘ2-DκmΘ1Dκκ(I-μ)]]×[Yallxallκ]=[0DallmmDκm]me-g+[HcvS-Dallmm-Dκm]Θ1ralldy,
whereHcv=[(Hv)T(H˙v)T]T,Θ1=[kc],Yallv=[(Yv)T(Y˙v)T]T,Dallmm=[(Dmm)T(D˙mm)T]T,Dallmκ=[(Dmκ)T(D˙mκ)T],Θ2=[kSTcST],xall=[xTx˙T]T,rall=[rTr˙T]T.The block matrices in (29) store all the necessary information about the dynamics of the systems and are independent of the damage. Thus, given the values of the moving vehicle parametersmi, vi, ciand the extent of the damageμi, the system (29) can be quickly assembled and then solved to obtain the moving excitations and virtual distortions, which can then be used in (24) to compute the response of the damaged structure to the moving vehicles.However, in application there are different vehicles passing over any given bridge, and it was necessary to set up the corresponding impulse frequency matrix of the vehicle(s),Hcvin (29), repeatedly. Based on the reanalysis concept of VDM, the problem can be improved. Denotem~i, J~i, k~i, c~ias the initial vehicle system parameters of theith moving vehicle, and denoteμim, μiJ, μik, μic, respectively, as the ratios between the actual, and initial, values; that is,
(30)mi=μimm~i,Ji=μiJJ~i,ki=μikk~i,ci=μicc~i.Substitute (30) into (2), and shift the modified counterpart compared to the initial vehicle value to obtain
(31)m~ivY¨iv(t)+c~ieikY˙iv(t)+k~ieikYiv(t)=μikk~iSi(xi(t)+ri(t))+μicc~iSi(x˙i(t)+r˙i(t))+(I-μim,J)m~ivY¨iv(t)+(1-μic)c~iekiY˙iv(t)+(1-μik)k~iekiYiv(t),
whereμim,J=diag[μimμiJ].Then, through the use of the impulse response matrix of the initial vehicle systemH~v, the vibration of the vehicle can be reformulated as follows:
(32)Yv=H~vS(μkk~(x+r)+μcc~(x˙+r˙))+H~v(I-μm,J)m~m,JY¨v+H~vSe((I-μk)k~Yv+(I-μc)c~Y˙v),
where matricesk~,c~,m~have a similar form to matricesk,c,m, as didμk,μc, andμm.m~m,J=∑i=1nmIi⊗m-iv, m-iv=m~iv⊗Int,nt, Se=∑i=1nmIi⊗Sei, Sei=eik⊗Intnt, μm,J=∑i=1nmIi⊗μ-im,J.Meanwhile, ifH~v, in (32), is replaced byH~˙v, it yields the expression forY˙vand using the reformulated expression for the vehicle vibrationYvandY˙v, (29) becomes
(33)[I-[H~allvI2Sem(I-μkcm)H~cvSΘ1μkc0DallmmΘ2μkcm-DallmmΘ1μkcDallmκ(I-μ)DκmΘ2μkcm-DκmΘ1μkcDκκ(I-μ)]]×[Yallvxallκ]=[0DallmmDκm]μmm~e-g+[H~cvS-Dallmm-Dκm]Θ1μkcrall,
whereH~allmm=diag[H~vH~˙v],  Sem=diag[Sek~Sec~m~m,J],  μkcm=diag[μkμcμm,J],μkc=diag[μkμc],I2=one(2)⊗I2nmnt.Equation (33), like (29), is a large linear system with full rank, and its solution is well conditioned. There, the block matrices, such as the related influence matrix or initial system matrices, are generated by the FE model, which is noiseless and only needs to compute once for each set of moving vehicle and damage parameters. Therefore, it provides a fast reanalysis approach to the estimation of the dynamic responses and the moving excitation, which makes the optimisation of the structural damage identification and vehicle parameters feasible.


## 3.1. Dynamic Response of the Damaged Structure under Known Excitation

Denoted byμithe damage extent of theith finite element, which is quantified by the ratio between its original stiffness matrixKiand the damaged stiffness matrixK~i,
(9)K~i=μiKi.Then, the stiffness matrix of the damaged bridge can be expressed by(10)K~=∑iμiKi.Substitute (10) into the equation of motion of the damaged beam to obtain
(11)Mu¨b(t)+Cu˙b(t)+∑iμiKiub(t)=∑i=1nmbi(t)pi(t).Compared with the intact structure, remove the modified part of the stiffness matrix on the right-hand side of (11); then the equation is transformed into the equation of motion of the distorted structure(in the terminology of VDM); that is, the intact structure was subjected to the same external excitation and to a certain response coupled pseudoloadp0(t),
(12)Mu¨b(t)+Cu˙b(t)+Kub(t)=∑i=1nmbi(t)pi(t)+p0(t),
where the pseudo-loadp0(t)is related to the extent of the damage by
(13)p0(t)=∑i(1-μi)Kiu(t).In VDM, virtual distortion is more often adopted to simulate stiffness-degrading damage, which is related to the pseudo-load as follows:(14)pe,i0(t)=Ke,i∑jκij0(t)φij,
wherepe,i0is the local pseudo-load applied on theith element expressed in the local DOFs,Ke,iis the stiffness matrix of theith element in its local coordinate,κij0(t)φijis thejth virtual distortion of the  ith element,φijis thejth nonzero eigenvector to thejth nonzero positive eigenvalueλijof the stiffness matrixKe,i=∑jλijφijφijT, andφijrepresents thejth local unit distortion,κij0(t)=(1-μi)φijTue,i(t)is the time-varying combination coefficient of the correspondingjth local unit distortionφij.ue,i(t)is the nodal displacement of theith element in its local coordinate system, which shows thatκij0(t)is coupled with the system vibration.Denoted bype,iis the local nodal load applied on theith element; then it can be similarly expressed by
(15)pe,i(t)=Ke,i∑jκij(t)φij,
whereκij(t)φijis thejth actual distortion of theith damaged element and the combination coefficientκij(t)=φijTue,i(t)has the relation to the coefficientκij0(t)as follows [23]:
(16)κij0(t)=(1-μi)κij(t).Based on (12), by virtual distortions, the response of the damaged bridge under moving vehicles can be expressed by
(17)yα(t)=∑i∫0tDαim(t,τ)pi(τ)dτ+∑i,j∫0tDαijκ(t-τ)κij0(τ)dτ.Here, it assumed zero value initial conditions andyα(t)was the response of theαth measurement, which is modelled as the sum of the linear responses of the intact structure to a known excitationpi(t)and to certain virtual distortions.Dαim(t,τ)denoted the impulse response of the intact structure at time t to unit impulsive excitation applied at timeτat the respective locations of the excitationpi(t). It is worth noticing that the excitation location is changing with the vehicles.Dαijκ(t)denotes the impulse response (in the scope of the VDM, it is called the dynamic influence matrix) to an impulse unit distortionφijof theith element which is equivalent to a local impulsive loadKe,iφij.Similarly to (17), the distortion responseκij(t)of the damaged structure can be formulated as
(18)κij(t)=∑k∫0tDijkκm(t,τ)pk(τ)dτ+∑k,l∫0tDijklκκ(t-τ)κkl0(τ)dτ,
whereDijkκm(t,τ)is thejth distortion of theith damaged element at timetto unit impulsive excitation applied at timeτat the respective location of the excitationpk(t), andDijklκκ(t)is thejth distortion of theith damaged element to an impulse unit distortionφklof thekth damaged element.Substitute (16) into (18) to obtain
(19)(1-μi)∑k∫0tDijkκm(t,τ)pk(τ)dτ=κij0(t)-(1-μi)∑k,l∫0tDijklκκ(t-τ)κkl0(τ)dτ.If combined, all the damaged elementsiand the distortionsj form a system of Volterra integral equations of the second kind, which is always well posed and thus uniquely solvable (Kress 1989 [27]). Then, the responseyα(t)can be computed by (17) with the obtained virtual distortions.


## 3.2. Response of the Damaged Structure under Moving Vehicles

In practice, the interaction between the moving vehicle and the bridge, (see (6)), is time varying and coupled with the system vibration, which is usually hard to measure in advance for monitoring and identification of the extent of damage in a bridge. While the vehicle parameters, such as mass, stiffness, and damping, can be obtained more easily through traffic count statistics on the bridge, in this case, the gain of the moving excitation, besides the virtual distortion, is the key procedure, by which the response of the damaged structure in (17) is computed. Equation (6) shows that the interaction was related to the vibration of the moving vehicle; that is,uiv(t), θiv(t), and the vertical response of the bridgexi(t)at the contact location of the vehicle wheels, as well as their corresponding vibration velocity: the solutions are discussed in the following paragraphs.For the responsexi(t), it can be formulated in a manner similar to
(20)xi(t)=∑β∫0tDiβmm(t,τ)pβ(τ)dτ+∑j,k∫0tDijkmκ(t-τ)κjk0(τ)dτ,
whereβ=2l+lw-2is the number of the contact point of thelwth wheel(lw=1,2)of thelth moving vehicle.Diβmm(t,τ)is the impulse response at the location of contact pointβto unit moving impulsive excitation applied at timeτat the location ofpβ(τ)and to an impulse unit distortionφjkof thejth damaged element. It is worth noticing that position ofxi(t)changes with vehicle move. MatrixDiβmm consisting ofDiβmm(t,τ)is thus called moving dynamic influence matrix.Dijkmκ(t)is the impulse response corresponding to an impulse unit distortionφjkof thejth damaged element.Similarly to (20), the velocity response of the bridgex˙i(t)can be expressed as follows:
(21)x˙i(t)=∑β∫0tD˙iβmm(t,τ)pβ(τ)dτ+∑j,k∫0tD˙ijkmκ(t-τ)κjk0(τ)dτ,
whereD˙iβmm(t,τ)andD˙ijkmκ(t)are the respective velocity impulse responses.Using the equation of motion of the moving vehicle (2), the dynamic response of the moving vehicle can be expressed by the convolution integral between the impulse response and the excitation (the right-hand side of (2))
(22)Yiv(t)=∫0tHiv(t-τ)si(ki(xi(τ)+ri(τ))+ci(x˙i(τ)+r˙i(τ)))dτ,
whereHiv(t)=[huu,i(t)huθ,i(t)hθu,i(t)hθθ,i(t)]consists of the impulse response with respect to each DOF of the vehicle. For example,huθ,i(t)is the impulse response along the vertical direction to unit impulsive excitation along the rotational direction: the velocity of the vehicle vibrationYiv(t)is formulated as
(23)Y˙iv(t)=∫0tH˙iv(t-τ)si×(ki(xi(τ)+ri(τ))+ci(x˙i(τ)+r˙i(τ)))dτ,
where matrixH˙iv(t)consists of the corresponding velocity impulse response.Then, given damage extents and the moving vehicle parameters, the moving excitationpi(t), and virtual distortionsκij0(t), can be computed by solving the combined equations (6), (19), (20), (21), (22), and (23), where the relevant impulse matrices of the intact bridge and of the vehicles are required. Then, substitute excitationpi(t)and virtual distortionsκij0(t)into (17); the responses of the damaged structure can be rapidly computed.


## 3.3. Discretization

In applications, the responses are usually measured or obtained by numerical simulation and are thus discrete. Collect the dynamic responsesyα(t)for all the time steps and arrange them in a certain sequence, and the discrete matrix form of (17) becomes
(24)y=Dmp+Dκκ0,
where vectorsy,  p,  κ0, respectively, consist of the discrete responses (of all the considered sensors), discrete excitations, discrete virtual distortions (of all the damaged elements) for all time steps, and thus they are of respective lengthsnant, 2nmnt, and ndnt, wherentdenotes the measured time step, andnadenotes the number of sensors, andnd denotes the considered number of virtual distortions. MatricesDmandDκare block matrices and, respectively, consist of the corresponding discrete integral operators in (17) with lower-triangularnt×ntblocks. Thereinto, the submatrix ofDκis Toeplitz matrix.Similarly, the discrete responsex,x˙and the discrete distortionsκdepend on the discrete excitationpand discrete virtual distortionsκ0in the following way:
(25)[xx˙κ]=[DmmDmκD˙mmD˙mκDκmDκκ][pκ0],
which is an aggregated discrete version of (20), (21), and (18).The discrete version of vehicle vibration in (22) and (23) is as follows:
(26)Yv=HvS(k(x+r)+c(x˙+r˙)),Y˙v=H˙vS(k(x+r)+c(x˙+r˙)),
whereHv=∑i=1nmIi⊗Hiv, S=∑i=1nmIi⊗S-i, S-i=Si⊗Int×nt,k=∑i=1nmkiI^i,c=∑i=1nmciI^i,I^i=Ii⊗I2nt,2nt,I2nt,2ntis a unit matrix with dimension2nt.Iiis a square matrix with dimensionnm, of which only elementIi(i,i)has value 1, and all other elements have zero values.Equation (27) shows the expression of discrete excitationp:
(27)p=kSTYv+cSTY˙v-k(x+r)-c(x˙+r˙)+me-g,
wherem=∑i=1nmmiI^i, andgis the vector of Earth’s gravitiesg.e-=[e1T,…enmT]T⊗1nt×1, ei=[eir/(eir+eil)eil/(eir+eil)]T.Further, the discrete version of (16) is as follows:
(28)κ0=(I-μ)κ,
whereμis a block diagonal matrix of respective dimensions with diagonal blocksμiI2nt,2nt.Finally, the previous discrete equations (25), (26), and (27) are combined and generate the following large linear system:
(29)[I-[0HcvSΘ10DallmmΘ2-DallmmΘ1Dallmκ(I-μ)DκmΘ2-DκmΘ1Dκκ(I-μ)]]×[Yallxallκ]=[0DallmmDκm]me-g+[HcvS-Dallmm-Dκm]Θ1ralldy,
whereHcv=[(Hv)T(H˙v)T]T,Θ1=[kc],Yallv=[(Yv)T(Y˙v)T]T,Dallmm=[(Dmm)T(D˙mm)T]T,Dallmκ=[(Dmκ)T(D˙mκ)T],Θ2=[kSTcST],xall=[xTx˙T]T,rall=[rTr˙T]T.The block matrices in (29) store all the necessary information about the dynamics of the systems and are independent of the damage. Thus, given the values of the moving vehicle parametersmi, vi, ciand the extent of the damageμi, the system (29) can be quickly assembled and then solved to obtain the moving excitations and virtual distortions, which can then be used in (24) to compute the response of the damaged structure to the moving vehicles.However, in application there are different vehicles passing over any given bridge, and it was necessary to set up the corresponding impulse frequency matrix of the vehicle(s),Hcvin (29), repeatedly. Based on the reanalysis concept of VDM, the problem can be improved. Denotem~i, J~i, k~i, c~ias the initial vehicle system parameters of theith moving vehicle, and denoteμim, μiJ, μik, μic, respectively, as the ratios between the actual, and initial, values; that is,
(30)mi=μimm~i,Ji=μiJJ~i,ki=μikk~i,ci=μicc~i.Substitute (30) into (2), and shift the modified counterpart compared to the initial vehicle value to obtain
(31)m~ivY¨iv(t)+c~ieikY˙iv(t)+k~ieikYiv(t)=μikk~iSi(xi(t)+ri(t))+μicc~iSi(x˙i(t)+r˙i(t))+(I-μim,J)m~ivY¨iv(t)+(1-μic)c~iekiY˙iv(t)+(1-μik)k~iekiYiv(t),
whereμim,J=diag[μimμiJ].Then, through the use of the impulse response matrix of the initial vehicle systemH~v, the vibration of the vehicle can be reformulated as follows:
(32)Yv=H~vS(μkk~(x+r)+μcc~(x˙+r˙))+H~v(I-μm,J)m~m,JY¨v+H~vSe((I-μk)k~Yv+(I-μc)c~Y˙v),
where matricesk~,c~,m~have a similar form to matricesk,c,m, as didμk,μc, andμm.m~m,J=∑i=1nmIi⊗m-iv, m-iv=m~iv⊗Int,nt, Se=∑i=1nmIi⊗Sei, Sei=eik⊗Intnt, μm,J=∑i=1nmIi⊗μ-im,J.Meanwhile, ifH~v, in (32), is replaced byH~˙v, it yields the expression forY˙vand using the reformulated expression for the vehicle vibrationYvandY˙v, (29) becomes
(33)[I-[H~allvI2Sem(I-μkcm)H~cvSΘ1μkc0DallmmΘ2μkcm-DallmmΘ1μkcDallmκ(I-μ)DκmΘ2μkcm-DκmΘ1μkcDκκ(I-μ)]]×[Yallvxallκ]=[0DallmmDκm]μmm~e-g+[H~cvS-Dallmm-Dκm]Θ1μkcrall,
whereH~allmm=diag[H~vH~˙v],  Sem=diag[Sek~Sec~m~m,J],  μkcm=diag[μkμcμm,J],μkc=diag[μkμc],I2=one(2)⊗I2nmnt.Equation (33), like (29), is a large linear system with full rank, and its solution is well conditioned. There, the block matrices, such as the related influence matrix or initial system matrices, are generated by the FE model, which is noiseless and only needs to compute once for each set of moving vehicle and damage parameters. Therefore, it provides a fast reanalysis approach to the estimation of the dynamic responses and the moving excitation, which makes the optimisation of the structural damage identification and vehicle parameters feasible.


## 4. Identification of Moving Vehicles and Damages

The analysis in Section2.3 shows that the coupled vehicle-bridge system is time variant, and the response is nonstationary. Time-domain analysis is the proper identification of unknown structural systems and moving vehicles.


### 4.1. Objective Function and the Optimization Variables

Equation (24) shows that the dynamic response of the damaged bridge was a function of the moving excitation and virtual distortions, which were in a one-to-one correspondence to the moving vehicle parameters and structural damage. Thus, here the inverse problem of identification of unknown vehicle and damage extent was stated as an optimisation problem involving the minimisation of the normalised mean-square distance between the measured structural responseyMand the computed response y. The optimisation variables are the modifying factors of vehicle parametersμim, μiJ, μik, μic, and damage extentμiwhich takes advantage of two aspects: on one hand, the magnitudes of the adopted optimisation variables were much less than the magnitudes between the vehicle parameters and bridge damage extents, for the latter can seriously impair the accuracy of many optimisation procedures: the damage extentsμibelong to the interval[0,1], but the vehicle parameters might be as large as tens of tonnes, and so forth. Moreover, this makes the direct adoption of (33) to obtain the moving excitation and dynamic response of the damaged bridge to given optimisation variables with the prerequisite of the initial system matrices which only are computed once.Therefore, the initial vehicle parameters were firstly computed using the measured responses of the damaged structure while assuming that the bridge was intact and was under the action of moving excitations which equal the weight forces of the vehicles with initial massm~i,  for which the corresponding responses equal the measured responses, therefore giving rise to the following expression:
(34)yM=Dmm~g,
whereyMandgare known, and matrixDmcan be constructed using the intact bridge system, matrix vectorm~consisting ofm~iis the unknown to be solved. For the number of the required sensors is bigger than the number of vehicles, (34) is an overdetermined system and the initial massm~ican be computed via the least-square method.Then, the rest of the initial parameters for the vehicles are determined by their empirical relationships to the vehicles’ masses. Given the initial parameters, the optimisation variables can be represented by the following dimensionless variablesμi⋆(i=1,…,4nm+ne), neis the number of the potential damaged elements. Denotei=4(j-1)+jm, there is
(35)μi⋆={μjmifj≤nm,jm=1μjJifj≤nm,jm=2μjkifj≤nm,jm=3μjcifj≤nm,jm=4μi-4nmifj>nm.All the optimization variablesμi⋆are of the same magnitude and have the natural initial value of 1 for the optimization. The objective function is thus built as
(36)minimizef(μ1⋆,…,μ4nm+ne⋆)=12∥yM-y∥2∥yM∥2subjecttoμi⋆≥0,i=1,…,4nm+ne,
whereyis the computed response of the structure to the given optimization variables by (35), (33), (27), and (28).


### 4.2. Sensitivity Analysis

For the optimization of the objective function (36), it can be performed quickly using the gradient-based algorithms provided that the gradient can be computed at a reasonable cost. The formulation based on (24), (27), (28), and (33) allows the discrete adjoint method to be used, which is quicker by one order of magnitude compared with the standard differentiation method [28].DenoteD=[DmmDκ],Λ=[PT(κ0)T]T; then the first derivative of (36) with respect to the variableμi⋆is stated as
(37)fi=-(yM-y)T∥yM∥2DΛi.Here and below, the subscriptiin each variable stands for its first derivative to the variableμi⋆.In order to obtain the expression ofΛi, (27) is firstly reformulated as follows:
(38)p=Γ1z1+w,
whereΓ1=[kSTcST-k-c],  z1=[(Yv)T(Y˙v)TxTx˙T]T, w=-kr-cr˙+me-g.Then(39)pi=Γ1z1i+Γ1iz1+wi.Equation (28) is differentiated to obtain
(40)(κ0)i=Γ2κi+Γ2iκ,
whereΓ2=I-μ.The combination of (39) and (40) generates the following expression:
(41)Λi=Γzi+Ωi,
wherez=[(Yallv)TxallTκT]T,Γ=diag[Γ1Γ2],Ωi=[Γ1iz1][(Γ1iz1+wi)T(Γ2iκ)T]T.Vectorz  represents the parameters in (33) required to be solved. By reformulating (33) as the formAz=b, the first derivativezisatisfies the following relationship:
(42)Azi=bi-Aiz.Substitute (41) into (37); there is
(43)fi=-(yM-y)T∥yM∥2DΓzi-(yM-y)T∥yM∥2DΩi
which includes the first derivativeszi. The discrete adjoint method adds the scalar product of the adjoint vectorλwith (42) to (43), to obtain
(44)fi=λT(bi-Aiz)-(yM-y)T∥yM∥2DΩi+(λTA-(yM-y)T∥yM∥2DΓ)zi.In this way, the gradient of the objective function can be expressed as(45)fi=λT(bi-Aiz)-(yM-y)T∥yM∥2DΩi,
where the adjoint vectorλis computed only at the cost of a single solution of the adjoint equation
(46)ATλ=ΓTDT(yM-y)∥yM∥2.


### 4.3. Remarks

If a small number of time steps are used, then the system matrix in (33) can be computed and used directly. However, for offline identification, in the case of a dense time discretization or a longer sampling time, the system can become prohibitively large and computationally hardly manageable. To reduce the numerical costs, the system matrix, which is a block matrix composed of lower triangular matrices, can be rearranged into a lower triangular block form; then the system can be solved by a specialized linear solver, like block forward substitution [29]. Due to the fact that all the data in the system (33) are computed based on the ideal FEM model and are away from any measurement noise pollution, the unknown vectors can be computed stepwise precisely.If the bridge is undamaged, the virtual distortions in (24) vanish; that is, the value of damage extent equals one in related Equations, and the proposed method can be also used for identification of moving vehicles only.


## 4.1. Objective Function and the Optimization Variables

Equation (24) shows that the dynamic response of the damaged bridge was a function of the moving excitation and virtual distortions, which were in a one-to-one correspondence to the moving vehicle parameters and structural damage. Thus, here the inverse problem of identification of unknown vehicle and damage extent was stated as an optimisation problem involving the minimisation of the normalised mean-square distance between the measured structural responseyMand the computed response y. The optimisation variables are the modifying factors of vehicle parametersμim, μiJ, μik, μic, and damage extentμiwhich takes advantage of two aspects: on one hand, the magnitudes of the adopted optimisation variables were much less than the magnitudes between the vehicle parameters and bridge damage extents, for the latter can seriously impair the accuracy of many optimisation procedures: the damage extentsμibelong to the interval[0,1], but the vehicle parameters might be as large as tens of tonnes, and so forth. Moreover, this makes the direct adoption of (33) to obtain the moving excitation and dynamic response of the damaged bridge to given optimisation variables with the prerequisite of the initial system matrices which only are computed once.Therefore, the initial vehicle parameters were firstly computed using the measured responses of the damaged structure while assuming that the bridge was intact and was under the action of moving excitations which equal the weight forces of the vehicles with initial massm~i,  for which the corresponding responses equal the measured responses, therefore giving rise to the following expression:
(34)yM=Dmm~g,
whereyMandgare known, and matrixDmcan be constructed using the intact bridge system, matrix vectorm~consisting ofm~iis the unknown to be solved. For the number of the required sensors is bigger than the number of vehicles, (34) is an overdetermined system and the initial massm~ican be computed via the least-square method.Then, the rest of the initial parameters for the vehicles are determined by their empirical relationships to the vehicles’ masses. Given the initial parameters, the optimisation variables can be represented by the following dimensionless variablesμi⋆(i=1,…,4nm+ne), neis the number of the potential damaged elements. Denotei=4(j-1)+jm, there is
(35)μi⋆={μjmifj≤nm,jm=1μjJifj≤nm,jm=2μjkifj≤nm,jm=3μjcifj≤nm,jm=4μi-4nmifj>nm.All the optimization variablesμi⋆are of the same magnitude and have the natural initial value of 1 for the optimization. The objective function is thus built as
(36)minimizef(μ1⋆,…,μ4nm+ne⋆)=12∥yM-y∥2∥yM∥2subjecttoμi⋆≥0,i=1,…,4nm+ne,
whereyis the computed response of the structure to the given optimization variables by (35), (33), (27), and (28).


## 4.2. Sensitivity Analysis

For the optimization of the objective function (36), it can be performed quickly using the gradient-based algorithms provided that the gradient can be computed at a reasonable cost. The formulation based on (24), (27), (28), and (33) allows the discrete adjoint method to be used, which is quicker by one order of magnitude compared with the standard differentiation method [28].DenoteD=[DmmDκ],Λ=[PT(κ0)T]T; then the first derivative of (36) with respect to the variableμi⋆is stated as
(37)fi=-(yM-y)T∥yM∥2DΛi.Here and below, the subscriptiin each variable stands for its first derivative to the variableμi⋆.In order to obtain the expression ofΛi, (27) is firstly reformulated as follows:
(38)p=Γ1z1+w,
whereΓ1=[kSTcST-k-c],  z1=[(Yv)T(Y˙v)TxTx˙T]T, w=-kr-cr˙+me-g.Then(39)pi=Γ1z1i+Γ1iz1+wi.Equation (28) is differentiated to obtain
(40)(κ0)i=Γ2κi+Γ2iκ,
whereΓ2=I-μ.The combination of (39) and (40) generates the following expression:
(41)Λi=Γzi+Ωi,
wherez=[(Yallv)TxallTκT]T,Γ=diag[Γ1Γ2],Ωi=[Γ1iz1][(Γ1iz1+wi)T(Γ2iκ)T]T.Vectorz  represents the parameters in (33) required to be solved. By reformulating (33) as the formAz=b, the first derivativezisatisfies the following relationship:
(42)Azi=bi-Aiz.Substitute (41) into (37); there is
(43)fi=-(yM-y)T∥yM∥2DΓzi-(yM-y)T∥yM∥2DΩi
which includes the first derivativeszi. The discrete adjoint method adds the scalar product of the adjoint vectorλwith (42) to (43), to obtain
(44)fi=λT(bi-Aiz)-(yM-y)T∥yM∥2DΩi+(λTA-(yM-y)T∥yM∥2DΓ)zi.In this way, the gradient of the objective function can be expressed as(45)fi=λT(bi-Aiz)-(yM-y)T∥yM∥2DΩi,
where the adjoint vectorλis computed only at the cost of a single solution of the adjoint equation
(46)ATλ=ΓTDT(yM-y)∥yM∥2.


## 4.3. Remarks

If a small number of time steps are used, then the system matrix in (33) can be computed and used directly. However, for offline identification, in the case of a dense time discretization or a longer sampling time, the system can become prohibitively large and computationally hardly manageable. To reduce the numerical costs, the system matrix, which is a block matrix composed of lower triangular matrices, can be rearranged into a lower triangular block form; then the system can be solved by a specialized linear solver, like block forward substitution [29]. Due to the fact that all the data in the system (33) are computed based on the ideal FEM model and are away from any measurement noise pollution, the unknown vectors can be computed stepwise precisely.If the bridge is undamaged, the virtual distortions in (24) vanish; that is, the value of damage extent equals one in related Equations, and the proposed method can be also used for identification of moving vehicles only.


## 5. Numerical Example

A multispan frame was used to illustrate and verify the proposed methods for simultaneous identification of moving vehicles and damage. Measurement error and model error were taken into account to test the robustness of the method.


### 5.1. Structural and Moving Vehicles

The frame, shown in Figure3, is made of steel with Young’s modulus 2.15 × 1011 Pa, and a density of 7.8 × 103 kg/m3. It has a uniform mass distribution of 15.3 × 103 kg/m3. The beam is 200 m long with a simplified rectangular cross-section ofb×h=0.89 m × 2.21 m, so that the inertial moment of the area is 0.8 m4; each of the two side spans is 50 m. Each pier is 20 m high with inertial moment of 0.16 m4.Figure 3
The vehicle-bridge coupled system.The road surface condition is considered to be good, and the roughness coefficientGd(n0) is 16 × 10−6m3. According to the vehicle velocity,nl=0.0221 m−1,nu=1.4142 m−1; then the length of the shortest effective pavement is1/nl=45.25 m, so the effective pavement can be the whole length of the bridge. In this way,Δn=1/L=0.01m−1,NT=fix((nu-nl)/Δn)= 278. Figure 4 shows the simulated road surface roughness via (3).Figure 4
The road surface roughness.Two moving vehicles passed through the bridge from opposite direction, with respective constant velocitiesv1=34 m/s,v2=-30 m/s. A two-axle, two-DOF, vehicle model was adopted. The vehicle masses werem1=61.2×103 kg, m2=53×103 kg. The other vehicle parameters such as the rotary inertia, stiffness, and damping are provided by multiplying the vehicle mass by a certain coefficient according to the literature [30]. The rest of the parameters of the two vehicles were defined, respectively, as k=396mN/m, c=1.4m N·s/m; the rotary inertiaJ=9mkg·m2. The distance between the centre of mass and either of the two axles was 3 m. Denoted byx1,0=-3 m,x2,0=204 m, respectively, were the initial centroidal points of each vehicle’s mass.Two sensors are employed:S1 at the location of 65.2 m and S2 at the location 145.2 m as shown in Figure 3. The sensors are placed at bottom surface of the beam.


### 5.2. Identification Cases

For the identification purpose, the beam is divided into 20 elements, and each of the piers is divided into two elements. Measurements are simulated with an uncorrelated Gaussian noise. Besides the measurement error, the influence of the model error is tested by using a different finite element (FE) model of the structure; that is, the theoretical FE model of the intact structure is modified as the actual intact structure.Two following identification cases are discussed.Case 1.
Measurement error is simulated at 5% RMS level. No model error is assumed. Assume that two pier elements (Nos. 21 and 23) are damaged with the damage extentsμ21=0.4andμ23=0.7.Case 2.
As in Case1, but the model error is additionally simulated. It considers the modification of the stiffness of all the elements. More precisely, uncorrelated Gaussian modifications with mean −2% and standard deviation 5% are used; see Figure 5, for the stiffness modifications. Then, the actual damage extents are slightly different from that listed in case 1, which include the model error besides the real damage.Figure 5
Stiffness reduction levels of the elements.In both of the two cases, it was assumed that the damage location was limited to the four pier elements. Two moving vehicles were identified simultaneously as being causally linked to the extent of the damage. Four stiffness modification coefficients were used in the optimisation, besides the eight variables related to the masses, rotary moments, stiffness, and damping of the two vehicles. In this way, the damage extent and its location, as well as vehicle parameters, were treated as unknowns and thus identified. Meanwhile, the dynamic interaction of the coupled vehicle-bridge system was established.


### 5.3. Moving Vehicle and the Damage Identification

Responses of the two sensors are used for vehicle and damage identification. In each case, the dynamic responses are calculated using the discrete FE model and the Newmark integration method with the parametersα=0.25andβ=0.5. The first two damping ratios are 0.01 and 0.015. The integration time step equals 0.01 s (100 Hz sampling frequency). A total of 100 time steps is used; hence, the measured time is 1 s. The simulated responses with measurement noise in Case 1 are show in Figure 6. Table 1 shows the natural frequencies of the original and damaged bridge model.Table 1
The first eight natural frequencies of the original and damaged bridge model (Hz).
ω
1
ω
2
ω
3
ω
4
ω
5
ω
6
ω
7
ω
8
Original model
5.023
5.707
14.432
17.279
20.6
35.324
53.397
59.437
Damaged model
4.552
5.526
13.958
17.054
20.56
35.036
52.782
58.606Figure 6
Simulated strain responses of the damaged and the intact structure in Case1 with measurement noise at 5% RMS level.The damage was limited to the two piers, that is, to the four pier elements. Together with two known vehicles (eight modification coefficients), there were twelve variables to be optimised by minimising the objective function (36). Responses from two sensors were used for this purpose. The initial trial mass values were estimated using (34). Then, the other initial vehicle parameters were provided based on their masses.The identification results for the two vehicle parameters are listed, respectively, in Table2 (Case 1) and Table 3 (Case 2).Table 2
Identified parameters of moving vehicles in Case1.
m
1 (kg)
J
1 (kg·m2)
k
1 (N/m)
c
1 (N·s/m)
m
2 (kg)
J
2 (kg·m2)
k
2 (N/m)
c
2 (N·s/m)
Actual value
6.12
×
10
4
5.51
×
10
5
2.42
×
10
7
8.57
×
10
4
5.30
×
10
4
4.77
×
10
5
2.10
×
10
7
7.42
×
10
4
Trial value
7.00
×
10
4
4.20
×
10
5
3.50
×
10
7
7.00
×
10
4
6.33
×
10
4
3.80
×
10
5
3.16
×
10
7
6.33
×
10
4
Identifiedμ
0.88
1.33
0.7
1.32
0.84
1.38
0.67
1.24
Identified value
6.18
×
10
4
5.58
×
10
5
2.44
×
10
7
9.25
×
10
4
5.35
×
10
4
5.23
×
10
5
2.12
×
10
7
7.84
×
10
4
Identified error (%)
0.97
1.39
0.59
8.01
0.86
9.62
0.79
5.59Table 3
Identified parameters of moving vehicles in Case2.
m
1 (kg)
J
1 (kg·m2)
k
1 (N/m)
c
1 (N·s/m)
m
2 (kg)
J
2 (kg·m2)
k
2 (N/m)
c
2 (N·s/m)
Actual value
6.12
×
10
4
5.51
×
10
5
2.42
×
10
7
8.57
×
10
4
5.30
×
10
4
4.77
×
10
5
2.10
×
10
7
7.42
×
10
4
Trial value
6.96
×
10
4
4.18
×
10
5
3.48
×
10
7
6.96
×
10
4
6.38
×
10
4
3.83
×
10
5
3.19
×
10
7
6.38
×
10
4
Identifiedμ
0.89
1.3
0.7
1.35
0.85
1.37
0.67
1.26
Identified value
6.16
×
10
4
5.44
×
10
5
2.44
×
10
7
9.41
×
10
4
5.45
×
10
4
5.23
×
10
5
2.14
×
10
7
8.02
×
10
4
Identified error (%)
0.7
1.29
0.79
9.79
2.82
9.66
2.05
8.07The identified extents of the structural damages in the two cases are listed, respectively, in Table4 (Case 1) and Table 5 (Case 2).Table 4
Identified extents of the damages in Case1.
μ
21
μ
22
μ
23
μ
24
Actual value
0.400
1.000
0.700
1.000
Identified value
0.454
0.930
0.709
0.965
Identified error (%)
13.50
7.00
1.29
3.5Table 5
Identified extents of the damages in Case2.
μ
21
μ
22
μ
23
μ
24
Actual value
0.398
1.004
0.666
0.977
Identified value
0.453
0.908
0.755
0.990
Identified error (%)
13.82
9.56
13.36
1.33The parameter identification results were assessed by their relative accuracy, while the damage identification results were more naturally assessed in terms of their absolute accuracy. It showed that, even with all the simulated errors, the parameters of the two vehicles and four potential damage identifiers could be identified acceptably. Since only two damages were actually assumed, the optimisation allowed their number and location (limited to the four considered pier elements) to be identified as well.With the identified vehicles and structural damage, the moving loads can be computed meanwhile by using (27). Figure 7 shows the estimated results which have reasonable accuracy. The relative error was 2.6% even under both 5% RMS measurement noise and model error in Case 2.Identified moving loads.
(a)
Moving loads of vehicle 1
(b)
Moving loads of vehicle 2


## 5.1. Structural and Moving Vehicles

The frame, shown in Figure3, is made of steel with Young’s modulus 2.15 × 1011 Pa, and a density of 7.8 × 103 kg/m3. It has a uniform mass distribution of 15.3 × 103 kg/m3. The beam is 200 m long with a simplified rectangular cross-section ofb×h=0.89 m × 2.21 m, so that the inertial moment of the area is 0.8 m4; each of the two side spans is 50 m. Each pier is 20 m high with inertial moment of 0.16 m4.Figure 3
The vehicle-bridge coupled system.The road surface condition is considered to be good, and the roughness coefficientGd(n0) is 16 × 10−6m3. According to the vehicle velocity,nl=0.0221 m−1,nu=1.4142 m−1; then the length of the shortest effective pavement is1/nl=45.25 m, so the effective pavement can be the whole length of the bridge. In this way,Δn=1/L=0.01m−1,NT=fix((nu-nl)/Δn)= 278. Figure 4 shows the simulated road surface roughness via (3).Figure 4
The road surface roughness.Two moving vehicles passed through the bridge from opposite direction, with respective constant velocitiesv1=34 m/s,v2=-30 m/s. A two-axle, two-DOF, vehicle model was adopted. The vehicle masses werem1=61.2×103 kg, m2=53×103 kg. The other vehicle parameters such as the rotary inertia, stiffness, and damping are provided by multiplying the vehicle mass by a certain coefficient according to the literature [30]. The rest of the parameters of the two vehicles were defined, respectively, as k=396mN/m, c=1.4m N·s/m; the rotary inertiaJ=9mkg·m2. The distance between the centre of mass and either of the two axles was 3 m. Denoted byx1,0=-3 m,x2,0=204 m, respectively, were the initial centroidal points of each vehicle’s mass.Two sensors are employed:S1 at the location of 65.2 m and S2 at the location 145.2 m as shown in Figure 3. The sensors are placed at bottom surface of the beam.


## 5.2. Identification Cases

For the identification purpose, the beam is divided into 20 elements, and each of the piers is divided into two elements. Measurements are simulated with an uncorrelated Gaussian noise. Besides the measurement error, the influence of the model error is tested by using a different finite element (FE) model of the structure; that is, the theoretical FE model of the intact structure is modified as the actual intact structure.Two following identification cases are discussed.Case 1.
Measurement error is simulated at 5% RMS level. No model error is assumed. Assume that two pier elements (Nos. 21 and 23) are damaged with the damage extentsμ21=0.4andμ23=0.7.Case 2.
As in Case1, but the model error is additionally simulated. It considers the modification of the stiffness of all the elements. More precisely, uncorrelated Gaussian modifications with mean −2% and standard deviation 5% are used; see Figure 5, for the stiffness modifications. Then, the actual damage extents are slightly different from that listed in case 1, which include the model error besides the real damage.Figure 5
Stiffness reduction levels of the elements.In both of the two cases, it was assumed that the damage location was limited to the four pier elements. Two moving vehicles were identified simultaneously as being causally linked to the extent of the damage. Four stiffness modification coefficients were used in the optimisation, besides the eight variables related to the masses, rotary moments, stiffness, and damping of the two vehicles. In this way, the damage extent and its location, as well as vehicle parameters, were treated as unknowns and thus identified. Meanwhile, the dynamic interaction of the coupled vehicle-bridge system was established.


## 5.3. Moving Vehicle and the Damage Identification

Responses of the two sensors are used for vehicle and damage identification. In each case, the dynamic responses are calculated using the discrete FE model and the Newmark integration method with the parametersα=0.25andβ=0.5. The first two damping ratios are 0.01 and 0.015. The integration time step equals 0.01 s (100 Hz sampling frequency). A total of 100 time steps is used; hence, the measured time is 1 s. The simulated responses with measurement noise in Case 1 are show in Figure 6. Table 1 shows the natural frequencies of the original and damaged bridge model.Table 1
The first eight natural frequencies of the original and damaged bridge model (Hz).
ω
1
ω
2
ω
3
ω
4
ω
5
ω
6
ω
7
ω
8
Original model
5.023
5.707
14.432
17.279
20.6
35.324
53.397
59.437
Damaged model
4.552
5.526
13.958
17.054
20.56
35.036
52.782
58.606Figure 6
Simulated strain responses of the damaged and the intact structure in Case1 with measurement noise at 5% RMS level.The damage was limited to the two piers, that is, to the four pier elements. Together with two known vehicles (eight modification coefficients), there were twelve variables to be optimised by minimising the objective function (36). Responses from two sensors were used for this purpose. The initial trial mass values were estimated using (34). Then, the other initial vehicle parameters were provided based on their masses.The identification results for the two vehicle parameters are listed, respectively, in Table2 (Case 1) and Table 3 (Case 2).Table 2
Identified parameters of moving vehicles in Case1.
m
1 (kg)
J
1 (kg·m2)
k
1 (N/m)
c
1 (N·s/m)
m
2 (kg)
J
2 (kg·m2)
k
2 (N/m)
c
2 (N·s/m)
Actual value
6.12
×
10
4
5.51
×
10
5
2.42
×
10
7
8.57
×
10
4
5.30
×
10
4
4.77
×
10
5
2.10
×
10
7
7.42
×
10
4
Trial value
7.00
×
10
4
4.20
×
10
5
3.50
×
10
7
7.00
×
10
4
6.33
×
10
4
3.80
×
10
5
3.16
×
10
7
6.33
×
10
4
Identifiedμ
0.88
1.33
0.7
1.32
0.84
1.38
0.67
1.24
Identified value
6.18
×
10
4
5.58
×
10
5
2.44
×
10
7
9.25
×
10
4
5.35
×
10
4
5.23
×
10
5
2.12
×
10
7
7.84
×
10
4
Identified error (%)
0.97
1.39
0.59
8.01
0.86
9.62
0.79
5.59Table 3
Identified parameters of moving vehicles in Case2.
m
1 (kg)
J
1 (kg·m2)
k
1 (N/m)
c
1 (N·s/m)
m
2 (kg)
J
2 (kg·m2)
k
2 (N/m)
c
2 (N·s/m)
Actual value
6.12
×
10
4
5.51
×
10
5
2.42
×
10
7
8.57
×
10
4
5.30
×
10
4
4.77
×
10
5
2.10
×
10
7
7.42
×
10
4
Trial value
6.96
×
10
4
4.18
×
10
5
3.48
×
10
7
6.96
×
10
4
6.38
×
10
4
3.83
×
10
5
3.19
×
10
7
6.38
×
10
4
Identifiedμ
0.89
1.3
0.7
1.35
0.85
1.37
0.67
1.26
Identified value
6.16
×
10
4
5.44
×
10
5
2.44
×
10
7
9.41
×
10
4
5.45
×
10
4
5.23
×
10
5
2.14
×
10
7
8.02
×
10
4
Identified error (%)
0.7
1.29
0.79
9.79
2.82
9.66
2.05
8.07The identified extents of the structural damages in the two cases are listed, respectively, in Table4 (Case 1) and Table 5 (Case 2).Table 4
Identified extents of the damages in Case1.
μ
21
μ
22
μ
23
μ
24
Actual value
0.400
1.000
0.700
1.000
Identified value
0.454
0.930
0.709
0.965
Identified error (%)
13.50
7.00
1.29
3.5Table 5
Identified extents of the damages in Case2.
μ
21
μ
22
μ
23
μ
24
Actual value
0.398
1.004
0.666
0.977
Identified value
0.453
0.908
0.755
0.990
Identified error (%)
13.82
9.56
13.36
1.33The parameter identification results were assessed by their relative accuracy, while the damage identification results were more naturally assessed in terms of their absolute accuracy. It showed that, even with all the simulated errors, the parameters of the two vehicles and four potential damage identifiers could be identified acceptably. Since only two damages were actually assumed, the optimisation allowed their number and location (limited to the four considered pier elements) to be identified as well.With the identified vehicles and structural damage, the moving loads can be computed meanwhile by using (27). Figure 7 shows the estimated results which have reasonable accuracy. The relative error was 2.6% even under both 5% RMS measurement noise and model error in Case 2.Identified moving loads.
(a)
Moving loads of vehicle 1
(b)
Moving loads of vehicle 2


## 6. Conclusion

This paper presented an effective method for the simultaneous identification of moving vehicles and structural damage based on a virtual distortion method. The vehicle parameters and the structural damage were treated as optimisation variables. Through the use of a moving dynamic influence matrix, the response of the damaged system to given optimisation variables could be computed quickly without the need of a numerical simulation and the repeated assembly of the time-variant system parameter matrix at each time-step. This reduced the numerical cost significantly. The adjoint method was adopted for the fast sensitivity analysis of the objective function. A numerical example of a three-span beam was presented to verify the proposed method, where two vehicles (eight variables) and two damaged piers (four damage instances) were identified using only two sensors. Even with both the pollution of 5% RMS measurement error and a certain model error, the maximum identification error was 13.82%. Moreover, the moving excitation effects were estimated accurately.



---

*Source: 963424-2013-12-08.xml*