Monitoring subsidence rates along road network by persistent scatterer SAR interferometry with high-resolution TerraSAR-X imagery

Ground subsidence is one of the key factors damaging transportation facilities, e.g., road networks consisting of highways and railways. In this paper, we propose to apply the persistent scatterer synthetic aperture radar interferometry (PS-InSAR) approach that uses high-resolution TerraSAR-X (TSX) imagery to extract the regional scale subsidence rates (i.e., average annual subsidence in mm/year) along road networks. The primary procedures involve interferometric pair selection, interferogram generation, persistent scatterer (PS) detection, PS networking, phase parameterization, and subsidence rate estimation. The Xiqing District in southwest Tianjin (China) is selected as the study area. This district contains one railway line and several highway lines. A total of 15 TSX images covering this area between April 2009 and June 2010 are utilized to obtain the subsidence rates by using the PS-InSAR (PSI) approach. The subsidence rates derived from PSI range from −68.7 to −1.3 mm/year. These findings show a significantly uneven subsidence pattern along the road network. Comparison between the PSI-derived subsidence rates and the leveling data obtained along the highways shows that the mean and standard deviation (SD) of the discrepancies between the two types of subsidence rates are 0.1 and ±3.2 mm/year, respectively. The results indicate that the high-resolution TSX PSI is capable of providing comprehensive and detailed subsidence information regarding road networks with millimeter-level accuracy. Further inspections under geological conditions and land-use categories in the study area indicate that the observed subsidence is highly related to aquifer compression due to groundwater pumping. Therefore, measures should be taken to mitigate groundwater extraction for the study area.


Introduction
As the primary mode of transport, the road network, which consists of highways and railways, handles the majority of traffic volume [1,2]. Maintaining the sustainability and stability of these transportation infrastructures is crucial to maintaining traffic safety [3]. Previous investigations indicate that ground subsidence caused by tectonic movements or anthropic activities, such as groundwater overuse, is a major concern in land-use planning, infrastructure sustainability evaluation, and engineering construction [4]. Subsidence (especially subsidence troughs) along highways and railways impairs the sustainability and stability of these transportation infrastructures, thus resulting in great risk in traffic safety [3]. Therefore, the effective and accurate monitoring of subsidence along road networks is necessary in preventing these negative effects. This is particularly important in urban areas where highway and railway lines are highly integrated and dense.
In the last decades, the subsidence of highways and railways has generally been measured by conventional pointbased surveying techniques, e.g., leveling, wire-flex extensometer, and global positioning system [4]. Nowadays, the differential synthetic aperture radar interferometry (DIn-SAR) [5,6] has exhibited great potential as a newly developed geodetic technique to map ground subsidence. The spatial resolution and coverage of DInSAR are advantageous compared with conventional methods. However, the DIn-SAR technique is inevitably affected by spatiotemporal decorrelation [7], topographic errors [5], and atmospheric delay [8], which significantly reduce the accuracy of DIn-SAR measurements. To overcome these drawbacks, Ferretti et al. [9] proposed the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) method. Numerous PS-InSAR (PSI) approaches have been developed since then and have been tested by using moderate spatiotemporal resolution synthetic aperture radar (SAR) data, such as ERS-1/2 and ENVISAT ASAR C band SAR images [10][11][12]. In spite of the differences of data processing strategies among the different PSI approaches, PSI techniques generally track the deformations of stable point-like targets, i.e., the persistent scatterers (PSs) [9]. The PSI method employs strategies of specialized modeling, and estimates the deformation rate and topographic error at a given target on the basis of multitemporal InSAR phase observations of the target. Therefore, it is possible to separate the deformation information from other phase components (e.g., topographic errors and atmospheric delay). The PSI method has been proved very useful for monitoring subsidence in urban areas [11][12][13][14]. More advantages of PSI approach can be found in [9][10][11][12][13][14].
In recent years, the InSAR techniques (both DInSAR and PSI) have been used to monitor deformations related to highways and railways. In 2006, the U.S. Department of Transportation reported the results of InSAR applications for deformation monitoring of slopes nearby the highway transportation projects connecting different states [15]. Shan et al. [16] used the InSAR technique to extract deformations of an expressway and road area in an isolated permafrost area of China in 2012. Froese et al. [17] used the InSAR technique to manage risks associated with ground movements along the railway corridors. Ge et al. [18] mapped the ground subsidence along the Beijing-Tianjin high-speed railway using InSAR technique, and studied the impact of the subsidence on the railway. These previous studies of using InSAR technique to monitor deformations related to highways and railways utilized the moderate-resolution SAR images as mentioned before, and they analyzed the capability of InSAR applications in these monitoring activities.
The recent X-band (that has a wavelength of 3.1 cm) radar sensor onboard the German satellite TerraSAR-X (TSX) is capable of producing SAR images at high spatial and temporal resolutions of approximately 1-3 m and at 11 days (i.e., satellite repeat cycle), respectively [13,14]. Thus, data availability for subsidence detection by PSI is extended [13,14,20]. First of all, the high spatial resolution of TSX can remarkably increase PS density, and thus leading to more detailed monitoring of ground movements. Furthermore, the short repeat cycle and short radar wavelength of TSX make it highly sensitive to ground movements, which means higher capability of tracking the slowly-accumulated ground subsidence. In this paper, we propose to apply PSI approach with the use of high-resolution TSX images to monitor the comprehensive subsidence rates of the road networks in urban areas. The PSI approach includes interferometric pair selection, interferogram generation, PS detection, PS networking and neighborhood differencing, phase parameterization, and subsidence estimation strategies. The Xiqing district in southwest Tianjin (China), which contains multiple highways and one railway line, is selected as the study area. A total of 15 TSX images obtained from this area between April 29, 2009 and June 2, 2010 are utilized to extract subsidence rates along the road network by using the TSX PSI approach. To assess the accuracy of the subsidence derived by using the PSI approach, the subsidence of the road network as measured by leveling is used as reference datum for comparative analysis.
2 Study area and data source

Study area
We select the Xiqing district in Tianjin (China) as the study area to investigate the applicability and potential of the TSX PSI approach in subsidence inversion along road networks. Figure 1 shows that Tianjin is located in the northeastern part of China, bordering Beijing at northwest, the Hebei province at northeast, and the Bohai Bay at east [19,20]. Being one of the largest and most important industrial cities in north China, Tianjin suffers from water shortage due to its natural geographic condition and semiarid climate [19]. Thus, a large amount of groundwater (especially the deep phreatic water) has been exploited to meet industrial and agricultural needs, resulting in severe land subsidence in several areas of Tianjin [19].
As shown in Fig. 1, the study area (Xiqing) is located in the southwestern part of Tianjin, bordering the Tianjin urban area at northeast. Figure 2 shows the averaged SAR amplitude image of the study area. In recent years, several new industrial parks have been constructed and launched for daily production in Xiqing. To meet the transportation needs from the increased industrial production, several new highways were constructed in the recent years. As illustrated in Fig. 2, multiple highways and one railway line are located in the study area, forming a road network. As a result of the progress in industrial and agricultural production, much of the groundwater in this area is being exploited to meet the increasing water needs for production activities. This excessive exploitation of groundwater may result in uneven subsidence.
Uneven subsidence can cause deformation on the highway surface and the railway tracks [1,2]. The dynamic and static loads from heavy vehicles (carrier trucks and trains) should also be considered. These additional burdens lead to roadbed compression (i.e., subsidence), resulting in the distortion of the road surface and railway tracks.

Data source
To obtain the subsidence rates of the road network in the study area, the 15 TSX SAR images acquired between March 27, 2009 and June 2, 2010 are used in PSI processing. All images are provided in single look complex format, with a pixel spacing of 1.36 m in slant range (i.e., 2.07 m in ground range) and 1.90 m in azimuth. Figure 2 shows the averaged TSX SAR amplitude image of the study area, which is approximately 11 9 12 km 2 . The roads, highways, railway line, buildings, fishponds, and crop parcels can easily be identified from the high-resolution TSX SAR amplitude image. As annotated in Fig. 2, 20 leveling points (LPs) were deployed along four highways, namely, Xiqing road (AA 0 ), Zhongbei road (BB 0 ), Jinjing road (CC 0 ), and South Haitai road (DD 0 ). Three leveling campaign epochs were conducted on all of the LPs. The acquired leveling data were used to validate the subsidence rates derived using TSX PSI.
The image acquired on September 30, 2009 was selected as the reference image and all other images were coregistered and resampled into the same grid space as the reference image. In this study, we first generate 105 interferometric pairs through a full combination of the 15 TSX images. The spatial (perpendicular) and temporal baselines (SB and TB) of the interferometric pairs range from 3.8 to 320 m and from 11 to 429 days, respectively. To reduce the spatial decorrelation and residual topographic errors (both of them are positively related to spatial baselines of the interferometric pairs) in the interferometric phases, only interferometric pairs with SBs shorter than 150 m are considered for processing. A total of 53 interferograms are generated. Table 1 lists the general information (acquisition dates of master and slave images, SBs, and TBs) of the interferograms. We generate 53 differential interferograms by using the two-pass differential InSAR method with the use of the digital elevation model (DEM) derived by the shuttle radar topography mission (SRTM). One potential disadvantage of differential InSAR with high-resolution TSX SAR images is lacking of external DEMs with comparable resolution, which is possible to induce extra topographic errors. Taking into consideration of the difference between resolution of the external DEM and the TSX SAR imagery, the external DEM was first oversampled to the same grid space as the reference TSX SAR image. Moreover, we rejected the interferometric pairs with longer perpendicular baselines to further confine the related topographic errors.

Detection of PSs
We detect PSs by following method proposed by Ferretti et al. [9]. A pixel is a PS if it satisfies the following empirical criteria:

Map of China
where D amp is the amplitude dispersion index (ADI) [8]; a and r amp are the mean and standard deviation (SD), respectively, of the amplitude time series at a pixel; and A and r A are the overall mean and SD of the amplitude values of all the pixels in the mean amplitude image. The first criterion in Eq. (1) indicates that the pixels with smaller ADI are more temporally stable in radar backscattering than those with higher ADI, whereas the second criterion postulates that pixels with higher amplitude values tend to be more temporally coherent [12].

PS networking, phase parameterization, and subsidence estimation
After PS detection, all PSs are connected to form a Delaunay triangulation network, which is taken as the subsidence observation network [11,20]. Phase modeling is based on the concept of neighborhood differencing [11,12,20] applied to each of the links in the Delaunay triangulation network. Given N differential interferograms, the phase values at two neighboring PSs (e.g., p and q) extracted from the ith differential interferogram can be expressed by [12]: where U Ã i is the wrapped phase at the PSs; v Ã and e Ã are the subsidence rates and the elevation residuals (due to uncertainties in the SRTM DEM used) at the PSs, respectively; B T i is the TB of the ith interferogram; B ? i;Ã is the SB of the PSs in the ith interferogram; k is the radar wavelength (3.1 cm for the TSX system); R Ã and h Ã are the sensor-to-target range and the radar incident angle at the PSs, respectively; / _ Ã i is the residual phase consisting of the nonlinear subsidence, atmospheric artifacts, orbit errors,  [12]. N differential equations can be obtained by using N differential interferograms. In each link, t and n can be estimated by maximizing the following objective function [12]: where R, h, and B ? i are the mean sensor-to-target range, the mean radar incident angle, and the mean perpendicular baseline between two PSs, respectively; c is the model coherence (MC) of the link; and j ¼ ffiffiffiffiffiffi ffi À1 p . t and n can be derived by searching within a given solution space (e.g., -5 to 5 mm/year for t and -20 to 20 m for n when high-resolution TSX images are used) to maximize the MC. In this study, the searching procedure was carried out with step values of 0.01 mm/year and 0.02 m for t and n, respectively.
Once the subsidence rate and elevation residual increments of all the links are estimated by using Eq. (4), the Delaunay triangulation network can be treated by the weighted least squares (LS) adjustment to estimate the subsidence rates and the elevation residuals of all the PSs. The square of the maximized MC value of each link is considered the weight [12]. An LP with subsidence rate obtained through leveling measurements can be considered as a reference point for the LS adjustment. In this paper, we focus on analyzing the PSI-retrieved subsidence along road networks. A detailed discussion on PSI approach is beyond the scope of this work and can be found in [11][12][13].

Subsidence rate map and interpretation
The subsidence rates at all PSs were extracted from 53 differential interferograms by using the PSI approach presented in Sect. 3. One corner reflector with leveling measurements known was taken as the reference point (RP) while carrying out the LS adjustment. Figure 3 shows the generated subsidence rate map of the entire study area and the position of the RP (the black square). The map comprises two layers: the color-coded layer; and the base map layer (the averaged amplitude image). The color-coded layer, which includes a scale bar, presents the magnitude and distribution of the estimated subsidence rates of all the PSs. Very dense PSs were detected in the study area due to the high-resolution of TSX imagery. The total number of PSs is 686,598, while the averaged density of PSs is 5,201  Fig. 3 reveal that the subsidence pattern in the study area is highly associated with local land-use categories. Four typical areas are selected as examples for further analysis. The four selected areas are identified by the dashed rectangles, namely, P 1 , P 2 , P 3 , and P 4 . Figure 4 shows the relevant optical images. To provide a good presentation of subsidence related to P 1 , P 2 , P 3 , and P 4 , we list their maximum, minimum, and mean subsidence rates in Table 2. The highest subsidence rate is observed in P 1 , with the maximum and mean magnitudes of -68.8 and -50.6 mm/year, respectively. P 3 and P 4 have relatively higher subsidence rates with the mean values of -33.9 and -44.8 mm/year, respectively. P 2 has relatively lower subsidence rates as compared with P 1 , P 3 , and P 4 . The maximum and mean subsidence rates of P 2 are -40.7 and -26.2 mm/year, respectively.
According to our in situ investigations, P 1 contains a thermal power plant (marked by the dashed oval), multiple industrial factories (marked by the dashed rectangle) and some residential quarters (marked by the dashed polygon). P 3 and P 4 mainly comprise industrial parks, whereas P 2 consists of residential quarters. The observed subsidence in these areas is caused by groundwater extraction for productive and domestic water needs. The aquifer in the study area was formed in the Neogene Period. This layer mainly belongs to the Neogene period and the Quaternary sedimentary formation period [4,19]. Geologically, the aquifer is a thick clayey stratum consisting of grits and loose, or semiloose argillaceous sediments [19]. Excessive exploitation of groundwater can result in severe water table depression. This depression further leads to interstitial water runoff and the consequence is aquifer compression which accounts for the subsidence. The heterogeneous subsidence observed in the study area (as shown in Fig. 3) is attributed to the various groundwater pumping intensities in different regions.

Closer inspections of the subsidence related to the road network
We select four highways and one railway line to further analyze the subsidence phenomena related to the road network (see Fig. 2). The selected highways and railway line, namely, Xiqing road, Zhongbei road, Jinjing road, South Haitai road, and Jinpu railway, are annotated by AA 0 , BB 0 , CC 0 , DD 0 , and EE 0 , respectively. We first estimated the full resolution (i.e., the original resolution as per the TSX image) subsidence rates along the highways and the railway line by using the Kriging interpolation method. We then extracted the subsidence rate profiles from the centerline of each highway and from the railway line. Figure 5a-e shows the subsidence rate profiles related to each highway and the railway line. Figure 5 exhibits the uneven subsidence distribution along the transportation facilities being studied. Table 3 lists the maximum, minimum, and mean subsidence rate values of each highway and the railway line. The Jinpu railway and the Xiqing highway have subsidence rates that are significantly higher than those of other highways. Further inspections of Figs. 2 and 4 indicate that dense industrial parks and residential quarters are located along both sides of the selected highways and railway line. The local land-use category of each highway and railway line is annotated in Fig. 5. Figure 5 shows that the subsidence troughs revealed by the profiles are mostly located near industrial areas, whereas the sections with relatively lower subsidence rates are near residential areas, agricultural areas, and commercial and storage areas. A closer inspection of Figs. 2, 4, and 5a, e shows that the maximum subsidence rates along the Xiqing road and the Jinpu railway are observed near the thermal power plant. As discussed in Sect. 4.1, groundwater pumping due to Fig. 4 Optical images related to P 1 , P 2 , P 3 , and P 4 productive and domestic water needs in the residential and industrial regions is accountable for the subsidence in and around these areas. Water pumping has indeed affected the neighboring highways and the railway line. Therefore, measures should be taken to optimize groundwater use planning and to restrict groundwater pumping. However, groundwater pumping is not the only factor affecting the transportation facilities being studied. Another factor responsible for the uneven subsidence of the highways and the railway line is the dynamic load from heavy vehicles,    For more comprehensive understanding of the uneven subsidence, the maximum slope changes of subsidence rates (mm/year) along different sections (being divided by the black vertical lines in Fig. 5) of the selected highways and the railway line are estimated from the subsidence rate profiles and are shown in Table 4. The maximum slope changes in industrial areas show the highest magnitudes ranging from 14.3 to 23.8 mm/year, and those in residential areas have the lowest magnitudes, while the subsidence rates in agricultural and commercial areas reveal medium maximum slope changes. Such circumstances demonstrate that the subsidence in industrial areas is highly dominated by the heterogeneous pattern, while the subsidence in areas with other land-use categories (e.g., agricultural land, residential land, and commercial land) is relatively less heterogeneous.
It is worth noting that dense PSs were identified along the studied highway and railway lines in our study area (see Fig. 3). One reasonable explanation is that the street lamps, stones, and fences distributed along the highway and railway lines can be easily and individually identified by the high-resolution TSX imagery (i.e., TSX SAR images have small pixel size). This characteristic is advantageous of providing more detailed subsidence information in PSI analysis, which makes the subsidence troughs along the road network be clearly observed (see Fig. 3). Moreover, the TSX imagery has wavelength of 3.1 cm, which is shorter than other regular used SAR systems such as ALOS PALSAR L band (23.6 cm) and ENVISAT ASAR and ERS1/2 C band (5.6 cm) SAR imagery. Shorter wavelength corresponds to relatively higher observation sensibility to ground movements, resulting in more precise monitoring of the land subsidence. The accuracy of the subsidence rates derived by TSX PSI will be validated by comparing with leveling measurements (see Sect. 4.3).

Subsidence validation with leveling data
For validation purposes, we used the in situ leveling data to assess the accuracy and reliability of the subsidence rate measurements derived by using the TSX PSI. We compared the subsidence rates derived by TSX PSI with those derived through leveling, i.e., subsidence rates obtained at the LPs (as marked in Fig. 2) by three precise leveling campaign epochs. Figure 6 shows the comparison between the two types of subsidence rates at the LPs. Statistical analysis shows that the discrepancies between the two types of subsidence rates at all LPs range from -5.1 to 4.2 mm/year, while those at the LPs on each highway range from -3.0 to 3.2 mm/year, -3.2 to 4.1 mm/year, -5.1 to 2.7 mm/year, and -2.4 to 4.2 mm/year, respectively. The mean and SD of the discrepancies at all LPs are 0.1 and ±3.2 mm/year, respectively. The results show that the high-resolution TSX PSI is useful in monitoring subsidence related to highways and railways at millimeterlevel accuracy.

Conclusions
We propose to apply PSI approach using high-resolution TSX imagery to extract the regional scale subsidence related to the road network. The Xiqing district in southwest Tianjin (China) is selected as the study area, and the subsidence rates along the highways and the railway line in this area are obtained by using the PSI approach utilizing 15 TSX SAR images acquired between April 29, 2009 and June 2, 2010. The quality of the subsidence rates derived using the PSI approach is assessed with use of ground truth data obtained by precise leveling.
The subsidence rate map shows that the study area is highly affected by uneven subsidence. The maximum, minimum, and mean subsidence rates are -68.7, -1.3, and -33.7 mm/year, respectively. The subsidence rate map, the optical images, and our in situ investigations in the study area demonstrate that the subsidence pattern is highly related to land-use categories (i.e., industrial lands, residential lands, agricultural lands, commercial lands, etc.). The uneven subsidence in this area is attributed to aquifer compression caused by the unbalanced pumping of groundwater in different types of lands. Further inspections of the subsidence rate profiles of the transportation facilities and the land-use categories around them indicate that groundwater pumping has indeed affected the road network being studied. Therefore, measures should be taken to optimize groundwater use planning and to mitigate groundwater extraction. The comparison between the subsidence rates derived using PSI and those obtained by leveling shows that the discrepancies between the two types of subsidence rates at all the LPs are within -5.1 to 4.2 mm/year. The mean and SD of the discrepancies are 0.1 and ±3.2 mm/year, respectively.
The results demonstrate that the PSI method based on high-resolution TSX imagery is useful in detecting road network subsidence. The subsidence measurements derived using the TSX PSI method have millimeter-level accuracy. The high-resolution and regional-scale characteristics of TSX PSI are helpful in detecting subsidence troughs that may damage the road network. The TSX PSI method is also advantageous in revealing the mechanism and origins of the subsidence phenomena.