Detecting land subsidence near metro lines in the Baoshan district of Shanghai with multi-temporal interferometric synthetic aperture radar

Land subsidence is a major factor that affects metro line (ML) stability. In this study, an improved multi-temporal interferometric synthetic aperture radar (InSAR) (MTI) method to detect land subsidence near MLs is presented. In particular, our multi-temporal InSAR method provides surface subsidence measurements with high observation density. The MTI method tracks both point-like targets and distributed targets with temporal radar backscattering steadiness. First, subsidence rates at the point targets with low-amplitude dispersion index (ADI) values are extracted by applying a least-squared estimator on an optimized freely connected network. Second, to reduce error propagation, the pixels with high-ADI values are classified into several groups according to ADI intervals and processed using a Pearson correlation coefficient and hierarchical analysis strategy to obtain the distributed targets. Then, nonlinear subsidence components at all point-like and distributed targets are estimated using phase unwrapping and spatiotemporal filtering on the phase residuals. The proposed MTI method was applied to detect land subsidence near MLs of No. 1 and 3 in the Baoshan district of Shanghai using 18 TerraSAR-X images acquired between April 21, 2008 and October 30, 2010. The results show that the mean subsidence rates of the stations distributed along the two MLs are −12.9 and −14.0 mm/year. Furthermore, three subsidence funnels near the MLs are discovered through the hierarchical analysis. The testing results demonstrate the satisfactory capacity of the proposed MTI method in providing detailed subsidence information near MLs.


Introduction
Due to the potential effects of land subsidence geo-hazards on metro lines (MLs) in terms of human lives and economic losses, great effort has been put forth to develop sustainable solutions for these hazards. To date, several methods have been investigated, including leveling [1], global positioning systems [2,3], interferometric synthetic aperture radar (In-SAR) [4], and multi-temporal InSAR (MTI) [5]. The MTI is a newly developed technique and provides a sound tool to assess land subsidence on Earth's surface. It is capable of providing millimetric-precision subsidence rates. Moreover, the method offers a synoptic view of subsidence funnels. Therefore, many methods using MTI to detect slow-moving land subsidence have been presented.
To increase the observation density of surface subsidence measurements, an improved MTI method is proposed to jointly process both PTs and DTs. The MTI method uses both low-amplitude dispersion index (ADI) [6,14] and high Pearson correlation coefficient [15] to select all valid pixels corresponding to PTs and DTs from the image series. Furthermore, to control error propagation, a hierarchical analysis strategy is presented. The method is first introduced at a high level. It consists of three major steps. They are briefly described below.
First, the subset of potential pixels (SPP) of the PTs is selected by considering pixels with low-ADI values. The SPP is processed using a least-squared estimator on an optimized freely connected network to obtain the PTs. Second, the SPP of the DTs is selected by considering pixels with high-ADI values, which are classified according to the ADI intervals. To reduce error propagation, the DTs and their subsidence rates are detected using the hierarchical analysis strategy and quality assessments on each SPP subgroup. Finally, the nonlinear subsidence components for all valid pixels are estimated using phase unwrapping and spatiotemporally filtering of the phase residuals.
The proposed MTI method is used to detect land subsidence near MLs in the Baoshan district of Shanghai using 18 TerraSAR-X images acquired between April 21, 2008 and October 30, 2010. Due to the high-spatial resolution (approximately 2 m in both ground range and azimuth directions), the TerraSAR-X images are capable of providing detailed subsidence information when used in the newly modified MTI method.
The remainder of the paper is organized as follows. In Sect. 2, the basic MTI method theories are explained. In Sect. 3, the data and experimental results are discussed. Finally, several conclusions are shown in Sect. 4.

Method
In this section, the core models of the three major steps in the MTI method are provided. Unlike the aforementioned MTI methods [6][7][8][9][10][11][12][13], both temporal and spatial phase data correlations are jointly considered using ADI and Pearson correlation coefficient. First, the equations for the subsidence rate extraction at the PTs with low-ADI values are discussed. Then, to extract the subsidence rates at the DTs, the pixels with high-ADI values are processed using the hierarchical analysis strategy. Finally, the nonlinear subsidence components are obtained through spatiotemporally filtering all valid pixels of the phase residuals.

Estimating subsidence rates at the PTs
The SPP with ADI values \0.4 [9] is selected to estimate the subsidence rates at the PTs. The SPP is connected using an optimized freely connected network [12], which provides basic observations for estimating linear subsidence rates. A least-squared estimator [13,16] is then applied on each arc of the network to calculate the relative subsidence parameters (i.e., the subsidence rate and DEM error).
It is assumed that N interferograms are obtained from N ? 1 images. For each interferogram k, the differential model is applied to each arc of the freely connected network: where Du k ðx i ;y i ;x j ;y j Þ is the differential phase between two pixels, k is the wavelength of the sensor, Dv ðx i ;y i ;x j ;y j Þ is the relative subsidence rate between two pixels, T k is the temporal baseline, B k ? is the local perpendicular baseline, R is the slant range distance from the reference pixel to the sensor in the master image, h is the local look angle, De ðx i ;y i ;x j ;y j Þ is the relative DEM error, and Du k RE x i ;y i ;x j ;y j ð Þ is the differential phase residual. We assume that there are no phase ambiguities in the arcs. The assumption allows applying the least-squared estimation method. And the phase residuals calculated from least-squared estimation can be used to locate the arcs affected by phase ambiguities. The problematic arcs determined should be discarded for further analysis. The least-squared estimation is used by jointly considering N interferograms [13,16] where The most probable values of Dv and De are retrieved according to whereÁ denotes the estimated value. The vector of estimated residual phases can be represented by After the least-squared estimation, the spatiotemporal phase data correlation between two adjacent pixels connected by the freely connected network is measured using the modeling coherence factor c ðx i ;y i ; x j ;y j Þ [6]: If c ðx i ;y i ; x j ;y j Þ is \0.8 [14], the arc is discarded from the freely connected network, which reasonably assures that low-quality pixels are eliminated.
Besides, an easy and efficient assessment used to detect outlier was applied by Zhang et al. [16] to discard the arcs affected by phase ambiguities: where Max(Á) means the maximum value in a vector or matrix. If the threshold value in Eq. (7) is reached, the tested arc is considered to contain an outlier at 95 % confidence level [16]. The subsequent processing estimates the absolute subsidence parameters for all PTs using the relative linking parameters in the freely connected network. The most probable estimates of the two subsidence parameters are calculated using the following formulas provided by Liu et al. [12]: where B is the design matrix for weighted LS estimation; the elements of B are either 1 or -1. Moreover, L is the vector of Dv, R is the vector of residuals, X is the vector of unknown PT linear subsidence rates (or DEM errors), Q is the number of arcs, and M is the number of PTs.

Estimating subsidence rates at the DTs
To increase the observation density of the subsidence rate map, the pixels with high-ADI values are further treated on a group-by-group basis. All pixels with ADI values [0.4 are classified into G groups according to the ADI intervals, i.e., (0.4, 0.5), …, (0.
The PTs are treated as Group 0. After classification, the subsidence rates in the radar line-of-sight direction are hierarchically analyzed for each group.
Group i is used here as an example. The phase quality of each pixel in Group i is assessed through Pearson correlation coefficient [15] to determine whether the POI phase is consistent with any neighboring pixels already accepted in the previous i groups (i.e., Groups 0 to i-1). The Pearson correlation coefficient is described as follows: When the Pearson correlation coefficient between two vectors is [0.75, the two vectors are considered to be highly correlated [15]. Therefore, we introduce 0.75 as a threshold to discard the point pairs whose phase time series are uncorrelated. The relative subsidence parameters between these pixels are derived using Eqs. (1)-(6) and consequently added to those at the neighboring pixel to determine the absolute parameters. Besides Pearson correlation coefficient, two other constraints are also introduced to assess the pixel quality when the subsidence parameters are obtained through the aforementioned processing. They are spatial autocorrelation of the two subsidence parameters for the given pixel. The spatial autocorrelation is denoted by the localized standard deviations. For example, if there are L valid pixels near a pixel (x i , y i ) in a given window size (e.g., 50 9 50 pixels), the standard deviations of the two subsidence parameters are where (x j , y j ) are the valid pixels that are maintained from the previous i groups. The two constraints ensure the quality of the pixels around the investigated pixel. If the pixels are totally stable in radar backscattering, and Dv ðx i ;y i ;x j ;y j Þ can be estimated without error, the value of v ðx i ;y i Þ À v ðx j ;y j Þ À Dv ðx i ;y i ;x j ;y j Þ will be zero. A high value of dv indicates a greater disturbance introduced by the Detecting land subsidence near metro lines in the Baoshan district of Shanghai 139 noise. Meanwhile, de is used on the same purpose.
Empirically, dv and de should not be greater than 5 mm/ year and 10 m, respectively [12]. The pixel is discarded for further analysis if either constraint is not met. All pixels in Group i are processed in the same manner. The reserved pixels are used to analyze the subsequent groups. After all G groups are treated according to this group-by-group basis, all valid pixels are obtained. Moreover, the DTs are affirmed. A detailed layout of the algorithm is shown in Fig. 1 to clearly present the processing sequences for calculating the two subsidence parameters.

Estimating nonlinear subsidence components
To determine the complete subsidence evolution, the phase residuals on the valid pixels must be unwrapped and the nonlinear subsidence components must be extracted and subsequently added to the linear components. According to previous studies [9,17,18], the nonlinear subsidence components exhibit different characteristics from noise and the atmospheric phase screen in both space and time domains. Therefore, spatiotemporally filtering [9,17,18] on the unwrapped phase residuals is applied to estimate the nonlinear subsidence values. The phase unwrapping technique is applied on the phase residuals embodied in Eq. (5). Because three phase data dimensions are provided, i.e., two in space and one in time, the three-dimensional phase unwrapping method described by Hooper and Zebker [19] is used. The unwrapping results can be summarized as where the three terms on the right-hand side denote the unwrapped phase of the nonlinear subsidence, atmospheric phase screen, and decorrelation noise, respectively. After phase unwrapping, spatial low-pass filtering is applied for each interferogram given the condition that the decorrelation noise is spatially uncorrelated, and the other two terms are spatially correlated [9]. The filtered phase data contain only the nonlinear subsidence components and atmospheric phase screen, i.e., where (Á) LP S is the spatial low-pass filter operator, typically the convolution with a two-dimensional Gaussian function [9]. The width of the Gaussian is as narrow as 50 m to include the entire useful signal except the noise.
To obtain the nonlinear subsidence components, which are expected to be temporally correlated, the results from  Eq. (12) are filtered using a temporal low-pass filter to obtain where (Á) LP T is the temporal low-pass filter operator, which is typically a Kaiser temporal filter [18]. The cut-off frequency of the filter is empirically set to 25 %.
Finally, for the given pixel (x i , y i ) of the k-th image, the full resolution subsidence value S k f ðx i ;y i Þ is the sum of the linear components S k lðx i ;y i Þ and nonlinear components S k nlðx i ;y i Þ : 3 Experimental results and discussion

Study area and data source
The experiment is conducted using images acquired over the Baoshan district, Shanghai, to detect land subsidence near the No. 1 and 3 MLs. As shown in Fig. 2a, Shanghai is located on the deltaic deposit of the Yangtze river. Consolidation of the soft layer causes land subsidence in Shanghai, where land subsidence was first observed in China [20]. It is reported that the largest accumulated subsidence value from 1921 to 1965 was 2.63 m [20]. More than 400 km 2 of land was affected by this geological hazard. However, since 1965, Shanghai has been controlled by restricting groundwater extraction and several subsidence funnels can still be observed [21,22]. Land subsidence directly damages sewerage systems, roads, and buildings, etc. Moreover, land subsidence is a major factor affecting the stability of the Shanghai Metro. The Shanghai Metro is one of the fastest-growing rapid transit systems in the world [23]. The first line was opened in 1993. On October 16, 2013, the operating route length reached 468 km. Moreover, the planned route length is 970 km, which would make the Shanghai Metro the world's longest metro system. Shanghai Metro has become one of the most popular means of travel. On March 9, 2013, the Shanghai Metro set a daily ridership record of 8.486 million passengers. The social position of the MLs is important for maintaining their sustainability and stability.
The proposed MTI method is applied to detect land subsidence near the No. 1 and 3 MLs in the Baoshan district of Shanghai. Eighteen TerraSAR-X images are used in the experiment (See Table 1). The TerraSAR-X images were provided by Infoterra GmbH. The observation period encompassed April 21, 2008, to October 30, 2010. All the images were collected with an incidence angle of 26 degrees in HH polarization mode. The original datasets were provided as single look complex images with slant range pixel spacing of 0.91 m (2.04 m in ground range) and azimuthal pixel spacing of 1.97 m.  The September 09, 2009 image is chosen as the unique master image. All other images are co-registered to the master image to produce 17 valid interferograms. Furthermore, a DEM provided by the Shuttle Radar Topography Mission is introduced to eliminate the topographic phase in each interferogram. The differential InSAR-related functions are all accomplished using the GAMMA DIFF module [24].
The amplitude image (4,688 9 7,208 pixels, equivalent to 9.6 9 14.2 km 2 ) of the test area is shown in Fig. 2b. The stations distributed along the No. 1 and 3 MLs are depicted with red squares and green squares, respectively. Line 1 originates at Fujin Road station; nine stations (red squares) along Line 1 are studied. Line 3 begins at North Jiangyang road station; 12 stations (green squares) along Line 3 are examined. The two MLs are distributed beneath the main roads in this district. Subsidence of the two MLs may directly damage roads and buildings in the surrounding areas.

Results and discussion
To obtain the subsidence fields in the test area, the three major steps of the aforementioned MTI method are followed. First, the SPP with ADI B 0.4 are processed. The number of PT pixels extracted from the SPP is 303,219. Afterward, the pixels with ADI [ 0.4 are analyzed on a group-by-group. Seven other groups are analyzed in this experiment. All eight groups are shown in Fig. 3. Each sub-figure i is plotted with the valid pixels extracted from all previous i groups (i.e., Groups 0 to i -1). For example, the first figure is plotted with Group 0; the eighth figure is plotted with Groups 0 to 7. Table 2 shows the number of valid pixels for the eight groups. The number of useful pixels decreases as the ADI values increase. Groups 0 to 2 contribute 93 % of the useful information to the total valid pixels in the study area. Group 7 which contains ADI values greater than 1.0 provides very little useful information. Therefore, the hierarchical analysis strategy stops at Group 7. After hierarchical analysis, the total number of valid pixels is 1,245,526, which is 311 % more than provided by the initial group. This result indicates a large increase in the observation density of the subsidence rate map. The figure series shows that three subsidence funnels emerge when the hierarchical analysis is performed. The three areas are located in the northwestern, northeastern, and east-central parts of the test area. These subsidence funnels are clearly depicted in Fig. 3h. However, the subsidence funnels cannot be identified in Fig. 3a.
For further analysis of the subsidence field, close-ups of Fig. 3a, h are shown in Fig. 4a, b respectively. The metro stations along the No. 1 and 3 MLs are depicted with red and green squares. Moreover, S 1 , S 2 , and S 3 are the aforementioned three subsidence funnels. S 1 is located at the No. 3 ML origination point. North Jiangyang road is located in this area. S 2 is in a coastal region; compression of the silty marine clay, which is approximately 40 m thick and less than 70 m deep may be the primary reason for the subsidence [20]. S 3 is a subsidence funnel presented in many previous studies [22,25]. The three subsidence funnels are not clearly depicted in Fig. 4a, which is plotted with valid pixels from Group 0. However, the subsidence funnels are clearly observed in Fig. 4b, which is plotted with valid pixels from Groups 0 to 7. Pixels extracted from Group 0 are primarily PTs, such as rocks, iron fences, and corner reflectors. However, pixels extracted from Groups 1 to 7 are mainly DTs, e.g., asphalt roads, building tops, and bare lands. The results indicate that the DTs are stable in radar backscattering time series. Therefore, the DTs can provide more subsidence details than the PTs.
If no hierarchical analysis strategy is applied, meaning that all of the points with ADI \ 1.1 are processed simultaneously. 1,891,004 point candidates are selected and finally 1,557,626 points are maintained after using the method presented in Sect. 2.1. We display the result in Fig. 5. In this figure, the pixel number is greater than that of the hierarchical analysis result. That is because the point quality is too low to be controlled simply by considering  6) and (7). Therefore, some invalid points are included in the final subsidence rate map. Meanwhile, all of the three subsidence areas detected in Fig. 4b are not exhibited in Fig. 5. It is concluded that the useful information on the subsidence funnels is totally overwhelmed by the noise information. The subsidence rate map in Fig. 5 might be biased or misleading. However, the in situ ground truth leveling data are needed for indepth validation.
The subsidence rates on the stations are plotted in Fig. 6. Mean subsidence rates for the No. 1 and 3 MLs are -12.9 and -14.0 mm/year, respectively. Along Line 1, Fujin Road station has the maximum subsidence rate and Tonghe Xincun has the minimum subsidence rate. The differential subsidence rate between the two stations is approximately 6.4 mm/year. The subsidence rate vibration (r dv ) along Line 1 is 2.0 mm/year, which indicates that homogeneous subsidence exists along Line 1. A similar result is obtained for Line 3, along which the maximum subsidence rate difference between the Shuichan Road and South Changjiang Road stations is 8.4 mm/year. r dv is 2.9 mm/year along this ML. The first station of Line 3 is at North Jiangyang Road, a metro station with one of the largest subsidence rates in S 1 . The following four stations are also affected by small subsidence: Tieli Road, Youyi Road, Baoyang Road, and Shuichan Road. However, subsidence rates change slowly among these four stations. Specifically, r dv is 0.5 mm/year for the four stations, indicating that the station subsidence is uniform and that the stations are relatively stable. The last three stations also exhibit similar subsidence rates, i.e., r dv for the West Yinggao Road, Jiangwan Town, and Dabaishu stations is 0.6 mm/year.
After obtaining the linear subsidence rate map, the nonlinear components are calculated to provide the complete subsidence evolution for each pixel. For definitiveness and to avoid loss of generality, the nonlinear subsidence components are shown for three stations, i.e., Jiangwan Town, Fujin Road, and Hulan Road, in Fig. 7. Figure 7 shows that the first two SAR images (with acquisition dates of April 21, 2008, andAugust 20, 2008) are isolated from the image cluster. However, the two images provide useful information for the experiment. In Fig. 7, the nonlinear subsidence components indicate a slow vibration during the first three observation dates. Afterward, the nonlinear components suggest that the ground surface rises to a maximum in August and decreases to a minimum in January. The Shanghai climate data [26] show that precipitation might be the primary reason for the weather-related trend. The climate data from 1971 to 2000 show that the average precipitation amounts in July and December are 169.6 and 37.1 mm, respectively, which also correspond to the maximum and minimum rain capacities. The precipitation indicates that the pore-fluid pressure in the Shanghai soft clay is recharged in summer; the pore-fluid pressure provided by groundwater is transferred to the granular skeleton in winter.

Conclusions
An improved MTI method that includes linear subsidence rate extraction, hierarchical analysis strategy, and nonlinear subsidence extraction is presented to analyze the stations distributed along No. 1 and 3 MLs in the Baoshan district, Shanghai. The technique provides a method to derive useful information from a set of SAR images by tracking both PTs and DTs, therefore, increasing the subsidence information observation density. The linear subsidence rates are first obtained from the PTs. Thereafter, the hierarchical analysis strategy is introduced to F u j i n R o a d W e s t Y o u y i R o a d B a o 'a n H i g h w a y G o n g f u X i n c u n H u l a n R o a d T o n g h e X i n c u n G o n g k a n g R o a d P e n g p u X i n c u n W e n s h u i R o a d Subsidence rate (mm/yr) Metro Line 1 N o r t h J i a n g y a n g R o a d T i e l i R o a d Y o u y i R o a d B a o y a n g R o a d S h u i c h a n R o a d S o n g b i n R o a d Z h a n g h u a b a n g R o a d S o n g f a R o a d S o u t h C h a n g j i a n g R o a d W e s t Y i n g g a o R o a d J i a n g w a n T o w n Detecting land subsidence near metro lines in the Baoshan district of Shanghai 145 extract the linear subsidence rate from the DTs. The nonlinear subsidence values are estimated using 3D phase unwrapping and spatiotemporally filtering on the phase residuals. Furthermore, 18 TerraSAR-X images are used to determine the subsidence rate map for the test site and stations along the two MLs. The number of valid pixels increases from 303,219 to 1,245,526 after performing the hierarchical analysis. Groups 0 to 2 contribute 93 % of the useful information to the total valid pixels in the study area. The number of valid pixels decreases as the ADI increases. Group 7 provides very little useful information. The valid pixels extracted from Group 0 are primarily PTs, and those extracted from the subsequent groups are mostly DTs. The experiment shows that DTs are largely capable of providing high observation density subsidence information. Furthermore, three subsidence funnels are identified after the hierarchical analysis, indicating that our MTI method is  In-depth analysis is performed on the No. 1 and 3 MLs. The mean subsidence rates for these MLs are -12.9 and -14.0 mm/year, respectively. Moreover, vibrations of the relevant subsidence rate are 2.0 and 2.9 mm/year, respectively, indicating a homogenous subsidence pattern along the MLs. This finding may help in the management of metro line stations and in hazard prediction and prevention. Subsidence rate maps with high observation density are expected to facilitate metro line stability assessment for the transportation department in Shanghai.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.