Characterization of the road profile and the rotational stiffness of supports in a bridge based on axle accelerations of a crossing vehicle

In an effort to find more cost‐effective and proactive ways to keep bridges in good condition, the use of instrumented vehicles has gained great interest in the last decade. Two bridge components that can wear rapidly are the bearings and the road surface. However, past research on drive‐by monitoring has placed focus mostly on detecting losses of bending stiffness in the bridge deck, while assuming ideal support conditions that may differ from real cases significantly, and ignoring the characterization of the road profile. Even further, the need for specialized vehicles equipped with high‐tech instrumentation, low speeds, or very good road profiles has been a major obstacle preventing its practical implementation. This paper investigates the use of axle accelerations from an ordinary two‐axle vehicle crossing the bridge to quantify the rotational stiffness of the supports and the height of the road irregularities while overcoming the limitations exposed above. In contrast to previous research where the response of the contact point has been derived from other vehicular locations based on complex differential equations of motion, transfer functions are employed here. The key advantage of transfer functions is their simple algebraic form that can be easily calibrated on the field. The road profile is then obtained by subtracting the displacement of the bridge under each axle from the displacement of the contact point. There is one prediction of the road profile per axle but only a unique value of rotational stiffness at each support that will yield the same prediction by both axles. The algorithm is successfully tested with a half‐car traveling at 5, 10, 15, and 20 m/s, over a 15‐m bridge beam model with ISO road classes “A,” “B,” and “C,” for boundary conditions ranging from simply supported to fixed. The solution's robustness to modeling inaccuracies and noisy data is also investigated.


INTRODUCTION
Bridge managers dedicate a significant portion of their budgets to emergency repairs, leaving limited funds for inspection, monitoring, and keeping bridges healthy. The periodicity of the inspections varies between countries, that is, a minimum of one visual inspection takes place every year in France, every 2 years in the United States and the United Kingdom, and every 3 years in Germany unless an issue is detected in the bridge. This inspection has the purpose of recording the condition of the bridge, determining any damage, recommending repairs or other services as needed, and rating the bridge regarding its load-carrying capacity (Waheed & Adeli, 2005). The gap between inspections can be relatively long, and resources for continuous structural health monitoring (SHM) are costly and mostly allocated to long-span bridges (Less & Adeli, 2010). If, for instance, a cable-stayed bridge is observed to vibrate considerably, then, active, passive or semi-active control systems can be adopted to prevent a potential failure of the cable connections (Ghaedi et al., 2017). SHM data are limited to a number of measurement points, whereas the number of parameters needed to build a computer model of the bridge can become enormous. The latter leads to illconditioned problems that may hinder the identification of damage. Bayesian frequency-domain substructure identification (Yuen & Huang, 2018) and convolutional neural networks (Youqi Zhang et al., 2019) are among the SHM approaches proposed to deal with these issues in large structures.
Over the past two decades, drive-by monitoring of bridges, pioneered by Y.B. Yang and C.W. Lin Yang et al., 2004), has attracted significant attention due to its low-cost implementation and high-energy efficiency (Malekjafarian et al., 2015;. Drive-by, or indirect monitoring, relies on sensors placed on vehicles to provide an alternative to the traditional approach of direct monitoring, in which the sensors are glued or attached to the bridge. The feasibility of the drive-by approach to extract the modal properties of the bridge has been widely studied in the literature, that is, frequencies (Gomez et al., 2011;McGetrick et al., 2009), mode shapes (Oshima et al., 2014;Yao Zhang et al., 2012), and damping Keenahan et al., 2014). Furthermore, changes in the extracted modal properties or their derivatives can be used to obtain information about damage on the bridge (Cerda et al., 2014;Hou & Xia, 2021;Kim & Kawatani, 2008). However, there are still some challenges associated with drive-by considering the fact that vehicular responses usually contain not only (1) the bridge component but also (2) the road and (3) vehicular components. The last two are the dominant components in the majority of cases, especially in the presence of rough road profiles and high vehicular speeds. Therefore, it is essential to develop techniques aimed at eliminating or mitigating components (2) and (3) from the vehicular response to isolate the bridge component before a high-level assessment can be achieved.
For the removal of the road component, several approaches have been tested based on the subtraction of measurements from two axles of the same vehicle (Keenahan & OBrien, 2018;Keenahan et al., 2014), bandpass filters (Liu et al., 2021), or signal separation techniques, that is, blind signal separation (Li et al., 2019). However, these works often adopt assumptions that limit their applicability, for example, the subtraction-based algorithm by Keenahan et al. requires a specialized vehicle with identical axle properties, whereas the filtering and separating techniques need a priori knowledge of the range where to look for the bridge frequency. Regarding the removal of the vehicular component, Yang et al. have recently proposed to use the response of the contact point between the vehicle and the bridge as a reference for assessment purposes . They obtain the theoretical acceleration of the contact point for a vehicle bridge interaction (VBI) system consisting of a sprung mass crossing a simply supported bridge. The bridge frequencies and mode shapes are identified from this acceleration response for a constant vehicular speed of 2 m/s. Yang et al. extend their previous work and successfully estimate the bridge damping by applying the Hilbert transform to the acceleration response of the contact point (Yang et al., 2019). These acceleration responses are computed by taking a double derivative of the displacement responses, which can be measured via a set of laser sensors mounted on the vehicle body. The latter are simulated using a VBI model consisting of a simplified two-axle vehicle with two degrees of freedom (DoFs) moving over a simply supported bridge at 10 m/s. Following this concept, other authors have used the contact point response to achieve modal identification in drive-by monitoring. For instance, Nayek and Narasimhan propose a Gaussian process latent force model to recover the contact point response (Nayek & Narasimhan, 2020). Then, the bridge frequency is extracted from the recovered contact point acceleration. The authors conduct numerical simulations in a VBI model consisting of a quartercar (QC) moving at low speed (2 m/s). Corbally and Malekjafarian derive the contact point acceleration response from the equations of motion for a QC model by differentiating the displacement twice with respect to the sampling time increment (Corbally & Malekjafarian, 2021). Further on, the same authors build upon their work to develop an artificial neural network-based algorithm that monitors the bridge behavior by using the Fourier spectrum of the contact point acceleration as the input (Corbally & Malekjafarian, 2022). Before these contributions are able to calculate the contact point response, they require to know values that may be difficult to gather in a real-life situation, such as the vehicular accelerations of all DoFs, including both the vehicle body and axle measurements for a given sampling rate, in addition to the vehicle properties. Li et al. put forward a novel two-step approach to determine the contact point response for a VBI system consisting of two uncoupled QC models and a simply supported beam model (Li et al., 2022). First, the Newmark-beta-based dual Kalman filter is utilized for computing the contact point forces, and then the corresponding contact point response is derived based on the contact forces and the vehicle properties. In summary, most of the research presented above implements the concept of the contact point response in a simple vehicle model, that is, sprung mass or QC. Moreover, some of the proposed methodologies require laser sensors (Yang et al., 2019) for the measurement of the contact point displacement, which is expensive and challenging to manage for general practice. Even further, a few studies rely on the ability to measure the body mass response, that is, body mass displacement and rotations (Corbally & Malekjafarian, 2021. Measuring the body mass response in a practical scenario can present challenges, such as determining the center of gravity or being able to place the sensors at its location, that is, the center of gravity may be located outside the vehicle. Ultimately, the implementation of some of these approaches creates immense difficulties, and, not surprisingly, most of them have not been tested in a real-life situation yet. This paper proposes a novel drive-by approach to SHM that uses as input only the axle accelerations of the moving vehicle to achieve two targets: (1) characterization of the road profile on the bridge, and (2) assessment of the bridge boundary conditions, that is, quantifying the rotational stiffness of the bridge supports. Target (1) is very relevant to traffic safety as well as bridge health given that the road profile can be responsible for significant impact forces on a bridge, especially when the road deteriorates and its roughness increases. Nonetheless, the surface of the bridge is rarely monitored or reported in the literature, except for a few publications on the presence of bumps in the approach to the bridge as a result of settlements or expansion joints (Du et al., 2015;Henderson et al., 2016;Marques Lima & de Brito, 2009). Most of the work conducted over recent years concentrates on the estimation of road profiles without a bridge (Nguyen et al., 2019), and again, some of these algorithms require dynamic responses measured at the body mass, that is, sprung mass bounce accelerations and pitch velocities , and body acceleration and angular velocity (Zhao et al., 2017(Zhao et al., , 2019).
An alternative method for detecting road damage, that is, potholes, is the use of machine learning techniques in which the original training set composed of images taken by an in-vehicle camera is enlarged by including artificially generated samples (Maeda et al., 2021). The assessment of boundary conditions, target (2), is also a significant problem in bridges and deserves attention. According to the design manual for roads and bridges published by Highways England (DMRB, 2020), the design working life for bridge structures is typically over 120 years, but it can be reduced to 50 years for the bearings. Some natural hazards involving scour and flooding can also shorten the working life of the bridge bearings. From modal analysis of field measurements before and after bearing replacement, it becomes clear that the rotational stiffness of the support changes significantly (Khan et al., 2022). The removal and replacement of bearings is a costly operation that requires road closures and the use of hydraulic jacks to lift the bridge. Ideally, bridge managers would like to delay this replacement beyond the design life once the bearing did not exhibit signs of deterioration. Clearly, a continuous monitoring of the boundary conditions is deemed to be necessary for ensuring bridge safety as the design life becomes closer or exceeded. This paper pursues to achieve this goal by means of a drive-by algorithm capable of quantifying the rotational stiffness of the bridge bearings.
Section 2 presents the VBI model, consisting of a four-DoF half-car and a 15 m simply supported discretized beam. Different rotational stiffness is assumed at both end supports. Former drive-by approaches commonly establish the contact point response in time domain. Instead, a frequency domain approach based on the transfer function (TF) associated with the axle displacement is proposed here and derived for a half-car model. Section 3 describes the methodology followed to meet targets (1) and (2), including a novel calibration of the initial condition (IC) of the axles. Section 4 presents detailed results for different values of the rotational stiffness, vehicle speed, and road roughness. The impact of modeling and measurement errors on the performance of the algorithm is also analyzed. Finally, Section 5 summarizes the main findings and makes suggestions for future work.

2.1
Definition of the VBI system Figure 1 illustrates the four-DoF half-car model with two connected axles used to represent the vehicle in this paper.
The half-car model consists of two axle masses ( , ), each of them connected to the body mass ( , ) by a spring-dashpot system ( , , , ) simulating the suspension. The axle masses are linked to the road profile and the bridge through springs representing the tires ( , ). It is worth noting the different mechanical properties adopted for the front and rear axles of this ordinary vehicle, in contrast to the strict values of axle properties of specialized vehicles required by other algorithms. The four DoFs of the model are , , , and , which represent the body mass displacement, body mass rotation, front axle displacement, and rear axle displacement, respectively. Moreover, and are the contact point displacements of the front and rear axles, respectively. The gravitational acceleration is taken as 9.806 m/s 2 . Table 1 shows the typical values of vehicle properties adopted for the VBI system, based on OBrien et al. (2014).
The bridge of the VBI system is discretized as a series of Euler-Bernoulli beam elements 1.0-m long each. gives the main parameters of the concrete bridge (Feng et al., 2021). The boundary conditions of the bridge are modeled using rotational springs at the left ( 1 ) and right ( 2 ) supports. Other researchers use the hysteretic Bouc-Wen model to represent the translational stiffness of the bearing and to predict the response of an existing bridge (Gutierrez Soto & Adeli, 2019) or the parameters of the bearing model in real time from the forces directly applied to the bearing using the Hilbert-Huang transform (González & Aied, 2015) or a Bayesian framework (Yuen et al., 2019). This paper assumes that a linear approximation of the rotational response will suffice to assess a change in support conditions in the presence of small traffic loads. For the analyses conducted to illustrate the methodology in Section 3, 1 = 2 × 10 6 kN⋅m/rad ( = 10 6.3 ) and 2 = 8 × 10 6 kN⋅m/rad ( = 10 6.9 ) are employed. In Section 4, values of 1 and 2 are varied to cover a wide range from simply supported to fixed support conditions. The first three undamped bridge frequencies of longitudinal vibration are 5. 66, 22.62, and 50.90 Hz in the case of simply supported support conditions. Bridge damping is typically small, that is, around 3% (Castellanos-Toro et al., 2018), and its effect on forced displacements is neglected in this paper. In free vibration, damping plays an important role in bringing the amplitude of dynamic vibrations due to the inertial forces down to zero after a number of cycles. However, in forced vibration, the bridge response is made of static and dynamic components. Damping only affects the dynamic component and it does so to a small extent due to the tiny differences between damped and undamped frequencies and the short duration of the vehicle crossing. The bridge is assumed to be at rest when the vehicle enters the bridge, that is, environmental vibrations are neglected, and a sufficiently large gap needs to be ensured with the preceding vehicle to prevent non-zero initial vibrations. The last condition is relatively easy to meet for a short-span bridge where the amplitudes and periods of vibration are comparatively small. For safety purposes, a driver should stay at least 2 s behind any vehicle that is directly in front of them, and the displacements in free vibration are reduced by half after approximately four periods of vibration for only 2%-3% bridge damping.
A dynamic transient simulation is used to replicate the response that would be measured in the field. This simulation involves the integration of the coupled differential equations of motion of the VBI system based on the properties of Tables 1 and 2 to obtain a time series of the axle accelerations due to the passing of the vehicle at a constant speed. In this paper, the Newmark-beta method (Newmark, 1959) is applied with = 0.5 and = 0.25 and a time step of 0.001 s. Further details of the integration methods for a coupled VBI system can be referred to the comprehensive review by González (2010). Additionally, a road profile is simulated and extended for a certain length before and after the bridge. The road profile is modeled as the sum of a series of sinusoids characterized by the power spectral density of its displacements. The numerical approach for the generation of the road profile can be found in the paper by Agostinacchio et al. (2014). According to the International Organization for Standardization, ISO 8608 (ISO, 2016), roads are classified into classes "A," "B," "C," "D," "E," "F," "G," and "H" based on their degree of roughness. Geometric means ( 0 ) of 16 × 10 −6 , 64 × 10 −6 , and 256 × 10 −6 m 3 /cycle, corresponding to classes "A," "B," and "C" respectively, are adopted in the following sections. In practice, it must be noted that accelerations will be measured via accelerometers mounted on each axle, and it is not necessary to know all the vehicle properties listed in Table 1 in order to implement the proposed algorithm.

Derivation of the response of the contact point using a TF
In spite of the latest progress, an accurate definition of the displacement of the contact point is still very challenging to obtain, especially in scenarios with high vehicular speeds or complex vehicle models. For this purpose, this paper employs the TF, and axle accelerations,̈and̈, as the only input. These can be easily measured in a real vehicle by mounting accelerometers on the axles. The latest research shows that smartphone technology can be an interesting alternative to wired accelerometers in terms of cost savings and implementation (McGetrick et al., 2017;Zhao et al., 2019).
First, the equations of motion of the vehicle can be expressed as where superscripts . and .. are the first and second derivatives of a variable, representing the corresponding velocity and acceleration responses. Variableṡ,̇, , and are related to the DoFs of the vehicle (Figure 1) as per Equation (5): Next, by replacing Equation (5) into Equations (1)-(4) and applying the Fourier transform to these equations, Equations (6)-(9) are obtained: where , Θ, , , , and denote the Fourier spectra of the body mass displacement, , body mass rotation, , front axle displacement, , rear axle displacement, , and the contact point displacements of the front and rear axles, and . The symbols and represent the imaginary unit and circular frequency, respectively.
The Appendix demonstrates that the Fourier spectra of the contact point displacements of the front and rear axles, and , can be determined from the Fourier spectra of the axle responses alone, and , without the need for any knowledge of the Fourier spectra of the DoFs associated to the motion of the body, and Θ. These equations are the so-called theoretical TFs and require the vehicle properties to be known. The time series of the contact point displacements of the front ( ( )) and rear ( ( )) axles can then be re-constructed by utilizing the inverse fast Fourier transform (FFT) of the frequency spectra ( and ). In practice, the TF can be obtained experimentally by driving the vehicle over a measured road profile and comparing the Fourier spectra of the measured axle response F I G U R E 2 Comparison of true and estimated displacements of the contact point for the (a) front axle and (b) rear axle for an initial zero velocity, (c) front axle and (d) rear axle for a non-zero initial velocity to the Fourier spectra of the road profile. Once the TFs of the instrumented vehicle are calibrated periodically, changes in the vehicle properties over time should not affect the performance of the proposed methodology. It is worth highlighting that the instantaneous axle forces, whether measured or calculated from displacements, are the basis for the application of the methodology. The time history of these forces can then be exploited to warn about a possible change in the TF by comparing an average value of the applied forces to the static axle weights predicted by a weigh-in-motion system in the route of the vehicle. Figure 2 gives an example of the extraction of the contact point displacements for a vehicle moving at 5 m/s. The axle accelerations,̈and̈, are simulated using the VBI system defined in Section 2.1 and integrated twice to obtain the axle displacements, and . Then, the FFT is applied to the displacements and the resulting Fourier spectra, and , are inputted into Equations (A5a) and (A5b) defined in the Appendix. Finally, the inverse FFT of and is calculated to obtain the contact point displacements, and . The true and estimated contact point responses of the front and rear axles are compared for two different scenarios: initial zero velocity (constant drift as in Figure 2a,b) and non-zero initial velocity (linear drift as in Figure 2c,d). Here, the true response refers to the exact displacement of the contact point between the vehicle and the roadway in the theoretical transient simulation, and the estimated displacement is obtained indirectly from the integration of the simulated axle accelerations and the application of the TF. In reality, the exact displacement of the contact point is very difficult to measure, except for high-tech instrumented vehicles such as the traffic speed deflectometer. The latter incorporates a set of laser Doppler vibrometers on a straight beam to measure the relative velocity between the beam and the pavement surface (Keenahan & OBrien, 2018).
The patterns of the contact point displacements are successfully obtained through the proposed TFs, except for the first and last points of the record. However, a significant shift can be observed in the predicted displacements. This problem arises from the ICs being unknown and assumed to be zero when conducting the double integration at the start of the process. If the initial displacement is different from zero, the true displacement and the displacement estimated from accelerations will be shifted by an amount that remains constant for an initial zero velocity but that varies linearly for an initial velocity different from zero. Previous research obviates this issue by assuming ICs to be known or measurable through sophisticated instrumentation mounted on the vehicle (Harris et al., 2010;. Figure 3 distinguishes two stages as the vehicle moves along the road profile. The vehicle is assumed to travel at constant speed ( ) over an approach road ( ) before reaching the 15-m long bridge. The road profile is extended for an additional 10 m ( ) after the front axle of the vehicle leaves the bridge. In the absence of a specialized vehicle able to provide instantaneous measurements of the contact point displacements and the forces applied to the bridge, the following inputs are required:

METHODOLOGY
1. The acceleration of all axles and their location with respect to the bridge. Here, the acceleration responses that would be measured for a two-axle vehicle in the field are simulated as per Section 2.1.

F I G U R E 3
Stages of the algorithm 2. Static axle weights, stiffness of the tires, and TFs. The latter are derived theoretically for a half-car in Section 2.2. 3. Bending stiffness and mass distribution of the bridge deck.
The axle weights can be measured using a static weighing scale. The impact of noise on the measurements or inaccurate input values for the stiffness of the tires and the bridge deck are considered in Sections 4.4 and 4.5. The ICs are unknown for each axle at the start of the simulations. Section 3.1 describes the methodology for calibrating the ICs in order to solve the drifting problem identified in Figure 2. For the calibration stage, only part of the acceleration record is used, namely, the segment of the approach road preceding the bridge. Once the ICs have been established, Section 3.2 defines the remaining steps of the methodology before the bridge supports and the road profile can be assessed. The proposed algorithm relies on knowledge of the forces applied to the bridge during the vehicle crossing. Hence, the algorithm is mainly aimed at short-and medium-span bridges in which a single traffic event consisting of the test vehicle can be implemented without interruption to the traffic flow. Moreover, the impact of wind forces should be minimal for such bridges.

Calibration of the ICs
The determination of the displacement from a measured acceleration has been widely studied in the literature (Gindy et al., 2008;K. T. Park et al., 2005; J. W. Park et al., 2014). However, if the constants of integration are unknown, the double-integral will cause a drifting problem as shown in Figure 2. This subsection aims to address this problem through the joined optimization of two objective functions. The scenario with initial zero velocity resulting in the constant drift shown in Figure 2a,b is tack-F I G U R E 4 Calibration of the initial conditions (ICs). CP, contact point led first. Figure 4 shows a schematic of the step-by-step procedure for calibrating the ICs. For demonstration purposes, the following parameters are selected: = 80 m, = 5 m/s, the road profile is a class "A" with ( 0 ) = 16 × 10 −6 m 3 /cycle, and the portion of acceleration selected to establish the ICs corresponds to the vehicle traveling over the approach from 40 to 80 m, that is, before entering the bridge. The four main steps in Figure 4 are then implemented as follows: Step 1: Accelerations are measured at the front and rear axles and utilized as the inputs (̈,̈).
Step 2: First, the axle displacements ( , ) are obtained by applying double integration to the two inputs (̈,̈) for a wide range of combinations of the ICs. In this case, the initial displacements of the front and rear axles are varied taking the road roughness as a reference from −2 × 10 −2 m to 2 × 10 −2 m in increments of 5 × 10 −4 m. Therefore, = 81 × 81 = 6561 combinations of ( , ) are considered potential solutions. A larger solution space would need to be explored in case of a rougher or unknown profile. Second, the corresponding contact point displacements ( , ) of the front and rear axles are calculated following the inverse FFT of Equations (A5a) and (A5b) (the Appendix), and the associated contact point forces ( , ) are determined according to Equations (10a) and (10b): where , represent the forces applied to the bridge by the front and rear axles, respectively. Finally, two objective functions, 1 and 2 , are established separately as Equations (11a) and (11b): where = 1, 2, 3, . . . , represents the number of the combinations of ICs, refers to the longitudinal location at which the road height is evaluated (i.e., time-shifted so that the location is the same for both axles), refers to the instant in time at which the gross vehicle weight (GVW) is evaluated, and and are the total numbers of equally spaced samples being compared. In this example, the objective functions are assessed every 0.005 m for 1 and 0.001 s for 2 . The objective function 1 (Equation 11a) represents the sum of the least squares errors in the time-shifted difference between the contact point displacements of the front and rear wheels. The displacement of the contact point under the front wheel estimated at a time is compared to the displacement of the contact point under the rear wheel at a time + Δ . The time shift Δ is equal to the axle spacing divided by the vehicular speed. Given that the vehicle has not yet reached the bridge in this analysis, the contact point displacement is expected to be identical to the road profile regardless of the axle being employed. Therefore, the aforementioned time-shifted contact point displacement of the rear wheel should be identical to the front wheel for those scenarios where the difference in the ICs between the two axles is equal to the real one. For illustration purposes, Figure 5a shows the result of the time shift for the true ICs, which are = −0.859 × 10 −2 m and = −0.711 × 10 −2 m. The objective function 2 (Equation 11b) compares the difference between the true and estimated GVWs. The estimated GVW is calculated as the sum of the contact point forces ( , ) at each time step. The mean value of the estimated GVW should be very close to the true GVW as illustrated by Figure 5b corresponding to the same combination of ICs used in Figure 5a. In this case, the mean value of the estimated GVW is 176.522 kN, with a relative error of 0.0079%, compared to the true GVW given in Table 1.
Step 3: Based on the computations of the objective functions, the contour plots shown in Figure 5c,d are obtained. Then, two sets of potential solutions, namely, 1 and 2 , are identified, one from each objective function. These sets of solutions correspond to the ridge of each contour plot and are signaled as dash-dotted lines in the figure.
Step 4: If a linear regression is fitted to S 1 and S 2 , the following equations are obtained: = − 0.0015 (S 1 ) and = −1.9978 × − 0.0227 (S 2 ). Therefore, the intersection point that satisfies both equations is: = −0.858 × 10 −2 m and = −0.708 × 10 −2 m, for which the relative errors are −0.0724% and −0.3830%, compared to the true ICs of −0.859 × 10 −2 m (front) and −0.711 × 10 −2 m (rear). Figure 5e,f presents the comparison between the true and estimated road profiles using the optimal ICs for the front and rear axles, respectively. The root mean square errors (RMSEs) between the true and estimated road profiles in the [40, 80] m range used for calibration are 6.696 × 10 −6 m and 0.273 × 10 −6 m for the front and rear axles, respectively.
Figure 5e,f shows that the undesired drifting highlighted in Figure 2a,b has been removed successfully. A good agreement is found between the true and estimated contact point responses, although a zoom into the beginning of the record allows appreciation of the inaccuracies within the 0-10 m range framed in Figure 5e,f. Thus, a long approach would be desirable in a practical scenario to guarantee high accuracy by the time the vehicle reaches the bridge. In this paper, a compromise between computational time and accuracy is sought by considering = 80, 160, 240, 320 m for = 5, 10, 15, 20 m/s, respectively.
A more challenging scenario involving non-zero initial velocity leads to significant linear drifting as shown in Figure 2c,d due to initial velocities of 0.0044 and −0.0240 m/s for the front and rear axles, respectively. The number of unknowns increases from two ( , ) to four

F I G U R E 5
Step 2: (a) axle displacements and (b) gross vehicle weights (GVWs); Step 3: (c) 1 and (d) 2 ; Step 4: (e) front axle and (f) rear axle for an initial zero velocity; (g) displacement optimization and (h) velocity optimization for a non-zero initial velocity ( , ,̇,̇) and a structured optimization procedure, that is, cross-entropy (Dowling et al., 2012), is applied instead of the least-squares approach described before. Again, it is sought to find an optimal solution that will minimize the objective functions defined by Equations (11a) and (11b). Rather than using a discrete range of potential values for the unknowns, cross-entropy relies on treating the unknowns as random variables for which the statistical distributions are updated through successive iterations. In the first iteration, values of 0 m for the mean ( ) and 0.0067 m for the standard deviation ( ) are assumed for the distribution of the unknown initial displacements of the front and rear axles. For the distribution of unknown initial velocities of both axles, = 0 m/s and = 0.0500 m/s are considered. Figure 5g,h shows the evolution of the mean values together with the true value of each unknown. The shaded area surrounding the mean values in the figure corresponds to three times the standard deviation, that is, [ − 3 , + 3 ]. Small relative errors of 0.5222% and −4.8164% for the initial displacements of the front and rear axles, respectively, and 0.4229% and 0.2260% for their respective initial velocities are obtained.

Assessment of rotational stiffness at end supports and road profile of the bridge
Once the ICs have been obtained, the next stage is to assess the bridge boundary conditions and estimate the road profile over the bridge. Here, it is assumed that = 320 m, = 10 m, = 20 m/s and the road profile is class "A" with ( 0 ) = 1 × 10 −6 m 3 /cycle. Unlike Section 3.1, the 15-m road surface of the bridge is considered here. Figure 6a presents the estimated time-shifted contact point displacements when the two axles are traveling over the bridge. The differences between the two records of estimated contact point displacements can be attributed to the contribution of the bridge deflection that is not identical for each axle. After the contact point displacements and corresponding contact forces ( , ) are calculated according to the methodology described in Section 3.1, the associated bridge deflections ( , ) can be computed using Equations (12a) and (12b).

[ ] {̈} + [ ] {̇} + [ ] { } = [ ]
(13) where = 1, 2, 3, . . . , Q, and it represents each combination of a total of Q values of 1 and 2 that are considered F I G U R E 6 Comparison of estimations by front and rear axles of (a) displacement of the contact point between vehicle and bridge, (b) contact point forces, and (c) bridge deflections as a potential solution, y refers to a spatial increment in sequential order, and Y denotes the length of the compared vectors. Again, the displacements of the contact point under the rear wheel must be shifted in the time domain to refer to the same point of the road profile where the front wheel is located at a given point in time.
The pair of ( 1 , 2 ) values that leads to a minimum value among those combinations is identified, that is, = , and the associated rotational stiffness is taken as the solution. The bridge displacements, and , under the contact point of the front and rear axles corresponding to the solution are noted. Moreover, the road profiles estimated by each axle can be averaged to characterize the road profile on the bridge as defined in Equation (14).
where denotes the minimum value of the objective function, and the vector { } contains the estimated height of the road irregularities along the bridge span. For clarity, Figure 7 summarizes the steps described above in the form of a flowchart.

Static analysis
The rotational stiffnesses of the bridge, 1 and 2 , can take a wide range of values, with the lower bound ( 1 = 2 = 0) corresponding to a simply supported beam and the upper bound to a fixed-fixed beam ( 1 = 2 = ∞). However, low values of ( 1 , 2 ) lead to a behavior of the bridge that is barely distinguishable from that of a simply supported beam, and the same is true for high values of the rotational stiffness and a fixed-fixed beam. For that reason, this subsection carries out a static analysis of the impact of the rotational stiffness of the two supports on the maximum bending moment at the midspan of the bridge. The applied load is concentrated on the midspan with a value equal to the GVW. The rotational stiffness of both the left support, 1 , and the right support, 2 , is varied from 10 3.5 kN⋅m/rad to 10 8.5 kN⋅m/rad with an increment of 10 5.5 kN⋅m/rad, which results in 1000 samples being generated for each support and a total of 10 6 combinations. The relative difference between the maximum bending moment at the midspan of a simply supported beam and the scenario is calculated for each ( 1 , 2 ) combination and plotted in the contour plot of Figure 8a using a base-10 logarithmic scale. In particular, the relative difference is 0% for a scenario with 1 = 2 = 0 (moment = × 4 ) and that value increases to 50% for a scenario with 1 = 2 = ∞ (moment = × 8 ).

F I G U R E 7
Drive-by algorithm for assessing condition of road and supports. FEM, finite element model; FFT, fast Fourier transform; TF, transfer function For the bridge under investigation, the relative difference is smaller than 1% when 1 = 2 < 10 4.7 kN⋅m/rad, and larger than 49% for scenarios with 1 = 2 > 10 8.1 kN⋅m/rad, which are highlighted by dashdotted lines in the lower left-hand side corner and the upper right-hand side corner of Figure 8a, respectively. Figure 8b highlights the scenarios for which 1 = 2 (i.e., along the diagonal direction of Figure 8a). The solid curve represents the relative variation of the bending moment at midspan with the same change in rotational stiffness at both supports. Additionally, two horizontal lines at the bottom of Figure 8b indicate relative differences within the range 0%-1% for 1 = 2 < 10 4.7 kN⋅m/rad, and another two horizontal lines at the top of the figure highlight the differences ranged from 49% to 50% for F I G U R E 8 (a) 1 versus 2 versus relative difference of the bending moment at midspan, (b) subset of cases for which 1 = 2 , (c) true and estimated rotational stiffness, and accuracy versus vehicular speed and road class, for the (d) left support, and (e) right support 1 = 2 > 10 8.1 kN⋅m/rad. For computational efficiency, a range of ( 1 , 2 ) values between 10 5 and 10 7 kN⋅m/rad seems reasonable and is selected for investigation in the following subsections.

Accuracy in assessing boundary conditions
The performance of the methodology described in Section 3.2 is illustrated first using the following scenario: = 20 m/s, class "A" road with ( 0 ) = 16 × 10 −6 m 3 /cycle, 1 = 2 × 10 6 kN⋅m/rad ( = 10 6.3 ), and 2 = 8 × 10 6 kN⋅m/rad ( = 10 6.9 ). Figure 8c presents the contour plot of the objective function (Equation 13) versus combinations of the rotational stiffness of the two supports, where both 1 and 2 are varied from 10 5 to 10 7 kN⋅m/rad with a constant increment of 1.98 × 10 5 kN⋅m/rad. The solution space to be explored consists of a total of 51 × 51 = 2601 combinations ( 1 , 2 ). The final solution proposed by the algorithm corresponds to the combination giving the minimum value of Equation (13), namely, K 1 = 2.08 × 10 6 kN⋅m/rad and K 2 = 8.02 × 10 6 kN⋅m/rad. The estimated solution and its corresponding exact values can be visualized as a solid triangle and dash-dotted lines in Figure 8c. Clearly, the location of the minimum least squares error in the contour plot is a very close match to the real solution, representing relative errors of 4% for 1 and 0.25% for 2 .
A good performance of the algorithm is noticed in Figure 8d,e. For instance, the average accuracy for both supports is over 95% for scenarios with 5 m/s, even for a class "C" road. For 10, 15, and 20 m/s, the average accuracies are over 90% for all scenarios involving classes "A" and "B" but decrease for scenarios with a class "C." For a given speed, accuracy tends to worsen as road roughness increases. For example, the average accuracies at 20 m/s take values of 96.5% ("A"), 95% ("B"), 81.5% ("C") for 1 , and 97.4% ("A"), 94.8% ("B"), 83.7% ("C") for 2 . A similar trend applies to the corresponding values of standard deviation: 4.6% ("A"), 3.9% ("B"), and 8.7% ("C") for 1 and 2.7% ("A"), 5.6% ("B"), and 13.3% ("C") for 2 . The sensitivity of the methodology to road roughness is more strongly felt for the highest speeds over the roughest class "C." Relative errors in estimations of rotational stiffness are less than 20% (accuracies ≥ 80%) considering all scenarios, including 20 m/s and class "C." Nevertheless, it must be stated that there is an evident source of inaccuracy associated with the discretization level chosen for the solution space. Given that solutions are sought every 1.98 × 10 5 kN⋅m/rad to save computational time, true values of falling in the middle of two consecutive candidate solutions will lead to a maximum relative error (%) = 1.98 × 10 5 2 × K × 100. For values near the minimum ( 1 = 1.5 × 10 6 kN⋅m/rad) and near the maximum tested ( 2 = 10 × 10 6 kN⋅m/rad), maximum errors due to the discretization of ≈ 6.6% and ≈ 0.99%, respectively, are possible depending on the distance of the true value to the nearest candidate solution.

Accuracy in characterizing the road profile
In addition to the end rotational stiffness of the bridge, the algorithm can also assess its surface (Equation 14). This output is illustrated using the scenarios in which the vehicle travels at 20 m/s. For each of the road classes "A," "B," and "C," there are 10 estimations of road profiles, corresponding to the 10 different boundary conditions defined in Section 4.2. Figure 9a compares the differences between the true and the mean ( ) of the estimated road profiles for a class "A" road, which are indistinguishable at this scale. Figure 9b zooms in a short segment of the road, where the solid line and dash-dotted line represent the true and mean of the estimated road profiles, respectively. Additionally, the shaded area surrounding the dash-dotted line (mean values) corresponds to [ − 3 , + 3 ], which means that 99.7% of the computed road profiles are located within the highlighted region. The differences are sub-millimetric. Following ISO 8608 (ISO, 2016), a set of parallel straight dotted lines are plotted in Figure 9c to define the thresholds between different road classes in the spatial frequency domain. The power spectral density (PSD) for the true and estimated road profiles are included in the figure. The PSD of the estimated road resembles the true PSD, and correctly classifies it as class "A," except for the lowest spatial frequencies below 0.51 cycle/m due to the limited duration of the data sample.
Similar results are obtained for scenarios with road classes "B" and "C" as shown by Figure 9d,e,f and 9g,h,i, respectively. Unlike Figure 9c, Figure 9f,i locates the PSD of the estimated profiles within the first-second and secondthird straight parallel dotted lines, respectively, which successfully identify the classes "B" and "C." The accuracy in predicting the height of road irregularities at each point in time is assessed via the RMSEs, which compare the difference between the time series of true and estimated road profiles presented in Figure 9a,d,g. The RMSEs are 0.0364, 0.0616, and 0.1234 mm for classes "A," "B," and "C," respectively. The RMSE is approximately doubled when employing the immediately above rougher road class, but even so, it is very small and allows for providing a highly accurate longitudinal profile.

Accuracy allowing for modeling errors
The algorithm inputs concerning the bending stiffness of the deck ( ) or the tire stiffness of the vehicle ( , ) may be difficult to accurately know sometimes. In this subsection, these input parameters are added to the unknowns to be sought by the algorithm together with the rotational stiffness at the supports ( 1 , 2 ). As a result, the solution space expands considerably, and it makes necessary to adopt a more efficient optimization procedure than in previous sections. Cross-entropy is selected as the optimization procedure to resolve Steps 2 and 3 of Figure 7. This procedure selects samples of potential solutions based on initial distributions for the unknowns that are updated using the elite samples and improved in successive iterations. The cross-entropy procedure is implemented here following four steps: 1. Initial normal distributions are assumed by establishing values of the mean ( ) and standard deviation ( ) for all unknowns;  Figure 7. Then, the objective function is calculated according to Equation (13) for each trial beam . The objective function will tend to zero as the timeshifted road profiles predicted at the front ( − ) and rear axles ( − ) get closer; 3. Results of the trial FEMs generated in Step 2 are ranked according to the objective function . Then, the FEMs yielding the top 10% minimum values of the objective function are used as the basis for updating the mean and standard deviation of the normal distributions in the next iteration; 4. Steps 2 and 3 are repeated until convergence occurs, that is, the iterative process will be stopped when the change in the objective function is found to remain below 1% for 10 consecutive iterations. Once this stopping criterion is achieved, the mean values of the normal distributions in the last iteration are proposed as the real solution.
Cross-entropy-Steps 1 to 4-is applied five times to account for the randomness inherent to the procedure and to detect solutions driven into a local minimum. Therefore, the results in Sections 4.4 and 4.5 correspond to the average of five independent applications of the cross-entropy procedure.
Modeling errors regarding the bending stiffness of the bridge are considered first. Here, a more realistic, nonuniform stiffness profile is adopted for the VBI simulations rather than the constant stiffness used in previous sections. The bending stiffness of each beam element is randomly assigned from normal distributions assuming a spatial correlation between them. The mean stiffness value is taken from Table 2, the standard deviation is set at 2% of the mean stiffness, and the correlation coefficient is given by Equation (15): where 0 = 0.1 and = 1 m represent the constant component of correlation and the scale of fluctuation, respectively; corresponds to the distance between the centers of elements (Casero et al., 2020). The resulting stiffness profile exhibits a smooth variation and has an average value of 18.26 × 10 6 kN⋅m 2 .
The values of standard deviation are smaller than those observed in Figure 10d and suggest that the degree of uncertainty surrounding the tire stiffness is smaller than the one around the bending stiffness of the deck. Indeed, the FEM of the bridge model employed in the optimization procedure has been based on an average constant inertia (the varying stiffness profile has been modeled as one unknown), while the input acceleration was obtained from a transient simulation using the stiffness profile defined by Equation (15). Such a discrepancy between the model used to generate the simulated measurements and the built-in model part of the algorithm does not exist for Figure 10h.

Accuracy allowing for noisy input data
This subsection investigates the impact of noise on the performance of the algorithm by corrupting the axle displacements with the additive noise model defined by Equation (16).
F I G U R E 1 0 Estimation at 20 m/s of (a) rotational stiffness of the left support vs iteration (unknown bridge stiffness), (b) rotational stiffness of the right support vs iteration (unknown bridge stiffness), (c) bridge stiffness vs iteration on a class "A" road, (d) overall results for rotational stiffness allowing for unknown bridge stiffness, (e) rotational stiffness of the left support vs iteration (unknown tire stiffness), (f) rotational stiffness of the right support vs iteration (unknown tire stiffness), (g) tire stiffness vs iteration on a class "A" road; (h) overall results for rotational stiffness allowing for unknown tire stiffness where { } represents the noisy signal; { } and represent the noise-free signal and its associated standard deviation, respectively; and is the signal to noise ratio. (0, 1) is a normal distribution with zero mean and a unit standard deviation. The cross-entropy procedure is applied based on = 500 FEMs and two unknowns, 1 and 2 , with the same initial distributions assumed in Section 4.4. Figure 11a,b shows the impact of = 100, 50, 20, and 10 (i.e., relative errors in the measurements of 1%, 2%, 5%, and 10%) when the vehicle in Table 1 travels at 20 m/s over the bridge in Table 2 for the 10 combinations of 1 and 2 tested before. The impact of noise is hardly felt for small measurement errors such as = 100, when the accuracy is found to be 89.54% ("A"), 84.33% ("B"), and 77.44% ("C") for 1 and 92.12% ("A"), 88.29% ("B"), and 82.76% ("C") for 2 . The decrease in accuracy is not too significant for = 50 with a good road profile. However, the accuracy worsens as decreases from 50 to 20 and 10, more rapidly the rougher the road class. For instance, in the case of a class "A" profile, accuracy drops to 89.04% ( 1 ) and F I G U R E 1 1 Accuracy of rotational stiffness versus four noise levels when vehicle travels at 20 m/s for the (a) left support, and (b) right support. SNR, signal-to-noise ratio 90.45% ( 2 ) for SNR = 50, 72.11% ( 1 ) and 86.71% ( 2 ) for SNR = 20, and 40.24% ( 1 ) and 75.52% ( 2 ) for SNR = 10. When the road class is "C," accuracies of 56.52% ( 1 ) and 78.43% ( 2 ) for SNR = 50, 25.92% ( 1 ) and 68.62% ( 2 ) for SNR = 20, and 26.91% ( 1 ) and 67.40% ( 2 ) for SNR = 10 are obtained.

CONCLUSION
The scope of drive-by methods is often limited to very low speeds, very smooth profiles, or the use of special-ized vehicles. This paper has proposed an algorithm able to characterize the road profile and to quantify the rotational stiffness at the bridge supports by minimizing the difference in the road irregularities estimated from the displacements or accelerations measured by two axles of an ordinary vehicle at speeds of 5, 10, 15, and 20 m/s. Using noise-free data, the accuracies of the estimation of rotational stiffness have been over 90% for class "A" and "B" road scenarios. There has not been a noticeable deterioration in accuracy with higher speeds, except for class "C." Even so, accuracies using the class "C" road have remained over 80% for all tested speeds. Furthermore, a very good agreement has been found between the true and estimated road profiles on the bridge. The testing has been extended to scenarios considering modeling errors without a substantial loss of accuracy. Results have then worsened as increasing levels of noise were added, more significantly for the left support and for the class "C" road. Still, for the vehicle traveling at 20 m/s, accuracies for both supports have been kept over 77.44% with SNR = 100 considering the three road classes, over 78.01% with SNR = 50 and classes "A" and "B," and over 72.11% with SNR = 20 and class "A." The proposed methodology can also be applied to 3D VBI models by applying wheel forces instead of axle forces to the bridge, and by optimizing two road profiles under the wheels on each side of the vehicle, that is, the wheels on the left-and right-hand sides of the vehicle would be traveling under their own distinct road profile. Naturally, further uncertainties should be expected linked to the increased number of unknowns and complexity of 3D models. Similarly, future work is needed for validating the algorithm using real data from laboratory and field tests.

A C K N O W L E D G M E N T S
This research has received funding from Science Foundation Ireland (SFI)'s US-Ireland R&D partnership program under the proposal ID 16/US/I3277 titled MARS-Fly.
Open access funding provided by IReL.

R E F E R E N C E S
How to cite this article: Feng, K., Casero, M., & González, A. (2023). Characterization of the road profile and the rotational stiffness of supports in a bridge based on axle accelerations of a crossing vehicle. Computer-Aided Civil andInfrastructure Engineering, 38, 1935-1954 Additionally, by taking the third and fourth rows of Equation (A1), the following expressions are derived for the Fourier spectra of the contact point displacements of the front ( ) and rear ( ) axles: Substituting Equation (A3) into Equations (A4a) and (A4b), the last two equations can be re-written as Equations (A5a) and (A5b)