Moving Principal Component Analysis Based Structural Damage Detection for Highway Bridges Under Moving Vehicles

1. Introduction

Bridges are key component of transportation infrastructure that is crucial for a society to function well. They are under increasing pressure from continuing deterioration due to the ageing and operational and environmental loading as results of population growth and climate change. The real-time assessment of bridge's status is urgently required for structural safety. Structural health monitoring (SHM) provides a practical tool to assess and predict the structural performance of bridges. SHM is a multi-discipline field involving data collection by sensor networks and the diagnosis of structural health based on the collected data [1]. The collected data are processed to extract features that can be analysed through model-based or data-driven techniques to enhance the decision making for structural condition assessment [2-4]. The key factor that hinders a reliable bridge condition assessment is the operational and environmental variety, especially the traffic excitation and environmental temperature [5]. The bridge condition assessment under operational environments has drawn enormous attentions to researchers and practical engineers [6].

There are two main research trends to deal with the uncertain operational environments, e.g. elimination and utilization. Different methods have been developed to separate the temperature induced response from the structural response and predict the temperature-induced response [7]. Modal parameters are widely used for structural damage detection. The change of modal parameters due to the temperature variations is differentiated and eliminated during the process of damage detection [8,9]. An artificial neural network model for the temperature induced response is trained by large amount of data and the model is used to compensate the temperature effect [10,11]. The Kalman filter is adopted to eliminate the change due to environmental variations [12]. A two-stage procedure based on adaptive Mahalanobis-squared distance and one-class K-nearest neighbors (KNN) is used to remove the environmental variability [13]. A long period monitored data (e.g., for months) is used to construct the prediction model of structural response under operational environments for structural damage detection [14,15]. The above data-driven methods under changing temperature conditions reveal great potential for practical use with high robustness. These methods can identify, eliminate, or utilize the influence of the environmental temperature on the bridge response. The biggest obstacle to the performance of these methods comes from the availability of the data for all possible operational environments in bridge health monitoring. The performance of the regression models or neural networks has been affected by the high-quality long-term monitored data.

The utilization, known as ambient vibration survey (AVS), aims use dynamic responses under natural excitations such as traffic loads, wind and micro-tremors [16]. The AVS is an economical, convenient and time-saving approach as it does not require a special excitation on structures. For bridge engineering, the vehicle-bridge interaction is a common problem that will affect the accuracy of dynamic analysis. A great effort has been made to fully understand its mechanism for practical applications [17]. Based on these studies, the bridge information could be extracted from the traffic induced response using the vehicle-bridge interaction model. Recently, the moving load driven bridge damage detection has drawn much attention. The moving vehicle could be an exciter and a moving sensor for the bridge real-time monitoring [18,19]. The duration of a vehicle passing the bridge is short and the change due to the temperature effect is relatively small within this short-time period. The traditional moving load driven methods have a good interpretability as these methods are based on the accurate vehicle-bridge interaction model [20]. These methods are suitable to the controlled environmental conditions, and it is laborious and time-consuming to ensure their performances due to their limited anti-noise ability [21].

The principal component analysis (PCA) is one of most widely used data-driven methods. It is based on an orthogonal decomposition of the covariance matrix of the process variables along the direction that explain the maximum variation of the data. The PCA has been used to separate the structural damage features from that by environmental changes and then the environmental effects are eliminated during the process of bridge condition assessment [22]. The PCA has been used to reduce the size of feature vectors [23], eliminate the operational and environmental effects [24,25] and reduce the noise effect [26] for structural damage detection. With the traditional PCA, the whole time series of measured responses is analysed and it cannot reflect the instant status of the structure [27]. Also the computational cost of the covariance matrix will be increased with the number of measurements and the length of the time series [28]. Posenato et al. [29] proposed a data-driven method named moving principal component analysis (MPCA) for long-term structural monitoring. MPCA calculates the covariance matrix in a fixed length window. Lanata et al. [28] used MPCA to capture the correlation in each small cluster obtained by the K-means method for structural damage detection. Cavadas et al. [30] compared the performance of structural damage detection by MPCA and robust regression analysis (RRA). Zhu et al. [31] proposed a temperature-driven damage detection approach for bridges considering the temperature variations and traffic loads using MPCA. Jin et al. [32] compared the performance of the modal analysis-based damage detection method under changing temperatures using PCA or MPCA. The results revealed that MPCA has more anti-noise ability with lower false alarm rates than PCA. Zhang et al. [33] deployed MPCA for damage detection on a rigid frame bridge under seasonal temperature variations with space and time windows. Jin and Jung [34] compared modified MPCA with the static linear principal component analysis (SPCA) and incremental linear principal component analysis (IPCA) on Z24 bridge dataset using the k-means clustering with Linde-Buzo-Gray algorithm (KMC-LBG) and Bayesian information criterion (BIC) to choose the window size of MPCA. Nie et al. [27] developed a narrow moving window for MPCA and successfully detected the real-time change of an actual suspension bridge. The window's size is decided based on the cumulative contribution ratio with a convergent spectrum. As the above, MPCA has a great potential for early warning and damage detection as it can reveal the data’s inherent correlation and structure in detail. The length of the moving window significantly affects the effectiveness of the MPCA-based method. How to select the appropriate window length is still a big challenge. The loading forms and the causes of bridge anomalies (e.g., damage) are often inconsistent. This inconsistency causes the results of MPCA are not interpretable.

This paper is to develop a new method for bridge condition assessment under operational environments. A new damage sensitive feature (DSF) based on moving principal component analysis (MPCA) has been proposed. Two main environmental factors, e.g. the temperature variations and traffic loads, are considered. The numerical and experimental studies have been conducted to verify the performance of the proposed method. In the numerical study, the vehicle-bridge interaction system is modelled as a simply supported Euler–Bernoulli concrete beam subjected to a moving load. The bridge damage is simulated by a breathing crack. The temperature influence is considered as the variation of the beam's parameters and the thermal stress induced by the vertical temperature difference. The results show that the vehicle's mass and the temperature do not have the effect on the accuracy of structural damage identification. The changing pattern of the proposed DSF reflects different damage locations. The experimental study of a T-section reinforced concrete beam subjected to moving vehicles with different weights has been conducted in the laboratory. A new type of windows for MPCA is proposed to filter out the effect of measurement noise and vehicle-bridge interaction. The results show that the proposed method is robust and accurate to detect the crack damage of the bridge under operational environments.

2. Theory

Principal component analysis (PCA) and moving principal component analysis (MPCA) are briefly introduced in this section. The detail information could be found the textbook [35]. Then the main target for current existed data processing algorithm is discussed.

2.1. Principal Component Analysis (PCA)

As the foundation of MPCA, PCA is a statistical learning method that decomposes the original data into linearly uncorrelated vectors-principal components (PCs) according to the maximum variance's direction. At the same time, the new coordinate axis's direction should be orthogonal to all previous coordinate axis's directions. This transformation ensures that the selection of the coordinate axis's direction can make each PC contain as much information as possible. It is widely used in data compression and feature extraction.

In this study, the singular value decomposition (SVD) method is adopted for PCA. Considering a data (signal) matrix ${\mathbf{X}}_{m\times n}$ with the k-order PCs, suppose the matrix's rank is r that is greater than or equal to k. The matrix ${\mathbf{X}}_{m\times n}$ can be factorized according to truncated SVD as

$${\mathbf{X}}_{m\times n}\approx {\mathbf{U}}_k{\mathsf{\Sigma }}_k{\mathbf{V}}_k^T$$

(1)

where ${\mathbf{U}}_k$ is $m\times k$ matrix, ${\mathbf{V}}_k$ is $n\times k$ matrix and ${\mathsf{\Sigma }}_k$ is the diagonal matrix of the order $k$. ${\mathbf{U}}_k$, ${\mathbf{V}}_k$ are respectively taken from the first $k$ columns of matrices $\mathbf{U},\mathbf{V}$ that are singular vectors of the matrix ${\mathbf{X}}_{m\times n}$. ${\mathsf{\Sigma }}_k$ is obtained from the first $k$ diagonal elements of matrix $\mathsf{\Sigma }$ which is the singular value matrix of the matrix ${\mathbf{X}}_{m\times n}$.

Before using PCA, the data in the matrix ${\mathbf{X}}_{m\times n}$ needs to be standardized to obtain the matrix ${\mathbf{X}}^{\prime }$ as

$$x_{ij}^{\mathrm{*}}=\frac{x_{ij}-\overline{x_i}}{\sqrt{s_{ii}}}$$

(2)

where $\overline{x_i}=\frac{1}{n}\sum _{\mathrm{j}=1}^{\mathrm{n}}{x}_{ij},i=\mathrm{1,2},...,m$ and $s_{ii}=\frac{1}{n-1}\sum _{\mathrm{j}=1}^{\mathrm{n}}(x_{ij}-\overline{x_i}{)}^2,i=\mathrm{1,2},...,m$. $x_{ij}$ is the ${\mathbf{X}}_{m\times n}$’s element at row $i$ and column $j$. $x_{ij}^{\mathrm{*}}$ is the standardized matrix ${\mathbf{X}}^{\prime }$’s element at row $i$ and column $j$. $\overline{x_i}$ is the average value in row $i$ of the matrix ${\mathbf{X}}^{\prime }$.

After obtaining the standardized data matrix ${\mathbf{X}}^{\prime }$, the traditional PCA uses the eigenvalue decomposition of the ${\mathbf{X}}^{\prime }$'s correlation matrix or covariance matrix to calculate the principal component matrix. According to characteristic of SVD, the principal component matrix can be obtained.

$${\mathbf{X}}^{″}=\frac{1}{\sqrt{n-1}}{\mathbf{X}}^{\prime T}$$

(3)

$${\mathbf{X}}^{″}=\mathbf{U}\mathsf{\Sigma }{\mathbf{V}}^T$$

(4)

$${\mathbf{Y}}_{k\times n}={\mathbf{V}}^T{\mathbf{X}}^{\mathbf{\text{'}}}$$

(5)

where ${\mathbf{X}}^{″}$ is constructed for the truncated SVD. The row of ${\mathbf{V}}^T$ is the eigenvector of the $\mathbf{X}\mathbf{\text{'}}$'s covariance matrix. ${\mathbf{Y}}_{k\times n}$ is the principal component or score matrix. ${\mathbf{V}}^T$ is the weight or coefficient matrix.

2.2. Moving Principal Component Analysis (MPCA)

For time series signals, the analysis within a sliding window can better excavate the feature of the data. MPCA is a method that deploys PCA on the signal truncated into the window length instead of the full signal. The window is like a filter that slides and decomposes the original signals along the time axis to different PCs. These PCs obtained from MPCA have more significant features than the PCs obtained from PCA. Likewise, the window expands the eigenvalues into a series of changing curves along the time axis. These curves can better reflect the inner variations of the data over time. The main difference between MPCA and PCA is that Eq. (1) in MPCA is calculated over the signal of the window length, not over the entire signal. When the window moves from the initial time to the end, the eigenvalue vector and the PCs are calculated simultaneously at each corresponding time $t$, and saved consecutively in the eigenvalue matrix and the principal component matrix respectively according to the order of the time axis. The coordinate of the eigenvalue vector and the PCs obtained at each time $t$ from Eqs. (4) to (5) is corresponding to the location of the window's center on the time axis in the window's each movement. Figure 1 shows the schema of MPCA.

3. Numerical Study

In this section, the vehicle-bridge interaction model with the temperature change and the crack damage is built firstly. Dynamic responses of the bridge with different damage scenarios subject to a moving vehicle are simulated. The dynamic responses are analysed by PCA and MPCA for comparison. Then the damage sensitive feature is constructed based on the analysis. The influence of the crack’s location and the breathing crack mechanism on the damage sensitive feature is discussed.

3.1. Vehicle-Bridge Interaction Model

3.1.1. Finite Element Model of a Beam Bridge

The bridge is discretized as $N-1$ beam elements and $\mathrm{N}$ is the number of nodes. The element mass matrix and stiffness matrix of a beam element can be obtained as

$${\mathbf{M}}_e=\frac{\rho Al}{420}\left[\begin{array}{cccc}156& 22l& 54& -13l\\ 22l& 4l^2& 13l& -3l^2\\ 54& 13l& 156& -22l\\ -13l& -3l^2& -22l& 4l^2\end{array}\right]{\mathbf{K}}_e=\frac{EI}{l^3}\left[\begin{array}{cccc}12& 6l& -12& 6l\\ 6l& 4l^2& -6l& 2l^2\\ -12& -6l& 12& -6l\\ 6l& 2l^2& -6l& 4l^2\end{array}\right]$$

(6)

where $\rho ,A,l$ are the density, section's area and length of the beam element respectively.

Figure 2 shows the i^th beam element. The response at point $x$ and time $t$ can be obtained by the Hermite interpolation $H\left(x\right)$ from the node responses. The Hermite interpolation $H\left(x\right)$ is as follows,

$$\begin{array}{c}{H}_1\left(x\right)=1-3{\left(\frac{x}{l}\right)}^2+2{\left(\frac{x}{l}\right)}^3\\ H_2\left(x\right)=x{\left(1-\frac{x}{l}\right)}^2\\ H_3\left(x\right)=3{\left(\frac{x}{l}\right)}^2-2{\left(\frac{x}{l}\right)}^3\\ H_4\left(x\right)=\frac{x^2}{l}\left(\frac{x}{l}-1\right)\end{array}$$

(7)

The displacement at point $x$ and time $t$ is obtained as,

$${\omega }_e\left(x,t\right)={\mathbf{H}\left(x\right)}^T\mathbf{R}\left(t\right)=\left\{\begin{array}{cccc}{H}_1\left(x\right)& H_2\left(x\right)& H_3\left(x\right)& H_4\left(x\right)\end{array}\right\}\left\{\begin{array}{c}{\omega }_i\left(t\right)\\ {\theta }_i\left(t\right)\\ {\omega }_j\left(t\right)\\ {\theta }_j\left(t\right)\end{array}\right\}$$

(8)

where ${\omega }_e\left(x,t\right)$ is the displacement at point x and time $t$. ${\omega }_i\left(t\right),{\theta }_i\left(t\right)$ are the node displacement and rotation at the i^th node and time $t$ of the beam element. $x$ is the position and $l$ is the length of the beam element. $\mathbf{R}\left(t\right)$ a the node responses.

The strain at point $x$ and time $t$ can also be obtained [5] as

$$\mathsf{\epsilon }\left(\mathrm{x},\mathrm{t}\right)=-\mathrm{z}\frac{\begin{array}{c}\partial \mathbf{H}\left(x\right)\mathbf{R}\left(t\right)\end{array}}{\partial x^2}$$

(9)

where $z$ is the distance from the bottom to the neutral axis.

3.1.2. Equation of Motion for the Bridge Subjected to a Moving Vehicle

The bridge is modelled as a simply supported beam and the vehicle is modelled as a mass $m$, as shown in Figure 3. The bridge length is $L$. The vehicle is moving along the bridge at a constant speed $v$. The crack damage is considered in this study. $l{c}_d$ is the crack's location from the left support. The beam bridge is discretized into $N-1$ elements and $N$ is the number of nodes. Considering Rayleigh damping for the bridge, the motion of equation for the bridge subjected to a moving vehicle can be obtained as

$${\mathbf{M}}_b\ddot{\mathbf{R}}\left(t\right)+{\mathbf{C}}_b\dot{\mathbf{R}}\left(t\right)+{\mathbf{K}}_b\mathbf{R}\left(t\right)={\mathbf{H}}_m\mathbf{P}$$

(10)

where ${\mathbf{M}}_b$, ${\mathbf{C}}_b$ and ${\mathbf{K}}_b$ are the mass, damping and stiffness matrices of the bridge respectively. The $\ddot{\mathbf{R}}\left(t\right)$, $\dot{\mathbf{R}}\left(t\right)$ and $\mathbf{R}\left(t\right)$ are the node's acceleration, velocity and displacement response vectors respectively. ${\mathbf{H}}_m\mathbf{P}$ is the node's equivalent force vector induced by the moving mass. $\mathbf{P}$ is the equivalent resultant force vector induced by the moving mass. The interaction force between the bridge and the mass is $P(x\left(t\right),t)$ that can be obtained as

$$P\left(x\left(t\right),t\right)=m\left\{g-\frac{{\mathrm{d}}^2\left(\omega \left(x\left(t\right),t\right)\right)}{\mathrm{d}{t}^2}\right\}$$

(11)

where $m$ is the mass of a moving vehicle. ${\mathbf{H}}_m={\left\{\begin{array}{cccccc}0& \cdots & \mathbf{H}\left(\eta \right(t){)}_j^T& 0& \cdots & 0\end{array}\right\}}^T$ when $\left(j-1\right)l⩽\eta \left(t\right)⩽jl$. $H\left(x\right)$ is the Hermite interpolation as the beam element's shape function. The beam's deflection at point $x$ and time $t$ can be written as

$$\omega \left(x,t\right)=\mathbf{H}\left(x\right)\mathbf{R}\left(t\right)$$

(12)

where $\mathbf{H}\left(x\right)={\left\{\begin{array}{cccccc}0& \cdots & H(x{)}_j^T& 0& \cdots & 0\end{array}\right\}}^T$ when $(j-1)l⩽\eta \left(t\right)⩽jl$. Combining Eqs. (10)-(12), the equation of motion can be written as

$$\mathbf{M}\left(t\right)\ddot{\mathbf{R}}\left(t\right)+\mathbf{C}\left(t\right)\dot{\mathbf{R}}\left(t\right)+\mathbf{K}\left(t\right)\mathbf{R}\left(t\right)=mg{\mathbf{H}}_m$$

(13)

where $\mathbf{M}\left(t\right)={\mathbf{M}}_b+m{\mathbf{H}}_m\mathbf{H}\left(x\right)$, $\mathbf{C}\left(t\right)={\mathbf{C}}_b+2mv{\mathbf{H}}_m\mathbf{H}\prime \left(x\right)$,$\mathbf{K}\left(t\right)={\mathbf{K}}_b+m{v}^2{\mathbf{H}}_m{\mathbf{H}}^{″}\left(x\right)$. $\mathbf{H}\prime \left(x\right)$, ${\mathbf{H}}^{″}\left(x\right)$ are the first and second derivatives of the Hermite interpolation vector $\mathbf{H}\left(x\right)$. Eq. (13) can be solved using the Newmark-$\beta$ method. The parameters are: $\alpha =0.5,\beta =0.25$. The time step is 0.01 s. Then, the bridge response at point $x$ can be obtained by Eqs. (7) to (9).

3.1.3. Temperature Influence

The temperature's influence can be divided into two parts. The first one is the variation of the beam's parameters which will directly influence the beam's dynamic property. In the model, the thermal coefficients of the temperature influence on each parameter are listed in Table 1 [16]. $A_0$, $E_0$ and $I_0$ are the cross-section area, Young’s modulus and second moment of inertia at the reference temperature $T_0$ respectively. The influence of the temperature on the section's Young's modulus is linearly weights according to the steel and concrete's bearing ratio of the section's ultimate bending moment at the reference temperature. The section's expansion only considers the concrete's area growth.

$\mathrm{Steel}\left({/}^{\circ }C\right)$$\mathrm{Concrete}\left({/}^{\circ }C\right)$$A=A_0(1+{\alpha }_L\mathsf{\Delta }T{)}^2$${\alpha }_L$$1.2\times {10}^{-5}$$1.3\times {10}^{-5}$$E=E_0(1+{\alpha }_E\mathsf{\Delta }T)$${\alpha }_E$$-3.2\times {10}^{-4}$$-7.2\times {10}^{-3}$$I=I_0(1+{\alpha }_L\mathsf{\Delta }T{)}^4$${\alpha }_L$$1.2\times {10}^{-5}$$1.3\times {10}^{-5}$The second influence is induced by vertical temperature gradient. This temperature gradient is to influence the bending curvature of the beam. For a simply supported beam, the influence of the vertical temperature gradient can be taken as [36]:

$$r_1=-\frac{\alpha l_0\cdot \mathsf{\Delta }T}{2\sqrt{3}\cdot h}\cdot \frac{\lambda +1}{\lambda -1}$$

(14)

$${\mathrm{r}}_2=-{\mathrm{r}}_1$$

(15)

$$M_1=\left[1-\frac{1}{\sqrt{3}}\cdot \frac{\lambda +1}{\lambda -1}\right]\cdot \frac{\alpha EI\cdot \mathsf{\Delta }T}{h}$$

(16)

$$M_2=M_1$$

(17)

where $r_1$ and $M_1$ are the vertical temperature gradient induced rotation and moment at the beam's left support. $r_2$ and $M_2$ are the vertical temperature gradient induced rotation and moment at the beam's right support. $\mathsf{\Delta }T$ is the temperature difference between the top surface and bottom, and the positive reflects the top surface is warmer. $h$ is the cross-section's height. $EI$ is the beam's rigidity. $\alpha$ is the concrete's thermal expansion coefficient in Table 1. $\lambda$ is a constant and it is equal to $\frac{\sqrt{3}+1}{\sqrt{3}-1}$. $l_0$ is the beam's length. According to Eqs. (14) to (17), the vertical temperature gradient induced effect is obtained and considered as the boundary condition of the simulated beam.

3.1.4. Crack Damage Model

The beam's damage starts from initial microcracks and develops by many factors such as degradation, loads, temperature impact, etc. While expanding, these microcracks keep opening and closing due to external dynamic excitation. This phenomenon is known as the breathing crack and dominates the beam's crack behavior in incipient crack stage [37]. Since the prestress is widely applied, the crack in the prestressed concrete bridges will perform as the breathing crack.

In this study, the breathing crack is used to model the bridge damage [38]. The breathing crack is simulated as a rotational spring at the crack's location $l_c$. Figure 4 shows the beam element with a breathing crack. This element is considered as two undamaged beam segments connected by the proposed rotational spring. $EI$ is the undamaged beam's flexural rigidity and $l$ is the length of this element. The stiffness matrix of this element can be written as

$$\left\{\begin{array}{c}{Q}_i\\ M_i\\ Q_d^L\\ M_d^L\end{array}\right\}=\frac{EI}{l_c^3}\left[\begin{array}{cccc}12& 6l_c& -12& {6l}_c\\ 6l_c& 4l_c^2& -6l_c& 2l_c^2\\ -12& -6l_c& 12& -6_c\\ 6l_c& 2l_c^2& -6l_c& 4l_c^2\end{array}\right]\left\{\begin{array}{c}{w}_i\\ {\theta }_i\\ w_d^L\\ {\theta }_d^L\end{array}\right\}$$

(18)

$$\left\{\begin{array}{c}{Q}_d^R\\ M_d^R\\ Q_j\\ M_j\end{array}\right\}=\frac{EI}{{(l-l_c)}^3}\left[\begin{array}{cccc}12& 6(l-l_c)& -12& 6(l-l_c)\\ 6(l-l_c)& 4{(l-l_c)}^2& -6(l-l_c)& 2{(l-l_c)}^2\\ -12& -6(l-l_c)& 12& -6(l-l_c)\\ 6(l-l_c)& 2{(l-l_c)}^2& -6(l-l_c)& 4{(l-l_c)}^2\end{array}\right]\left\{\begin{array}{c}{w}_d^R\\ {\theta }_d^R\\ w_j\\ {\theta }_j\end{array}\right\}$$

(19)

where $w_i$, $w_j$, ${\theta }_i$, ${\theta }_j$ are the displacements and rotations at the i^th and j^th nodes respectively. $Q_i$, $Q_j$, $M_i$, $M_j$ are their corresponding transverse shear forces and moments. $w_d^L$, ${\theta }_d^L$, $w_d^R$, ${\theta }_d^R$ are the spring's displacements and rotations at the left and right hand sides of the joint, and $Q_d^L$, $M_d^L$, $Q_d^R$, $M_d^R$ are their corresponding shear forces and moments.

By the equilibrium and compatibility condition at the crack's location, the cracked beam's element stiffness matrix can be obtained.

$${\mathbf{K}}_d={\mathbf{K}}_1+{\mathbf{K}}_2{\mathbf{K}}_3^{-1}{\mathbf{K}}_4$$

(20)

Where

$${\mathbf{K}}_1=\frac{EI}{l^3}\left[\begin{array}{cccc}\frac{12}{{\delta }^3}& \frac{6l}{{\delta }^2}& 0& 0\\ \frac{6l}{{\delta }^2}& \frac{4l^2}{\delta }& 0& 0\\ 0& 0& \frac{12}{(1-\delta {)}^3}& -\frac{6l}{(1-\delta {)}^2}\\ 0& 0& -\frac{6l}{(1-\delta {)}^2}& \frac{4l^2}{(1-\delta )}\end{array}\right]{\mathbf{K}}_2=\frac{EI}{l^3}\left[\begin{array}{cccc}-\frac{12}{{\delta }^3}& \frac{6l}{{\delta }^2}& 0& 0\\ -\frac{6l}{{\delta }^2}& \frac{2l^2}{\delta }& 0& 0\\ 0& 0& -\frac{12}{(1-\delta {)}^3}& -\frac{6l}{(1-\delta {)}^2}\\ 0& 0& \frac{6l}{(1-\delta {)}^2}& \frac{2l^2}{(1-\delta )}\end{array}\right]$$

$${\mathbf{K}}_3=\left[\begin{array}{cccc}-\frac{12}{{\delta }^3}& \frac{6l}{{\delta }^2}& 0& 0\\ \frac{6l}{{\delta }^2}& S-\frac{4l^2}{\delta }& 0& 0\\ 0& 0& \frac{12}{(1-\delta {)}^3}& \frac{6l}{(1-\delta {)}^2}\\ 0& -S& -\frac{6l}{(1-\delta {)}^2}& S-\frac{4l^2}{(1-\delta )}\end{array}\right]{\mathbf{K}}_4=\left[\begin{array}{cccc}-\frac{12}{{\delta }^3}& -\frac{6l}{{\delta }^2}& 0& 0\\ \frac{6l}{{\delta }^2}& \frac{2l^2}{\delta }& 0& 0\\ 0& 0& -\frac{12}{(1-\delta {)}^3}& \frac{6l}{(1-\delta {)}^2}\\ 0& 0& -\frac{6l}{(1-\delta {)}^2}& \frac{2l^2}{(1-\delta )}\end{array}\right]$$

$$\delta =\frac{l_c}{l},S=\frac{K_{rd}{l}^3}{EI}$$

where $K_{rd}$ is the tangent stiffness which reveals the spring's instant rigidity. As shown in Figure 5, there is a crack opening at the edge of a rectangular section. 2$h$ is the cracked element's length and $b$ is the element's height. The crack is started at the center of the element's long edge. The rotational displacement due to the crack opening at the edge can be obtained using linear-elastic fracture mechanics as [39]

$${\theta }_{crack}=\frac{4\sigma }{E}S\left(\frac{a}{b}\right)$$

(21)

where $\mathsf{\sigma }$ is the applied stress of the whole cracked element induced by bending moment. The $\frac{\mathrm{a}}{\mathrm{b}}$ is the ratio of the crack's depth. For $\frac{h}{b}>2$, the $S\left(\frac{a}{b}\right)$ can be written as [39]

$$S\left(\frac{a}{b}\right)={\left(\frac{\frac{a}{b}}{1-\frac{a}{b}}\right)}^2\left\{5.93-19.69\left(\frac{a}{b}\right)+37.14{\left(\frac{a}{b}\right)}^2-35.84{\left(\frac{a}{b}\right)}^3+13.12{\left(\frac{a}{b}\right)}^4\right\}$$

(22)

The relationship of $\sigma$ and $M$ for the pure bending element can be written as

$$\sigma =\frac{6M}{{\mathrm{b}}^2}$$

(23)

where $M$ is the cracked element's bending moment. $b$ is the length of the section's short edge. Combining Eqs. (21) to (23), the instant tangent stiffness of the virtual rotational spring $K_{rd}$ is obtained.

$$K_{rd}=\frac{M}{{\theta }_{crack}}=\frac{E{b}^2}{24S\left(\frac{a}{b}\right)}$$

(24)

3.2. Results and Discussions

3.2.1. Numerical Simulation

The numerical model is validated by comparing with the results by Zhu and Law [40]. A simply supported beam with 50 m long, 0.5 m wide and 1 m high is used. The elastic modulus of the beam is $2.1\times {10}^{11}$ Pa and the density is $7860\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$. The moving force is 10 kN. The first six natural frequencies of the beam are listed in Table 2.

Figure 6 shows the normalized deflection at the mid-span. The sampling rate is 100 Hz. The deflection is normalized by $\frac{{\mathrm{F}}_0{\mathrm{L}}^3}{48\mathrm{E}\mathrm{I}}$, which is the static deflection when the force is at the mid-span. The scattered points are the analytical solution obtained by Zhu and Law [40] and the line curves are the numerical solution obtained by the proposed model in this study. $a/h$ is the crack depth ratio at the mid-span. For $v$=5 m/s, the number of elements is 13. For $v$=40 m/s, the number of elements is 7. The number of elements is consistent with the compatibility condition in Eq. (22). The results by the proposed model are close to those in the reference paper [40]. This numerical model will be used in the study.

3.2.2. Comparison of Results by PCA and MPCA

A comparison study is conducted in this section. The moving load is 10 kN. The velocity of the moving load is 5 m/s and the number of elements is 13. Other parameters are the same as that in Section 3.2.1. The beam is intact. A sudden and slight change of the moving mass is simulated to illustrate the sensitivity and reliability of MPCA. The acceleration responses are analysed by PCA and MPCA. As listed in Table 3, three cases are studied.

Figure 7 shows the first PCs for those three cases using PCA. There are no obvious changes for those cases in Figure 7. The results show that the PCA could not indicate out the changes of moving mass. As PCA analyses the whole time series, the instant status of the vehicle-bridge interaction system cannot be detected.

The eigenvalues are more sensitive to abrupt changes than PCs [27]. The first eigenvalue is taken for comparison in this study. The size of the window is 50 times the sampling interval. Figure 8 shows the result using MPCA. From Figure 8, there are clear changes at 5 s in Cases 2 and 3. The results show that the 1% mass variation can be detected immediately in the first eigenvalue curve using MPCA. For Case 2, the magnitude of the first eigenvalue is the same with Case 1 after 6 s as the mass returns back the original value. For Case 3, the first eigenvalue keeps approximate the same magnitude after 5 s since the mass of the moving load does not change after that.

3.2.3. The Effect of Damage Patterns

This section is to study the effect of the crack depth. The crack occurred at the mid-span. The crack depth is increased from 0% to 50% of the thickness. Other parameters are the same as Case 1 in Section 3.2.2. Figure 9 shows the pattern change of the first eigenvalue induced by the growing crack depth. The figure below is an enlarged view of the area in the red box of the figure above. From Figure 9, the results show that changes of the crack depth mainly affect the distance between two adjacent peaks in the first eigenvalue curve. When the damage occurs, the distance between two adjacent peaks is increased with the crack depth.

3.2.4. Orthogonality

MPCA can decompose the data into different coordinate axes in which each axis is orthogonal with other axes. For bridge SHM, the changes induced by structural damage need to be orthogonal with changes induced by other factors such as the vehicle's mass, the temperature, road surface roughness, etc. In this section, the orthogonality between the vehicle's mass and damage is studied. As listed in Table 4, six cases have been studied. The parameters Cases 1, 2 and 3 are the same as those in Section 3.2.2. Except for the crack depth, other parameters for Cases 4, 5 and 6 are the same as that in Section 3.2.2. Table 4 shows all six simulated scenarios.

Figure 10 shows results of six cases. The figure below is an enlarged view of the area in the red box of the figure above. Comparing Cases 1 to 3 with Cases 4 to 6, the results show that the influence of the damage is orthogonal with the influence of the moving mass. The influence of the damage is embodied in the distance between two adjacent peaks. Cases 4 to 6 show similar patterns which are different with the intact beam's patterns in Cases 1 to 3. The influence of the moving load's mass is embodied in the magnitude of each peak. The first eigenvalue represents the amount of information at the corresponding time in acceleration responses. For the intact beam, the increase of the moving mass leads to the decrease of the first eigenvalue. For the damaged beam, the increase of the moving mass leads to the increase of the first eigenvalue. For Cases 1, 2, 4 and 5, the peak magnitude after 6 s is the same and the results shows that there is correlation between the peak magnitude and the moving mass. It could be used for estimating the moving mass.

3.2.5. Temperature Influence

This section is to investigate the temperature effect using MPCA. As listed in Table 5, there are six cases. Two scenarios are simulated: the change for the uniform temperature of the whole beam and the change for the temperature gradient at the cross-section. The temperature gradient is slightly changed when the vehicle is passing the bridge. All the considered temperature changes are linear and the moving vehicle passes the bridge in a short time. There is no damage for Cases 1, 3, and 5 and the crack depth is 50% of the height for Cases 2, 4, and 6. The uniform temperature difference between Cases 3, 4 and Cases 5, 6 to simulate the temperature difference between day and night in mid-summer, and the corresponding temperature gradient is to simulate strong sunlight conditions. The other parameters are the same as Case 1 in Section 3.2.4.

Figure 11 shows the first eigenvalues of MPCA for different cases. The figure below is an enlarged view of the area in the red box of the figure above. From the figure, the first eigenvalue curve is changed with both the uniform temperature and the temperature gradient. Compared the results of Cases 1, 3 and 5, the time interval between two adjacent peaks for a case is increased with the uniform temperature and the amplitude of the first eigenvalue curve is increased with the temperature gradient. Compared the results of Cases 1, 3 and 5, and that of Cases 2, 4 and 6, the deviation caused by the temperature change is concentrated in a specific small range, and the change by the damage is much large as curve patterns.

3.3. Damage Sensitive Features

3.3.1. Observation

According to previous discussions, the crack depth mainly affects the time interval between two adjacent peaks in the first eigenvalue curve and the time interval is increased with the crack depth. The first eigenvalue curve has an obvious and identifiable pattern. There is a small sawtooth-interference information around the peak and the damage detection accuracy will be affected by these oscillations. To avoid the influence of these perturbation, a smoothing treatment is adopted. Following the PCA, the first eigenvalue curve of the principal component is needed. Figure 12 shows that the steadiest part which can reflect the main trend of the first eigenvalue is the two limbs of each peak. Thus, the mean line of the first eigenvalue curve is taken from the first eigenvalue curve. The midpoints of each pair of intersections in each peak are taken as the foundation of the DSF’s construction.

3.3.2. Construction

Figure 13 shows that the growth trend of the x-axis location of each midpoint is linear. The thumbnail in this figure is an enlarged view of the red boxed area in Figure 12. The crack depth's change will influence the inclination of each midpoint's x-axis location's growth trend line. In this section all cracks are occurred at the mid-span.

Thus, the gradient of the line in Figure 13 is used as the DSF. The numerical derivatives of each pair of discrete midpoints are obtained. The mean of all numerical derivatives in each line is calculated corresponding to its crack depth. The angle of each line is obtained by the mean's arctangent. Figure 14 shows the angle of each line corresponding to its crack depth.

3.3.3. Influence by Crack's Location

Figure 15 shows the crack location's influence. The result reveals that for each crack depth the crack located around the midpoint of the beam has a larger influence than the crack located near the beam's end. This phenomenon will be more obvious when the crack grows deeper. This result complies with the beam's dynamic analysis theory. In PCA, the first eigenvalue represents the variance of the first dimension. In this study, the first dimension is dominated by the acceleration changes caused by the moving loads. The vertical fluctuation and slow decline of the first eigenvalue reveal that the amount of the information (energy) brought by the moving load is gradually dissipated in the beam's response over time due to the beam's vibration. The existence of cracks will change the rate of this dissipation. Therefore, the increases in the distance between each midpoint represent these crack induced changes. Additionally, this DSF reflects the damage of the beam from the overall perspective, which is the so-called equivalent crack depth. Since cracks in actual structure are distributed near the damaged area of the beam, this DSF can provide a more realistic beam damage situation. When the crack depth is greater than 30%, the Maxwell-Betti reciprocal theorem is no longer valid due to the nonlinearity caused by the cracks. Thus, although the two damage locations with a distance of the 1/10 L on both sides of the beam’s midpoint are symmetrical in space, the damage extent reflected by these two DSFs is no longer consistent due to the directionality of the moving load on the time axis. In this case, the first eigenvalue is dominated by the moving load and the breathing crack. Within this range, the moving load passing through the crack earlier means that the beam has more time to dissipate, so the reflected degree of the damage will be slightly higher than the other places that the moving load passes the crack later. Therefore, the proposed DSF describes the damage extent of the beam from a dynamic perspective. In other words, it depicts beam's "rhythm". The traditional modal analysis describes the beam’s vibration from a static perspective. Since the analysis is based on the simply supported beam model, the magnitude of the first eigenvalue represents a measure of the maximum variance direction dimension. For a simply supported beam, the direction of the maximum variance of the measuring points' acceleration is in the direction of the gravity axis when a uniaxial moving load is passing. Therefore, MPCA can capture a continuous peak fluctuation of the same magnitude in the first eigenvalue curve.

4. Experimental Investigation

A laboratory study has been conducted in this study. The experimental data from Zhu and Law [40] have been used. The first eigenvalue curve has been obtained using the proposed method in Section 3. The Gaussian window is adopted to reduce the measurement noise and the vehicle-bridge interaction. The selection of the window parameters is also discussed. The results of the obtained damage sensitive feature on this laboratory data using the proposed method in Section 3 are presented and discussed.

4.1. Experimental Setup

Figure 16 shows the experimental setup. The cross-section of the concrete beam is shown in Figure 16 (a). As shown in Figure 16 (b), the whole experimental beam is composed of three T-section reinforcement concrete beams: the front beam, the main beam and the tail beam. The front and tail beams are 4.5 m long each. The main beam is 5.0 m long. The gaps between these three beams are 10 mm. An electric motor is used to pull the vehicle along the beam at a speed of approximately 0.5 m/s. The vehicle's axle spacing is 0.8 m, and its wheel spacing is 0.39 m. There are two vehicle models with different weights in this study. The whole weight of the first vehicle model (without elastic spring) is 10.60 kN with the front axle load 5.58 kN and the rear axle load 5.02 kN. The whole weight of the second vehicle model (with elastic spring) is 15.00 kN with the front axle load 6.20 kN and the rear axle load 9.00 kN. Since the mass of the whole concrete beam is 1050 kg, the weight ratios between the vehicle and the beam bridge for these two vehicle models are 1.01 and 1.43 respectively. Figure 16 (c) and (d) are the photos taken during this experiment. Figure 16 (c) shows the vehicle model passing through the beam and Figure 16 (d) shows the large damage case being generated.

As shown in Figure 16 (b), seven accelerometers are evenly installed along the beam at the bottom surface. Thirteen photo-electric sensors are distributed on the lead and main beams with 0.56 m spacing to measure the vehicle's moving velocity. The third and thirteenth photo-electric sensors are installed at the entry and exit points of the main beam respectively. The INV300E data acquisition system is used to obtain the response data. The duration for each test is 30 seconds and the sampling frequency is 2024.292 Hz.

A three-point load system is used to create the damage. The small damage is created by deploying the load at the 1/3 L from the beam's right support as marked in Figure 17 (a). The load is gradually added with a 2 kN increment. Several tensile cracks are obviously appeared on the beam rib when the load reaches 36 kN. When the load is 50 kN, the largest crack at the beam's bottom is measured as a 0.10 mm width. This crack is located close to the loading point but on the span inside with a 213 mm depth and a 760 mm wide crack zone visually. The beam is unloaded after the load is kept for 30 minutes. Then the crack at the beam's bottom decreases to a 0.025 mm width and closes partly. These descriptions are referred to as the small damage case.

For the large damage case, a 50 kN load is first loaded at the 2/3 L of the beam from the right support by the three-point load system. This produces a crack pattern similar in the extent and magnitude to the existing crack zone at the 1/3 L. After that, a four-point load system is used for a further loading as marked in Figure 17 (b). The final total load was 105 kN without the main reinforcement yielding and the largest crack is located near the beam's midpoint with a 281 mm depth. This crack is a 0.1 mm width at the beam's bottom when the load is 105 kN. When the beam is unloaded after keeping the 105 kN static load on the top surface for 30 minutes, this crack width reduces to 0.038 mm. The crack zone is 2371 mm long.

4.2. The Gaussian Window

Figure 18 shows the first eigenvalue curves under different moving vehicle models. Figure 18 (a) shows the results for the undamaged beam, and the beams with the small and large damage under a 10.6 kN moving vehicle. Figure 18 (b) shows the results under a 15 kN moving vehicle. The results show that the pattern of the first eigenvalue curve is severely affected due to the existence of measurement noise and the vehicle-bridge interaction. Measurement noise and the vehicle-bridge interaction are two main influencing factors for the first eigenvalue curve. Comparing Figure 18 (a) with Figure 18 (b), since the latter uses a 30% heavier vehicle, its curve is less affected by those influencing factors. The window length in these two cases is randomly selected with 54 sampling intervals for a comparison and illustration. Therefore, following this idea, if a smaller window can be used, more components caused by the moving load can be extracted. However, due to the limitation of PCA algorithm, the window's length cannot be less than the number of input signal's channels. Additionally, a small window length may lose the important information caused by the moving loads. Thus, the Gaussian window is proposed.

The Gaussian window draws on the idea of regularization. Figure 19 shows a Gaussian window (𝜎 (0,5)) when the normal window's length is 50 sampling intervals. It is equivalent to add a penalty term: the farther the time is from the current moment, the smaller the impact on the current moment. The greater the weight in the middle, the deeper the consideration of instantaneous effects. The influencing factors' impact on the time axis is diffuse, so their effects can be significantly reduced by the Gaussian window and the proportion of the vehicle excitation in this window is magnified at the same time. In other words, PCA is a multi-channel data processing method and the additional normal moving window expands the eigenvalue along the time axis. Through the expansion of eigenvalues, not only the same overall movement trend of each measuring point on the bridge is displayed along the time axis, but also each interference factor is expanded on the time axis. Therefore, the Gaussian window realigns this expansion on the time axis again. This minimizes the effect of interference factors on the detection results.

4.3. Parametric Study

4.3.1. The Effect of the Window Length

This section will discuss the selection of the window's parameters. The first is to investigate the influence of the window length on the first eigenvalue curve. In this section, the responses of the undamaged experimental beam under the 10.6 kN moving load are used. Figure 21 shows the first eigenvalue curves obtained under three commonly used lengths of the window. The image below is an enlarged view of the area in the red box of the image above. It is not recommended to use a length more than 200 times the sampling interval since a large window will greatly increase the MPCA's computational cost. From the first eigenvalue curves under different window lengths, we can see they have the same shape. Their only difference is the magnitude. The window length will not influence the damage detection since the damage influences the distance between each pair of adjacent peaks instead of the magnitude. Therefore, the second contribution of the Gaussian window is that it simplifies the problem of the window length’s selection. The Gaussian window transforms the selection of the window length into the hyperparameter 𝜎's selection of itself and allows the use of a small window to reduce the amount of calculation.

4.3.2. Hyperparameter of the Gaussian Window

The choice of the hyperparameter 𝜎 will affect the Gaussian window's attention degree to the current moment and the tolerance of the influencing factors. If the value is too small, the anticipated effect cannot be achieved. If the value is too large, all the information carried in the data will be destroyed (only taking the value at the current time t into account). The numerical model is used to find the optimal value. The parameters of the simulated beam are same as Case 1 in Section 3.2.2 and the length for both normal and Gaussian windows is 50 sampling intervals.

Figure 22 (a) shows the first eigenvalue curve processed by the Gaussian window (𝜎 (0,5)). Comparing with the curves in Figure 12, the curve obtained by the Gaussian window is much smooth. The Figure 22 (b) is the enlarged view of the area in the red box of Figure 12 and Figure 22 (a) correspondingly. Figure 22 (b) reveals that the Gaussian window can erase the distortions near the peak and trough areas of the first eigenvalue curve. Except for the window's type, the damage conditions and other parameters for obtaining the DSFs in Figure 12 and Figure 22 (a) are the same. Figure 22 (c) shows that the DSFs obtained by the normal window and the Gaussian window are almost the same. They both have the same gradient which could be an indicator of the beam's damage extent. Thus, the 𝜎 (0,5) can be used as the optimal hyperparameter of this study. The Gaussian window (𝜎 (0,5)) can smooth the first eigenvalue curve without affecting the accuracy of the DSF for eliminating the influence of interference factors.

4.4. Experimental Results and Discussions

Figure 23 shows that the proposed DSF can distinguish the beam's damage extent well on the experimental dataset. Since the heavier vehicle can better excite the response due to the crack, the result in Figure 23 (b) is better than Figure 23 (a) in distinguishing the damage extent. Due to the existence of the elastic spring, the wheel can better maintain the contact with the concrete surface. The result in Figure 23 (b) is smoother and more continuous than that in Figure 23 (a). The size of the gradient reflects the "rhythm" of the beam. The cracks weaken the effective cross-sectional area of the beam, thereby hindering the transmission of information in the beam. In this case, the speed at which the first eigenvalue reaches each local extreme value will slow down and the gradient will become larger. Due to the existence of the elastic spring, part of energy in the response is transferred in the beam using the vehicle as the transmission path. This leads to an increase in the information transmission bandwidth of the beam. When a visible crack zone occurs in the beam, the load in this area will be more borne by the steel bars. Since the vehicle in Figure 23 (b) is 30% heavier than the one in Figure 23 (a), this phenomenon will become deeper as the load increases. Since the load carried by the steel bars increases as the cracks deepen, the overall "rhythm" of the beam becomes faster and the gradient becomes smaller. Thus, these two points will reverse the change pattern of the DSF when the cracks grow but require a further study.

5. Conclusions

A new structural damage detection method based on moving principal component analysis (MPCA) has been developed for condition assessment of the highway bridge under moving vehicles. Numerical and experimental results show that the proposed method is effective and accurate to detect the structural damage. The following points are concluded in this paper,

• The gradient of the first eigenvalue curve obtained from raw acceleration signals by MPCA is used as the damage sensitive feature (DSF) of the highway bridge. The DSF can clearly reflect the existence of the breathing crack on the bridge. The change pattern of the first eigenvalue curve induced by the different vehicle's mass, temperature fluctuations, different damage depths and locations has been studied and the results show the robustness, accuracy and practicality of the proposed DSF.
• The DSF is not limited to several few pre-considered parameters but reflects the beam's damage extent from a dynamic perspective. As the damage in the concrete structures is a crack zone in the actual situation, the equivalent crack depth indicated by this DSF could reflect the damage extent of the beam.
• The experimental results show that the Gaussian window is useful to improve the performance of MPCA on actual datasets. This window can filter out the impact of effects like measurement noise and the vehicle-bridge interaction. The experimental result also shows that the DSF can detect and distinguish the crack damage of different extents under the different vehicle's weight.