Monitoring the health of bridges using accelerations from a fleet of vehicles without knowing individual axle weights

ABSTRACT

This paper proposes a new indirect bridge structural health monitoring concept that uses acceleration data from a fleet of different vehicles with unknown weights. When a vehicle passes the bridge, the vertical displacement under its axles can be inferred from its vertical accelerations. This displacement, termed the “apparent profile”, contains two components: bridge profile elevations and bridge deflections under the axle. The two deflection component can be used to find the moving reference influence function (MRIF), defined as the deflection at a (moving) reference point due to a unit load at another point, moving at the same speed. The MRIF can be found when all axle weights are known. In this paper, a new method is proposed to obtain road profile and bridge health condition from the vehicle acceleration, without knowing individual axle weights. Numerical simulation results show that the inferred bridge profile changes when the bridge health condition changes. The difference can be used as an indicator of bridge damage and is illustrated here through an example of bearing damage.

KEYWORDS:

• Acceleration
• bridge
• structural health monitoring
• drive-by
• fleet monitoring
• moving reference influence function

Introduction

Bridge collapse subjected to non-extreme loading conditions, e.g. I-35W Mississippi River bridge (Hao, 2010), Wuxi highway bridge (Hao, 2010; Peng et al., 2020), are of significant concern to owners, operators and users of bridges. Bridge structures deteriorate due to aging, environmental conditions and traffic loading and such deterioration can in part lead to collapse due to non-extreme loading conditions. To prevent unanticipated failures, structural health monitoring (SHM) techniques are sometimes used to monitor the health conditions of bridges. Depending on the type of measurement, bridge SHM methods can be divided into two groups: direct methods with sensors on the structure itself (Chang et al., 2003) and indirect or drive-by methods, with sensors in vehicles passing over the bridge (Yang & Yang, 2018). Direct methods can use, for example, images, strain, deflection, rotation or acceleration (Carden & Fanning, 2004; Graybeal et al., 2002; Li, 2010; Martinez et al., 2020; McGeown et al., 2021). Indirect methods use inertial measurements, such as acceleration or rotation, in passing vehicles. The implicit assumption is that the signals measured in the passing vehicle are affected by the bridge and contain bridge condition information. When the vehicle passes over the bridge, its vibration is affected by both the vehicle properties (speed, mass, etc.) and the bridge properties (flexural rigidity, boundary conditions, etc.). Hence, the bridge condition can be obtained by analysing the signals from the vehicles. Compared with direct methods, indirect methods have the potential to easily monitor a large stock of bridges. Once the measurement system is installed on the vehicles and the system is calibrated, the vehicles can monitor every bridge they cross. In contrast, direct monitoring requires each bridge of interest to be separately monitored which has significant cost implications.

From vehicles’ vibration signals, many bridge dynamic parameters can be obtained: natural frequencies, mode shapes and damping ratio (Malekjafarian et al., 2015; Yang et al., 2019). Yang et al. (2004) first introduce a way to find bridge frequency and successfully identify four frequencies from passing vehicle acceleration signals. Lin and Yang (2005) validated the theory of Yang et al. (2004) through an experimental test. Siringoringo and Fujino (2012) found the frequency of a bridge being crossed by a train (30 km/h speed) in a field test. Siringoringo and Fujino (2012) observed that the vehicle vibration has frequency components that are associated with the vehicle frequencies, vehicle speed, bridge frequencies, and bridge surface profile. Unfortunately, the amplitudes of the bridge frequencies tend to be relatively small. Zhang et al. (2012) firstly extract bridge mode shape from vehicle vibrations and use it to detect damage. In a different strategy to assess bridge condition, bridge acceleration can be decomposed in the frequency domain to produce mode shapes (Malekjafarian & O’brien, 2014). Another approach is to monitor bridge damping according to the magnitude of Power Spectral Density at peaks of bridge and vehicle frequencies (McGetrick et al., 2009). González et al. (2012) found the value of damping ratio through minimising the error between inferred profiles under each axle. In general, damping is difficult to measure in practice as it is influenced by environmental conditions and is not well correlated with structural condition (Cunha et al., 2013).

In recent years, some researchers have focused on axle impact forces or the displacements under the vehicle’s axles (El-Hattab et al., 2015). When the vehicle passes over the bridge, the displacement under the axle is the sum of the bridge profile elevation and the bridge deflection at that point, most of which is static. The location of this bridge deflection changes as the vehicle progresses across the bridge, so the deflection under the axle is termed the moving reference deflection (MRD). The displacement under the axle is termed as the “apparent profile” as it is the profile experienced by the vehicle axle. To infer apparent profile from vehicle vibrations, different approaches have been proposed. O’brien and Keenahan (2015) use a simulated traffic speed deflectometer to determine the apparent profile. They suggest that the differences between the apparent profile from two axles can be used as a damage indicator. However, the measurement accuracy needs to be very high and the concept is sensitive to small differences in axle properties. O’brien et al. (2014) propose a Moving Force Identification method to obtain the dynamic force acting between the bridge and the vehicle. This force is equal to the apparent profile multiplied by the stiffness of the tyre. Yang et al. (2020) propose a way to infer apparent profile when both axle accelerations and the bridge influence line are known. Keenahan et al. (2020) developed an inverse integration method to obtain apparent profile from vehicle accelerations. As all vehicle properties were required in the study by Keenahan et al. (2020), they also propose a fleet monitoring concept to obtain the road profile and vehicle properties from accelerations measured on a fleet of vehicles. A similar approach has been used for railway bridges by Quirke et al. (2017) but with a forward integration method. Yang et al. (2018) introduced the idea of the “contact point response”. Similar to the apparent profile concept, this is the response at the point of contact between the axle and the bridge, which is the same position as the apparent profile. From “contact point” response, bridge frequency can be inferred accurately (Xu et al., 2022). Corbally and Malekjafarian (2022) found that operating deflection shape ratios are correlated with contact point acceleration responses.

Ren et al. (2022) successfully inferred the influence line for MRD from apparent profiles under multiple axles. This influence line is a special case of the moving reference influence function (MRIF) introduced by Wang (2023). MRIF is proposed to describe the relationship between MRD and the axle weights that cause it. By using the apparent profiles under the first axle of each vehicle in the fleet, MRIF is accurately estimated. Wang (2023) also shows that MRIF can be a good indicator of bridge bearing damage. However, this system requires a knowledge of all axle weights on the bridge, which poses a significant practical challenge when dealing with uncontrolled vehicles from the general traffic. In this paper, a similar fleet monitoring concept is developed to monitor bridge bearing damage, but without the need to know the individual axle weights. On the assumption that acceleration data can be procured from vehicles in the general traffic, this will make it possible to monitoring a bridge’s health using that data.

Obtaining bridge profile without knowing individual axle weights

Relationship between bridge profile and apparent profile

By using a similar method to that proposed by Keenahan et al. (2020), an apparent profile calculated from vehicle acceleration is shown in Figure 1 (the downward direction is defined as negative for both bridge profile and MRD). The clear difference between the apparent profile and real bridge profile can be observed in Figure 1, which is termed as moving reference deflection (MRD). Equation (1)(1) $R_x^a=R_x^p+{\delta }_x$(1) gives the relationship between apparent profile and MRD at a point, x:

Figure 1. Apparent profile (blue dash line) and real bridge profile (black solid line).

(1) $R_x^a=R_x^p+{\delta }_x$(1)

where $R_x^a$ denotes the apparent profile at that point, $R_x^p$ denotes bridge profile and ${\delta }_x$ denotes the MRD. An example of the MRD under the second axle of a multi-axle vehicle is illustrated in Figure 2.

Figure 2. Moving reference deflection concept.

The MRIF, $J_{x,d}$, is defined as the bridge deflection response at a moving point (distance x from the bridge start) due to a unit load at distance, d from that location, x (the right direction is defined as positive). According to this definition, the MRD is:

(2) ${\delta }_x=W_1{J}_{x,d_1}+\dots +W_i{J}_{x,d_i}+\dots$(2)

where W_i is the weight of the i^th axle. Combining Equations (1)(1) $R_x^a=R_x^p+{\delta }_x$(1) and (2(2) ${\delta }_x=W_1{J}_{x,d_1}+\dots +W_i{J}_{x,d_i}+\dots$(2) ) gives,

(3) $R_x^p=R_x^a-\left(W_1{J}_{x,d_1}+W_2{J}_{x,d_2}+\dots \right)$(3)

In Equation (3)(3) $R_x^p=R_x^a-\left(W_1{J}_{x,d_1}+W_2{J}_{x,d_2}+\dots \right)$(3) , the apparent profile ($R_x^a$) can be inferred from vehicle accelerations by using an inverse Newmark-Beta method as shown in Appendix. The relative vehicle properties are obtained by using the optimisation method proposed by Keenahan et al. (2020) and more details are going to be shown in Section 3. In this paper, only one vehicle is considered to pass the bridge at a time. In practice, it could be other vehicles pass at the same time. The additional vehicles could be considered as the additional contribution axles in (Equation 2(2) ${\delta }_x=W_1{J}_{x,d_1}+\dots +W_i{J}_{x,d_i}+\dots$(2) ), but their impact has not been assessed in this paper. As the MRIF ($J_{x,d}$) of the bridge is difficult to obtain in practice, a theoretical MRIF can be used as the base line. The difference between the real MRIF and theoretical MRIF will be identifiable in the inferred bridge profile ($R_x^p$). For a simply supported beam bridge, the deflection due to a unit load (distance d from location x) at location x can be estimated through (Equation 4(4) $J_{x,d}=\int \frac{M_x\times {M^v}_x}{EI}dx$(4) ) (Udoeyo, 2020):

(4) $J_{x,d}=\int \frac{M_x\times {M^v}_x}{EI}dx$(4)

where M_x, which is expressed in terms of the x, is the internal moment in the beam caused by the unit load (acting at location x+d). M^v_x is the internal virtual moment in the beam caused by the virtual unit load (acting at location x).

By using a theoretical MRIF, the bridge profile ($R_x^p$) can be found through (Equation 3(3) $R_x^p=R_x^a-\left(W_1{J}_{x,d_1}+W_2{J}_{x,d_2}+\dots \right)$(3) ) when the axle weights are known for each vehicle. In practice, while acceleration data may be provided, the axle weights in uncontrolled vehicles from a fleet of general traffic will not be known. For a fleet of vehicles, all vehicles are subject to the same profile. This feature is exploited here to obtain bridge profile by using the cross entropy optimisation method (Rubinstein & Kroese, 2004). During the cross entropy optimisation process, the mean of Gross Vehicle Weight (GVW) and the mean of first axle weight ratio to GVW are required to generate initial information. It is well known that the mean of gross vehicle weights at a given site is repeatable, especially for vehicles of a given class (Guerson et al., 2016; Van Loo & Lees, 2015). This is illustrated in Figures 3–5 with example of Weigh-in-Motion statistical data of 2-axle vehicles from three different sites in the U.S.A: a) Pennsylvania, b) New Mexico, and c) Minnesota (Walker & Cebon, 2012). Figures 3–5 show the mean of GVW on each month at three sites.

Figure 3. The mean GVWs at the site a in Pennsylvania, U.S.A.

Figure 4. The mean GVWs at the site b in New Mexico, USA.

Figure 5. The mean GVWs at the site c in Minnesota, U.S.A.

From Figures 3–5, although the mean GVWs vary from site to site, the mean GVWs at within each site are close to each other and only vary slightly from month to month. The percentage difference between the mean GVW and the mean of mean GVW is calculated in Table 1. The maximum differences of mean GVW at different sites are +7.1% (site a), −3.0% (site b), and 2.8% (site c).

Table 1. Percentage difference between the mean GVW and the mean of mean GVW.

	Site a in Pennsylvania Mean of mean: 9105 kg		Site b in New Mexico Mean of mean: 9361 kg		Site c in Minnesota Mean of mean: 8786 kg
Month	Mean GVW (kg)	Percentage difference Mean GVW & Mean of means (%)	Mean GVW (kg)	Percentage difference Mean GVW & Mean of means (%)	Mean GVW (kg)	Percentage difference Mean GVW & Mean of means (%)
1	9049	−0.6	9539	1.9	9030	2.8
2	9032	−0.8	9353	−0.1	8902	1.3
3	9055	−0.5	9310	−0.5	8819	0.4
4	8970	−1.5	9076	−3.0	8965	2
5	8950	−1.7	9392	0.3	8873	1
6	8946	−1.8	9285	−0.8	8734	−0.6
7	8959	−1.6	9349	−0.1	8566	−2.5
8	9499	4.3	9587	2.4	8722	−0.7
9	9162	0.6	9380	0.2	8599	−2.1
10	8796	−3.4	9245	−1.2	8781	−0.1
11	9756	7.1	9340	−0.2	8850	0.7
12	9090	−0.2	9473	1.2	8593	−2.2

Similarly, the percentage differences between the mean of the first axle weight ratios to GVW and the mean of mean are also calculated in Table 2. The maximum differences of mean ratio at different sites are+1.3% (site a), −2.3% (site b), 1.9% (site c). This suggests that the mean of GVW and the mean of first axle weight ratio to GVW at a specific site can be assumed to be known with reasonable accuracy.

Table 2. Percentage difference between the mean of first axle weight ratio and the mean of mean.

	Site a in Pennsylvania Mean of mean: 0.388		Site b in New Mexico Mean of mean: 0.374		Site c in Minnesota Mean of mean: 0.418
Month	Mean of 1st axle weight ratio to GVW	Percentage difference 1st axle weight ratio to GVW and the mean of mean (%)	Mean of 1st axle weight ratio to GVW	Percentage difference 1st axle weight ratio to GVW and the mean of mean (%)	Mean of 1st axle weight ratio to GVW	Percentage difference 1st axle weight ratio to GVW and the mean of mean (%)
1	0.385	−0.8	0.371	−0.9	0.422	1.0
2	0.393	1.1	0.370	−1.2	0.424	1.5
3	0.389	0.2	0.365	−2.3	0.425	1.9
4	0.391	0.6	0.369	−1.3	0.419	0.2
5	0.385	−0.9	0.373	−0.3	0.421	0.7
6	0.389	0.2	0.378	1.0	0.415	−0.7
7	0.387	−0.3	0.377	0.8	0.412	−1.5
8	0.388	−0.1	0.379	1.2	0.414	−0.9
9	0.389	0.1	0.377	0.9	0.412	−1.5
10	0.385	−0.8	0.377	0.8	0.419	0.3
11	0.386	−0.5	0.376	0.4	0.413	−1.2
12	0.393	1.3	0.377	0.8	0.419	0.2

Optimisation method to find the bridge profile

The cross entropy optimisation method is an iterative process which first generates a population of GVWs and first axle weight ratio to GVW according to their mean value through random Monte Carlo simulation. The root mean square (RMS)s are used to select “elite sets”. The “elite sets”, which contain the lowest RMSs between inferred bridge profiles and best fit profile. Those “elite sets” are then used to generate a new population. This iterative process will converge to a minimum error between inferred bridge profile and the best fit profile. The flow chart of the process is shown in Figure 6.

Figure 6. Flow chart of the cross entropy optimisation process.

Step 1: For each set, randomly generate GVW and first axle weight ratio for each vehicle according to the mean of GVWs and the mean of first axle weight ratio to GVW. Repeat this process H (= 50 in this paper) times to give H sets of vehicle axle weights for each vehicle.

Step 2: For each vehicle in each set, the profile can be calculated through the (Equation 3).(3) $R_x^p=R_x^a-\left(W_1{J}_{x,d_1}+W_2{J}_{x,d_2}+\dots \right)$(3) This process will give H × N profiles (N is the total number of vehicles in the fleet).

Step 3: Find the best fit bridge profile by using the error minimisation process introduced in Section 2.3.

Step 4: Calculate the RMS between each inferred bridge profile and the best fit bridge profile. The objective function, O, is calculated as the sum of the RMS for each vehicle in each set.

Step 5: For each vehicle, find the 5 “elite sets” with the lowest RMS.

Step 6: According to the “elite sets”, calculate the mean and standard deviation of the axle weights for each vehicle.

Step 7: For each vehicle, regenerate H sets axle weight according to the mean and standard deviation inferred from Step 6.

Repeat Steps 2–7 until the objective function is lower than a threshold. Once the objective function threshold is satisfied, the best fit profile is the output inferred profile.

Best fit profile calculation

In Step 2 in Figure 6, a total of H × N bridge profiles will be calculated. Forcing all calculated profiles to be equal provides a basis for an objective function to be minimised:

(5) $E=\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^K{\left(R_{x_k}^{p,best}-R_{x_k,n,h}^p\right)}^2$(5)

where H is the total number of sets and N is the total number of vehicles in the fleet, $R_{x_k}^{p,best}$ denotes the best fit bridge profile at scan k and location x_k, $R_{x_k,n,h}^p$ denotes the inferred bridge profile at scan k for vehicle n, set h. To reduce the number of unknowns required to describe the best fit bridge profile, Kernel Density Estimation is used as:

(6) $R_{x_k}^{p,best}=\sum _{a=1}^A{r}_{a}exp\left\{\frac{-{\left(x_k-{x^r}_a\right)}^2}{2{\sigma }^2}\right\}$(6)

The best fit profile is formulated as a sum of series Gaussian functions ($r_{a}exp\left\{\frac{-{\left(x_k-{x^r}_a\right)}^2}{2{\sigma }^2}\right\}$). The total number of functions is denoted as “A”. For the Gaussian function, r_a is the elevation of the peak and x^r_a denotes the position of the peak. The term, σ represents the standard deviation or spread of each Kernel function and must be chosen by Engineering judgement – a small value will give a close fit to the representative points while a larger value will smooth out local variability. Using the Kernel representation of the best fit bridge profile, the error function of (Equation 5(5) $E=\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^K{\left(R_{x_k}^{p,best}-R_{x_k,n,h}^p\right)}^2$(5) ) becomes:

(7) $E=\sum _{a=1}^A\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^K{\left(r_{a}exp\left\{\frac{-{\left(x_k-{x^r}_a\right)}^2}{2{\sigma }^2}\right\}-R_{x_k,n,h}^p\right)}^2$(7)

Taking the partial derivative of E with respect to r₁ gives:

(8) $\frac{\mathrm{\partial }E}{\mathrm{\partial }{r}_1}=2\sum _{a=1}^A\left(\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^{K}exp\left\{\frac{-{\left(x_k-{x^r}_1\right)}^2}{2{\sigma }^2}\right\}exp\left\{\frac{-{\left(x_k-{x^r}_a\right)}^2}{2{\sigma }^2}\right\}\right)r_a-2\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^K{R}_{x_k,n,h}^p\times exp\left\{\frac{-{\left(x_k-{x^r}_1\right)}^2}{2{\sigma }^2}\right\}$(8)

To minimise the error function, the partial derivative of (Equation 8(8) $\frac{\mathrm{\partial }E}{\mathrm{\partial }{r}_1}=2\sum _{a=1}^A\left(\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^{K}exp\left\{\frac{-{\left(x_k-{x^r}_1\right)}^2}{2{\sigma }^2}\right\}exp\left\{\frac{-{\left(x_k-{x^r}_a\right)}^2}{2{\sigma }^2}\right\}\right)r_a-2\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^K{R}_{x_k,n,h}^p\times exp\left\{\frac{-{\left(x_k-{x^r}_1\right)}^2}{2{\sigma }^2}\right\}$(8) ) is set to zero, giving:

(9) $\sum _{a=1}^A\left(\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^{K}exp\left\{\frac{-{\left(x_k-{x^r}_1\right)}^2}{2{\sigma }^2}\right\}exp\left\{\frac{-{\left(x_k-{x^r}_a\right)}^2}{2{\sigma }^2}\right\}\right)r_a=\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^K{R}_{x_k,n,h}^p\times exp\left\{\frac{-{\left(x_k-{x^r}_1\right)}^2}{2{\sigma }^2}\right\}$(9)

Taking partial derivatives of E with respect to each r_a gives the system of equations:

(10) $\left\{\begin{array}{c}{B}_1\\ ⋮\\ B_A\end{array}\right\}=\left[\begin{array}{ccc}{C}_{1,1}& \cdots & C_{1,\mathrm{A}}\\ ⋮& \hspace{0.17em}& ⋮\\ C_{A,1}& \cdots & C_{A,\mathrm{A}}\end{array}\right]\left\{\begin{array}{c}{r}_1\\ ⋮\\ r_A\end{array}\right\}$(10)

The components of the matrix, C can be expressed as:

(11) $C_{q1,q2}=\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^{K}exp\left\{\frac{-{\left(x_k-{x^r}_{q1}\right)}^2}{2{\sigma }^2}\right\}exp\left\{\frac{-{\left(x_k-{x^r}_{q2}\right)}^2}{2{\sigma }^2}\right\}$(11)

where q₁ denotes the row number and q₂ the column number. The components of vector, B can be expressed as:

(12) $B_{q1}=\sum _{n=1}^N\sum _{h=1}^H\sum _{k=1}^K{R}_{x_k,n,h}^p\times exp\left\{\frac{-{\left(x_k-{x^r}_{q1}\right)}^2}{2{\sigma }^2}\right\}$(12)

The kernel value matrix can be calculated through:

(13) $\left\{\begin{array}{c}{r}_1\\ ⋮\\ r_A\end{array}\right\}={\left[\begin{array}{ccc}{C}_{1,1}& \cdots & C_{1,\mathrm{A}}\\ ⋮& \hspace{0.17em}& ⋮\\ C_{A,1}& \cdots & C_{A,\mathrm{A}}\end{array}\right]}^{-1}\left\{\begin{array}{c}{B}_1\\ ⋮\\ B_A\end{array}\right\}$(13)

Once the kernel value matrix is obtained, the profile can be found from (Equation 6).(6) $R_{x_k}^{p,best}=\sum _{a=1}^A{r}_{a}exp\left\{\frac{-{\left(x_k-{x^r}_a\right)}^2}{2{\sigma }^2}\right\}$(6)

Numerical model and simulation results

To test the concept numerically, a finite element beam model interacting dynamically with a vehicle model was used to generate vehicle axle accelerations for a fleet of 100 vehicles. The bridge model was simply supported, 4 m wide and 20 m long, made up of three prestressed concrete Y3 beams, topped by a 210 mm deep slab. The bridge was represented numerically by ten beam finite elements and 11 nodes. The properties are summarised in Table 3.

Table 3. Properties of bridge model.

Span	20 m
Number of finite elements	10
Total degrees of freedom	21
Young’s Modulus, E	35 × 10^9 N/m^2
Cross sectional area	2.09 m^2
Second moment of area, I	0.301 m^4
Damping ratio, ξ	3%

To generate vehicle accelerations, the vehicle is considered as a half-car model, as shown in Figure 7. The sprung mass, m_s represent the vehicle body and the sprung mass moment of inertia is denoted by I_s. Two axle components, which are represented by m_u,1 and m_u,2, are separately connected with the body mass by a spring (stiffness is K_s,i) and a damper (coefficients of viscous damper is C_s,i). Each unsprung mass connects to the bridge surface via a spring with stiffness K_t,i. The four independent degree of freedoms correspond to unsprung masses vertical displacement (y_u,1 and y_u,2), body mass vertical displacement y_s and rotation θ_s. The acceleration is collected from the unsprung mass of the first axle. The bridge profile is generated according to the ISO 8608:2016 (ISO, 2016).

Figure 7. Half-car vehicle and bridge model.

To make this numerical study more realistic, we have considered several aspects to make the simulations more close match the physical reality: (1) we have added 2% noise in the accelerations; (2) We have varied some vehicle properties in the fleet to make sure each vehicle has different properties; (3) we used some real vehicle data from physical sites, e.g. speed, axle spacing, gross vehicle weight, ratio between two axle weights. In this study, typical vehicle properties, such as speed, axle spacing and axle weights, are taken from the Weigh-in-Motion database at site b in New Mexico (Walker & Cebon, 2012). Some of the vehicle properties for the fleet are generated randomly by Monte Carlo simulation using the population means and standard deviations shown in Table 4, taken from the literature (Martinez et al., 2020). The vehicle speed, axle spacing, axle location (start/end of the bridge) are assumed to be known in this paper. In practice, those information can be collected by using a axle detector installed near the bridge. Such a detector has been widely used in bridge weigh in motion system (Huseynov et al., 2022; Lydon et al., 2017).

Table 4. Vehicle fleet properties.

Property	Symbol/Units	Mean value	Standard deviation
Unsprung mass	Mu,1: kg	1000	100
	Mu,2: kg	1000	100
Damping	Cs,1: Ns/m	1×10^4	1×10^3
	Cs,2: Ns/m	2×10^4	2×10^3
Unsprung stiffness	Ks,1: N/m	4×10^5	4×10^4
	Ks,2: N/m	8×10^5	8×10^4
Tyre stiffness	Kt,1: N/m	1.75×10^6	1.75×10^5
	Kt,2: N/m	3.5×10^6	3.5×10^5

As only one accelerometer is assumed to be installed on the vehicle and the current methods cannot calculate the apparent profile from single measurements when half-car model is used, a quarter car model is used to represent the first axle of the vehicle to calculate the apparent profile. Quarter-car model shows in Figure 8. The quarter-car model represents the first axle of the Half car. The acceleration is assumed to be collected from the unsprung mass (M_u) of the quarter-car, as shown in Figure 8. For the quarter-car model, the sprung mass, M_s and an unsprung mass, M_u are connected through a viscous damper with coefficient, C_s and a spring with stiffness, K_s. The unsprung mass is connected to the bridge through a spring with linear stiffness, K_t.

Figure 8. Quarter-car vehicle and bridge model.

Three different fleets (each fleet contains 500 different vehicles) are used here to infer the bridge profile. In each fleet, vehicle speed, axle spacing, axle weight are randomly chosen from the database. Other properties are randomly generated according to Table 4. In this method, vehicles in the fleet are assumed to pass the same profile, i.e. the same transverse location on the bridge. In practice, conceptually it could seem that vehicles in a fleet may not pass the same transverse location. However, the distribution of the transverse location of the vehicles in the same traffic line normally follow the normal distribution or Gumbel distribution. Zhou et al. (2015) showed the standard deviation of the transverse location distribution less than 30 cm. Blab and Litzka (1995) showed that transverse position does not greatly affect profile. Hence, it is reasonable to assume that all vehicles pass the same or very similar profiles in the fleet.

As 2% noise is added into the vehicle accelerations, the inferred apparent profile for each vehicle in (Equation 3(3) $R_x^p=R_x^a-\left(W_1{J}_{x,d_1}+W_2{J}_{x,d_2}+\dots \right)$(3) ) contains some error. Although each individual apparent profile is not accurate, the “best fit profile” inferred from the entire fleet of vehicles is accurate, as shown in Figure 9. Figure 9 shows that the three inferred bridge profiles are repeatable and match well with the real bridge profile (black solid line). The difference between the real and the inferred bridge profile is caused by the fact that the different vehicle models are used between acceleration generation and apparent profile calculation. However, the three inferred bridge profiles are close to each other which suggests that good repeatability of the inferred bridge profiles.

Figure 9. Surface profile on bridge, real and inferred (Fleet 1, 2 & 3).

The standard deviation parameter, σ, used in the kernel estimation affects the inferred profile. A simple example of the impact which σ has on the inferred profile is provided to demonstrate how an appropriate value of σ is selected. Firstly, different σ values are used to infer the profile. Figure 10 shows inferred profiles with different σ. When a small σ value is used, the inferred profile is not smooth and contains some isolation arising from individual kernel peaks and troughs. For this example, the function becomes smooth when σ equals 0.23 and the value is adopted.

Figure 10. Inferred profile with different σ values: (a) 20 m of profile; (b) 4 m of profile (10–14 m).

During the profile calculation process, the bridge and its MRIF are assumed to behave perfectly according to the equations of statics. Deviations in the true MRIF from the theory may occur due to 3D effects or non-structural elements (e.g. stiffness of parapet or footpath concrete). Such deviations in MRIF will result in errors in the calculated surface profile. However, these deviations will be repeatable, provided they are consistent between vehicle passes. To simulate the difference between true MRIF and theoretical MRIF in a simple way, 1.05×EI (true) are used as stiffness value to calculate the MRIF used in bridge profile calculation process in the next section.

Bridge bearing damage detection using inferred profile

When the bridge is damaged, the response is affected and hence the MRIF changes. It is challenging to predict how the function will change with damage and no attempt is made here to determine this. A healthy profile is defined as the inferred profile from a healthy bridge and using the theorical MRIF of a healthy bridge. When the bridge is damaged, the deviations of the MRIF between the healthy version and the damaged version will generate further inaccuracies in the inferred profile. A damaged profile is defined as the inferred profile from a damaged bridge and using a theorical MRIF of the healthy bridge. The differences between the healthy profile and damaged profile are proposed here as damage indicators, evidence of bridge damage. The bearing damage model as shown in Figure 11 (Khan et al., 2022) is used here as an example of a change in bridge behaviour that results in a change of MRIF. An increase of rotational stiffness of 10⁹ Nm is assumed as the bearing damage in the first example, i.e. deterioration of the bearing is assumed to be equivalent to this increase in rotational spring stiffness (Khan et al., 2022).

Figure 11. Bridge bearing damage model.

The differences in inferred bridge profiles on the healthy and damaged bridge are shown in Figure 12 (healthy minus damaged). Most of the differences are positive. Damage to the right bearing (red chain line) causes a skew to the right side while damage to the left bearing causes a relative skew to the left. The degree of the difference is greater when both bearings are damaged (yellow solid line).

Figure 12. Difference between inferred bridge profile with healthy bridge and inferred bridge profile with damaged bridge.

Focusing on the case of both bearings being damaged, this calculation is repeated for different damage scenarios. The RMS difference between the inferred profile (damaged bridge) and inferred profile (healthy and theoretical bridge) is given by:

(14) $\mathrm{R}\mathrm{M}\mathrm{S}=\sqrt{\frac{{\sum }_{k=1}^K{\left(R_{x_k}^p\left(healthy\right)-R_{x_k}^p\left(damaged\right)\right)}^2}{K}}$(14)

Figure 13 provides a RMS plot with different increases of rotational stiffness. As would be expected, the RMS increases with rotational stiffness as the bridge support behaviour changes from being predominantly pinned towards predominantly fixed. As has been noted by Kliewer and Glisic (2017), most of the changes occur in a small log-scale range, in this case when K_rot is between 10⁸ Nm and 10¹¹ Nm. The RMS difference does not change significantly outside of this interval, giving the characteristic S-shaped behaviour of Figure 13. Clearly, early detection of bearing failure should focus on the time when the increase in RMS starts to become significant (around 10⁸ to 10⁹ Nm in this case). In practice, the difference due to damage might not be observable for small levels of damage and the effect of noise might make it difficult to detect damage. A threshold is proposed that is sufficient to overcome the natural variation due to noise and different vehicle properties.

Figure 13. Variation of RMS with increase of rotational stiffness.

Conclusion

This study proposes a new concept for drive-by fleet monitoring that uses accelerations from single axles in a fleet of two-axle vehicles to identify bridge damage. For the first time, the bridge profile can be reliably predicted without knowing the weight of each axle in advance. When bridge is damaged, the inferred bridge profiles are influenced by the deviation of the damaged moving reference influence function from the corresponding healthy function. In simulations, a non-zero difference between inferred (healthy) bridge profile and inferred (damaged) bridge profile shows when bridge bearings are damaged and gives an indication of which one is damaged. For a simply supported bridge, the results show that the damage indicator follows an S-shaped pattern, not changing significantly for very low or very high stiffness values.

Acknowledgments

The authors also gratefully acknowledge the Federal Highway Administration and the Long-Term Pavement Performance Program for access to WIM data.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The author, Wang Shuo, wishes to acknowledge the financial support received from University College Dublin and the Chinese Scholarship Council, grant number: 201908300011.

Appendix

The equations of motion of the vehicle (Quarter-car in Figure 8) can be written as:

(A1) $\left[M_v\right]\{{\ddot{u}}_v\}+\left[C_v\right]\left\{{\dot{u}}_v\right\}+\left[K_v\right]\left\{u_v\right\}=\left\{F_v\right\}$(A1)

where [M_v], [C_v] and [K_v] are mass damping and stiffness matrix of the vehicle model, respectively:

(A2) $\left[M_v\right]=\left[\begin{array}{cc}{m}_s& 0\\ 0& m_u\end{array}\right]$(A2)

(A3) $\left[C_v\right]=\left[\begin{array}{cc}{c}_s& -c_s\\ -c_s& c_s\end{array}\right]$(A3)

(A4) $\left[K_v\right]=\left[\begin{array}{cc}{k}_s& -k_s\\ -k_s& k_t+k_u\end{array}\right]$(A4)

The vectors, $\left\{{\stackrel{\cdot \cdot }{u}}_v\right\}$, $\left\{{\dot{u}}_v\right\}$, $\left\{u_v\right\}$ are vehicle accelerations, velocities and displacements, respectively. The displacement vector of the vehicle is:

(A5) $\left\{u_v\right\}=\left\{\begin{array}{c}{u}_s\\ u_u\end{array}\right\}$(A5)

The force vector is:

(A6) $\left\{F_v\right\}=\left\{\begin{array}{c}0\\ f_v\end{array}\right\}$(A6)

where f_v denotes the force acting on the vehicle wheel. In the Newmark-Beta method, the integration constant is listed:

Time step, $\Delta t=0.001\gamma =\mathrm{0.8.}\beta =0.25\times {\left(0.5+\gamma \right)}^2,a0=1/\left(\beta \times \Delta t^2\right),a_1=\gamma /\left(\beta \times \Delta t\right),a_2=1/\left(\beta \times \Delta t\right),a_3=1/\left(\beta \times 2\right)-1,a_4=\gamma /\beta -1,a_5=\Delta t/2\times \left(\gamma /\beta -2\right),a_6=\left(1-\gamma \right)\times \Delta t,a_7=\gamma \times \Delta t.$.

The displacement and velocity of the unsprung mass at each time step can be calculated using the Newmark-Beta method:

(A7) $u_{u,t+1}=uu,t+1+a_2{\dot{u}}_{u,t}+a_3{\dot{u}}_{u,t}{a}_0+u_{u,t}$(A7)

(A8) ${\dot{u}}_{u,t+1}={\dot{u}}_{u,t}+a_6{\ddot{u}}_{u,t}+a_7{\ddot{u}}_{u,t+1}$(A8)

According to the (Equation A1)(A1) $\left[M_v\right]\{{\ddot{u}}_v\}+\left[C_v\right]\left\{{\dot{u}}_v\right\}+\left[K_v\right]\left\{u_v\right\}=\left\{F_v\right\}$(A1) , the equilibrium of sprung mass gives:

(A9) $k_s{u}_{s,t+1}-k_s{u}_{u,t+1}+c_s{\dot{u}}_{s,t+1}-c_s{\dot{u}}_{u,t+1}+m_s{\ddot{u}}_{s,t+1}=0$(A9)

In the Newmark-Beta method:

(A10) ${\dot{u}}_{s,t+1}=a_1\left(u_{s,t+1}-u_{s,t}\right)-a_4{\dot{u}}_{s,t}-a_5{\ddot{u}}_{s,t}$(A10)

(A11) ${\ddot{u}}_{s,t+1}=a_0\left(u_{s,t+1}-u_{s,t}\right)-a_2{\dot{u}}_{s,t}-a_3\hspace{0.17em}{\ddot{u}}_{s,t}$(A11)

Substitute Equations (A10)(A10) ${\dot{u}}_{s,t+1}=a_1\left(u_{s,t+1}-u_{s,t}\right)-a_4{\dot{u}}_{s,t}-a_5{\ddot{u}}_{s,t}$(A10) and (A11)(A11) ${\ddot{u}}_{s,t+1}=a_0\left(u_{s,t+1}-u_{s,t}\right)-a_2{\dot{u}}_{s,t}-a_3\hspace{0.17em}{\ddot{u}}_{s,t}$(A11) into (A9)(A9) $k_s{u}_{s,t+1}-k_s{u}_{u,t+1}+c_s{\dot{u}}_{s,t+1}-c_s{\dot{u}}_{u,t+1}+m_s{\ddot{u}}_{s,t+1}=0$(A9) gives:

(A12) $k_s{u}_{s,t+1}-k_s{u}_{u,t+1}+c_s{a}_1\left(u_{s,t+1}-u_{s,t}\right)-a_4{\dot{u}}_{s,t}-a_{5}us,t-c_s{\dot{u}}_{u,t+1}+m_s{a}_0\left(u_{s,t+1}-u_{s,t}\right)-a_2{\dot{u}}_{s,t}-a_{3}us,t=0$(A12)

Rearrange (Equation A12)(A12) $k_s{u}_{s,t+1}-k_s{u}_{u,t+1}+c_s{a}_1\left(u_{s,t+1}-u_{s,t}\right)-a_4{\dot{u}}_{s,t}-a_{5}us,t-c_s{\dot{u}}_{u,t+1}+m_s{a}_0\left(u_{s,t+1}-u_{s,t}\right)-a_2{\dot{u}}_{s,t}-a_{3}us,t=0$(A12) gives:

(A13) $u_{s,t+1}=\frac{\left(c_s{\dot{u}}_{u,t+1}+k_s{u}_{u,t+1}\right)+\left[\left(\begin{array}{c}{c}_s{a}_1+\\ m_s{a}_0\end{array}\right)u_{s,t}+\left(c_s{a}_4+m_s{a}_2\right){\dot{u}}_{s,t}+\left(c_s{a}_5+m_s{a}_3\right)u_{s,t}\right]}{\left(k_s+c_s{a}_1+m_s{a}_0\right)}$(A13)

The effective stiffness matrix is:

(A14) $\left[\bar{K}\right]=\left[K_v\right]+a_0\times \left[M_v\right]+a_1\times \left[C_v\right]$(A14)

The effective force matrix is:

(A15) $\left\{{\bar{F}}_{t+1}\right\}=\left[\bar{K}\right]\times \left\{u_v\right\}$(A15)

The relationship between effective force and the force vector is:

(A16) $\left\{F_{v,t+1}\right\}=\left\{{\overline{F}}_{t+1}\right\}-\left[M_v\right]a_0\left\{u_{v,t}\right\}+a_2\left\{{\dot{u}}_{v,t}\right\}+a_{3}uv,t-\left[C_v\right]a_1\left\{u_{v,t}\right\}+a_4\left\{{\dot{u}}_{v,t}\right\}+a_{5}uv,t$(A16)

According to (Equation A16)(A16) $\left\{F_{v,t+1}\right\}=\left\{{\overline{F}}_{t+1}\right\}-\left[M_v\right]a_0\left\{u_{v,t}\right\}+a_2\left\{{\dot{u}}_{v,t}\right\}+a_{3}uv,t-\left[C_v\right]a_1\left\{u_{v,t}\right\}+a_4\left\{{\dot{u}}_{v,t}\right\}+a_{5}uv,t$(A16) , the force (f_v) can be inferred. As the relation between the force and the displacement under the wheel (u_w) is:

(A17) $f_v=k_t\times u_w$(A17)

Hence, the displacement under the wheel can be calculated by:

(A18) $u_w=f_v/k_u$(A18)