A Most-Unfavorable-Condition Method for Bridge-Damage Detection and Analysis Using PSP-InSAR

Abstract

The main contribution of this study is to provide a new idea to detect bridge damage by using PSP-InSAR technology. A most-unfavorable-condition method is proposed for bridge-damage detection and analysis. The method can determine the specific damaged location and occurrence time by using the differential deformation values of persistent scatterer (PS) points on bridge piers. Taking Beijing Suzhou Bridge as an experimental area, 96 COSMO-SkyMed time-series SAR images were used from September 2011 to November 2017. Deformation values of PS points around Suzhou Bridge were acquired and analyzed. Experimental results show that in July 2017, the unusual maximum differential deformation value was 25.73 mm. It occurred between piers D₃ and D₄ of Suzhou Bridge, and it was deduced that the main girder between piers D₃ and D₄ may have been damaged in July 2017. As a validation, taking the differential deformation value between piers D₃ and D₄ as an input, the maximum tensile stress, and the maximum compressive stress were calculated as 2.1 MPa and 8.4 MPa, respectively, through a finite element model. The tensile stress exceeded the design value of the concrete, further confirming the damage of the girder between piers D₃ and D₄. Moreover, all results are consistent with the Suzhou Bridge damage information shown in existing records, which verify the accuracy and reliability of the proposed method.

1. Introduction

Bridges are an important part of highway and railway networks. They play a fundamental role in ensuring and promoting social and economic development. However, factors such as traffic congestion, vehicle overload, frequent wind load, environmental degradation and climate disasters inevitably exacerbate the degradation of bridges. With an increase of bridge deformation, if not properly maintained, a bridge may collapse suddenly [1,2]. Thus, in order to carry out effective maintenance and ensure the operational safety of bridges, it is very important to acquire the damaged-location and time information accurately, and timely, based on deformation monitoring information [3], especially for potentially damaged bridges.

At present, the commonly used methods for bridge monitoring generally include on-site measurement and remote sensing measurement [4]. On-site measurement methods such as infrared thermography or spectroscopy [5], laser scanning [6], ground-penetrating radar [7], GPS receiver [8], etc., have a high accuracy presupposing that enough instruments or sensors are currently equipped on the bridges. Otherwise, the spatial sampling rate of point-by-point measurements would be low, which further affects the accuracy of bridge monitoring. Remote sensing measurement methods include ground-based radar [6,9], airborne/laser radar [10], and Interferometric Synthetic Aperture Radar (InSAR) [11], etc. Ground-based radar is always expensive and difficult to deploy in a fixed location for the collection of continuous observations. Airborne/laser radar is greatly affected by weather conditions, which are difficult to calculate and eliminate. InSAR technology has been widely used in bridge micro-deformation monitoring due to its advantageous all-weather, high precision and high resolution (mm) features [2,12]. However, InSAR monitoring accuracy is related to phase and is easily restricted by atmospheric delay. The Persistent Scatterers Interferometric Synthetic Aperture Radar (PS-InSAR) technique can overcome atmospheric delay and improve deformation monitoring accuracy [2]. The method selects the relative stable points as Persistent Scatter (PS) points. It uses SAR images with multiple views (generally more than 25 views) to statistically analyze the amplitude information, and to remove the phase values of the atmospheric effect. It has been found that the density of high-quality PS points has a great influence on monitoring accuracy. In order to overcome the PS-point-density limitation of PS-InSAR technology, Costantini M et al. proposed the Persistent Scatterer Pairs InSAR (PSP-InSAR) method and utilized a high-resolution Cosmo-SkyMed image to obtain the PS point deformation rate of a bridge across the Sochi River [13,14]. The principle of the PSP-InSAR technique is to establish a PS grid by connecting pixel pairs with the same phase characteristics, to expand the grid, and finally to select the PS point-sets. The PSP-InSAR technique can not only overcome the influence of atmospheric delay on the InSAR technique, but also improve the accuracy of PS point deformation measurements. As there have been few scientific studies to analyze bridge situations using PSP-InSAR, we aim to provide a new idea and method in this study to obtain the PS point-sets and to detect bridge damage based on PSP-InSAR technology.

Studies have shown that the uneven deformation phenomenon between bridge piers is directly related to the damage of bridges [15,16]. Therefore, it is necessary to monitor and analyze the deformation of bridge piers. In practical problems, the distribution of PS points obtained by PSP-InSAR technology is uneven. Furthermore, there are not enough PS points on bridge piers to provide data for research. In order to solve this problem, we plan to use the curved-surface interpolation fitting method to obtain the approximate position and deformation information of bridge piers based on the existing PS points that are available. The commonly used surface-interpolation techniques mainly include the polynomial surface-interpolation method [17], the thin-plate spline-interpolation method [18], and Green’s function-based interpolation method [19,20], etc. The polynomial surface fitting method is a surface fitting method based on the principle of least squares. The thin-plate spline interpolation fitting method is based on the minimum curvature method to ensure the smoothness of the interpolated surface. These two methods both have high requirements for data quality, for example, the control points should be evenly distributed in the grid. Additionally, they both require large amounts of analysis and calculation [21]. Green’s function-based interpolation method can avoid the strict demands on the number and distribution of control points. It can not only achieve consistency between the values obtained for the fitting point and the real value at the control point, but also achieve a good level of accuracy for the overall fitting results. In this study, Green’s function-based interpolation method is used to obtain the three-dimensional (3D) deformation field and time-series deformation values of bridge piers for the first time in the literature.

As stated above, in order to acquire the damaged location and occurrence time accurately for carrying out effective maintenance, and ensuring the operational safety of bridges, the most-unfavorable-condition method for bridge-damage detection and analysis was proposed, and three critical items were investigated in this study: (i) The 3D time-series deformation model based on the existing bridge-surrounding, and encrypted interpolated PS points, which were established using Green’s function-based interpolation method. The deficiency of PS points on bridge piers can be solved and the deformation information can be obtained; (ii) The most-unfavorable-condition method containing the two most unfavorable working condition modes (I and II), which was proposed based on the deformation values of PS points on bridge piers from constructed 3D time-series deformation models. The damaged location and occurrence time for the bridge can then be obtained; (iii) Taking Beijing Suzhou Bridge as an experimental area, 96 issued COSMO-SkyMed time-series images from September 2011 to November 2017 were used. In addition, the finite element method, and the existing record were taken as references to further verify the accuracy and reliability of the results. The remainder of the paper is organized as follows. In Section 2, the most-unfavorable-condition method for bridge-damage detection and analysis is explained, and the experimental research introduced, including a description of the research area and materials used, as well as the experimental design. In Section 3, the applied experiments and results for the proposed method on the Suzhou Bridge are described. Section 4 presents the conclusions.

2. Materials and Methods

2.1. Study Area and Materials Description

The Beijing Suzhou Bridge, built in 1994, is composed of a viaduct on the north-west of the Third Ring road and three ramp bridges. Existing records have shown that one of the ramps was damaged, and the concrete cracked, around the end of 2017. This damaged ramp is referred to as D-ramp later in the study. The D-ramp is a continuous box-girder section and consists of 5 bridge piers. The position of the D-ramp and its 5 piers are shown in Figure 1.

In this study, 96 SAR images acquired from 3 m high-resolution Cosmo-SkyMed images in strip mode, which were developed in Italy, are used. The data covered the area around Suzhou Bridge from 30 September 2011, to 6 November 2017. The basic parameters are shown in

Table 1. Basic parameters of InSAR data used in the study area.

Item	Value
Satellite	COSMO-SkyMed
Imaging mode	Stripmap
Ground resolution	3 m × 3 m
Incident angle	~20.07°
Polarization mode	HH
Amount of data	96
Time span	Sep. 2011~Nov. 2017

.

PSP-InSAR technology was used to process a total of 96 images from September 2011 to November 2017. Based on the deformation rate of PS points obtained by geocoding, the annual average deformation rate in the line-of-sight direction was calculated and drawn by using 96 deformation values, as shown in Figure 2. It was found that the deformation rate is −7.5~3.1 mm/y in the line-of-sight direction of the study area during 2011~2017. Moreover, the deformation rate in the southwest direction of Suzhou Bridge was relatively high. It was also found that the quantity of PS points in the D-ramp bridge area was low. Furthermore, they were unevenly distributed, especially in the vicinity of bridge piers D₄ and D₅. This makes it impossible to directly obtain the deformation data of pier positions.

2.2. Method

Figure 3 shows the entire technical framework of the most-unfavorable-condition method for bridge-damage detection and analysis, which consists of the following two steps: (i) The temporal 3D deformation model is built by using Green’s function-based interpolation method, in order to solve the problem of the PS points’ uneven distribution and deficiency on the bridge piers. Then the deformation values of PS points on the bridge piers can be acquired as inputs for the most-unfavorable-condition method; (ii) Regarding the deformation situation of each PS point on the piers from each SAR image as one working condition, the most-unfavorable-condition method, containing the two most unfavorable working condition modes (I and II), was proposed, so as to determine the specific damage location and occurrence time for the bridge.

2.2.1. Three-Dimensional Deformation Model Construction

Green’s function-based interpolation method is used primarily to solve the problem of the PS points’ uneven distribution and deficiency on bridge piers. Green’s function-based interpolation method has the advantages of a high calculation accuracy, a stable display of local abnormal changes, and a suppression of false changes in areas lacking data [20]. The surface generated by interpolation in this study is a linear combination of Green functions centered on the corresponding slope values of each point. As long as the interpolated surface satisfies Green’s function, the interpolated surface has the characteristic of minimum curvature. Green’s function-based interpolation algorithm can be described as follows [20].

Suppose that there are n points $q_i\left(i=1,2,\cdots n\right)$ scattered in m-dimensional space, the m-dimensional interpolation problem of Green’s spline can then be simplified to the solution of Equation (1).

$$\{\begin{array}{c}{\nabla }^{4}w\left(Q\right)={\displaystyle \sum _{i=1}^n{\alpha }_i\delta \left(Q-Q_i\right)}\hfill \\ w\left(Q_i\right)=w_i\hfill \end{array}$$

(1)

where ${\nabla }^4$ is Green’s spline interpolation operator, $\delta \left(\right)$ is a unit impulse function, $Q_i$ represents the position information of point $q_i$ in the m-dimensional space, $w\left(Q_i\right)$ is the corresponding attribute value (a deformation value in this study) of point $q_i$.

The general solution satisfying Equation (1) is:

$$w\left(Q\right)={\displaystyle \sum _{i=1}^n{\alpha }_i{\phi }_m\left(Q-Q_i\right)}$$

(2)

The coefficient $a_i$ satisfies the linear Equation (3):

$$w_i={\displaystyle \sum _{j=1}^n{\alpha }_j{\phi }_m\left(Q_i-Q_j\right)}$$

(3)

where ${\phi }_m\left(Q\right)$ is the m-dimensional Green function. The m-dimensional Green function expressions are shown in

Table 2. m-dimensional Green function.

Dimension	$\mathbf{Green} \mathbf{Function} {\mathit{\phi }}_{\mathit{m}}\left(\mathit{Q}\right)$
1	${\left|Q\right|}^3$
2	${\left|Q\right|}^2\left(\ln \left|Q\right|-1\right)$
3	$\left|Q\right|$
4	$\ln \left|Q\right|$
5	${\left|Q\right|}^{-1}$
6	${\left|Q\right|}^{-2}$
m	${\left|Q\right|}^{4-m}$

. In order to obtain the 3D deformation model, the two-dimensional (m = 2) Green function is adopted in this study [20].

Green’s function-based interpolation method can interpolate data points that are randomly distributed in an m-dimensional plane. The known data points are called definition points. In this paper, Green’s function-based interpolation method is applied to produce a large number of deformation values for PS points on bridge piers by using the known deformation values of all known PS points shown in Figure 2. The specific Green’s function-based interpolation algorithm that is applied in this study to obtain the deformation values on bridge piers is described as follows:

(1) Coordinates of definition points $q_i\left(i=1,2,\cdots n\right)$ and piers’ points (needed to be calculated) are expressed as $Q_i\left(x_i,y_i\right)$ and $Q_u\left(x_u,y_u\right)$, respectively, $z_i$ and $z_u$ are the deformation values of corresponding points, $z_u$ was unknown and needed to be acquired through Green’s function-based interpolation.

(2) In order to obtain the coefficients $a_i$, a system of linear equations can be constituted by using Equation (4):

$$\{\begin{array}{c}{\alpha }_1{\phi }_2\left(Q_1-Q_1\right)+{\alpha }_2{\phi }_2\left(Q_1-Q_2\right)+{\alpha }_3{\phi }_2\left(Q_1-Q_3\right)+\cdots {\alpha }_n{\phi }_2\left(Q_1-Q_n\right)=z_1\hfill \\ {\alpha }_1{\phi }_2\left(Q_2-Q_1\right)+{\alpha }_2{\phi }_2\left(Q_2-Q_2\right)+{\alpha }_3{\phi }_2\left(Q_2-Q_3\right)+\cdots {\alpha }_n{\phi }_2\left(Q_2-Q_n\right)=z_2\hfill \\ {\alpha }_1{\phi }_2\left(Q_3-Q_1\right)+{\alpha }_2{\phi }_2\left(Q_3-Q_2\right)+{\alpha }_3{\phi }_2\left(Q_3-Q_3\right)+\cdots {\alpha }_n{\phi }_2\left(Q_3-Q_n\right)=z_3\hfill \\ \cdots \hfill \\ {\alpha }_1{\phi }_2\left(Q_n-Q_1\right)+{\alpha }_2{\phi }_2\left(Q_n-Q_2\right)+{\alpha }_3{\phi }_2\left(Q_n-Q_3\right)+\cdots {\alpha }_n{\phi }_2\left(Q_n-Q_n\right)=z_n\hfill \end{array}$$

(4)

where $Q_i-Q_j\left(i=1,2\cdots ,n;j=1,2\cdots ,n\right)$ is the horizontal distance between two PS points, that is, $Q_i-Q_j=\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}=S_{ij}$, and thus $\phi \left(Q_i-Q_j\right)=\phi \left(S_{ij}\right)=\left|S_{ij}\right|\left(\ln \left|S_{ij}\right|-1\right)$. It should be noted that when $i=j$, $\phi \left(Q_i-Q_j\right)=0$.

(3) The coefficients ${\alpha }_i\left(i=1,2,\cdots n\right)$ of linear equations are obtained by using the Gauss elimination method [21].

(4) The deformation values of unknown PS points on bridge piers can be acquired using Equation (5):

$$z_u={\alpha }_1{\phi }_2\left(Q_u-Q_1\right)+{\alpha }_2{\phi }_2\left(Q_u-Q_2\right)+{\alpha }_3{\phi }_2\left(Q_u-Q_3\right)+\cdots {\alpha }_n{\phi }_2\left(Q_u-Q_n\right)$$

(5)

Finally, according to the deformation values of original and interpolated PS points, a 3D deformation surface-field, namely the 3D deformation model, can be acquired. In the 3D deformation model, the probable site of damage on the bridge is the point with the maximum deformation value, that is, the center of the subsidence funnel. Therefore, it can be deduced that the information regarding the damage of the bridge can be obtained by further analyzing the subsidence-funnel-center of the built model.

2.2.2. The Most-Unfavorable-Condition Method

The deformation state of bridge piers can be treated as the working condition. In this study, each working condition for PS points on bridge piers can be acquired from each SAR image. Safe working conditions refer to the bridge in normal operation, whereas unfavorable working conditions mean unfavorable differential deformation has occured. Once the bridges are prone to deterioration, or are even damaged, it is called the most-unfavorable working condition. By analyzing the most unfavorable working conditions of the bridge, the state of the damaged bridge can be understood.

When the deformation values or the deformation rate of PS points on bridge piers increase too much, it is most likely the stage for bridge damage to occur. It is very important to acquire the damaged-location and occurrence time information of bridges in a timely manner, to carry out effective maintenance and maintain normal traffic and transportation. In order to obtain the damaged location and time, local analysis of each bridge pier under each working condition is needed, which is time-consuming and labor-intensive. Therefore, the most-unfavorable-condition method including the two most unfavorable working condition modes is established. First, in the most-unfavorable working condition I, the damaged location and occurrence time interval can be determined by using the deformation values and deformation rate of PS points on bridge piers. Then in the most-unfavorable working condition II, the specific location of damage and its corresponding occurrence time can be acquired by analyzing the maximum differential deformation values between adjacent piers. The two most-unfavorable working-condition modes are described as follows:

1.

The most-unfavorable working-condition I:

Considering that 96 SAR Cosmo-SkyMed images were used in this study, the deformation situation of each PS point on each image is defined as a working-condition, and each PS point has a time-series deformation value with a length of 96. The deformation values of PS points on bridge piers can be obtained through the 3D deformation model. The time-series working conditions based on the deformation values of PS points on bridge piers can be expressed as follows:

$$\left[\begin{array}{c}{D}_1:P_{1,1},P_{1,2},P_{1,3},\cdots ,P_{1,j},\cdots ,P_{1,M}\hfill \\ D_2:P_{2,1},P_{2,2},P_{2,3},\cdots ,P_{2,j},\cdots ,P_{2,M}\hfill \\ ⋮\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdots \hfill \\ D_N:P_{N,1},P_{N,2},P_{N,3},\cdots ,P_{N,j},\cdots ,P_{N,M}\hfill \end{array}\right]$$

(6)

where $N$ represents the number of PS points on bridge piers, $M\left(M=96\right)$ represents the number of working conditions of each PS point.

(1) The maximum absolute deformation values from the data of 96 $N$ PS points are firstly acquired and expressed as:

$$\left\{\begin{array}{cccc}{P}_{1,j_1}& P_{2,j_2}& \begin{array}{ccc}\cdots & P_{i,j_i}& \cdots \end{array}& P_{N,j_N}\end{array}\right\}with\left\{\begin{array}{cccc}{j}_1& j_2& \begin{array}{cc}\cdots & \begin{array}{cc}{j}_i& \cdots \end{array}\end{array}& j_N\end{array}\right\}\in \left[1,\mathrm{M}\right]$$

(7)

(2) $j_{\max }=\max \left(\begin{array}{cccc}{j}_1& j_2& \begin{array}{cc}\cdots & \begin{array}{cc}{j}_i& \cdots \end{array}\end{array}& j_N\end{array}\right)$. Taking $P_{i,j_{\max }}\left(i=1,2,\cdots N\right)$ as the cut-off conditions of the right side, a new time-series working condition of the PS points can be obtained as:

$$\left[\begin{array}{c}{D}_1:P_{1,1},P_{1,2},P_{1,3},\cdots ,P_{1,j_{\max }}\hfill \\ D_2:P_{2,1},P_{2,2},P_{2,3},\cdots ,P_{2,j_{\max }}\hfill \\ ⋮\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdots \hfill \\ D_N:P_{N,1},P_{N,2},P_{N,3},\cdots ,P_{N,j_{\max }}\hfill \end{array}\right]$$

(8)

In order to make it easier to understand, Figure 4 shows the process that is taken to decide $j_{\max }$. Taking time-series $\left(M=10\right)$ deformation values for three PS points $\left(N=3\right)$ as an example (only for illustration, not real data), it is easy to obtain the maximum absolute deformation values of 3 points, that is, $j_1=j_2=7,j_3=5$. Thus, $j_{\max }=\max \left(j_1,j_2,j_3\right)=7$. The rebuilt time-series working conditions will be:

$$\left[\begin{array}{c}{D}_1:P_{1,1},P_{1,2},P_{1,3},\cdots ,P_{1,7}\hfill \\ D_2:P_{2,1},P_{2,2},P_{2,3},\cdots ,P_{2,7}\hfill \\ D_3:P_{3,1},P_{3,2},P_{3,3},\cdots ,P_{3,7}\hfill \end{array}\right]$$

(9)

(3) Based on the maximum absolute deformation values from the data concerning 96 $N$ PS points, that is, $\begin{array}{cccc}{P}_{1,j_1}& P_{2,j_2}& \begin{array}{ccc}\cdots & P_{i,j_i}& \cdots \end{array}& P_{N,j_N}\end{array}$, it is possible to calculate the deformation rate values from other instances of each PS point:

$$\{\begin{array}{l}{R}_{1,k}=\left|\frac{P_{1,j_1}-P_{1,k}}{j_1-k}\right|\\ R_{2,k}=\left|\frac{P_{2,j_2}-P_{2,k}}{j_2-k}\right|\\ ⋮\\ R_{N,k}=\left|\frac{P_{N,j_N}-P_{N,k}}{j_N-k}\right|\end{array}\left(k=1,2,\cdots ,j_{\max }\right)$$

(10)

(4) Find the maximum deformation-rate value of each PS point:

$$\{\begin{array}{l}\max \left(R_{1,k}\right)=R_{1,k_1}\\ \max \left(R_{2,k}\right)=R_{2,k_2}\\ ⋮\\ \max \left(R_{N,k}\right)=R_{N,k_N}\end{array}$$

(11)

(5) Define $k_{\min }=\min \left(\begin{array}{cccc}{k}_1& k_2& \cdots & k_N\end{array}\right)$. Taking $k_{\min }$ as the starting working conditions (left side), a new time-series working condition for PS points can be obtained as:

$$\left[\begin{array}{c}{D}_1:P_{1,k_{\min }},P_{1,k_{\min }+1},\cdots P_{1,h},\cdots ,P_{1,j_{\max }}\hfill \\ D_2:P_{2,k_{\min }},P_{2,k_{\min }+1},\cdots P_{2,h},\cdots ,P_{2,j_{\max }}\hfill \\ ⋮\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdots \hfill \\ D_N:P_{N,k_{\min }},P_{N,k_{\min }+1},\cdots P_{N,h},\cdots ,P_{N,j_{\max }}\hfill \end{array}\right],h=k_{\min },k_{\min }+1,\cdots j_{\max }$$

(12)

Using the above example as an illustration, Figure 5 is a simple schematic diagram, in which the black lines represent the maximum deformation rate of each point, and the black arrows represent the deformation value involved in calculating the rate. It is simple to obtain the maximum deformation rate of three PS points $k_1=4,k_2=6,k_3=3$, therefore $k_{\min }=\min \left(k_1,k_2,k_3\right)=3$.

2.

The most-unfavorable working condition II:

The simultaneous uneven deformation of PS points on adjacent piers may cause damage to the bridge. The uneven deformation of adjacent piers to obtain the specific damage location and occurrence time, based on the damage time interval, and determined from the most-unfavorable working condition I should be further analyzed.

First, taking the time-series working conditions of the first pier (D₁) as a reference, the relative deformation values of PS points on other piers can be calculated. A new time-series working-conditions expression is:

$$\left[\begin{array}{c}{D}_1:0,0,\begin{array}{ccc}\cdots & 0& \cdots \end{array},0\hfill \\ D_2:P_{2,k_{\min }}-P_{1,k_{\min }},P_{2,k_{\min }+1}-P_{1,k_{\min }+1},\begin{array}{ccc}\cdots & P_{2,w}-P_{1,w}& \cdots \end{array},P_{2,j_{\max }}-P_{1,j_{\max }}\hfill \\ ⋮\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdots \hfill \\ D_i:P_{i,k_{\min }}-P_{1,k_{\min }},P_{i,k_{\min }+1}-P_{1,k_{\min }+1},\begin{array}{ccc}\cdots & P_{i,w}-P_{1,w}& \cdots \end{array},P_{i,j_{\max }}-P_{1,j_{\max }}\hfill \\ ⋮\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdots \hfill \\ D_N:P_{N,k_{\min }}-P_{1,k_{\min }},P_{N,k_{\min }+1}-P_{1,k_{\min }+1},\begin{array}{ccc}\cdots & P_{N,w}-P_{1,w}& \cdots \end{array},P_{N,j_{\max }}-P_{1,j_{\max }}\hfill \end{array}\right]$$

(13)

Then, to calculate the differential deformation values between the adjacent piers from D₁ to D_N, pick out the maximum differential deformation values as:

$$diff\_D=\left({\left|P_{2,r}-P_{1,r}\right|}_{\max },{\left|P_{3,r}-P_{2,r}\right|}_{\max },\cdots {\left|P_{N,r}-P_{N-1,r}\right|}_{\max }\right),r=k_{\min },k_{\min }+1,\cdots j_{\max }$$

(14)

The PS point with the maximum differential-deformation value in Equation (14), that is, $D_{\max }=\max (diff\_D)$, is the probable site that may lead to damage. Moreover, the occurrence time also can be acquired from the corresponding working-condition information (that is the issue-time of the working condition).

2.3. Test Design

In order to analyze the damaged location and occurrence time of a damaged bridge, Suzhou Bridge in Beijing was taken as the experimental area. 96 issues of COSMO-SkyMed image data of this area were collected to verify the effectiveness of the proposed method. PSP-InSAR technology was used to process the time series data to obtain the overall land deformation information of the PS points. The PS point data for Suzhou Bridge and its surroundings were used for spatial interpolation to deal with the lack of PS points on Suzhou Bridge’s piers. In this study, Green’s function-based interpolation method was adopted to obtain PS point deformation values where deficient, which is the first of its kind. Then, a 3D time-series deformation model was built by interpolating the existing PS points in the study area. Regarding the deformation situation of each PS point on piers in each SAR image as a working condition, the time-series working conditions can be obtained. The proposed most-unfavorable-condition method is used for pier-local analysis based on the two most-unfavorable working-condition modes. Through the above process, the damaged position and occurrence time can be obtained. Results were then compared with information from the finite element method and the existing record to further verify the accuracy and reliability of the proposed method.

3. Results

3.1. 3D Deformation Model Construction

3.1.1. Model Construction Method Comparison and Analysis

In order to reflect the high precision, high reliability, and universal applicability of Green’s function-based interpolation method, different methods including the quadratic polynomial surface interpolation method, the thin-plate spline interpolation method, and Green’s function-based interpolation method were applied to fit the same set of PS points. Figure 6, Figure 7 and Figure 8 are 3D deformation model construction results obtained by using the three different methods, respectively. Figure 9 is a residual error diagram generated by Green’s function-based interpolation method.

Table 3. Residual error statistics of the three interpolation methods.

Fitting Method	Max_Err/mm	SSE/mm2	R2
Quadratic polynomial surface interpolation	11.296	4.237	0.0177
Thin plate splines interpolation	6.3394	1.655	0.6163
Green’s function-based interpolation	0.33	0.185	0.9182

shows the residual error statistics acquired from the three interpolation methods. Max_Err is the maximum absolute error value, SSE is the sum of squares due to error, and R² is the coefficient of determination. The lower the Max_Err and SSE values, the better the interpolation. The higher the R² value, the better the interpolation [22].

$$\{\begin{array}{c}Max\_Err=\max \left|d_i-{\stackrel{\wedge }{d}}_i\right|\hfill \\ \hfill \\ SSE={\displaystyle \sum _{i=1}^n{\left(d_i-\stackrel{\wedge }{d_i}\right)}^2}\hfill \\ R^2=1-\frac{{{\displaystyle \sum _{i=1}^n\left(d_i-\stackrel{\wedge }{d_i}\right)}}^2}{{\displaystyle \sum _{i-1}^n{\left(d_i-\overline{d}\right)}^2}}\hfill \end{array}$$

(15)

where $d_i$ is the deformation value of the defined PS point, ${\stackrel{\wedge }{d}}_i$ is the deformation value obtained from the constructed model, and $\overline{d}$ is the mean value of deformation.

It can be concluded from Figure 6, Figure 7, Figure 8 and Figure 9 and Table 3 that the fitting result obtained by the interpolation method based on Green’s function is closest to the true value. As shown in Table 3, the maximum absolute error of Green’s function-based interpolation fitting is 0.33 mm and its SSE is 0.185 mm², which are smaller than the fitting results from quadratic polynomial surface fitting and thin-plate spline interpolation. The fitting coefficient of Green’s function-based interpolation is 0.98182, which is much closer to 1 than the other two methods. It proves that the fitting effect of the interpolation method based on Green’s function is the best.

3.1.2. Three-Dimensional Deformation Model Construction

Taking the first-issue data from September 2011 as a reference, Green’s function-based interpolation method was performed for PS points within the demonstrated perimeter in Figure 2 (X∈[116.299°, 116.305°]; Y∈[39.958°, 39.964°]). The deformation values of PS points in the first issue are taken as zeros and the corresponding plane model is shown in Figure 10. For the sake of description and display, a total of 6 issued PS points’ deformation values from the same area were randomly selected from 2012 to 2017. Then, the Green’s function-based interpolation method was used to generate 3D deformation surface models, as shown in Figure 11. X- and Y-axes represent latitude and longitude respectively, and the Z-axis represents cumulative ground deformation values.

As can be seen from the 3D deformation models shown in Figure 11, ground subsidence is indicated the when the deformation value is negative, and ground uplift is indicated where the deformation value is positive. Since 2013, the subsidence of the ground where the Suzhou Bridge is located has become increasingly serious, and there is an obvious subsidence funnel in this area. In 2017, the maximum cumulative subsidence value of the subsidence funnel was close to 70 mm. Excessive ground subsidence will have a huge impact on the bridge. The coordinate values of the maximum subsidence funnel center can be obtained from the 3D deformation model, which is consistent with the position of the fourth pier (D₄). Then, the time-series deformation values of PS points on the other four piers can be obtained from the time-series 3D deformation model. The time-series deformation surfaces of the five piers are shown in Figure 12.

It can be seen from Figure 12 that the PS points of five piers have been declining slowly from 2013 to the beginning of 2017. The maximum deformation values all occurred on 2 August 2017. During the period from 1 July to 2 August 2017, the bridge deformation rate was the largest, as indicated by red lines in Figure 12. By the end of 2017, the maximum annual average deformation rate of PS points on pier D₄ was −14.5 mm/y.

3.2. The Most-Unfavorable Working-Condition Analysis

The 96 issued time-series deformation values of five piers were sorted into sequence working conditions, as shown in

Table 4. Time series working conditions of PS points on bridge piers.

No.	D1	D2	D3	D4	D5	No.	D1	D2	D3	D4	D5
1	0	0	0	0	0	49	−21.2	−14.3	−13	−26	−9.3
2	−1.3	0.1	2.68	1.28	1.3	50	−18.0	−14.5	−12.9	−24.5	−9
3	−2.1	−1.3	1.88	1.88	0.8	51	−23.0	−14.1	−13	−25	−10.2
4	−1	−1.5	2.18	1.48	1	52	−22.8	−16.7	−13.5	−26	−10.9
5	−3.24	−1.9	−1.65	0.98	0.7	53	−23.9	−16.3	−15.2	−27.1	−8.8
6	−2.98	−2.3	2.62	2.58	−1.1	54	−23.0	−16.1	−14.3	−28.1	−8.8
7	−4.56	−1.9	1.12	2.78	0.7	55	−22.7	−14.7	−19.2	−29.5	0.6
8	−4.8	−4.6	1.42	2.78	1.2	56	−19.0	−14.7	−16.8	−27.5	−10.7
9	−5.0	−4.6	2.32	1.98	1.5	57	−27.5	−16.5	−21.5	−28.3	−9.7
10	−4.4	−2.4	1.52	1.68	1.2	58	−24.8	−15.2	−18.8	−28.4	−10.5
11	−3.9	−2.2	4.22	3.37	3.1	59	−25.0	−16.2	−19.2	−32.8	−13.3
12	−7.2	−2.5	1.32	−1.63	−0.9	60	−27.4	−19.1	−22.9	−35.8	−13.1
13	−7.2	−2.9	−2.98	−4.43	−8.5	61	−26.1	−18.6	−23.6	−33	−13.2
14	−0.8	−2	2.67	1.67	6.9	62	−20.0	−21	−27.5	−37.6	−13.5
15	5.7	−3.2	9.61	11.07	9.7	63	−28.5	−22.3	−26.9	−37.2	−13.4
16	1.1	2.5	2	2.27	7.2	64	−28.8	−21.1	−25.8	−35.8	−12.2
17	−3.7	3.6	1.6	3.37	7.4	65	−29.7	−20.1	−29.5	−36.7	−14.9
18	−7.3	−5.1	−2.1	−5.83	4.9	66	−30.9	−23	−27.5	−38.5	−13.8
19	1.5	0.9	2.9	−4.43	6	67	−28.2	−20.9	−27.7	−37	−14.7
20	−7.1	−4.6	0.2	−10.1	4.8	68	−29.0	−23.3	−27.4	−39	−14.8
21	−5.9	0.6	2.1	−14.1	3.1	69	−29.9	−21.1	−28.2	−39.5	−14.8
22	−6.0	2.6	3.1	−16.5	4.6	70	−28.7	−19.6	−26.3	−39.2	−16
23	−8.2	3.1	2.2	−16.1	3.6	71	−30.6	−19.7	−28.5	−42	−16.3
24	−6.3	0	2.2	−16.3	3.7	72	−32.1	−21.3	−29	−46	−17.9
24	−5.4	3.8	3.3	−13.7	4.7	73	−32.5	−19.2	−28.2	−46.3	−19.5
26	−8.3	4.1	2.5	−17	4.6	74	−26.0	−20.2	−27.5	−46.6	−22.4
27	−2.3	3.7	−7.2	−15.6	3.8	75	−32.0	−19.4	−29.6	−49.5	−24.1
28	−13.2	−6.7	−7.8	−17.5	7.4	76	−28.0	−17.8	−27.9	−50.69	−23.7
29	−14.0	−6.2	−8.4	−20.9	−1.9	77	−26.0	−21.8	−32.2	−55.1	−24.2
30	−12.1	−6.2	−6.5	−18.4	4	78	−28.0	−21.3	−35.5	−53.7	−24.6
31	−10.8	−4.7	−7.6	−17.6	4.1	79	−25.0	−22.8	−36.1	−48.7	−24.9
32	−11.8	−9.3	−9.9	−20.6	2.7	80	−26.0	−33.6	−44.1	−57.1	−35
33	−13.0	−8.2	−8.6	−19.5	3.2	81	−25.9	−33	−43.8	−56.7	−36.6
34	−13.5	−10.3	−10.4	−19.9	2.3	82	−26.1	−33.9	−42.8	−60.1	−37
35	−13.3	−11.9	−8.7	−22	3.3	83	−35.7	−34.3	−41.9	−60.4	−37.7
36	−13.5	−11.5	−8.2	−20.5	2.8	84	−35.2	−35.1	−42.7	−61.7	−37.8
37	−10.9	−11.9	−7.6	−16.7	2.7	85	−36.5	−33.9	−44.1	−62.9	−38.9
38	−12.7	−9.5	−7.5	−19.8	−5.4	86	−36.8	−33.8	−43.6	−63.8	−39.3
39	−18.0	−15.1	−12	−22.3	−6.8	87	−38.0	−34.5	−44.8	−68.7	−45.5
40	−13.9	−9.6	−9.5	−23.8	−5.5	88	−43.0	−39.2	−47.2	−71	−47.1
41	−22.6	−17.5	−14.7	−22.9	−6.3	89	−40.8	−39.4	−47	−70.9	−46.2
42	−21.1	−13.7	−10.9	−21.4	−5.9	90	−45.7	−43.5	−49.8	−75.53	−50.3
43	−20.7	−15.1	−12.5	−23.4	−6.4	91	−48.9	−45.7	−51.5	−76.2	−51.2
44	−21.1	−13.9	−9.5	−22.9	−5.3	92	−44.6	−40.6	−48.5	−70.2	−47.6
45	−18.7	−9.5	−8.3	−23.5	−3.4	93	−41.4	−36.5	−44.3	−67.9	−44.6
46	−19.7	−8.6	−4.8	−20.1	−3.1	94	−38.1	−32.8	−42.2	−65.3	−42
47	−20.4	−12.4	−11.3	−23.3	−7.8	95	−35.7	−30.7	−41.6	−66.4	−41.9
48	−21.0	−14.9	−13.1	−23.2	−9	96	−35.3	−31.6	−40.8	−66.4	−40.2

. The deformation value is shown in mm. Working condition 1 corresponds to the data of PS points on piers from the first issue. Since it is the first time that data has been obtained in this area, the deformation values of condition 1 were all set to zero.

When the absolute differential deformation values and corresponding deformation rates of adjacent piers are at their maximum, Suzhou Bridge is most likely to be damaged. Therefore, the two most unfavorable working-condition modes were built to acquire the deformation values and to detect the bridge damage.

• Condition I:

According to the most-unfavorable condition I, the maximum absolute deformation values obtained from the five piers are all from working condition 91, as shown in Table 4. The corresponding date was 2 August, 2017. The maximum absolute deformation values from piers D₁ to D₅ are −48.9 mm, −45.7 mm, −51.5 mm, −76.2 mm, and −51.2 mm, respectively. Then, according to Equations (10) and (11), the maximum deformation rate values based on condition 91 are from condition 89. The maximum deformation rate values of all five piers are 4.05 mm/issue, 3.15 mm/issue, 2.25 mm/issue, 2.65 mm/issue and 2.5 mm/issue, respectively. Thus, under Equation (12), a new time-series working condition of PS points on piers can be built from conditions 89 to 91, as shown in Equation (16). The corresponding period is from 1 July to 2 August and it indicates that the Suzhou Bridge may have been damaged during this period.

2.

Condition II:

According to the most-unfavorable condition II, the deformation values of other piers relative to D₃ are also listed in Equation (16). Then, the differential deformation of adjacent piers was calculated. According to Equation (14), in working condition 90, the maximum value of D₄–D₃ was 25.73 mm, as shown below in Equation (17). The corresponding time of working condition 90 is 17 July 2017.

$$\left[\begin{array}{c}{\mathrm{D}}_1:-40.8,-45.7,-48.9\hfill \\ {\mathrm{D}}_2:-39.4,-43.5,-45.7\hfill \\ {\mathrm{D}}_3:-47,-49.8,-51.5\hfill \\ {\mathrm{D}}_4:-70.9,-75.53,-76.2\hfill \\ {\mathrm{D}}_5:-46.2,-50.3,-51.2\hfill \end{array}\right]\to \left[\begin{array}{c}{\mathrm{D}}_1^{’}:0,0,0\hfill \\ {\mathrm{D}}_2^{’}:1.4,2.2,3.2\hfill \\ {\mathrm{D}}_3^{’}:-6.2,-4.1,-2.6\hfill \\ {\mathrm{D}}_4^{’}:-30.1,-29.83,-27.3\hfill \\ {\mathrm{D}}_5^{’}:-5.4,-4.6,-2.3\hfill \end{array}\right]$$

(16)

$$\begin{array}{c}diff\_D=\left(3.2,7.6,25.73,25.23\right)\hfill \\ D_{\max }=25.73mm\hfill \end{array}$$

(17)

The two most-unfavorable working-condition results were plotted as shown in Figure 13. Working condition I is the deformation value corresponding to the maximum deformation rate of each bridge pier. Working condition II is the deformation value of other pier columns relative to pier D₁. It can be clearly seen from the figure that the maximum relative deformation occurs between piers D₃ and D₄, as shown in Figure 14. The results indicate that the damage occured between piers D₃ and D₄ in July 2017. The damage results from this proposed method are consistent with the recorded damage information. Besides this, it can also be seen that the cause of the damage to Suzhou Bridge is the unusual differential deformation between piers D₃ and D₄.

3.3. Damage Detection Using Finite Element Model

In this section, finite element software was used to explore the structural deformation and analyze the normal stress of the Suzhou bridge. The large-scale general structural calculation software MIDAS/Civil 2017 was used for simulation analysis. It was mainly applied to the linear-elastic small deformation beam element calculations [23,24]. The element type of the finite element model is the beam element. The internal forces of the bridge, such as shear force, bending force, and normal stress were calculated by inputting the deformation values of the bridge. The position where the maximum normal stress of the bridge appears indicates that the bridge structure is subjected to the maximum tension. When the tension exceeds the design value of the bridge, the bridge will be damaged. According to the structural design of Suzhou Bridge, the calculation of live load is based on the steam super–20 and hanging–120 automobile load classes, and the uneven deformation of the support is comprehensively considered based on the uneven deformation values shown in Equation (17). The results acquired from the finite element model are shown in Figure 15. The calculation model is from piers D₁ to D₅. The positions of piers are shown by black arrows in Figure 15. The maximum tensile stress is 2.1 MPa, and the maximum compressive stress is 8.4 MPa. In the finite element calculation results, the tensile stress of the fourth girder between piers D₃ and D₄ exceeds the design value of the concrete (the tensile stress design value of the D-ramp concrete of Suzhou Bridge is 1.8 MPa). In the results of the most-unfavorable-condition method, the main girder, caused by unusual differential deformation between pier D₄ and D₃, produces a large tensile force in the middle. The results of the two methods are consistent.

4. Conclusions

In this paper, the most-unfavorable-condition method for bridge-damage detection and analysis is proposed. Firstly, the 3D deformation model of the area was built using Green’s spline interpolation method to obtain the time-series deformation information of PS points on bridge piers. Then, the most unfavorable deformation position and occurrence time were obtained based on the most-unfavorable-condition-method using time-series deformation information from the 3D deformation model. Finally, the validity and reliability of the calculated damage location and occurrence time were verified by the existing records and finite element analysis. The following conclusions are drawn:

(1) The 3D deformation field obtained by using Green’s spline interpolation shows that the deformation of foundations in the Suzhou Bridge area is increasingly serious, and an obvious deformation funnel is formed at the pier position;

(2) The most unfavorable deformation of the pier is analyzed by using the most-unfavorable-condition method. The maximum deformation rate and the corresponding deformation values of PS points on bridge piers were obtained by using the algorithm of the most-unfavorable working conditions I and II. The damage location of the most-unfavorable differential deformation was between piers D₃ and D₄, and the occurrence time was 17 July 2017.

(3) Through the existing records and the finite element analysis of the damage at the fourth girder (between piers D₃ and D₄), the proposed most-unfavorable-condition method for bridge-damage detection and analysis is verified accurately and reliably.

(4) At the same time, it should also be noted that in Equation (17), the large relative deformation had occurred on piers D₄ and D₅. Although no damage has been found so far, the girder between piers D₄ and D₅ should be classified as a potential damage area and should be closely monitored.