A direct method to detect and localise damage using longitudinal data of ends-of-span rotations under live traffic loading

Existing work on rotation-based bridge monitoring has focused on indirect methods, such as bridge weigh-in-motion or influence line approaches. However, these approaches require increased instrumentation complexity, and require calibration, necessitating bridge closures. In this paper, we explore the potential of using rotation measurements to create a more practical and cost-effective monitoring system. To this end, we present a damage detection method which directly analyses bridge rotation data measured under live, free-flow traffic loading. We show how the Earth Mover’s Distance, typically used in statistics and image processing, can be applied directly on end-of-span rotation measurement data to achieve effective damage detection and localisation. Numerical simulation results demonstrate the approach’s robustness to the confounding effects of temperature variation and traffic diversity (vehicle type, loading, and velocity). The direct rotation measurement approach is applied to data from an in-service short-span bridge to demonstrate the technique’s capability with free-flow traffic loading.


Introduction
The challenge of managing bridges is critical in maintaining modern transport infrastructure. This is evidenced by the disruption and economic losses caused by a bridge closure such as London Hammersmith Flyover in 2011 [4]. The scale is significant as, for example, one third of 12,000 state-maintained bridges in France require repairs, and 840 are at risk of collapse [13]. Many bridge Structural Health Monitoring (SHM) approaches use bridge acceleration as the input signal to the analysis method [1,9]. With recent advances in non contact computer vision methods, SHM approaches that track bridge displacement by identifying a change in stiffness, have started to be proposed [12,16,33]. These approaches are perceived to offer advantages as information on the static response is more evident in the displacement signal and hence a change in stiffness may be easier to identify. Unfortunately, long term, bridge displacement measurement is logistically difficult and expensive.
In recent years, researchers have started to explore end rotation measurements as they contain similar information to displacement but rotation can be more easily measured and is seen as a more practical means of obtaining information on the bridge static response. The ease in measuring rotation is facilitated by modern, high-grade accelerometers, which can reliably detect the small rotations of the bridge/beam by measuring the corresponding rotation of the Earth's gravity acceleration vector. Existing work can be broadly categorised as rotation-based bridge weigh-in-motion (BWIM), and rotation influence line approaches. Rotation-based BWIM approaches [18] seek to detect and potentially localise damage by tracking the change in the statistical distribution of the calculated vehicle weights. Rotation influence line-based approaches [2,10] typically utilise load testing with vehicles of known weights, which when repeated over a long period of time can highlight damage in the structure.
Both approaches are based on a two step process, specifically: (1) A bridge closure to carry out load testing using a vehicle of known axle weights and spacings to calculate the rotation influence line. (2) For the WIM approach [18], an axle detection system needs to be permanently installed on the bridge and a complex optimisation procedure used to calculate the vehicle weights. For the influence line approaches [2,10], the approach carried in (1) needs to be carried out in similar environmental conditions, although this is not highlighted in [2,10].
Both approaches have major limitations from a practical engineering perspective either due to road closures or the need to perform measurements with the same environmental conditions. The underlying principle [2,10,18] is that if damage occurs in the structure, the rotation response to a given load will increase. This work proposes a novel direct method for damage detection which performs analysis directly on measured rotation signal data. This works on the observation that the distributions of axle weights and spacings are statistically repeatable. Thus, the distribution of maximum end-of-span rotation amplitude values for single truck events will also be repeatable, if considered over long enough time periods that environmental conditions become repeatable. If the distribution of maximum rotations changes significantly, it will thus indicate that the structure may have experienced damage.
A key challenge is to quantify this change for which we propose the use of a technique commonly used in image processing, namely the Earth Movers' Distance (EMD), which has been described in Sect. 2.2. Even when variation between measurement intervals due to temperature variation, traffic diversity, and sensor noise is allowed for, the increased EMD values post-damage are clearly evident. Damage localisation can be carried out as the increase in each end-of-span rotation post-damage is proportional to the proximity of the damage to a given sensor. Thus, the differential between both ends-of-span EMD values postdamage infers the damage location. Furthermore, we have shown our approach performs equivalently when only considering subsets of all vehicles by axle count.
Our main contributions are as follows: -Creation of a damage detection method based on direct analysis of maximum rotation signal amplitude datasets to detect shift in distributions that occur due to damage, thereby avoiding the need for calibration requirements and specific direct rotation measurement. -The approach is shown to be robust to confounding effects of vehicle diversity and temperature variation through Monte Carlo simulations including a realistic sensor noise model based on experimental measurements. -The work presents the first results to show that grouping rotation signals by axle count allows the same clar-ity of results as looking at the full population of trucks, resulting in significant reductions in dataset volume. -For the first time, the feasibility of the proposed approach for extracting rotation signals from direct measurements under live traffic loading for long term monitoring on a short-span bridge is shown. This is achieved by performing a pilot study on an in-service short-span bridge.
The structure of this paper is as follows: Sect. 2 introduces related work and the bridge finite-element (FE) model used is outlined in Sect. 3. Section 4 presents our accelerometer model based on experimental noise profiles, and the signals processing scheme used to extract maximum rotation amplitudes from noisy rotation signals via free-flow traffic loading. Our proposed damage detection approach is set out in Sect. 5. Section 6 describes our numerical simulations and results for the proposed damage detection approach. Section 7 outlines the case study using live traffic on an in-service bridge.

Related work
Conventionally, managing structural condition relies on periodic visual inspection by engineers, but its accuracy has been questioned, both generally [26] and specifically for bridges [7]. Therefore, an effective and reliable means of monitoring structural performance data in real-time is needed, to either remove or supplement the qualitative, subjective and unreliable human element of conventional approaches [3]. Measurements rely on physical variables consisting of strain, bridge natural frequency (typically determined from vertical acceleration measurements), displacement, or rotation. Whilst strain is a relatively localised effect, frequency can provide a global measure of a structure providing damage detection but not localisation. As in-field measurements of displacement are problematic due to technical and operational challenges, rotation or inclination measurements have been proposed as an alternative, more practical means of measuring bridge deflection [8,25].

Rotation-based Bridge SHM Work
To date, relatively little work has been carried out to look at the effectiveness of using rotation measurements for bridge SHM. Huseynov et al. in [10] proposed a damage detection approach based on the difference between rotation influence-lines obtained pre-and post-damage from direct rotation measurements using both numerical simulations and laboratory tests. Alamdari et al. [2] proposed a similar technique but used rotation influence lines obtained indirectly via strain measurements as opposed to directly measuring rotation. O'Brien et al. [18] presented a damage detection approach based on the overestimation of vehicle weights by a rotation-based BWIM system following damage. The work in [28,34], proposes an approach which uses rotation measurements in combination with strain measurements and a baseline FE model of the bridge to detect damage. It is demonstrated using numerical simulations of a bridge with simplified vehicle loading.

Earth Movers' Distance
Following the occurrence of damage, the resulting loss of stiffness leads to a systemic increase in rotation values, which will be evident in the statistical distribution of rotation values. Thus, by identifying and quantifying this shift in the distribution, the presence of damage can be established. The Earth Movers' Distance (EMD), which is a measure of the distance between two distributions [15], is used in this work as it offers a consistent measure of dissimilarity between two distributions [24]. The Earth Movers' Distance, which is also known as the first Mallows, or Wasserstein, distance in mathematics, is given by: where (u, v) are the set of probability distributions on ℝ × ℝ , with marginals, u and v. If U and V are the cumulative distribution functions of u and v, this also equals: with a fuller exposition of the above equations found in [22].
This work considers only one-dimensional, empirical distributions, which are a special case for which the EMD has an exact solution. If X and Y are empirical distributions of two datasets, i. e. X = {x (1) (1) ≤ ⋯ ≤ y (n) } , both of size n, then the Earth Movers' Distance can be calculated as described in [15] as follows: It is evident that to apply Eq.(3) requires that both distributions are the same size. However, in practice the distributions used in this application are relatively large, and thus, can be resampled or binned to meet this requirement [15].
To date, EMD has been used in image processing [24], biology [20], and environmental science [32]. Existing work has shown the potential of EMD for damage detection based on mode frequencies in a limited study using a numerical cantilever beam model [19].

Rotation signal generation
Numerical simulations based on an FE bridge model are used to test the capability of the proposed bridge damage detection and localisation approach. The key aspects of the numerical simulations carried out within this work are highlighted in Fig. 1. Ultimately, the FE model outputs the end-of-span rotation signals, which are affected by: (a) the structural properties of the modelled bridge; (b) a damage model allowing for a localised loss of stiffness due to a cracked element; (c) the confounding effect of temperature variation, i. e. stiffness variation with changing temperature; and (d) the diversity of vehicle axle weights, spacings, and velocities.

Finite-element model
A finite-element model developed within Queen's University Belfast [6,9] is used here and is presented in Fig. 2b. Responses used from the model include, displacement at midspan, indicated as v MID , and rotation at both left-and right-hand supports, indicated as LHS and RHS respectively. The numerical values of the model's structural properties have been chosen to approximate the 34 m span concrete beam and slab bridge shown in Fig. 2a. The simulated transverse width is slightly wider than one traffic lane to allow for some transverse load distribution into adjacent beams. This model has been chosen as it represents a very common bridge type, and therefore is characteristic of many in-service bridges. In addition, this bridge is used in a pilot study on the feasibility of long term monitoring of rotation values reported in Sect. 7. The bridge is modelled as a discretised simply supported finite-element beam as depicted in Fig. 2b. The axles of the trucks are modelled as a series of discreet point forces which move across the bridge. The image in the figure shows a three axle vehicle with axles spacing's of a1 and a2 respectively. However using a similar approach the model can simulate any number of axles. To allow for dynamics, the response of the beam to a moving force is solved using the second order matrix differential equation given as:

Structural Properties Environment
where (t) contains the vertical displacement and rotation of the degrees-of-freedom of the model at time, t; , , and denote the consistent global mass, damping and stiffness matrices of the bridge; is the vector of applied forces at each degree-of-freedom, which is calculated by distributing the moving force to the nodes of the underlying element using the Hermitian functions [23]; and , i. e. the global stiffness and mass matrices of the bridge, are assembled from the elementary stiffness and mass matrices, and [14]. As damping in bridges is typically low (1-3%) the effect of damping on the forced static response is minimal. Furthermore, as damping will reduce the magnitude of the dynamic rotation in the forced response, and later it will be shown that our approach relies on filtering these out to extract the static response, a no damping scenario actually presents a tougher challenge. Therefore, in the simulations the damping matrix, = .
When the matrices, and , and the forcing vector, , have been populated, Eq.(4) can be solved using a discrete time integration scheme [29] yielding the unknowns: ̈ , ̇ and , i. e. acceleration, velocity, and displacement respectively at each timestep.

Damage model
Despite a wide array of bridge damage mechanisms, one commonality amongst many of these is a resultant loss of (4) structural stiffness, either globally, or locally. Hence, damage is introduced in the FE model by allowing for a localised loss in stiffness due to a crack. The stiffness reduction used in this work is the one proposed by Sinha et al. [27] which quantifies a gradual loss of bending stiffness in the vicinity of the crack that extends 1.5 times the depth of the beam at both sides of the crack location. The latter is modelled by introducing a reduced moment of inertia in the elementary matrixes of those elements close to the crack. Previous work on damaged beams has focused on rectangular beams, where the ratio, = h∕d , of crack height, h, and beam depth, d, has often been used to characterise the severity of the damage. values of 0.1 and 0.2 represent a localised loss in stiffness of 27.1% and 48.8% respectively of the inertia of a healthy rectangular section. The same equivalency between and associated percentage stiffness loss is maintained for the beam sections and can be calculated using the equation, The effect of damage on the bridge deflection response is shown in Fig. 3a, b by plotting the midspan displacement and left end-of-span rotation responses respectively, for the same point load, P 1 = 300kN , at velocity, ̇b = 10m s −1 , with the bridge initially in a healthy condition, i. e. = 0.0 . Subsequently, a damage located at the third-span point is modelled at = 0.1 and = 0.2 , i. e. stiffness losses of 27.1% and 48.8% respectively. As expected, the midspan displacement and end-of-span rotation responses have similar profiles, as rotation is a spatial derivative of displacement, and the amplitude of both increase with increasing damage levels.

Vehicle loading parameters
The moving vehicle loading on the bridge is modelled as an ensemble of N point loads, P 1 , … , P N , representing each individual axle weight at separations of, a 1 , … , a N−1 , and this can be seen in Fig. 2b. The movement of the vehicle across the bridge is modelled as a constant velocity by incrementing b(t) at each timestep.   The effects on rotation signals due to changes in the inputs are demonstrated in Fig. 4a, b. Specifically, Fig. 4a plots the rotation responses to a 5-axle vehicle when (i) laden, 400kN , and (ii) unladen, 175kN . It demonstrates that when the vehicle is laden, it approximately doubles the maximum rotation response relative to the unladen vehicle, as expected given the approximate doubling in load magnitude. Figure 4b highlights the laden, 5-axle vehicle at a range of velocities, {15, 20, 25}m s −1 . Figure 4b shows that the response's dynamic component increases with increasing velocity.

Temperature-dependent stiffness model
As the Young's modulus of a structure's materials is temperature dependent, changes in a structure's temperature can result in significant shifts in response magnitudes. To model this effect, we use the negative sloped, linear relationship from [11,Eq. (10)]: where E c is the Young's modulus of concrete and T is the environmental temperature in • C . This allows the stiffness matrix, to be varied according to the desired environmental temperature of the model.
To illustrate the effect of temperature variation, a 5-axle truck weighing 120kN is modelled crossing for T = {0, 25} , with resultant rotation responses plotted in Fig. 4c. From Fig. 4c, it can be seen that between T = 0 and T = 25 that the maximum rotation value increases by ∼ 13%.

Rotation measurement and signal processing
Whilst our finite-element model simulates the rotation response of a bridge under vehicle loading at a particular temperature, this does not match the reality of noisy, measured signal data. Therefore, to allow the practical application of our damage detection approach, methods for rotation measurement with free-flow traffic loading are presented herein.

Accelerometer rotation model
Whilst the FE model shown in Sect. 3 provides end-of-span rotation signals LHS and RHS , these represent ideal measurements and do not consider the effect of the transducer's (sensor) transfer function on the actual measured signal.
Here, we present our signals model, derived from real-world measurements, with the aim of capturing key accelerometer errors such as bias instability and random error. The basic principle of rotation sensing using an accelerometer is illustrated by Fig. 5a, b. An accelerometer assumed fixed to the bridge deck is orientated with the x-axis co-linear to the bridge's longitudinal axis and the y-axis orthogonal to the bridge deck. As the bridge deflects, this causes the accelerometer to rotate by an equivalent angle, , as seen in Fig. 5a. In the reference frame of the accelerometer, this results in a shift in the resolved x and y components of the measured acceleration due to gravity, g. The vector diagram of this is shown for the x-axis component in Fig. 5b. Thus, the relationship between the ideal (noise-free) x-axis acceleration measurement component, x, and the sensor rotation, , is given by: Errors in inertial sensors, such as the accelerometers used in this work, encompass a number of effects, namely: quantisation noise, velocity random walk, bias instability, random walk, and accelerometer trend [21]. While adding Gaussian noise is an effective technique to approximate these errors, and has been used by others in bridge SHM, one of its' weaknesses is that it assumes the noise amplitude is the same across the full frequency spectrum. In practice, this assumption is not always true, for example in Fig. 6b, it can be seen that when the noise spectra (of the JA accelerometers used in this study) is measured empirically the noise spectrum shows peaks and troughs; thus, showing this assumption is not correct.
The presence of low frequency noise (e. g. below 1 Hz, i. e. 10 0 in Fig. 6b), which won't be removed by our low-pass filtering, will both cause a bias to the measured signal and introduce baseline drift.
To model this noise, we consider a Fourier-transformed surrogates method as described in [30], which ensures the simulated signals more closely match these real-world noise profiles, as seen in the equivalent frequency spectra between the measured noise in Fig. 6b, d. This procedure is represented by the block diagram in Fig. 5c. The top branch of the block diagram produces surrogate noise data, which have identical frequency spectra but randomly varied phases, . Thus, the noise component, n, in our model can be tuned to match experimental noise measurement data, by updating the Fourier spectrum, M, used to generate the surrogate noise signal, n, via an inverse discrete Fourier transform (IDFT). The surrogate noise signal, n is then added to the ideal rotation, x, to produce the equivalent noisy acceleration signal, a. Figure 6 illustrates the operation of our sensor model, wherein firstly, static noise measurements are carried out with the sensor placed on a solid concrete slab floor, with care taken to minimise any sources of external vibration; an example of this is shown in Fig. 6a using a JAE JA70-SA triaxial MEMS accelerometer. Figure 6b shows the corresponding frequency spectra. The components of this measured frequency spectrum, can then be assigned random phase shifts, and transformed back to the timedomain. This results in a surrogate for the original measured noise signal, whereby the surrogate signal has an equivalent frequency spectrum but a unique, uncorrelated time series, as seen in Fig. 6c

Rotation signal processing
In this section, we cover the signal processing required to first extract transient rotation signals from noisy accelerometer data, and carry out feature extraction on the filtered rotation data, in this instance via peak prominence extraction. The measurement of a transient rotation signal is challenging, especially on shorter/stiffer structures with lower magnitudes of deflection/rotation. The measured rotation signals have two main components, a relatively slowly varying static component associated with the deflected shape of the beam as it is loaded, and a dynamic component due to the vehicle-bridge interaction system dynamics. The appreciably lower frequency of the static component compared to the dynamic components, allows these to be reliably separated using a low-pass filter. The measured noise figure, is typically quite large relative to the rotation signals, as seen in Fig. 6e. However, the low pass filtering used to separate the static and dynamic components will also reduce the noise power, and allow the static component to be extracted.
The flow chart given in Fig. 7, shows the four stages of our signal processing pipeline. It should be noted that due to the small magnitude of the transient signals, the step of converting acceleration back to rotation using inverse sine can be avoided due to the small angle approximation of sine, i. e. sin( ) ≈ . The function of each stage shall now be described below: Windowing: As we are interested in the short portions of signal corresponding to vehicle crossing events, a window function is used to extract these events. The sliding rectangular window used herein achieves this by zeroing all bar the most recent L number of data points, where L is the length of the window, i. e. : The length, L, of the window can be determined by setting an upper limit on the rise time based on the slowest vehicle crossing to be captured. Detrending: To reduce the DC bias present in the accelerometer signal, the mean of the first m samples of the window, in this case chosen to be 10% of the window length, L, is subtracted from the signal, as follows: Denoising Filter: To remove the measured noise, and the dynamic rotation components, a fourth order low-pass Butterworth filter is used. The cut-off frequency, f c , for this filter can be selected based on three constraints: , is equal to its cut-off frequency, f c , the commonly used approximation for the relationship between bandwidth and rise time, B ⋅ r = 0.35 , can be used to obtain a minimum cut-off frequency, i. e. f c ≥ 0.975 Hz . Hence, the cut-off frequency, f c , for this implementation was set at 1 Hz. As predicted, with the cut-off frequency set close to the lower bound, in this case f c = 1 Hz , the peak value becomes consistent with the equivalent static model peak shown in Fig. 8. Whilst a small delay can be noted in Fig. 8, as our approach only uses the signal magnitudes, and does not rely upon processing any temporal or spectral information in the signal, this phase shift does not affect performance.
Peak Prominence: To extract the maximum rotation values, a peak prominence is calculated from the difference in the minimum and maximum values within the current window, i. e. a PK = max a f − min a f . Thus, a sequence of peak prominence values are obtained, and a subsequent filtering can be carried out to extract the peak value for each vehicle event.
To illustrate the application of peak prominences, three 5-axle lorries of varying load and speed and speed are modelled crossing in short succession in Fig. 9. The blue plot in Fig. 9 shows the static response, and it can be seen that each of the vehicle crossing events results in an approximately half sine shape with a response of zero between the events. The pink plot in Fig. 9 shows the signal obtained after filtering and while there are still three peaks evident the exact static response has not been recovered, specifically: (i) the times the peaks occur at are shifted due to the delay that happens in the filter (discussed in previous comment), and (ii) in the period between the peaks when there was no vehicle on the bridge the response is non-zero. This is due to the residual low frequency noise in the filtered signal.
The residual low frequency noise comes from the noise at less than 1 Hz shown in Fig. 6d, as this falls within the filter's passband, we cannot remove it without compromising the static component. The effect of this noise can be seen in the baseline wander/drift in the filtered signal in Fig. 9. From this time series, the peak prominence values for each vehicle were extracted and shown in Table 1 with the equivalent static response for comparison.

Proposed damage detection approach
Using the FE models (Sect. 3) and signal processing (Sect. 4), we set out our proposed approach here for damage detection and show it remains performant in spite of the presence of confounding effects, i. e. temperature  variation and vehicle diversity. Numerical simulations supporting this are presented in Sect. 6 and a field-trial to demonstrate the feasibility of the proposed front-end measurement system is shown in Sect. 7. We first show the basic principle of our approach in Sect. 5.1, wherein we note that as long as these confounding effects converge to repeatable distributions the systemic effect of damage on maximum rotation distributions may be unearthed. By exploiting the superposition principle of linear time-invariant (LTI) systems, which we show maximum rotation response obeys, we subsequently show the convergence of the confounding variables of vehicle diversity (Sect. 5.2) and temperature variation (Sect. 5.3).

Basic principle of earth Movers' distance-based damage detection
We consider a bridge with three conditions, a healthy ( = 0.0 ) state and two damaged states ( = 0.1 and = 0.2 , both at L/3). The confounding effects of temperature and vehicle variation are limited by using a fixed temperature, T = 25 • C and vehicle velocity, v = 15m s −1 , and restricted vehicle diversity to the most common axle geometry of five axle truck. Firstly, we randomly select N = 200 vehicles to cross the bridge for each condition. Observing the trend over the sampling index (Fig. 10a), the shift for even a relatively large damage in a very ideal scenario is not immediately apparent. When a cumulative distribution function (CDF) (Fig. 10c) is plotted for the data, it is still hard to draw any reliable conclusions (even with the correct state labels known), due to the dissimilarity of the distributions, e. g. the CDFs for When N is increased ( N = 2000 ) it is still not obvious from the sampling index plot (Fig. 10b) that damage has occurred. However, the damage is evident in the N = 2000 CDF plot (Fig. 10d). Specifically, the rotation CDFs converge, wherein it can be noted the shape of the CDFs becomes consistent, and hence, the systemic increase in rotation values due to damage is more reliably observed.
Quantifying this shift in CDFs post-damage then forms the basis of our proposed damage detection approach. The EMD which we propose to quantify this shift, can be visually interpreted as the area enclosed between the two CDF plots being compared. This is illustrated in Fig. 10c, d by the pale pink hatching shown between = 0.0 (green solid line) and = 0.1 (orange dashed line), and the pale grey shading between = 0.0 (green solid line) and = 0.2 (purple dotdash line). Table 2 provides details of this metric for pairs of condition states of the same bridge. S 1 represents the initial condition of the bridge (i. e. the value of 1 and the damage position, x d 1 ) and S 2 is the second condition. Once again, N = 1000 randomly chosen five axle trucks are simulated and the maximum rotation value sets are denoted, L | S i and R | S i , for the left and right ends-of-span respectively given condition, S i . EMD L | S 1 , L | S 2 denotes the EMD between the distribution of left hand rotations in condition S 1 and in condition S 2 , and EMD R | S 1 , R | S 2 for the equivalent right hand side rotations. The procedure is repeated for 100 iterations to obtain the average EMD values, to obtain the range of EMD values likely to be observed given this scenario. The first row of Table 2 considers when both conditions, S 1 and S 2 , are healthy. This indicates that the natural variation between successive healthy conditions is likely in the order of 2.86 × 10 −6 rad . Thus, any damage producing variations smaller than this would be unlikely to be detectable which fortunately, is the not the case as all EMD values for varying damage conditions as shown in rows 2-7 are larger than this, indicating that these damages would likely be identifiable.
From Table 2, it can also be noted that: (a) when damage occurs at midspan, L/2, the increases in EMD are equivalent at both left and right ends-of-span; (b) whenever the damage is located away from midspan, in this case L/3, the increase in EMD is larger at the end-of-span closest to the damage; (c) the magnitude of the increase in EMD is proportional to the increase in the damage severity, ; and (d) the EMD values between conditions 1 = 0.0 and 2 = 0.2 are approximately equal to the sum of the EMD values between conditions 1 = 0.0 to 2 = 0.1 and 1 = 0.1 to 2 = 0.2.
Effectively the confounding effects of vehicle payload, axle geometry, vehicle speed and temperatures can be thought of as inputs to the system, with rotation being the output. Table 2 demonstrates the approach using a simplified scenario with fixed vehicle speeds, axle geometry and temperatures, so only vehicle payload variation remains. Subsequently, as the other confounding effects such as vehicle speeds, axle geometry and temperature are introduced these will increase the variation of rotation values observed. However, as with vehicle payload these other confounding effects (vehicle speeds, axle geometry and temperatures) will still have statistically repeatable distributions given appropriate sampling.
The significance of each of the confounding effects having repeatable distributions is that this will ultimately result in repeatable rotation distributions. Hence, the shift in rotation distributions post-damage can be observed. To demonstrate this point mathematically, the system equations for an Euler-Bernoulli beam are presented below as the signals processing approach (Sect. 4.2) produces values equivalent to the maximum values of a static rotation response to within a small margin of error (Fig. 8).
for a simply supported beam of length, L, moment of inertia, I, and Young's Modulus, E. Wherein, the load, or input, to our system consists of a vehicle with C axles modelled as a series of point loads of individual magnitudes, P i , at positions, x = x i . The system outputs are in this case, the rotation at the left hand support, L , and at the right hand support, R .
Based on Eqs.(10) and (11) and Eq.(6) for a given structural condition, these equations are an LTI system, i. e. the output is a linear combination of the inputs and does not depend on the time these inputs are applied. Hence, any variation in the system's output can be related to the  system's input. For example, the results of Fig. 10 and Table 2 are based on inputs that only vary in axle weights, P i , and I, which is reduced by damage; thus, any change in I due to damage is relatively easy to observe. In a more realistic scenario, the axle geometry, x i , will also change with differing vehicles and likewise the stiffness, E, with varying temperatures. However, as long as the distributions of these effects are repeatable, then the systemic effect of I being reduced by damage may be unearthed. Due to the superposition principle of LTI systems, we can treat the vehicle axle weights and spacings, and temperature as independent inputs which will then combine linearly to give our system's output, rotation. Therefore, in Sects. 5.2 and 5.3 using the EMD we show the convergence of these input distributions, and thus, the convergence of our EMD damage metrics.

Confounding effect of vehicle diversity
In live traffic, individual vehicles exhibit a wide array of weights, axle spacings, and velocities, all of which influence the measured maximum rotation value. To ensure that these will not cause unwanted variation in the maximum rotation distributions, we show that by choosing an appropriately large sample size, these effects can be made repeatable to within a small margin of error.
However, the measurement scheme presented in Sect. 4 provides a maximum rotation value that is equivalent to the static response to within a small margin of error. Thus, vehicle velocity can be largely ignored and it is reasonable to assume that the overall distribution of vehicle speeds converges with increasing sample size as evidenced by the distribution in Fig. 11, taken from a vehicle speed compliance dataset available from UK Department for Transport [5].
Using a sample of 500,000 vehicles with axle counts, 2 ≤ C ≤ 6 from the weigh-in-motion dataset [31], the convergence of distributions of axle weight versus the axle location measured from the lead axle at x = 0 , is shown in Fig. 12a. The vertical bands represent the variation in the axle weights for a given axle count and they have a linear scaling effect on rotation response as seen in Eqs. (10) and (11). On the other hand, the variation in axle geometries will have a non-linear effect on rotation responses. However by using a large enough sample, this converges to a repeatable distribution.
The variation in axle geometry and axle weights for a given set of vehicles can be viewed using a CDF of each axle's location, expressed as a lag from the lead axle, a i , weighted by its load, P i , i. e. : This is plotted in Fig. 12b. The value of these curves at x = 0 represents the average proportion of a vehicle's total weight to be exerted by the lead axle. The steepness of the curve at a given location indicates both the number of axles observed at this location and their weight. For example, by examining the 5 and 6 axle curves, we can note the relatively flat region between approximately x = 4 and x = 8 , as the load for vehicles within these classes is largely observed at the tractor unit's axles at the front, i. e. x < 4 , and the trailer's axles at the rear, i. e. x > 8.
To show the convergence of these vehicle load distributions, a series of cumulative distribution plots of load proportion plotted against axle position are presented in Fig. 13, for varying sample sizes, N = {10, 100, 1000} , and axle counts, C = 5 , C = 6 , and 2 ≤ C ≤ 6 . The shared x and y axes are shown along the bottom and the right hand side respectively. The CDF plots in first row deals with 5 axle trucks , and the first plot is for three different samples, {S 1 , S 2 , S 3 } each with just N = 10 trucks in the sample. The inset allows the variation between the CDFs for the 3 different samples to be better visualised. The second plot in row 1 presents the CDF for three different samples of 5 axle trucks (12) with N = 100 trucks per sample. Comparing the insets in Fig. 13a, b, it is evident that by increasing the number of trucks in the sample particularly to over N = 1000 , we start to get a clear convergence. The 6 axle trucks in row 2 of Fig. 13 shows a similar pattern of convergence with increasing sample number. The third row deals with samples with mixed vehicle axle counts and it is clear from examining the insert in first CDF that there is significantly greater variation when the samples were made up of a trucks all with the same number of axles. However, once the sample number reaches N = 1000 , there is still convergence despite the mixed axle count. Thus, it can be seen that these sample distributions converge with increasing sample size, N; however, by considering only vehicles of a single axle count, e. g. C = 5 , these converge at smaller sample sizes.
The variation between these samples can be quantified by computing the EMD between pairs of these samples. We , and for axle counts, 2 ≤ C ≤ 6 , C = 5 , and C = 6 . For each pair of samples, the average EMD value is computed by repeating the computation 100 times. Table 3 shows that both the EMD mean and standard deviation decreases as the sample size, N, increases, and that taking samples from subsets of vehicles by axle count, e. g. C = 6 , further reduces these.

Confounding effect of temperature variation
As stiffness depends on temperature which varies over both diurnal and seasonal timescales, this can result in significant natural variation in structural responses. If unaccounted for, this can potentially mask the presence of other response trends due to damage. However, this can be overcome by choosing long enough measurement windows so that temperature variation converges towards its annual distribution. This is now demonstrated using the dataset of historical hourly air temperatures from 1949 to 2019 at Aldergrove, Northern Ireland [17]. In the top row, even with increasing sample sizes, no convergence is observed due to the large biases and differing variances of the samples taken from different 6 month windows. Whenever considering twelve month windows in the second row, convergence is observed as the bias introduced by seasonal temperature variation is mostly averaged out, particularly for larger values of N. The longer multi-annual windows as shown in the third row with M = 24 months, provide only limited additional convergence. This is confirmed by the summary EMD metrics in Table 4    increases. With M = 6 month windows, there is only a modest reduction in mean EMD and the standard deviation remaining proportionally high, thus showing the lack of convergence. With M = 12 months, a significant reduction can be observed in the EMD's mean and standard deviations. Moving to M = 24 months only provides a small reduction in EMD values for all values of N.
As a major source of variations between these random sample sets lies seasonal temperature trends, which have a period of 12 months, meaning window lengths less than 12 months length introduce a bias between the sets, e. g. consider the case of two 6-month windows covering April-September and October-March, hence why these show higher EMD values despite increasing the sample size. In aggregate, as each annual cycle is broadly the same as the previous, the advantage of moving from M = 12 to M = 24 , which span an integer number of annual cycles, lies mostly in the tails of the distributions. However, as the tails of the distribution have relatively low weights, these only make a marginal contribution to reducing the calculated EMD value. This suggests M = 12 month long windows provide a reasonable balance between reducing the confounding effect of temperature, and ensuring the measurement windows are practical.

Numerical simulations and results
To show the performance of our proposed approach with confounding effects of vehicle diversity and temperature variation, a multi-year, progressive damage simulation has been carried out using the FE model described in Sect. 3. Damage has been modelled as described in Sect. 3.2 with the damage located at x d = L∕3 . The largest damage level simulated ( = 0.1 ) for this localised damage case is in line with those commonly quoted in previous literature, e. g. [10,28], with the remaining damage levels being smaller or much smaller than those commonly reported.
As damage severity can be assumed to be a monotonically increasing function of time, this has been modelled as a series of step functions (Fig. 15) with this profile having the advantage of allowing the natural variation that occurs for a specific damage level, e. g. = 0.05 , to be observed.
The first year serves as the baseline with years 2-11 continuing to simulate the healthy bridge as this allows the natural variation that occurs withing the healthy case to be observed over a 10 year period. Subsequently, the damage progression includes ten years at each damage level, = {0.0, 0.01, 0.025, 0.05, 0.1} , where again the decade long simulation at each damage level allows the natural variation within a given damage level to be observed. It should be noted that Fig. 15 only presents simulation parameters, specifically the damage level to be simulated in each year of the simulation, with the output/results of these simulations reported in the later figures, e. g. Fig. 17.
A fifty-one year window covering the years 1969-2019 of hourly air temperature readings from Aldergrove, Northern Ireland [17], (see Sect. 5.3) is used to provide the temperature profile. The simulation models approximately one hundred thousand, single truck (vehicles with a gross weight ≥ 3.5t ) crossing events per annum, which are distributed based on the weekday and hour in proportion to the WIM dataset shown in Sect. 5.2. Vehicle speeds have been assigned using distributions by axle counts based on motorway speed compliance data [5], and these are shown in Fig. 11.
Thus, each sequential one hour interval consists of: 1. Updating our model stiffness, E, using the temperaturedependent stiffness relationship described in Sect. 3.4; 2. Getting the numbers of vehicles, N 2 … N 6 , for each axle count, 2 ≤ C ≤ 6 , by taking a random sample of M vehicles from the WIM dataset (equivalent to the expected weekly traffic volume, in this case, M = ⌊100000 ÷ 52⌉ = 1923 ) and counting the occurences matching the axle count, weekday, and hour; 3. Sampling N 2 … N 6 vehicles from the the WIM dataset subsetted by the corresponding axle count; 4. Assigning each vehicle a speed based on the distributions shown in Fig. 11; 5. Simulating each vehicle through our FE model, sensor model, and signal processing pipeline, and storing extracted values in the maximum rotation datasets.
A CDF plot of the maximum left and right end-of-span rotation values for year 1 is indicated as the black dashed line in Fig. 16a, b, and these form the baseline for future comparisons. The steps 1 to 5 are then repeated for years  Fig. 16a, b, and from these figures, it can be seen that there is some variation between years, largely due to the confounding effect of variations in seasonal temperatures. The remainder of the years are simulated similarly, with the damage level in the FE model updated as described in Sect. 3.2 according to the profile shown in Fig. 15. From  Fig. 16a, it is possible to see the overall trend of shifting towards higher rotation values as damage, albeit for the lowest level of damage ( = 0.01 ) the shift is quite small. From Fig. 16b, it can been seen that the same trend is present for the right end-of-span rotation values but the shift is smaller, as the damage is further from the sensor on the right hand end of the span .
We can use the EMD between a baseline and successive data to both automate the analysis of, and quantify, this shift in CDFs. For the left and right hand ends-of-span respectively, these damage indicators are computed as: wherein, L i and R i , denote the set of left and right hand endsof-span rotation values from the ith year; thus, L 0 and R 0 indicate the baseline data, i. e. the first year, as Sect. 5.3 shows that twelve month long windows can effectively minimise the variation between temperature distributions.  Fig. 17a, which was generated by taking 1969 as the healthy baseline year and the blue data markers (1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979) show 10 years where the bridge is still simulated as healthy. In this 'healthy decade', there is some variation present, e. g. for = 0.0 , these values range from 0.05 × 10 −6 and 0.6 × 10 −6 . This is to be expected, as there will still be some variation in the temperature profile of the baseline and successive periods. This variation is plotted as the blue bar chart on the right hand side of the figure, which can be thought of as the distribution plot based on 10 years of data. The periods 1980-89, 1990-1999, 2000-2009 and 2010-2019 simulate damage levels of = 0.01 (orange), = 0.025 (green), = 0.05 (pink), and = 0.1 (brown) respectively. As was the case for the 'healthy decade', natural variation occurs between successive years within a given damage level. This is to be expected, as there will still be some variation in the temperature profile of the baseline and successive periods.
Despite the natural variation which occurs, it can be seen that there are clear increases in DI LL which are approximately proportional to the relative increase in damage severity, . The separation between the different damage levels is evident in the bar charts for each of the damage level (right of the figure). From this, it could be difficult to reliably identify a damage level of = 0.01 due to how close its distribution sits with respect to the healthy (blue) one. However, for all higher levels of damage, there is clear separation, indicating that the proposed method could reliably detect them. Similarly, less noticeable trends can be observed in Fig. 17b  for DI RR ; however, this is to be expected as the damage is located further from the right end-of-span.
The significance is that the difference in DI LL and DI RR can be used to approximately locate of the damage. In isolation, the occurrence of damage is clear in Fig. 17a, but when viewed in conjunction with Fig. 17b, the observation that the value of DI LL is larger than the corresponding values of DI RR indicates that the damage is in left hand portion of the span. If the damage were in the right hand portion of the span, DI RR would show a greater level of increase. Finally, in a damage situation (i. e. DI LL and DI RR are both elevated with respect to their baseline values) where DI LL and DI RR are showing similar levels of increase, it can be concluded that the damage is close to midspan. With the the relative increases in DI LL and DI RR in Fig. 17a, b, it is possible to get a qualitative indication of the damage location. However, it is also possible to use the information contained in the relative magnitudes of these to obtain a numerical estimator for the damage location,  Figure 17a, b show the results when the rotation response of every truck crossing the bridge is recorded, which equates to approximately 100,000 vehicles per annum. However, by only processing rotation data for subsets of similar vehicles, such as by axle count, e. g. C = 5 or C = 6 , we can reduce the amount of data that needs to be collected and processed with only marginal reductions in damage detection performance. From the point of view of practically implementing this detection approach, the reduction the data volume is significant and important, as the power consumption and communications bandwidth of any embedded monitoring , otherwise.
. system, are typically constrained; hence, collecting, processing and transmitting less data helps alleviate these constraints. Figure 17c, d illustrate this by showing the same damage indicators for only trucks with axle counts, C = 5 , equating to approximately 33,000 vehicles per annum. The trends of DI LL and DI RR increasing with increasing damage levels can be observed; however, due to the smaller sample size, a larger variation is present for a given damage level, thus the shifts due to very small damages, e. g. = 0.00 to = 0.01 , are less distinct. By only considering 6 axle trucks, this reduces the vehicles per annum to approximately 12,500; however, as shown in Fig. 17e, f, the increased variation swamps the shift from = 0.0 to = 0.01 . Finally, Fig. 18b shows the damage location estimators calculated using the values of DI LL and DI RR from Fig. 17e, f and Eq. (15).
The results presented in Figs. 17 and 18 show the feasibility of our proposed damage detection and localisation approach based on identifying the shift in maximum rotation response distributions with post-damage, in spite of the presence of the confounding effects of vehicle dynamics, diversity and temperature variation.

Rotation measurement field trial
A field study was also undertaken to demonstrate the feasibility of recording direct rotation measurements from live traffic loading. The field trial was carried out on the bridge shown in Fig. 2a, which is reproduced and annotated in Fig. 19a. The bridge is a single 34 m span, pre-cast, reinforced concrete beam and slab bridge, constructed adjacent to an extant masonry arch bridge.
For the test, a single JAE JA70-SA triaxial MEMS accelerometer was used to measure the rotation at the southern end-of-span as seen in Fig. 19a and placed on the pedestrian  (Fig. 19b). The sensor was located approximately over the web of the outermost bridge beam. The camera was placed at the northern end of the bridge to record the vehicles crossing the bridge. Using the signal processing approach of Sect. 4.2, the end-of-span rotation response of the bridge for a number of single truck events was obtained and is presented in Fig. 19c, d. The first peak in Fig. 19c (622 seconds) is due to a five axle tanker truck and has an amplitude of approximately 0.7 × 10 −4 rad . A frame grab of the truck is included (labelled as truck 1), although somewhat obscured by water on the camera lens due to the heavy rain. The second peak (630 seconds) has an an amplitude of approximately 2.0 × 10 −4 rad and was generated when a six axle bulk feed truck crossed the bridge. The peaks in Fig. 19d at 740 seconds ( 1.0 × 10 −4 rad ), 759 seconds ( 0.5 × 10 −4 rad ) and 762 seconds ( 0.4 × 10 −4 rad ) are from a 4 axle car transporter lorry, a 4 axle tipper lorry and a 2 axle lorry respectively.
While the weights of the trucks are not known accurately, based on the engine tone of the truck 2 (six axle), the truck appeared to be heavily laden, and if so the amplitude of of rotation observed ( 2.0 × 10 −4 rad ) is close to the FE model prediction for a fully loaded six axle truck. In this context, the smaller rotation amplitudes observed for trucks 1, 3, 4, and 5, and their relative proportions are consistent with the values expected based on numerical simulations. The results from the pilot study indicate that that the proposed sensing and filtering approaches are working effectively. This demonstrates that the sensing and signal processing approach set out herein are feasible, and are giving credible values based on engineering judgement.

Discussion, limitations, and future work
While the results shown in Sect. 6 are encouraging, in the sense that low levels of damage are identified despite significant confounding effects, it is acknowledged that this was achieved via numerical simulations rather than field experiments. In reality, fully validating this approach would entail field experiments and a situation whereby damage could be introduced to the bridge in the second year of the experiment. Unfortunately, this is not logistically feasible. However, in future work, we hope to do the initial step of this validation to establish that the healthy rotation distributions obtained from long-term (i. e. multi-year) field measurements are in line with the distributions generated via numerical simulations.
In this paper, only a single-span bridge was considered with rotation measured at the ends of the span, and hence, further work is required to look at multi-span bridges. In the case of a continuous multi-span bridge deck subject to a moving load it is anticipated that sensors will be required over each support, i. e. not just at the ends of the deck. This extra instrumentation is based on the the fact that the change in rotation at the end of a continuous multi span deck due to damage in an interior span will likely be quite small, and hence difficult to detect. The numerical estimate for the damage location (Eq.(15)) shown herein would also not hold as it was developed based on a single simply supported span with a continuous I value. However, if rotation values at each support were known, some degree of damage localisation would likely be possible. In the case where the individual spans of a multi-span bridge are simply supported, it is anticipated that an accelerometer would be placed at each end of a given span, and Eq.(15) still holds.
This paper uses the commonly used model of a single damage location in a simply supported beam. However, in a situation where the damage occurs at several locations on the bridge, the proposed algorithm will still detect this combined damage's occurrence, as the stiffness losses will still cause a resultant increase in the end-of-span rotations. Identifying the approximate location of the different local damages would be more difficult. For example, in the case of multiple damage locations, the result from Eq.(15) would indicate an averaged damage location based on the weighted contribution of the various losses in I value across the span, rather than identifying the individual damage locations.
The bridge model used in this work assumed pin roller support conditions at the ends-of-span; however, this assumption may not be valid for all bridges. Prior to installation on-site, rotation responses will need to be checked over the bearings and at incremental distances away from the bearing to determine the degree of fixity in each bearing. Future simulation work will need to look at the significance of this effect. In particular different degrees of fixity at each end would need to be examined, as it is likely this would make damage localisation more challenging. Furthermore, the effect of changing bearing support conditions, e. g. due to bearings rusting, should be investigated to determine if this can be identified distinctly from damage in the bridge span.
The proposed damage detection approach uses two accelerometers, and while a signals model was used to model the accelerometer (Sect. 4.1), which considers sensor noise, bias, and baseline drift. However, this model did not consider the potential scenario of varying accuracy between the two sensors, for example, due to installation issues, such as axis misalignment, and damage or degradation during operation. If some of these effects occur to one accelerometer, that sensor will be less accurate than the other. Given the case that one sensor is less accurate, and depending on the exact nature of this error, this will likely make it harder to detect damage using that end-of-span's data and also impede successful damage localisation. Therefore, further study should be conducted to investigate the change in damage detection and localisation performance associated with this issue and if diagnostics can be implemented to preempt any subsequent performance loss.
In Sect. 6 ( Figs. 17 and 18), it is observed that by grouping the bridge rotation response by the type of vehicle that caused the response, a similar level of damage detection capability can be achieved by using a significantly smaller data set. Essentially, this indicates that if the loads causing the response in the healthy bridge are similar to the loads causing the response in the damaged bridge, it is easier to identify damage. More specifically, fewer vehicle runs are needed to expose the damage than with more diverse vehicle populations. For the moment, this has only been shown using numerical simulations of one bridge. Based on engineering judgement, it seems plausible that grouping responses by truck type on other bridges will achieve similar results. However, more work, both numerical simulation and experiment, is needed to verify this. As this was a numerical study, it was easy to group responses by the truck type causing it. Hence, an effective means of carrying out real-time vehicular classification on-site should be investigated and deployed for the practical implementation of this approach.

Conclusion
This paper presented a damage detection and localisation method based on using the EMD to quantify the shifts in maximum rotation response distributions due to damage. It was shown that the confounding effects of temperature and vehicle variation are related to rotation response by an LTI system. The convergence of the rotation distributions, which is required for repeatable distributions to measure shift postdamage is demonstrated via the convergence of the temperature distributions for year-long measurement windows, to reduce seasonal trends, and by the convergence of vehicle load distributions for sufficiently large sample sizes.
Numerical simulations demonstrate the performance of our proposed damage detection and localisation approach. From these simulation results, it can be seen that: 1. Despite the inclusion of confounding effect of temperature and vehicle variation, damages as small as = 0.01 may be detected. 2. An estimator of the damage location can be obtained based on the differential shifts in left and right end-ofspan rotation distributions post-damage. 3. Processing only subsets of rotation response based on a specific axle count, e. g. C = 6 , causes only a marginal reduction in performance but can be achieved with far less data, i. e. an 87% reduction in the amount of data for C = 6.
Finally, a field trial has been carried out which shows the effectiveness of the proposed sensing and filtering approach for extracting rotation signals from a bridge under live traffic loading.