Bridge pier scour level quantification based on output-only Kalman filtering

Abstract

Soil scour near a bridge pier foundation is one of the leading causes of bridge failures. Traditional vibration-based scour monitoring methods are nearly incapable of quantifying scour levels using a single acceleration response without knowledge of excitation information. In this paper, a new output-only scour level prediction method is introduced via the integration of an unscented Kalman filter (UKF), random decrement (RD), and newly derived continuous Euler–Bernoulli beam addressing river water, traffic loads, and the linear and nonlinear behavior of sediments around the pier as external effects. We conducted extensive simulation studies and applied the method to an existing medium-span bridge with a steel girder and concrete deck in service in the province of Manitoba, Canada. These studies show that our proposed method can accurately estimate scour levels using only one accelerometer, which was validated with an independent bathymetric survey of the soil level at the pier foundation. Furthermore, three different linear and nonlinear soil profiles representing the soil behavior around the pier were also investigated as case studies in the scour level estimation process. The results confirm that a cubic function exhibits the best performance in quantifying the scour level around bridge piers.

Introduction

Scour is one of the leading causes of bridge failures, which happen suddenly without any warning. ¹ Until 2010, around 80% of 614,387 bridges in the US were built over water streams, and all of these bridges experienced floods frequently during their life spans. ^2,3 Failure analysis of 1000 bridges in the US determined that foundation scour was the first cause of failure for about 60% of bridges. ³ Therefore, developing a scour monitoring system is critical to preventing catastrophic bridge failures in transportation systems. The failures of the Bonnybrook Bridge in Calgary in 2013 and the I-10 bridge in southern California in 2015 are two recent notable cases of bridge failure in North America attributed to scour around their piers and abutments.

Many approaches have been proposed to detect scours or quantify the severity of scours around the foundations of bridge piers and abutments. These approaches can be categorized into two groups. The first group focuses on directly measuring scour depth using instruments installed in critical places under the water. The most straightforward technology uses buried switches, which detect the scour when it reaches the buried sensors. ⁴ Different types of buried sensors, called float-out devices, are available. These devices are buried under the riverbed upstream and send a signal as soon as they are unburied and tilted in a predefined orientation. However, there is no way to approve their operations. ^5,6 Another type of sensor is a magnetic sliding collar sensor, which takes advantage of a manual or automated gravity-based physical probe resting on the streambed and moves downward as scour develops. A sufficiently large gravity rod should be selected to avoid its penetration into the streambed in static state. A remote sensing element is used to send the signal of a change in the depth of the gravity sensor as an indicator of scour. ⁷

However, these sensors are only for one-time use, and the sensor installation process is challenging. Other research has used electromagnetic signals to track changes in physical properties in areas where soil and water interact. ⁸ With this technique, replacing the sensor in cases of damage or dysfunctionality is almost impossible. Placing a gravity rod protected by tubes beside the bridge pier is another approach. The rod moves down when the soil level goes down because of scour. ⁹ However, the protection of these tubes during the flood season is a massive challenge given the processes involved in their reinstallation and calibration. ⁹ Some other techniques, such as the use of fiber sensors, ¹⁰ piezo sensors, ^11,12 and electrical conductivity devices, ¹³ are also used for scour monitoring.

Electromagnetic sensors are another type of probe that has been installed around bridge piers and piles to track changes in the dielectric permittivity of the surrounding bridge foundation. ¹⁴ These wireless probes send information about the scour depth and hole refill in a real-time manner. ¹⁵ Still, the installation of these probes, their protection during the flood season, and the calibration of sensor response after flood seasons are significant challenges associated with these sensors. Recently, a technology called BridgeCat ¹⁶ has been used in the UK; it puts a sonar scanner, a high-resolution camera, and a digital altimeter on a hydraulic arm. The whole system is installed on a vehicle, and the riverbed height is measured using these sensors. This technique is costly and requires closing the bridge, causing traffic issues. ¹⁷ In general, direct approaches are expensive, and most of them are for one-time use only. Measurement-related uncertainty given the nature of sensors installed underwater, is also a major challenge for these kinds of methods. ¹⁸ Therefore, structural vibration-based methods have been developed to overcome the disadvantages of direct methods. ¹⁹

The second group of approaches use changes in dynamic modal properties (e.g., natural frequencies and mode shapes) caused by scour in piers and foundations by analyzing vibration measurements. ^20,21 The stiffness of a bridge pier depends on the mechanical properties of its structural components and their boundary conditions. Scour phenomena change the boundary condition of a bridge pier by reducing the height of the soil that supports the pier and pier foundation—consequently, the stiffness of the pier changes, which also changes the natural frequencies and mode shapes. Many researchers have reported that the natural frequencies of bridges are reduced because of increased scour depth. ^22,23 Therefore, to facilitate the use of changes in natural frequencies and mode shapes caused by scour, researchers have developed a bridge–soil interaction model using the finite element method (FEM). ²⁴ The core concept of this method is to compare the natural frequencies calculated from actual measurements with those calculated from the developed tuned models. As the comparison index, a frequency response function with the Winkler model considering the mass of soil surrounding the pier was used, and modal updating using an iterative method was conducted. ²⁵ The limitation of this method is that it requires a specific vehicle FEM model in the developed bridge–soil interaction model; therefore, its application is limited to problems subjected to a specific vehicle vibration only.

Xiong et al. ²⁶ used a flexibility matrix based on FEM to predict the scour level in a bridge pier. The flexibility matrix is a combination of natural frequencies and their mode shapes. Some experimental studies have been conducted to investigate the effects of scour on the change in mode shapes. ²⁷ Choudhury et al. ²⁸ developed a FEM model and updated it with the actual natural frequencies of the pier; they then used the updated model to estimate scour depth with empirical equations and a modal updating process. To apply these methods to real structures, the mode shapes are required, but extracting accurate mode shapes from measured vibrations is very challenging; it also requires the installation of many sensors. Furthermore, these methods require precise modal updating, and which mode shapes and natural frequencies should be used is unclear; extensive simulations and case studies are thus required to determine sensitive modal properties.

Palanisamy and Sim ²⁹ used an extended Kalman filter (EKF) to estimate the change in foundation stiffness caused by scour. They assumed white noise as an input excitation for the bridge pier, which is not applicable to a full-scale bridge system that is usually subjected to different ambient, periodic, and impact excitations. These excitations are not easily measurable, so the development of an output-only scour quantification method is a must. Another issue is that the EKF works based on first-order Taylor series linearization, and its performance decreases by increasing the order of nonlinearity in the system. ^30,31 Therefore, the method cannot be used for highly nonlinear soil-structure interaction functions. Furthermore, Palanisamy and Sim developed a simple pier model, but this model cannot be generalized for different types of bridge piers.

The literature above shows that traditional approaches can detect whether a scour has occurred using natural frequencies or can localize the scour by comparing the mode shapes of bridge systems. Some of these approaches require specific excitations, such as impact or cyclic loadings, or specific loading, such as a specific vehicle type, while most can quantify the level of scour to some degree using ambient vibration data.

This study aims to develop a comprehensive mathematical model aims at introducing an output-only model-based scour quantification method considering various types of boundary conditions, and soil profiles around the pier of a bridge system. The output-only Kalman filtering-based scour quantification approach is proposed using a continuous Euler–Bernoulli beam model of a bridge pier that considers various types of predefined soil–bridge interaction functions.

The advantages of the proposed method are as follows. First, it does not require many sensors; rather, only a single accelerometer is attached to the top of the pier close to the girder. Second, the proposed method does not require information on excitation sources in that it can consider any type of excitation sources, such as ambient vibration, impact, or vehicle-induced vibration. Third, the method can quantify the scour level of the bridge pier. Finally, the different boundary conditions of bridge piers and different soil properties can be considered.

This paper is organized as follows. In Proposed method, the continuous governing equation of motion of a bridge pier is developed based on Euler–Bernoulli beam theory. This model is used for the state-space formulation of Kalman filtering; then, we introduce the unscented Kalman filter (UKF). ³¹ The original UKF needs corresponding measured excitation information for estimation, but it is not viable for measuring bridge excitations. ³² Therefore, the random decrement (RD) technique ^33,34 is introduced as a data-driven approach that extracts the free vibration of a bridge pier, which is used as an input of the UKF to predict the scour level. ³⁵ In Case studies, two illustrative systems are selected to demonstrate the efficacy of the proposed method. At first, simulation studies are done to show the capability of the proposed method. As a full-scale experiment, the proposed method is validated with vibration data measured from an in-service bridge pier located in the province of Manitoba, Canada. A survey rod and an acoustic Doppler profiler (ADP) system are used to measure the soil level around the pier of interest and verify the estimated results of the proposed method.

Proposed method

To predict the scour level, a physics model-based scour estimation method is proposed. This method is composed of the RD method and the UKF. The role of RD is to extract free vibration responses from the ambient acceleration vibrations measured by the accelerometer installed on a bridge pier. The extracted free vibration response is fed into the UKF to estimate the soil level, an indicator of the scour level. The UKF estimates the variables of the state-space equation defined with a continuous Euler–Bernoulli pier model. Therefore, the proposed scour monitoring method is an output-only approach because it does not use any bridge excitation information (Figure 1). The proposed method’s details, including the RD, UKF, and newly developed continuous Euler–Bernoulli pier model methods, are described in the next sections.

Figure 1.

Overall procedure of the proposed method.

Random decrement method

RD is a data-driven algorithm for extracting the free vibration response of a dynamic system from its measured response and subjected to random or ambient vibrations. The technique averages out to a finite number ( $\tilde{N}$ ) of segments, $x_j\left(t_r+\tau \right),$ with length $\tau$ from the response measured when the signal amplitude reaches a specific trigger threshold point ( $\theta$ ) for the first time, as shown in Figure 2.

Figure 2.

Schematic view of the RD method.

$x_j\left(t\right)$ is the $\mathit{j}\mathrm{th}$ measured response of the system (e.g., it can be the acceleration response of a bridge pier), and the trigger threshold point ( $\theta$ ) set on the signal usually equals $\sqrt{2}\hspace{0.17em}\sigma$ ( $\hspace{0.17em}\sigma$ is the standard deviation of the time history response), as proposed by Brincker and Rodrigues. ³² When the signal reaches a trigger point, a segment of the signal is stored. The free vibration response, ${\delta }_j\left(t\right)$ , is then calculated by averaging these segments, as formulated in equation (1). The initial idea of the method was developed for a system subjected to random excitations. However, it also works for impact and ambient excitations. More details of the algorithm with the mathematical background are available in ³⁴

$${\delta }_j\left(t\right)=\frac{1}{\tilde{N}}{\displaystyle \sum _{r=1}^{\tilde{N}}{x}_j\left(t_r+\tau \right)}$$

(1)

where $\tilde{N}$ is the number of extracted segments from each measurement. Increasing the number of segments improves the accuracy of the extracted free vibration response.

Eulerian beam model for a bridge pier

In this paper, the UKF uses the extracted free vibration, ${\delta }_j\left(t\right)$ , in equation (1) by the RD algorithm to estimate the scour level around the pier foundation. To determine the level of scour, we track the changes in the soil stiffness and damping values of the underlying pier. The changes in these values are estimated by the UKF using a state-space equation, which is driven by the dynamic equation of pier motion. Therefore, we establish a dynamic equation of motion by developing a continuous Euler–Bernoulli beam model for a bridge pier in this section. A schematic view of the bridge pier is shown in Figure 3. $\omega$ is the transverse displacement of a beam (i.e., pier) $P$ is an axial force attributed to the mass of the bridge deck ( $M_D$ ), $k_s$ is the stiffness, and $c_s$ is the damping of the soil over time, $t$ . $k_b$ and $c_b$ denote the stiffness and damping of the rocker bearing in this study, respectively.

Figure 3.

Schematic view of Euler–Bernoulli beam model of a bridge pier with rocker bearings.

The bridge pier experiences lateral and vertical forces because of the crossing vehicles and the river water beating the pier. The bridge deck is supported in the center of the cross-sectional area of the pier, limiting the lateral displacement (i.e., the maximum allowed lateral displacement during the year in this rocker beaning is approximately 8 cm) on the top of the bridge pier; therefore, assuming the Euler–Bernoulli beam model for the bridge pier and neglecting the shear deformations are reasonable. This very general and well-established modeling approach has been adopted in numerous published journal articles in the analysis area. ^36-38 Considering an Euler–Bernoulli beam with a modulus of elasticity E, moment inertia I, cross-section area and density of the bridge pier $\rho$ , the beam partial differential equation is written as

$$EI\frac{{\partial }^4\omega \left(x,t\right)}{\partial x^4}+\rho A\frac{{\partial }^2\omega \left(x,t\right)}{\partial t^2}+P\frac{{\partial }^2\omega \left(x,t\right)}{\partial x^2}=F\left(x,t\right),$$

(2)

where $t$ is the time, and $EI{\partial }^4\omega \left(x,t\right)/\partial x^4$ is the relationship between the bending moment and deflection based on Euler–Bernoulli beam theory: $M\left(x,t\right)=EI{\partial }^2\omega \left(x,t\right)/\partial x^2$ . $M$ is the bending moment through the centroidal axis of the beam. $\rho A{\partial }^2\omega \left(x,t\right)/\partial t^2$ is the inertia force of the beam, and $P{\partial }^2\omega \left(x,t\right)/\partial x^2$ is given by the dead weight of the bridge decks and girders. The external force $F\left(x,t\right)$ is applied as water flowing around the pier

$$\begin{array}{l}F\left(x,t\right)=f_{\mathrm{water}}\left(x,t\right)\left[{\rm H}\left(x-a\right)-{\rm H}\left(x-b\right)\right]\\ \hspace{1em}\hspace{0.5em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-k_s\omega \left(x,t\right)\left[{\rm H}\left(x\right)-{\rm H}\left(x-a\right)\right]\\ \hspace{1em}\hspace{0.5em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-c_s\frac{\partial \omega \left(x,t\right)}{\partial t}\left[{\rm H}\left(x\right)-{\rm H}\left(x-a\right)\right]\\ \hspace{1em}\hspace{0.5em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-k_b\omega \left(x,t\right)\delta \left(x-L\right)-c_b\frac{\partial \omega \left(x,t\right)}{\partial t}\delta \left(x-L\right)\\ \hspace{1em}\hspace{0.5em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-M_D\frac{{\partial }^2\omega \left(x,t\right)}{\partial t^2}\delta \left(x-L\right),\end{array}$$

(3)

while ${\delta }_{nm}$ is the Dirac delta function, and ${\rm H}\left(x-x_0\right)=\{\begin{array}{cc}1& x\ge x_0\\ 0& x<x_0\end{array}$ is the Heaviside function. For any arbitrary function, $g\left(x\right),$ ${\displaystyle {\int }_0^{L}g\left(x\right){\rm H}\left(x-x_0\right)dx}={\displaystyle {\int }_{x_0}^{L}g\left(x\right)dx}$ , where $x_0$ is a constant value in the domain $\left[0,\hspace{0.17em}\hspace{0.17em}L\right]$ , and $L$ is the length of the pier.

Based on the assumed mode approach, ³⁹ $\omega \left(x,t\right)={\displaystyle \sum _{n=1}^N{X}_n\left(x\right)q_n\left(t\right)}$ is defined, where $N$ is the number of modes, $X_n\left(x\right)$ is the mode shape function of the beam for each mode, and $q_n\left(t\right)$ is a time-dependent generalized coordinate function. ⁴⁰ $\omega \left(x,t\right)$ is substituted in equations (2) and (3), respectively. Considering the orthogonality of the mode shapes, ${\displaystyle \int X_n\left(x\right)X_m\left(x\right)dx={\delta }_{nm}},$ the equation is integrated over the length of the beam and it is given as below

$$\begin{array}{l}\left(\rho A\hspace{0.17em}{B}_{mn}+M_D\hspace{0.17em}{X}_m^2\left(x=L\right)\right){\ddot{q}}_n\left(t\right)\\ \hspace{1em}\hspace{1em}+\left(c_s\hspace{0.17em}{D}_{nm}+c_b\hspace{0.17em}{X}_m^2\left(x=L\right)\right)\hspace{0.17em}\hspace{0.17em}{\dot{q}}_n\left(t\right)\hspace{0.17em}\hspace{0.17em}\\ \hspace{1em}\hspace{1em}+\hspace{0.17em}\hspace{0.17em}\left[EI{V}_{mn}\hspace{0.17em}+P\hspace{0.17em}{Q}_{mn}+k_s\hspace{0.17em}{D}_{nm}\hspace{0.17em}+\hspace{0.17em}{k}_b\hspace{0.17em}{X}_m^2\left(x=L\right)\right]\hspace{0.17em}{q}_n\left(t\right)\\ \hspace{1em}={\displaystyle {\int }_a^b{f}_{\mathrm{water}}\left(x,t\right)X_m\left(x\right)}\hspace{0.17em}dx,\end{array}$$

(4)

where

$$\begin{array}{l}{V}_{mn}={\displaystyle {\int }_0^L{X}_m\left(x\right)X_n^4\left(x\right)dx\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\hspace{0.17em}{B}_{mn}={\displaystyle {\int }_0^L{X}_m\left(x\right)X_n\left(x\right)dx,\hspace{0.17em}}\\ Q_{mn}={\displaystyle {\int }_0^L{X}_m\left(x\right){X^{″}}_n\left(x\right)dx,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{D}_{nm}={\displaystyle {\int }_0^a{X}_n\left(x\right)X_m\left(x\right)\hspace{0.17em}\hspace{0.17em}dx}.\end{array}$$

(5)

The following equation derived from equation (4) has a standard dynamic equation of motion form as follows

$$\mathbf{M}\ddot{q}\left(t\right)+\mathbf{C}\dot{q}\left(t\right)+\mathbf{K}\mathbf{q}\left(t\right)=\mathbf{F}\left(t\right),$$

(6)

where $M_{mn}=\rho \hspace{0.17em}A{B}_{mn}+M_D\hspace{0.17em}{X}_m^2\left(x=L\right),$ $C_{mn}\hspace{0.17em}\hspace{0.17em}=c_s\hspace{0.17em}{D}_{nm}+c_b\hspace{0.17em}{X}_m^2\left(x=L\right)\hspace{0.17em},\hspace{0.17em}$ and $K_{mn}=\hspace{0.17em}\hspace{0.17em}EI\hspace{0.17em}{V}_{mn}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+P\hspace{0.17em}{Q}_{mn}\hspace{0.17em}\hspace{0.17em}+k_s\hspace{0.17em}{D}_{nm}\hspace{0.17em}+\hspace{0.17em}{k}_b\hspace{0.17em}{X}_m^2\left(x=L\right)\hspace{0.17em},\hspace{0.17em}{F}_m\left(t\right)={\displaystyle {\int }_a^b{f}_{\mathrm{water}}\left(x,t\right)X_m\left(x\right)}\hspace{0.17em}dx\hspace{0.17em}.$ This set of equations can be numerically simulated based on different boundary conditions. This study assumes a fix-free boundary condition in which the corresponding mode shape function is introduced in the following section.

Kalman filtering

Kalman filtering ^41,42 is a type of Bayesian recursive filter assuming a Gaussian distribution for measured data and a least square loss function in a Bayesian inference framework. A linear Kalman filter works for a limited number of states showing linear behavior. An EKF ⁴³ is defined to compensate for the linear Kalman filter’s incapability for nonlinear systems identification. The EKF works based on the first-order linearization of the Taylor series, but it still does not estimate an actual value for systems with a high level of nonlinearity. Besides, for the first-order Taylor series, the derivative of the transition and measurement function needs to be calculated, which is hardly viable in some cases. ³⁰ The UKF is a derivative-free filter developed to address the drawbacks of the EKF. The UKF tries to approximate the distribution of the states of the dynamic equation of motion (equation (6)) using carefully chosen sigma points instead of approximating nonlinear functions. There have been some applications and improvements of the UKF because of damage identification problems in structural systems. ^44-46 In this paper, we use the UKF for scour detection and quantification around a bridge pier. We need to convert equation (6) into a state-space form considering the desired states for estimation purposes.

The transition and measurement functions of a dynamic system in a general form of the state-space model in the discrete domain are defined as follows

$${\mathbf{X}}_{k+1|k}=\mathrm{F}\left({\mathbf{X}}_k,{\mathbf{u}}_k,{\mathbf{z}}_k\right),$$

(7)

$${\mathbf{Y}}_k=\mathrm{H}\left({\mathbf{X}}_{k+1|k},{\mathbf{u}}_k,{\mathbf{v}}_k\right),$$

(8)

where $\mathrm{F}$ is the transition function that presents the first-order Markov process between the states, and $\mathrm{H}$ is a measurement function that defines the relationship between the states and the measurements. The functions $\mathrm{F}$ and $\mathrm{H}$ can take the form of any arbitrary nonlinear functions in the identification process. $k$ is a time step, ${\mathbf{X}}_k$ is a state variable vector, and, in the simplest form, it could be equal to ${\mathbf{X}}_k={\left\{\begin{array}{cc}\mathbf{q}\left(k\mathrm{\Delta }t\right)& \dot{\mathbf{q}}\left(k\mathrm{\Delta }t\right)\end{array}\right\}}^{\mathrm{T}}$ in equation (6). It could also include any other parameters desired to be estimated, such as the stiffness and damping of a structure. ${\mathbf{z}}_k$ is a discrete process white noise, $\mathbf{Q}$ is its covariance matrix, ${\mathbf{Y}}_k$ is the simulation measurement (simulation output) vector, ${\mathbf{v}}_k$ is a discrete measurement process noise vector, and $\mathbf{R}$ is its covariance value. ${\mathbf{u}}_k$ is the excitation vector, which can have any form of white noise and stationary excitation, periodic excitation, ambient vibration, or any combination.

To begin the UKF algorithm, we take advantage of some predefined points, sigma points, to approximate the desired parameter distributions of each iteration. Let us assume ${\mathbf{X}}_0$ is an initial value for the desired state vector and ${\mathbf{P}}_0$ is its covariance matrix; then, the equations (9) to (11) need to be solved to find the sigma points

$${\displaystyle \sum _{i=0}^{2n_s}{w}_i}=1,$$

(9)

$${\displaystyle \sum _{i=0}^{2n_s}{w}_i{\mathbf{\chi }}_i}={\mathbf{X}}_0,$$

(10)

$${\displaystyle \sum _{i=0}^{2n_s}{w}_i\left({\mathbf{\chi }}_i-{\mathbf{X}}_0\hspace{0.17em}\right){\left({\mathbf{\chi }}_i-{\mathbf{X}}_0\hspace{0.17em}\right)}^{\mathrm{T}}}={\mathbf{P}}_0,$$

(11)

where $\mathbf{\chi }$ is a sigma point, $w$ is the corresponding weight, and $n_s$ is the number of states. There is no unique solution for the above system of equations. The most common solution, sigma points, of the above three equations is given below

$${\mathbf{\chi }}_{\mathrm{0,0}}={\mathbf{X}}_0,\hspace{0.17em}\hspace{0.17em}{\mathbf{\chi }}_{i,0}={\mathbf{X}}_0\hspace{0.17em}+\sqrt{\left(\left(n_s+\lambda \right){\mathbf{P}}_0\right)}{\hspace{0.17em}}_i,\hspace{0.17em}\hspace{0.17em}{\mathbf{\chi }}_{i+n,0}={\mathbf{X}}_0\hspace{0.17em}-\sqrt{\left(\left(n_s+\lambda \right){\mathbf{P}}_0\right)}{\hspace{0.17em}}_{i+n_s},\hspace{0.17em}i=1\hspace{0.17em},\hspace{0.17em}2,\hspace{0.17em}\mathrm{....},\hspace{0.17em}{n}_s,$$

(12)

where $\lambda$ is a scaling parameter and $i$ is the index of the column vector number. With the sigma points, the weighting functions can be calculated, as shown below

$$w_0^m=\frac{\lambda }{n_s+\lambda },\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{w}_0^c=\frac{\lambda }{n_s+\lambda }+\left(1-{\alpha }^2+\beta \right),w_i^m=w_i^c=\frac{1}{2\left(n_s+\lambda \right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\hspace{0.17em}\dots ,2n_s,$$

(13)

where $\alpha$ determines the spread of the sigma points and is a positive number $\left(0\le \alpha \le 1\right)$ and $\lambda ={\alpha }^2\left(n_s+\kappa \right)-n_s$ . $\kappa$ is a secondary scaling parameter, which usually equals $3-n_s$ . $\beta =2$ is optimal for Gaussian distributions.

For each time step $k,$ the UKF procedure is presented in equation (12) to (21), and more detail can be found in ³¹

$${\overline{\mathbf{\chi }}}_{i,k+1|k}=F\left({\mathbf{\chi }}_{i,k},{\mathbf{z}}_k\right),$$

(14)

$${\mathbf{P}}_{k+1}^{XX}={\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^c\left[{\overline{\mathbf{\chi }}}_{i,k+1|k}-{\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^m{\overline{\mathbf{\chi }}}_{i,k+1|k}}\right]}{\left[{\overline{\mathbf{\chi }}}_{i,k+1|k}-{\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^m{\overline{\chi }}_{i,k+1|k}}\right]}^{\mathrm{T}}+{\mathbf{Q}}_k,$$

(15)

$${\mathbf{\Psi }}_{i,k+1|k}=H\left({\overline{\mathbf{\chi }}}_{i,k+1|k},{\mathbf{v}}_k\right),$$

(16)

$${\mathbf{P}}_{k+1}^{YY}={\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^c}\left[{\mathbf{\Psi }}_{i,k+1|k}-{\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^m{\mathbf{\Psi }}_{i,k+1|k}}\right]{\left[{\mathbf{\Psi }}_{i,k+1|k}-{\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^m{\mathbf{\Psi }}_{i,k+1|k}}\right]}^{\mathrm{T}}+{\mathbf{R}}_k,$$

(17)

$${\mathbf{P}}_{k+1}^{XY}={\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^c}\left[{\overline{\mathbf{\chi }}}_{i,k+1|k}-{\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^m{\overline{\mathbf{\chi }}}_{i,k+1|k}}\right]{\left[{\mathbf{\Psi }}_{i,k+1|k}-{\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^m{\mathbf{\Psi }}_{i,k+1|k}}\right]}^{\mathrm{T}},$$

(18)

$${\mathbf{{\rm K}}}_{k+1}={\mathbf{P}}_{k+1}^{XY}{\left({\mathbf{P}}_{k+1}^{YY}\right)}^{-1},$$

(19)

$${\mathbf{X}}_{k+1}={\widehat{\mathbf{X}}}_{k+1}^{-}+{\mathbf{{\rm K}}}_{k+1}\left(\mathbf{\delta }\left(k+1\right)-{\displaystyle \sum _{i=0}^{2n_s}{\omega }_i^m{\mathbf{\Psi }}_{i,k+1|k}}\right),$$

(20)

$${\mathbf{P}}_{k+1}={\mathbf{P}}_{k+1}^{XX}-{\mathbf{{\rm K}}}_{k+1}{\mathbf{P}}_{k+1}^{YY}{\mathbf{{\rm K}}}_{k+1}^T.$$

(21)

The term $\mathbf{\delta }\left(k+1\right)$ in equation (20) refers to the real output of the system, which would be the extracted free vibration of the system. ${\mathbf{P}}_{k+1}$ and ${\mathbf{X}}_{k+1}$ are the estimated covariance matrix and desired state vector, respectively. They would be used in the next iteration to calculate the new sigma points.

In the Kalman filtering approach, we need to define a process error for the system, compensating for the mathematical model’s uncertainty and other assumptions involved in the estimation process. For this specific problem, measuring the external excitations of the bridge pier is almost impossible. To deal with this issue, we introduced the RD, which extracts the free vibration response of the vibrating system. This means that the reference measured data aimed for in the UKF would be independent of the excitation—whatever it is. ³⁴ As the RD eliminates the variables related to the input excitations in the transition and measurement functions, and the performance of the UKF decreases with the increasing size of the state vector, the combination of the RD and UKF improves the estimation accuracy of the UKF for large state vectors. Moreover, the uncertainty related to the input estimation does not affect the unknown state estimation, which reduces the total uncertainty in the estimation process. The robustness of the RD-UKF combination has been demonstrated using extensive numerical and experimental studies. ^34,35 The error associated with the free vibration extraction process will compensate for the white noise parameter (v _k ) in the measurement function. Eventually, the UKF estimates the system parameters (e.g., stiffness and damping) more accurately.

Case studies

In the previous sections, the Euler–Bernoulli beam theory was employed to determine the general dynamic equation of a bridge pier. This model could be adapted for various types of bridge piers with different boundary conditions and external excitations, water levels, and soil heights around the pier. This section simulated a bridge pier vibration using MATLAB ⁴⁷ software and used the output vibration response for scour quantification. We also implemented our proposed approach on a real bridge pier located in Morris, Manitoba (Canada). Details on each of these studies are presented in the following subsections.

Numerical simulations

The mathematical model of a virtual bridge pier, equation (6), was used for simulation studies via MATLAB software. We considered a concrete pier with a rectangular cross-sectional area, $A=13\times 1.7{\mathrm{m}}^2$ with a length ( $L$ ) equal to 20 m, and a modulus of elasticity equal to $E=3\times {10}^{10}\hspace{0.17em}\mathrm{Pa}$ with a density of $3000\hspace{0.17em}\hspace{0.17em}\mathrm{kg}/{\mathrm{m}}^{\mathrm{3}}.$ Soil stiffness and damping values were $k_s=8\times {10}^8\mathrm{N}/\mathrm{m}$ and $c_s=8\times {10}^6\hspace{0.17em}\mathrm{N}/\mathrm{m}\mathrm{.s}$ , respectively. The soil height was assumed to equal 5 m, and the water level was assumed to be 5 m from the top of the soil. $P$ , the pier axial force, is equal to $M_D\hspace{0.17em}g,\hspace{0.17em}$ and $M_D$ and $\hspace{0.17em}g$ are $42\times {10}^3\hspace{0.17em}\mathrm{kg}$ and $9.81\hspace{0.17em}\mathrm{m}/{\mathrm{s}}^2,$ respectively. We used a rocker bearing, and the stiffness and damping of the bearing were set to $k_b=56\times {10}^6\mathrm{N}/\mathrm{m}$ and $c_b=2.4\times {10}^4\hspace{0.17em}\mathrm{N}/\mathrm{m}\mathrm{.s}$ ., respectively. ⁴⁸

The first three vibrational mode functions of the Euler–Bernoulli beam are considered for analysis. We apply a combination of a triangle impulse force with a 0.5 s duration, a random vibration with an amplitude of 0.01 g, and an El-Centro earthquake time history (dividing the total signal amplitude by 10) as ambient vibration, all of which are external excitations applied to the bridge; see Figure 4(a). The total vibration response of the model, with the excitation introduced in Figure 4(a), is shown in Figure 4(b).

Figure 4.

(a) Excitation applied to the simulating structure; (b) Total acceleration vibration response of the structure.

We assume a fixed-free boundary condition for the Euler−Bernoulli beam with its corresponding mode shape function that is ${\mathrm{X}}_N\left(x\right)=\sqrt{1/L}\left[\left(\sin {\beta }_{N}x-\sinh {\beta }_{N}x\right)-{\alpha }_N\left(\cos {\beta }_{N}x-\cosh {\beta }_{N}x\right)\right],$ where ${\alpha }_N$ is equal to $\sin {\beta }_{N}L+\sinh {\beta }_{N}L/\cos {\beta }_{N}L+\cosh {\beta }_{N}L,$ $N$ is a mode number, and ${\beta }_{N}L$ is, respectively, equal to 1.87, 4.65, 7.85, and 10.99 for the first four modes. ³⁹

We define the desired states and transfer equation (6) to the state-space form. We chose to predict the soil parameters as desired states to eventually predict soil stiffness and height; this is more convenient for quantifying the scour level. It is assumed that no damage exists in the bridge pier, itself. The desired state vector consists of

$$\mathbf{X}={\left\{\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cc}{x}_1& x_2\end{array}& x_3& {\dot{x}}_1& {\dot{x}}_2\end{array}& {\dot{x}}_3& k_s& c_s\end{array}\right\}}^{\mathrm{T}}={\left\{\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cc}{x}_1& x_2\end{array}& x_3& x_4& x_5\end{array}& x_6& x_7& x_8\end{array}\right\}}^{\mathrm{T}}.$$

(22)

Due to the nature of our proposed output-only scour quantification method, we simply had to extract the free vibration response of the bridge pier from its ambient vibration (Figure 5). The RD was applied to the measured acceleration response in Figure 4(b), with a trigger point of $1.9\sqrt{\mathrm{var}\left(\ddot{\mathbf{x}}\right)}$ and 3000 (six-second) samples for each segment. The sampling rate was 500 Hz, and based on this trigger point, the algorithm picked about 420 segments for averaging. Depending on these RD tuning parameters, the first 5 seconds of the extracted free vibration response and its corresponding fast Fourier response (FFT) are shown in Figure 5.

Figure 5.

(a) Free vibration response of the system; (b) FFT of the free vibration response.

The free vibration response reduces the uncertainty of the mathematic model due to unknown excitations. Considering equation (6) and desired state vector equation (22), the free vibration response, the discrete domain state-space transition, and measurement functions of the model are expressed as

$${\mathbf{X}}_{k+1|k}=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_1+\mathrm{\Delta }t\hspace{0.17em}{x}_4,\\ x_2+\mathrm{\Delta }t\hspace{0.17em}{x}_5,\end{array}\\ x_3+\mathrm{\Delta }t\hspace{0.17em}{x}_6,\\ \begin{array}{l}{x}_4+\frac{\mathrm{\Delta }t}{M_{11}}[-\left(x_8{D}_{11}+{\tilde{C}}_{11}\right)x_4-\left(x_8{D}_{12}+{\tilde{C}}_{12}\right)x_5-\left(x_8{D}_{13}+{\tilde{C}}_{13}\right)x_6\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\left({\tilde{K}}_{11}+D_{11}{x}_7\right)x_1-\left({\tilde{K}}_{12}+D_{12}{x}_7\right)x_2-\left({\tilde{K}}_{13}+D_{13}{x}_7\right)x_3],\end{array}\\ \begin{array}{l}{x}_5+\frac{\mathrm{\Delta }t}{M_{22}}[-\left(x_8{D}_{21}+{\tilde{C}}_{21}\right)x_4-\left(x_8{D}_{22}+{\tilde{C}}_{22}\right)x_5-\left(x_8{D}_{23}+{\tilde{C}}_{23}\right)x_6\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\left({\tilde{K}}_{21}+D_{21}{x}_7\right)x_1-\left({\tilde{K}}_{22}+D_{22}{x}_7\right)x_2-\left({\tilde{K}}_{23}+D_{23}{x}_7\right)x_3],\end{array}\end{array}\\ \begin{array}{l}{x}_6+\frac{\mathrm{\Delta }t}{M_{33}}[-\left(x_8{D}_{31}+{\tilde{C}}_{31}\right)x_4-\left(x_8{D}_{32}+{\tilde{C}}_{32}\right)x_5-\left(x_8{D}_{33}+{\tilde{C}}_{33}\right)x_6\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\left({\tilde{K}}_{31}+D_{31}{x}_7\right)x_1-\left({\tilde{K}}_{32}+D_{32}{x}_7\right)x_2-\left({\tilde{K}}_{33}+D_{33}{x}_7\right)x_3],\end{array}\\ x_7,\\ x_8\end{array}\right]=\mathrm{F}\left({\mathbf{X}}_k,\mathrm{\Delta }t\right),$$

(23)

$${\mathbf{Y}}_k=\left[\begin{array}{c}\begin{array}{l}\frac{1}{M_{11}}[-\left(x_8{D}_{11}+{\tilde{C}}_{11}\right)x_4-\left(x_8{D}_{12}+{\tilde{C}}_{12}\right)x_5-\left(x_8{D}_{13}+{\tilde{C}}_{13}\right)x_6\\ -\left({\tilde{K}}_{11}+D_{11}{x}_7\right)x_1-\left({\tilde{K}}_{12}+D_{12}{x}_7\right)x_2-\left({\tilde{K}}_{13}+D_{13}{x}_7\right)x_3],\end{array}\\ \begin{array}{l}\frac{1}{M_{22}}[-\left(x_8{D}_{21}+{\tilde{C}}_{21}\right)x_4-\left(x_8{D}_{22}+{\tilde{C}}_{22}\right)x_5-\left(x_8{D}_{23}+{\tilde{C}}_{23}\right)x_6\\ -\left({\tilde{K}}_{21}+D_{21}{x}_7\right)x_1-\left({\tilde{K}}_{22}+D_{22}{x}_7\right)x_2-\left({\tilde{K}}_{23}+D_{23}{x}_7\right)x_3],\end{array}\\ \begin{array}{l}\frac{1}{M_{33}}[-\left(x_8{D}_{31}+{\tilde{C}}_{31}\right)x_4-\left(x_8{D}_{32}+{\tilde{C}}_{32}\right)x_5-\left(x_8{D}_{33}+{\tilde{C}}_{33}\right)x_6\\ -\left({\tilde{K}}_{31}+D_{31}{x}_7\right)x_1-\left({\tilde{K}}_{32}+D_{32}{x}_7\right)x_2-\left({\tilde{K}}_{33}+D_{33}{x}_7\right)x_3]\end{array}\end{array}\right]\hspace{0.17em}=\mathrm{H}\left({\mathbf{X}}_{k+1|k},\mathrm{\Delta }t\right).$$

(24)

where ${\tilde{C}}_{mn}=c_b\hspace{0.17em}{X}_m^2\left(x=L\right)\hspace{0.17em},$ and ${\tilde{K}}_{mn}=EI\hspace{0.17em}{V}_{mn}+P\hspace{0.17em}{Q}_{mn}+\hspace{0.17em}{k}_b\hspace{0.17em}{X}_m^2\left(x=L\right).$ Using these state-space equations, the stiffness and damping of the soil around the pier are estimated using the RD-UKF technique. The soil stiffness and damping values for the constant soil profile are estimated in Figure 6.

Figure 6.

(a) Estimated stiffness; (b) Damping value of the soil.

The initial displacement $\mathbf{x}$ , and velocity, $\dot{x}$ , are assumed to be zero, and the initial stiffness value, $x_{7,\hspace{0.17em}\mathrm{int}.}$ , and damping values, $x_{8,\hspace{0.17em}\mathrm{int}.}$ , are $1\times {10}^8\mathrm{N/m}$ and $1\times {10}^6\mathrm{N}\mathrm{.s/m,}$ respectively. The initial covariance matrix is defined as $\mathrm{diag}\hspace{0.17em}\left[\begin{array}{cccc}{\mathbf{P}}_x& {\mathbf{P}}_{\dot{x}}& {\mathbf{P}}_k& {\mathbf{P}}_c\end{array}\right]$ , where $\hspace{0.17em}{\mathbf{P}}_x={\mathbf{P}}_{\dot{x}}=\mathrm{ones}\left(\mathrm{3,1}\right)$ , ${\mathbf{P}}_k={10}^{12.5}\text{, }$ and $\hspace{0.17em}{\mathbf{P}}_c={10}^{12}.$ The tuning parameter of the UKF is set as $\lambda =1$ in equations (12) and (13). The process and measurement noise covariance matrices are ${10}^{-20}\times {\text{I}}_8$ and ${10}^{-5}\times {\text{I}}_8$ , respectively. To investigate the performance of the proposed method, we assumed a 20% scour (i.e., change in soil height) in the model and used the previously tuned RD-UKF algorithm for scour quantification without any change in the tuning parameters of the algorithm. The technique quantifies the scour within an acceptable range, as shown in Figure 7.

Figure 7.

Estimated stiffness (a) and damping (b) values of the soil with the constant profile affected by scour.

Based on Table 1, the error percentage of the stiffness estimation is constant and equals 2%, which is a reasonably small value. The error associated with the damping value estimation is more considerable than the stiffness value, and about 11% and 5% for no scour, and 20% scour, respectively. Notably, this poor estimation of damping has been found in numerous previous studies on Kalman filtering-based system identification since there is a huge difference between the two stiffness and damping absolute values. ^49-51

Table 1.

Estimation of the soil stiffness and damping values using the introduced method.

Parameter		Soil stiffness, %	Soil damping, %
Initial value		$1\times {10}^8\mathrm{N}/\mathrm{m}$	$6.4\times {10}^6\mathrm{N}/\mathrm{m}\mathrm{.s}\mathrm{,}$
No scour	Clean answer	$8\times {10}^8\mathrm{N}/\mathrm{m}$	$8\times {10}^6\hspace{0.17em}\mathrm{N}/\mathrm{m}\mathrm{.s}$
	Estimation	$8.21\times {10}^8\mathrm{N}/\mathrm{m}$	$7.1\times {10}^6\hspace{0.17em}\mathrm{N}/\mathrm{m}\mathrm{.s}$
	Error (%)	2	11
20% scour	Clean answer	$6.4\times {10}^8\mathrm{N}/\mathrm{m}$	$7.1\times {10}^6\hspace{0.17em}\mathrm{N}/\mathrm{m}\mathrm{.s}$
	Estimation	$6.54\times {10}^8\mathrm{N}/\mathrm{m}$	$7.5\times {10}^6\hspace{0.17em}\mathrm{N}/\mathrm{m}\mathrm{.s}$
	Error (%)	2	5

We estimated the soil stiffness around the pier accurately; however, to estimate the soil height around the pier, we had to define the soil profile as a function of the soil height. The soil reaction equation, equation (25), represents the lateral soil response around the pier (soil reaction, $p,$ lateral displacement curve, $\omega$ ). Details on the equation and constant values are available in American Petroleum Institute (API) design code ⁵²

$$p=\mathrm{\Lambda }{p}_u\tanh \left(\frac{kx}{\mathrm{\Lambda }{p}_u}\omega \right),$$

(25)

where $p_u$ is the ultimate resistance of the soil at depth $x$ . Detailed information about $\mathrm{\Lambda },$ and $k\hspace{0.17em}\left(\mathrm{kN}\hspace{0.17em}.{\mathrm{m}}^{-\mathrm{3}}\right)$ -constant values related to the soil specification -are available in the API standard. The first derivative of the soil reaction equation, equation (25), is approximately equal to the soil stiffness because the amplitude of the lateral displacement of the soil around the pier is small and remains in the linear region ⁵³

$${\overline{k}}_s=\frac{\mathbf{d}p}{\mathbf{d}\omega }{|}_{\omega \simeq 0}=\mathrm{\Lambda }{p}_u\frac{\frac{kx}{\mathrm{\Lambda }{p}_u}}{{\cosh }^2\left(\frac{kx\omega }{\mathrm{\Lambda }{p}_u}\right)}{|}_{\omega \simeq 0}=kx,$$

(26)

where ${\overline{k}}_s$ in equation (26) has a unit of $\mathrm{kN}\hspace{0.17em}\hspace{0.17em}{\mathrm{m}}^{\mathrm{-}2}$ . ⁵³ It is multiplied by the soil height and results in the stiffness unit. Considering this simplification, we define soil profiles around the pier to evaluate the performance of the proposed method to estimate soil height using linear, quadratic, and cubic profiles.

The initial values and details of the state estimation process are mentioned in Table 2 for the three profiles for the sake of reproducibility for future users. As we will estimate the soil height using different profiles, $h_s=x_7$ will be defined as the soil height, and we will have only seven states in each case

$$\mathbf{X}=\left\{\begin{array}{cccc}\begin{array}{ccc}\begin{array}{cc}{x}_1& x_2\end{array}& x_3& {\dot{x}}_1\end{array}& {\dot{x}}_2& {\dot{x}}_3& h_s\end{array}\right\}\hspace{0.17em}=\left\{\begin{array}{cccc}\begin{array}{ccc}\begin{array}{cc}{x}_1& x_2\end{array}& x_3& x_4\end{array}& x_5& x_6& x_7\end{array}\right\}.$$

(27)

Then, equations (23) and (24) are rearranged as follows in equations (28) and (29)

$${\mathbf{X}}_{k+1|k}=\left[\begin{array}{c}{x}_1+\mathrm{\Delta }t\hspace{0.17em}{x}_4,\\ x_2+\mathrm{\Delta }t\hspace{0.17em}{x}_5,\\ x_3+\mathrm{\Delta }t\hspace{0.17em}{x}_6,\\ \begin{array}{l}{x}_4+\frac{\mathrm{\Delta }t}{M_{11}}[-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{11}+{\tilde{C}}_{11}\right)x_4-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{12}+{\tilde{C}}_{12}\right)x_5-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{13}+{\tilde{C}}_{13}\right)x_6\\ -\left({\tilde{K}}_{11}+D_{11}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_1-\left({\tilde{K}}_{12}+D_{12}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_2-\left({\tilde{K}}_{13}+D_{13}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_3],\end{array}\\ \begin{array}{l}{x}_5+\frac{\mathrm{\Delta }t}{M_{22}}[-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{21}+{\tilde{C}}_{21}\right)x_4-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{22}+{\tilde{C}}_{22}\right)x_5-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{23}+{\tilde{C}}_{23}\right)x_6\\ -\left({\tilde{K}}_{21}+D_{21}\hspace{0.17em}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_1-\left({\tilde{K}}_{22}+D_{22}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_2-\left({\tilde{K}}_{23}+D_{23}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_3],\end{array}\\ \begin{array}{l}{x}_6+\frac{\mathrm{\Delta }t}{M_{33}}[-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{31}+{\tilde{C}}_{31}\right)x_4-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{32}+{\tilde{C}}_{32}\right)x_5-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{33}+{\tilde{C}}_{33}\right)x_6\\ -\left({\tilde{K}}_{31}+D_{31}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_1-\left({\tilde{K}}_{32}+D_{32}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_2-\left({\tilde{K}}_{33}+D_{33}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_3],\end{array}\\ x_7\end{array}\right]=\mathrm{F}\left({\mathbf{X}}_k,\mathrm{\Delta }t\right),$$

(28)

$${\mathbf{Y}}_k=\left[\begin{array}{c}\begin{array}{l}\frac{1}{M_{11}}[-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{11}+{\tilde{C}}_{11}\right)x_4-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{12}+{\tilde{C}}_{12}\right)x_5-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{13}+{\tilde{C}}_{13}\right)x_6\\ -\left({\tilde{K}}_{11}+D_{11}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_1-\left({\tilde{K}}_{12}+D_{12}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_2-\left({\tilde{K}}_{13}+D_{13}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_3],\end{array}\\ \begin{array}{l}\frac{1}{M_{22}}[-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{21}+{\tilde{C}}_{21}\right)x_4-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{22}+{\tilde{C}}_{22}\right)x_5-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{23}+{\tilde{C}}_{23}\right)x_6\\ -\left({\tilde{K}}_{21}+D_{21}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_1-\left({\tilde{K}}_{22}+D_{22}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_2-\left({\tilde{K}}_{23}+D_{23}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_3],\end{array}\\ \begin{array}{l}\frac{1}{M_{33}}[-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{31}+{\tilde{C}}_{31}\right)x_4-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{32}+{\tilde{C}}_{32}\right)x_5-\left({\mathrm{f}}_{c_s}\left(x_7\right)D_{33}+{\tilde{C}}_{33}\right)x_6\\ -\left({\tilde{K}}_{31}+D_{31}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_1-\left({\tilde{K}}_{32}+D_{32}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_2-\left({\tilde{K}}_{33}+D_{33}{\mathrm{f}}_{k_s}\left(x_7\right)\right)x_3]\end{array}\end{array}\right]=\mathrm{H}\left({\mathbf{X}}_{k+1|k},\mathrm{\Delta }t\right).$$

(29)

${\mathrm{f}}_{k_s}\left(x_7\right)$ and ${\mathrm{f}}_{c_s}\left(x_7\right)$ are defined as the soil stiffness and damping functions, respectively. Different forms of these functions are tabulated in Table 2 for different soil profiles.

Table 2.

Initial and tuning parameters for different soil profiles.

Soil profile	Equation	${\mathbf{m}}_0$	${\mathbf{P}}_0$	$\mathbf{Q}$	$\mathbf{R}$
Linear	${\mathrm{f}}_{k_s}\left({\mathit{x}}_7\right)\to 1.6\times {10}^8{x}_7$	$\mathbf{x}=\dot{\mathbf{x}}=0$	$\hspace{0.17em}{\mathbf{P}}_{\mathit{x}}={\mathbf{P}}_{\dot{\mathit{x}}}={10}^0$	$\hspace{0.17em}{10}^{-9.16}\times {\mathbf{I}}_{7\times 7}$	$\hspace{0.17em}{10}^{-2}\times {\mathbf{I}}_{3\times 3}$
	${\mathrm{f}}_{c_s}\left({\mathit{x}}_7\right)\to 1.6\times {10}^6{x}_7$	$\hspace{0.17em}{x}_{7,\hspace{0.17em}\mathrm{int}.}=0\mathrm{m}$	${\mathbf{P}}_{x_7}={10}^{0.5}$
Quadratic	${\mathrm{f}}_{k_s}\left(x_7\right)\to 3.2\times {10}^7{x}_7^2$	$\mathbf{x}=\dot{\mathbf{x}}=0$	$\hspace{0.17em}{\mathbf{P}}_x={\mathbf{P}}_{\dot{x}}={10}^0$	$\hspace{0.17em}{10}^{-9.34}\times {\mathbf{I}}_{7\times 7}$	$\hspace{0.17em}{10}^{-5}\times {\mathbf{I}}_{3\times 3}$
	${\mathrm{f}}_{c_s}\left(x_7\right)\to 3.2\times {10}^5{x}_7^2$	$\hspace{0.17em}{x}_{7,\hspace{0.17em}\mathrm{int}.}=0\mathrm{m}$	${\mathbf{P}}_{x_7}={10}^{0.6}$
Cubic	${\mathrm{f}}_{k_s}\left(x_7\right)\to 6.4\times {10}^6{x}_7^3$	$\mathbf{x}=\dot{\mathbf{x}}=0$	$\hspace{0.17em}{\mathbf{P}}_x={\mathbf{P}}_{\dot{x}}={10}^0$	$\hspace{0.17em}{10}^{-11}\times {\mathbf{I}}_{7\times 7}$	$\hspace{0.17em}{10}^{-15}\times {\mathbf{I}}_{3\times 3}$
	${\mathrm{f}}_{c_s}\left(x_7\right)\to 6.4\times {10}^4{x}_7^3$	$\hspace{0.17em}{x}_{7,\hspace{0.17em}\mathrm{int}.}=0\mathrm{m}$	${\mathbf{P}}_{x_7}={10}^{0.849}$

The estimated soil height is shown in Figure 8(a).

Figure 8.

Results of soil height estimations.

To show the proposed method’s performance in quantifying a scour level, we applied a one-meter scour to the model and applied the previously tuned RD-UKF algorithm to scour quantification with no change in the method’s tuning parameters. The method quantified the scour within an acceptable range. Based on Figure 8(b), the accuracy of the scour quantification is improved with an increase in the order of the stiffness function. The estimated soil heights for soil profiles under the no-scour and 20% scour situations are tabulated in Table 3. Although the cubic and square functions have the same level of error, our simulation studies and the estimation trends in Figure 8 show that the cubic function displays a more stable and steady trend to reach a clean response compared to the two other functions. Hence, we proceeded with this cubic function in the next study.

Table 3.

Estimated soil height for different soil profiles.

	No scour	20% scour
Actual soil height (m)	5	4
Soil profile	RD-UKF estimation (m)
Cubic function	5.4	3.9
Square function	5.3	3.8
Linear function	4.8	4.5

To compare the data-driven approach to this proposed method, we have illustrated the frequency response spectrum of the model under scour in Figure 9. The first natural frequency of the pier (Figure 5(b)) has been changed from 3 Hz to 2.83 Hz. While 20% of the scour was assumed for the soil around the pier, the natural frequency was decreased by about seven percent. This means that the natural frequency-based scour monitoring technique is not reliable in scour quantification. The data-based approach can only detect the existence of scour. ¹⁸ Therefore, the introduced method will be beneficial to compensate for the drawback of current data-based methods for scour monitoring purposes.

Figure 9.

(a) Free vibration response under scour; (b) FFT of the free vibration response of the system under scour.

This section introduced the capabilities of the proposed method for scour quantification in a simulation study. In the next section, the proposed method is examined using a real large-scale bridge system.

Application to a large-scale bridge pier

In the previous section, we demonstrated the capability of the proposed method in scour quantification of a bridge pier in simulation studies. To investigate the method’s performance in an actual structure, the ambient vibration responses of a multi-span bridge under normal operation (located in Morris, Manitoba, Canada) were measured over the duration of a year, starting in October 2019. This five-span continuous bridge system is about 300 m long and was supported on 10 piers. Approximately 7500 vehicles travel across the bridge each month. We installed a uniaxial accelerometer to pier five to measure ambient vibrations on the pier caused by the beating river water, crossing traffic, wind, and pedestrians. We selected pier five because it is surrounded by water all year, and it is vulnerable to scour, as shown in Figure 10. Images of the bridge and the installed sensor on pier five are shown in Figure 11.

Figure 10.

Side view of the large-scale bridge system drawing.

Figure 11.

Views of the bridge, installed sensor, and data acquisition system.

A MEMS-Piezoresistive BDI model A1316 accelerometer was used to measure the vibration response of the pier. It has an amplitude range of ±2 g, a frequency range (min) of 0–300 Hz, and a voltage sensitivity of 1000 mV/g. A four-channel wireless data acquisition system was used for data communication (BDI model, STS-PRIME-0337), as shown in Figure 11(c).

The river water always flows around piers five and six during a year, but the water level changes during the flood season in the early spring and Fall. We recorded the ambient vibration responses of pier five for four months (February, July, August, and October). The soil stiffness around the pier was estimated each time. The ambient vibrations of pier five in February and August (2020), both measured around 3pm, are shown in Figure 12.

Figure 12.

Ambient vibration responses of pier 5.

The sampling rate of the data acquisition system was 100 Hz, and the data were recorded for about 30 min each time. Considering that the bridge length is about 300 m, and the maximum speed limit on the bridge is $50\hspace{0.17em}\hspace{0.17em}\mathrm{km}/\mathrm{h},$ it takes at least 21.6 s for a vehicle to cross the bridge. This time is useful for selecting a reasonable length of segments for the RD to extract the free vibration of the pier using the recorded ambient vibration. For RD implementation, $0.01\sqrt{\mathrm{var}\left(\ddot{\mathbf{x}}\right)}$ is selected as a trigger point, with a segment length equal to 30 s. Figure 13 shows the first 20 s of the extracted free vibrations of pier 5 in February and August 2020.

Figure 13.

Free vibration responses of pier 5.

As in the simulation study, the fixed-free boundary condition was considered for the modeling of the Euler–Bernoulli beam, as shown in Figure 3. Four rocker bearings are at the top of the pier, and they are connected to the girders, as shown in Figure 14. The stiffness and coefficient of friction of this bearing are $k_b=14\times {10}^6\hspace{0.17em}\mathrm{N}/\mathrm{m}$ and $\mu =0.04,$ respectively. ⁴⁸ The equivalent damping corresponding to the $\mu =0.04$ is calculated based on $C_b=\left(4\mu N_p\right)/\left(\pi {\omega }_p{X}_p\right)$ , where $N_p$ is the normal force, ${\omega }_p$ is the fundamental natural frequency, and $X_p$ is the maximum displacement. ³⁹ As we have four bearings on the top of the pier 5, and they are installed in parallel, the total stiffness and damping of the bearings are four times the stiffness and damping of each bearing.

Figure 14.

The four high-type rocker bearings of the pier 5.

We consider pier 5 as an Euler–Bernoulli beam with a rectangular cross-section that has an area of $A=14\times 1.5{\mathrm{m}}^2$ and a length of 20 m. The concrete modulus of elasticity is set to $E=15\times {10}^9\hspace{0.17em}\mathrm{Pa}$ with a density of $3000\mathrm{kg}/{\mathrm{m}}^{\mathrm{3}}.$ The initial soil stiffness and damping values and their corresponding covariance values for the UKF tuning are $k_s^0=1\times {10}^7\hspace{0.17em}\mathrm{N}/\mathrm{m},$ $c_s^0=1\times {10}^3\hspace{0.17em}\mathrm{N}/{\mathrm{m}}^{\mathrm{2}}$ , $\hspace{0.17em}{\mathbf{P}}_{k_s}={10}^{17}$ , and $\hspace{0.17em}{\mathbf{P}}_{c_s}={10}^{15}$ , respectively. Also, the covariances of process noise and measurement noise are defined as $\hspace{0.17em}{10}^{-10.8}\times {\mathbf{I}}_{8\times 8}$ and $\hspace{0.17em}{10}^{-3.2},$ respectively. The estimated stiffness value at each measurement time is shown in Figure 15. The equivalent stiffness of the sediments around the pier 5 was estimated using the mathematical model developed in Proposed method and the geometrical specification of the pier 5.

Figure 15.

Estimated soil stiffness around the bridge pier in different months.

Figure 15 shows the estimated stiffness value around $k_s=1.5-2.4\times {10}^8\mathrm{N}/\mathrm{m}\hspace{0.17em}$ during all the seasons. The difference between February and the other months should be addressed because of the thick ice around the pier during wintertime. The ice thickness can be about 60 cm and it increases the stiffness of the pier. ⁵⁴ Gonzales et al. ⁵⁵ have shown that the natural frequency of structures increases by up to 35% during winter. Similar behavior was observed in the bridge considered in this study, with a 30% increase in the estimated stiffness value during winter.

The simulation study results in Numerical simulations show that soil stiffness behavior leads to a more accurate estimation of the soil height around the pier using a cubic function rather than the two other functions. Therefore, the soil height will only be estimated using the cubic function. The initial soil height is assumed to be zero with $\hspace{0.17em}{\mathbf{P}}_{h_s}={10}^{0.8195}.$ The covariances of process and measurement noises are $\hspace{0.17em}{10}^{-18.6}\times {\mathbf{I}}_{7\times 7}$ and $\hspace{0.17em}{10}^{-10.3},$ respectively. To estimate soil height, the average estimated stiffness for each month in Figure 16 was used. The estimated soil height was 4–5 m for all months, except in February (6.2 m), as shown in Figure 16. The thick ice increases the pier stiffness, and eventually, the proposed method overestimates the soil height value in February.

Figure 16.

Estimated soil height around the bridge pier in different months.

To verify the estimated soil height in Figure 16, we excavated the river ice in February 2020 and measured the height of the sediment around the pier using the Canadian Geodetic Vertical Datum of 2013 (CGVD2013). ⁵⁶ We measured the water height and water level using the surveying rod equipment shown in Figure 17(a). The pier-footing elevation is available from the bridge drawings. The soil height around the pier was extracted, as shown in Figure 17(b).

Figure 17.

(a) The survey used for the soil height around the pier 5; (b) and schematic of the height of the soil around the pier 5.

Figure 17(b) shows that the orange-colored points denote the downstream heights, and the black points denote the upstream heights. The estimated soil height is almost in agreement with Figure 16; both show that the soil height was 4–5 m, with an average of 4.7 m. This height can be used for future analysis to create an initial estimate of the sediment profile. The thick ice in February caused a difference between the estimated soil height and the measured height in February. As the model does not consider the temperature effect, and we also assumed that all other components had fixed stiffness values, the RD-UKF intelligently increased the soil stiffness (height) to compensate for this difference.

In addition to the February measurement, for further confirmation, we verified our estimated soil height in August 2020 as well. A highly accurate ADP system was used to measure the river soil height around the pier, as shown in Figure 18. Grounded in the second field measurement, the soil height around the pier was 4–5 m, with an average of 4.67 m. The RD-UKF algorithm then accurately estimated the soil height around the pier during summer and fall, but not in winter. It should be emphasized that the estimated soil heights in the full-scale study (over a year) did not change considerably (no scour happened). It was also validated by river bathymetry tests in both winter and summer.

Figure 18.

An acoustic Doppler profiler (ADP) boat and measured soil height around the pier 5 in August 2020.

In addition to the model-based scour quantification, we have shown the FFT of the extracted free vibration response for the 4 months considered in this study in Figure 19. The first three natural frequencies of the Euler–Bernoulli model estimated by the RD-UKF are 2 Hz, 5.3 Hz, and 11 Hz, respectively. The natural frequencies of the measured vibration responses, as shown in Figure 19, showed good agreement with the calculated values for all 4 months. The magnitude of the natural frequencies in February appeared to be larger than the other 3 months. The second natural frequency (5.8 Hz) was measured only in February, and the measured natural frequencies in February were slightly higher than those of the other months. The third natural frequency in February was about 7% greater than the corresponding values in the other 3 months. In addition to the three modes discussed above, modes of 1.57 Hz, 2.6 Hz, and 3.57 Hz were also measured. It is assumed that these modes are related to the modes of the entire bridge system. ^57-59

Figure 19.

FFT of the free responses of the pier 5 in different months.

It is challenging to find the dominant peak triggered by the soil effect from Figure 19 to track its change for the purpose of scour monitoring. Another difficulty is determining whether the scour affects all of these peaks or only some of them. Moreover, it is impossible to define a clear relationship between the rates of change in the natural frequencies and soil stiffness because of the scour (Figure 9). However, Figures 8 and 16 show that the Kalman filtering-based method can estimate the soil height around the pier. In a nutshell, the introduced Kalman filtering-based method could help overcome the limitations of traditional data-driven methods in estimating the scour level, despite its complex mathematical equations, difficulties in finding the tuning parameters, and the initial information for the UKF.

Conclusion

The scour phenomenon around bridge piers and foundation systems is one of the most critical problems threatening the safety of bridge systems. The current vibration-based scour monitoring methods focus on the extraction of the natural frequencies and mode shapes of bridge decks and piers. These features can only be used to detect the scour around the pier, but they are not enough for scour quantification. In this paper, a model-based scour monitoring method was introduced to estimate the level of scour. To realize this, a continuous mathematical model was developed for a bridge pier based on the Euler−Bernoulli beam theory. The beam boundary condition was defined as fixed-free, with the soil layer at the fixed end and the load from the superstructure above at the free end. The model accepts various external loadings or combinations of them; therefore, it can be customized for a broad range of bridge pier types with different boundary conditions and loading conditions. The RD technique was employed to extract the free vibration response of the bridge pier, which was fed into the UKF for the desired state estimation. The pier model parameters were updated using RD-UKF and acceleration responses were measured at the pier coping. The introduced method was capable of scour quantification in the different soil profiles around the pier. The experimental study’s estimated soil heights were verified using two different bathymetry tests in the winter and summer seasons. The technical contributions of this study are summarized as follows:

(1) A novel output-only physics model-based scour level quantification method was introduced.

(2) The Euler−Bernoulli beam theory was used to develop the pier model, considering loads due to traffic, water flow, wind, and sediments around the pier as external loadings.

(3) The introduced method requires only one measurement of uniaxial acceleration response.

(4) The introduced method can quantify the level of scour accurately.

(5) The introduced approach was verified using a comprehensive simulation and field study that considered the linear and nonlinear behavior of the soil around the pier.

(6) The model can be generalized to a broad range of bridges with various types of loadings, boundary conditions, and soil-structure interaction functions.

(7) The experimental data estimation results were verified by two independent bathymetric surveys, and the results were well in agreement with actual measurements.

Developing more accurate mathematical models with different types of predefined soil-structure interaction functions and considering the winter ice effect and early spring ice collision effect in the model is the next step in our future study.