Ground Subsidence Analysis in Tianjin (China) Based on Sentinel-1A Data Using MT-InSAR Methods

: Multi-temporal InSAR (MT-InSAR) methods have been widely used in remote sensing monitoring of ground subsidence, which occurs at many places around the world. Land subsidence, caused by excessive extraction of groundwater, has always been a problem to be solved in Tianjin, China. Although the subsidence in the urban area has been controlled at a low rate, the subsidence issue has not been e ﬀ ectively solved in the suburban area recently, which should be paid much attention. This paper aims to present two multi-temporal di ﬀ erential interferometry techniques, persistent scatterer (PS) and small baseline subset (SBAS), for monitoring the latest surface subsidence in a Tianjin study area on the basis of 20 Sentinel-1A images obtained from March 2017 to March 2019. Our research showed that the average velocity map obtained from the SBAS method closely followed the outcomes of the PS technique from the perspective of identifying similar subsidence patterns. Subsidence rate gradually increased from the urban area of Tianjin to the suburbs and high subsidence zones were mainly distributed at the junction of the Wuqing, Xiqing and Beichen districts. In the past two years, the annual average subsidence rate in the high settlement area mostly exceeded − 50 mm / year, which caused serious damage to local infrastructures. Besides, high-resolution remote sensing images combined with ﬁeld investigations further veriﬁed the successful application of MT-InSAR technology in Tianjin’s subsidence monitoring. E ﬀ ective ground subsidence control measures need to be taken as soon as possible to prevent the situation from getting worse.


Introduction
Land subsidence is a geological hazard due to natural and human factors [1,2], such as settlement of loose sediments, compression of unconsolidated strata and overexploitation of underground resources. It is usually a slow process and not easy to be noticed in the early stage. Once the subsidence continues to develop, it could cause the slow lowering of ground surface elevation [3], resulting in critical infrastructures damage [4]. Today, land subsidence has become a global problem and cases of subsidence have been detected in many places all over the world [1]. A number of cities like Mexico [5], Tianjin [6], Jakarta [7], Tokyo [8] and Venice [9] are facing the problem of land subsidence.
Traditional land subsidence monitoring methods include global navigation satellite system (GNSS), GPS, leveling, extensometer and geophysical investigation [10]. However, these means are Appl. Sci. 2020, 10 not only high-cost and laborious, but also generally restricted to a small spatial extent, failing to depict a spatial display of ground subsidence [10]. Interferometric synthetic aperture radar (InSAR) is a confirmed remote sensing method, making use of SAR images' phase information to monitor ground surface displacement. The standard differential interferometry SAR (D-InSAR) method can obtain qualitative information on the displacement from two interferograms. Nevertheless, there are intrinsic shortcomings that could limit its operational capability, such as aerial disturbance, temporal decorrelation and imprecision of the reference DEM [11]. To get over the limitations of D-InSAR, multi-temporal InSAR (MT-InSAR) was proposed in the early 21 st century for mapping ground subsidence with high spatial and temporal resolution [12]. Persistent scatterer InSAR (PS-InSAR) [13] and small baseline subset (SBAS) [14] are the two most widely utilized InSAR time series means.
In addition to them, other MT-InSAR approaches developed in recent years, including along track interferometry (ATI) [15], coherent pixels technique (CPT) [16], corner reflector interferometry SAR (CRInSAR) [17] and quasi-PS [18] to reduce decorrelations by combining numerous SAR images. Taking advantage of a set of pixels, the MT-InSAR technique can get the surface displacement history and the corresponding mean displacement rate [19]. By measuring the subsidence area and velocity, the MT-InSAR technique can not only get the early signal of geological hazards which is helpful for forecasting potential events, but can also monitor the trend of ground subsidence, allowing to detect the law of its change [20]. At present, the efficiency and precision of the MT-InSAR technique has been proven in many applications, especially in monitoring ground subsidence. Using the distributed scatterers InSAR method, Han et al. [21] analyzed the land subsidence based on COSMO-SkyMed images over the coast of Hangzhou Bay, with an area of 1400 km 2 , proving the MT-InSAR technique achieves an accuracy of ±5.0 mm/year. The study presented by Casu et al. [22] showed that the SBAS technique provides an estimate of the mean deformation velocity with a standard deviation of about 1 mm/year for the ERS SAR dataset. Similarly, Du et al. [4] found there is strong correlation between the MT-InSAR and leveling measurements and the root mean square error was about 6.1 mm/year, suggesting the reliability of InSAR-derived subsidence rates. At the same time, some studies focused on a comparative analysis of inversion results on land subsidence obtained by different MT-InSAR methods. For example, Yan et al. [23] applied two MT-InSAR approaches on 38 ENVISAT images to map Mexico City subsidence, finding that the subsidence information gained from the SBAS method portrays the local characteristics associated with urban constructions and basic facilities, while the PS method quantitatively characterizes the movement of single objects. Tiwari et al. [24] used PS and SBAS to estimate the surface deformation caused by an earthquake in Italy based on 14 ASAR images, indicating that both the methods successfully pick up measurement pixels and display a comparable displacement pattern. Huang et al. [25] investigated the ground subsidence of Tongzhou district in Beijing using 21 COSMO-SkyMed images, showing that the PS and SBAS methods have high consistency in the pattern of land subsidence, i.e., in the time series of corresponding points with the same geographical coordinates.
Tianjin is one of the four municipalities in China. Since the year of 2018, it has had a permanent population of about 15.6 million, with an urbanization rate of 83.15%. There is no doubt that human activities have increased the demand for geological and water resources. In the evolution of urbanization, Tianjin has experienced large scale groundwater extraction, inducing serious land subsidence [26]. Land subsidence has been affecting Tianjin for the past 60 years. From the year of 1959 to 2008, the maximum subsidence value reached 3.22 m, and the affected area was almost 8000 km 2 [27][28][29]. Nowadays, Tianjin is still the city with the most serious ground subsidence in China [12]. Since the 1980s, the ground subsidence rate in Tianjin's urban area has been controlled at a relatively low value (−10 mm/year) by importing the Luan River to Tianjin and limiting the extraction of groundwater. However, the region of groundwater exploitation has expanded to the periphery following the economic increase since the 1990s, forming new suburban subsidence centers [29]. Although it has attracted the attention of the local government, the subsidence problem in the suburbs has not been fundamentally solved due to various reasons such as insufficient control measures.
In order to understand the subsidence evolution in Tianjin in recent years, the goal of the paper is to (1) obtain and compare the remote sensing monitoring results of Tianjin ground subsidence using PS and SBAS methods; (2) analyze the characteristics of the latest spatial distribution of ground subsidence in urban and suburban areas and (3) put forward suggestions for the prevention and control of suburbs subsidence. This paper is organized as follows: Section 2 describes the materials and methods. It includes the description of the study area and the used datasets; the detailed steps of the PS and SBAS implementation; and the cross-comparison analysis method between the two MT-InSAR techniques. Section 3 discusses the main results of this study. Finally, the conclusions of this research are presented in Section 4.

Study Area
Tianjin is situated 130 km southeast of Beijing and is on the western coast of the Bohai Sea ( Figure 1). The city is located in the North China Plain (NCP), within the East Asia monsoon region, and is a typical city lacking surface and groundwater resources [3]. Tectonically, Tianjin is situated in a Meso-Cenozoic basin zone which was filled up by continental Tertiary and Quaternary deposits [4]. Considering the subsidence situation of Tianjin in recent years and the data availability, we choose the Tianjin urban area and its western suburbs as our study area. It covers an area of about 460 km 2 ( Figure 1), including parts of urban areas (Hongqiao, Hebei, Heping and Nankai districts) and suburbs (Xiqing, Wuqing and Beichen districts).

Data Description
To analyze the ground subsidence using MT-InSAR methods, Sentinel-1A single look complex (SLC) SAR images, SRTM DEM and precise orbit determination (POD) ephemerides are necessary in our research.
Sentinel-1A was launched on 3 April 2014 by the Europe Space Agency (ESA), which provides C-band imagery continuity following the retirement of ERS-2. Sentinel-1A includes four exclusive imaging modes with different resolutions and coverage. Its default acquisition pattern, known as terrain observation by progressive scan (TOPS), minimizing image scalloping and reducing noise, is designed for topographic mapping with a 12-day repeat cycle [30]. In this paper, for interferometry analysis to estimate the ground subsidence rate and subsidence time series of the Tianjin area, a dataset consisting of 20 Sentineal-1A SLC images in interferometric wide (IW) mode from 16 March 2017 to 18 March 2019 were selected. The datasets are provided by the ESA via Copernicus Open Access Hub. Specific parameters of the dataset are displayed in Table 1, and a detailed description about the images used in our study is provided in Table 2.  The 3-arc second (~90 m) SRTM3 version-4 DEM offered by NASA was applied to take out the topographic phases from the interferograms. Besides, the corresponding POD ephemerides released by ESA were adopted to carry out the orbital refinement and phase reflattening procedure.
In addition, for the purpose of explaining InSAR time series results, Google Earth images and the multi-spectral Sentinel-2A images with a 10 m spatial resolution were utilized. These highresolution and up-to-date images have an important impact in helping better interpret InSAR outcomes and comprehend the development of the ground subsidence in the study area.  The 3-arc second (~90 m) SRTM3 version-4 DEM offered by NASA was applied to take out the topographic phases from the interferograms. Besides, the corresponding POD ephemerides released by ESA were adopted to carry out the orbital refinement and phase reflattening procedure.
In addition, for the purpose of explaining InSAR time series results, Google Earth images and the multi-spectral Sentinel-2A images with a 10 m spatial resolution were utilized. These high-resolution and up-to-date images have an important impact in helping better interpret InSAR outcomes and comprehend the development of the ground subsidence in the study area.

In-SAR Techniques
In this study, the 20 Sentinel-1A SLC images were processed through the SARscape 5.2 module based on the ENVI platform. Before MT-InSAR processing, the whole Sentinel-1A images were cut to match the study area. Figure 2

In-SAR Techniques
In this study, the 20 Sentinel-1A SLC images were processed through the SARscape 5.2 module based on the ENVI platform. Before MT-InSAR processing, the whole Sentinel-1A images were cut to match the study area.

PS-InSAR
The method of PS-InSAR is based on identifying point targets, named permanent scatters (PS), in the multi-temporal interferometric SAR images [13]. The PS are natural or man-made targets that give out a stable signal phase, such as artifacts, parts of buildings and rock outcrops. They usually remain consistent for years of SAR images collected over the region with the same sensor, displaying dominant scattering by an object generally smaller in dimension than the resolution cell and they are exploited for obtaining millimeter land deformation and improving DEM accuracy [11,31]. The fundamental supposition of the PS method is that point objectives are all related in land scope and height across the SAR dataset, expecting to have an equal strength when processing images of different looks.
The PS-InSAR main processing steps are as follows [32]: (1) PS connection graph The super master image was selected automatically using the optimal selection algorithm and then 19 interferograms were generated by phase differencing between each slave image and the super master image accordingly ( Figure 3a). The perpendicular baseline was from 8.88 to 98.64 m, while the temporal baseline ranged from 24 to 732 d.
influence of the decorrelation effect of the temporal-spatial baseline on the interferogram (maximum temporal baseline: 750 d and maximum normal baseline: 50%) and each slave image was coregistered to the super master image. In this study, interferometric pairs with poor unwrapping and low coherence were eliminated. At last, a combination with 173 differential interferograms was produced ( Figure 3b).
(2) Phase unwrapping for the wrapped phase signal The minimum cost flow (MCF) network [34,35] was used to unwrap the phase with an unwrapping coherence threshold of 0.25. A few interferogram pairs with unwrapped phase error and low coherence were not taken into account in the next step.
(3) Refinement and re-flattening By choosing and refining stable ground control points (GCPs), the unwrapped phase was corrected using the method of the calculation for residual phase content and phase slope. The criteria for selecting GCPs is described by Abir et al. [36] and 25 GCPs were selected in the study area. Meanwhile, the topographic phase was removed based on the SRTM DEM and the nonlinear subsidence part was moved apart by means of spatial and temporal bandpass filtering [13]. Next, time series deformation in the study area was retrieved by subtracting the orbital ramps and atmospheric artifacts.
(4) Geocoding and raster to shape conversion Geocoding was done in the LOS direction using SRTM DEM. The parameter of velocity precision threshold was set to 1 mm/year. Finally, in order to compare with PS results, the time series displacement of SBAS was converted to a shape file with the WGS-84 coordinate system.

Cross-Comparison Analysis Between PS and SBAS Methods
The subsidence velocity maps obtained by PS and SBAS were classified with the same grade, so that they can be compared and analyzed on the same scale. According to the different subsidence velocities, the study area was divided into different levels of subsidence. Then, a statistical analysis of the area, maximum subsidence velocity and average subsidence velocity in different subsidence zones was performed. Besides, corresponding points were selected to further explore the crosscomparison between PS and SBAS based on regression analysis and time series analysis, respectively.
In the regression analysis, the subsidence velocity difference of the corresponding points was compared using the correlation coefficient (R). The calculation formula of R is as follows [37]: (2) PS interferometric process Oversampling was done by a factor of 4 to detect a larger number of pixels and avoid aliasing effects, which could help to judge the subpixel location of the main PS in a resolution cell [24]. Then, the phase due to topography was removed using STRM DEM from every interferogram following the theory of D-InSAR. Next, preliminary PS candidates were selected on the basis of the amplitude dispersion index, which is the ratio of the standard deviation and the average of amplitude values for every pixel in the image dataset. (

3) PS inversion
The phase obscurity was settled by the application of a linear temporal inversion method [5]. The parameters of displacement sampling and minimum displacement velocity were set to 1 mm/year and −200 mm/year, respectively. Then, displacement and time measures for phase and time monitoring at each point were acquired. In this research, we only considered the PS points with a coherence value greater than 0.75.

(4) PS geocoding
Based on SRTM DEM, geocoding was executed in the original satellite line of sight (LOS) orientation with a max number of points of 100,000 in each shape file. As a result, the time series displacement of PS was converted from the SAR coordinate system to the WGS-84.

SBAS-InSAR
The SBAS technique retrieves displacement rate maps and subsidence time series of coherent targets by efficiently implementing an appropriate combination of numerous interferograms with the relevant small baselines [12,14,31]. The resulted small baseline interferograms are fulfilled in a linear model. The combination of small baselines is determined by the vector including values of the unwrapped differential interferometric phase. The sampling rate will be increased by the solution of the linear model. Meanwhile, the backscatter is preserved and phase noise is reduced by setting a coherence threshold [31]. This method, devised especially to settle the unwrapping issue in fields of large and comparatively stable subsidence, is described by [33] in detail.
There are four main steps for SBAS-InSAR processing: (1) Generation of multiple differential interferograms Similar to PS-InSAR, the super master image was selected automatically considering the influence of the decorrelation effect of the temporal-spatial baseline on the interferogram (maximum temporal baseline: 750 d and maximum normal baseline: 50%) and each slave image was co-registered to the super master image. In this study, interferometric pairs with poor unwrapping and low coherence were eliminated. At last, a combination with 173 differential interferograms was produced ( Figure 3b).
(2) Phase unwrapping for the wrapped phase signal The minimum cost flow (MCF) network [34,35] was used to unwrap the phase with an unwrapping coherence threshold of 0.25. A few interferogram pairs with unwrapped phase error and low coherence were not taken into account in the next step.
(3) Refinement and re-flattening By choosing and refining stable ground control points (GCPs), the unwrapped phase was corrected using the method of the calculation for residual phase content and phase slope. The criteria for selecting GCPs is described by Abir et al. [36] and 25 GCPs were selected in the study area. Meanwhile, the topographic phase was removed based on the SRTM DEM and the nonlinear subsidence part was moved apart by means of spatial and temporal bandpass filtering [13]. Next, time series deformation in the study area was retrieved by subtracting the orbital ramps and atmospheric artifacts.
(4) Geocoding and raster to shape conversion Geocoding was done in the LOS direction using SRTM DEM. The parameter of velocity precision threshold was set to 1 mm/year. Finally, in order to compare with PS results, the time series displacement of SBAS was converted to a shape file with the WGS-84 coordinate system.

Cross-Comparison Analysis between PS and SBAS Methods
The subsidence velocity maps obtained by PS and SBAS were classified with the same grade, so that they can be compared and analyzed on the same scale. According to the different subsidence velocities, the study area was divided into different levels of subsidence. Then, a statistical analysis of the area, maximum subsidence velocity and average subsidence velocity in different subsidence zones was performed. Besides, corresponding points were selected to further explore the cross-comparison between PS and SBAS based on regression analysis and time series analysis, respectively.
In the regression analysis, the subsidence velocity difference of the corresponding points was compared using the correlation coefficient (R). The calculation formula of R is as follows [37]: where n is the total number of corresponding points, and x and y represent the values of the points taken from the corresponding points. When the value of R is 0, it means there is no linear relationship between the subsidence velocities of the corresponding points. On the contrary, when R is 1, it indicates a perfect linear relationship. The standard deviation (SD) value was used to evaluate the precision of mean subsidence rate derived from the MT-InSAR technique [12].
where n is the total number of pixels involved in the MT-InSAR method, x i stands for each temporal subsidence rate of the pixel, and µ represents the mean value of all the temporal subsidence rata. A small standard deviation indicates that the dispersion is low, and the corresponding accuracy is higher.

Comparison of Results Obtained by the Two MT-InSAR Techniques and Field Validation
After the processing of 20 Sentinel-1A SLC images based on PS and SBAS methods, the annual mean subsidence velocity maps from March 2017 to March 2019 were obtained, respectively (Figure 4a,b). The InSAR technique measures the projection of 3D ground displacement in the LOS direction [38]. Negative values correspond to an increase of the sensor-to-object slant range distance [12], revealing the ground is sinking, whereas positive values indicate surface uplift. As shown in Figure 4a,b, the surface displacement results derived from SBAS have a relatively consistent spatial distribution with PS, and the Tianjin study area has exhibited prominent uneven subsidence characteristics in the analyzed two years. The severe settlement area is mainly distributed in the northwest of the study area, concentrated at the junction of the Wuqing, Xiqing and Beichen districts. The subsidence velocity gradually increases from the urban district of Tianjin to the suburbs, roughly in the SE-NW direction. The extension of the total area affected by high subsidence (in red) is 40.84 km 2 according to PS analysis and 72.13 km 2 according to SBAS analysis (Figure 4a,b).
According to the surface deformation results of InSAR inversion, the study area can be divided into three levels of subsidence, namely low subsidence area, medium subsidence area and high subsidence area (Figure 4a,b). Based on Google Earth images and field investigations, it is found that high subsidence areas (24.75% of the study area) are mostly concentrated in emerging suburban industry areas. There are many industrial parks in the high subsidence area, where the mean subsidence rates are −52.08 mm/year (PS) and −56.97 mm/year (SBAS).
In sites TJ10 and TJ11 near Wangqingtuo town, in the high subsidence area, some old buildings have experienced serious sinking (Figure 5b,c), compromising the security of local residents. Remote sensing images and a field survey showed that a number of factories existed near the town and subsidence caused the uneven deformation of the ground, even inducing the inclination of the lampposts (Figure 5d). Low subsidence areas (35.61% of the study area) are mostly distributed in the urban area of Tianjin, and the surface subsidence rate is at a low level, with mean rates of −3.23 mm/year (PS) and −1.87 mm/year (SBAS). The medium subsidence area (39.64% of the study area) is distributed in the "buffer" zone between the Tianjin urban area and suburban industrial zone. There are basically no industrial parks or high-rise buildings in this area, while shopping malls, playgrounds, hospitals and schools are widely distributed. The mean subsidence rates are −25.40 mm/year (PS) and −19.76 mm/year (SBAS). In the medium subsidence area, we investigated the settlement of Hebei University of Technology and its surrounding area (TJ05-TJ08). There are obvious ground subsidence phenomena in the dormitory and teaching buildings on the campus (Figure 5f). As shown in Figure 5g,h, visible road sink and wall cracks have appeared around the university.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 22 ground subsidence phenomena in the dormitory and teaching buildings on the campus (Figure 5f). As shown in Figure 5g,h, visible road sink and wall cracks have appeared around the university.

Comparison of Corresponding Points
To further compare the displacement results derived from the two MT-InSAR techniques and analyze the accuracy of them, 1684 corresponding points ( Figure 6, Table 3), with the same geographical coordinates, were selected in the study area for statistical analysis.

Comparison of Corresponding Points
To further compare the displacement results derived from the two MT-InSAR techniques and analyze the accuracy of them, 1684 corresponding points ( Figure 6, Table 3), with the same geographical coordinates, were selected in the study area for statistical analysis.

Comparison of Corresponding Points
To further compare the displacement results derived from the two MT-InSAR techniques and analyze the accuracy of them, 1684 corresponding points ( Figure 6, Table 3), with the same geographical coordinates, were selected in the study area for statistical analysis.  The densities of the corresponding points in the high, medium and low subsidence areas increase sequentially, which are 3.26, 3.39 and 4.22 units/km 2 , respectively. The maximum subsidence rate in the corresponding points is −107.47 mm/year, derived from the method of SBAS. For all the corresponding points in the whole study area, the average subsidence rate obtained by PS is −22.99 mm/year, while the average subsidence rate of SBAS is −20.21 mm/year, slightly less than PS. The subsidence rates of PS and SBAS in the high subsidence areas are −53.25 mm/year and −53.92 mm/year, respectively, which are significantly higher than those in the medium and low settlement areas. The histogram of the PS-SBAS velocity difference is shown in Figure 7. The proportion of the corresponding points whose PS subsidence rate minus SBAS is less than 0 accounted for 87.11% (Figure 7), suggesting that the deformation rate obtained by PS is higher than that of SBAS in the study area. In addition, 64% of the corresponding points' velocity difference is between −5 and 0 mm/year. mm/year, respectively, which are significantly higher than those in the medium and low settlement areas. The histogram of the PS-SBAS velocity difference is shown in Figure 7. The proportion of the corresponding points whose PS subsidence rate minus SBAS is less than 0 accounted for 87.11% (Figure 7), suggesting that the deformation rate obtained by PS is higher than that of SBAS in the study area. In addition, 64% of the corresponding points' velocity difference is between −5 and 0 mm/year. The good agreement between these two different MT-InSAR-based measurements can be checked by plotting the PS subsidence rate versus SBAS, as shown in Figure 8. On the whole, the linear correlation coefficient R of all the corresponding points reached 0.974, indicating the annual subsidence rate acquired by the two methods was highly consistent. Separately, the correlation coefficients of the high, medium and low subsidence areas are 0.707, 0.948 and 0.962, (Figure 8), showing the existence of more discrete points in the high subsidence area than in the medium and low subsidence areas. The average absolute value of the corresponding points' velocity differences in the high, medium and low subsidence areas are 5.18, 4.14 and 3.62 mm/year, and the points with rate differences more than 10 mm/year account for 13.71%, 1.45% and 0%, respectively. Therefore, the consistency of the subsidence velocity derived from the two MT-InSAR techniques increases from the high subsidence area to the low subsidence area. The good agreement between these two different MT-InSAR-based measurements can be checked by plotting the PS subsidence rate versus SBAS, as shown in Figure 8. On the whole, the linear correlation coefficient R of all the corresponding points reached 0.974, indicating the annual subsidence rate acquired by the two methods was highly consistent. Separately, the correlation coefficients of the high, medium and low subsidence areas are 0.707, 0.948 and 0.962, (Figure 8), showing the existence of more discrete points in the high subsidence area than in the medium and low subsidence areas. The average absolute value of the corresponding points' velocity differences in the high, medium and low subsidence areas are 5.18, 4.14 and 3.62 mm/year, and the points with rate differences more than 10 mm/year account for 13.71%, 1.45% and 0%, respectively. Therefore, the consistency of the subsidence velocity derived from the two MT-InSAR techniques increases from the high subsidence area to the low subsidence area. In order to further explore the spatial distribution of the corresponding points with different velocities in the study area, the points were divided into three levels (0-5; 5-10; > 10 mm/year) according to the magnitude of their velocity difference ( Figure 9A). Statistics for the absolute value of velocity difference suggested that 75.12% of the corresponding points' difference is smaller than 5 mm/year, while only 3.56% is larger than 10 mm/year. The corresponding points with a velocity difference greater than 10 mm/year are mainly distributed in the high subsidence area. Besides, combined with remote sensing optical images, it was found that almost all these corresponding points appeared in the emerging industrial parks in the suburbs ( Figure 9B). Newly built factories are densely distributed in the industrial park, requiring a large amount of freshwater, mostly extracted from groundwater. Due to the excessive exploitation of groundwater in the industrial park, uneven ground subsidence occurs. Considering the different sensitivity of PS and SBAS methods to such unstable subsidence, it is likely that the subsidence rates gained at the same location will eventually be slightly different. Comparatively, the corresponding points' velocity difference is less than 5 mm/year in areas where the subsidence rate is relatively stable, such as high-rise buildings in urban area and transportation facilities (highways and railways, for example). In order to further explore the spatial distribution of the corresponding points with different velocities in the study area, the points were divided into three levels (0-5; 5-10; > 10 mm/year) according to the magnitude of their velocity difference ( Figure 9A). Statistics for the absolute value of velocity difference suggested that 75.12% of the corresponding points' difference is smaller than 5 mm/year, while only 3.56% is larger than 10 mm/year. The corresponding points with a velocity difference greater than 10 mm/year are mainly distributed in the high subsidence area. Besides, combined with remote sensing optical images, it was found that almost all these corresponding points appeared in the emerging industrial parks in the suburbs ( Figure 9B). Newly built factories are densely distributed in the industrial park, requiring a large amount of freshwater, mostly extracted from groundwater. Due to the excessive exploitation of groundwater in the industrial park, uneven ground subsidence occurs. Considering the different sensitivity of PS and SBAS methods to such unstable subsidence, it is likely that the subsidence rates gained at the same location will eventually be slightly different. Comparatively, the corresponding points' velocity difference is less than 5 mm/year in areas where the subsidence rate is relatively stable, such as high-rise buildings in urban area and transportation facilities (highways and railways, for example). It is necessary to analyze the cumulative displacement of the corresponding points with different velocities for knowing the different performance of the two MT-InSAR techniques. In our research, to compare the subsidence results of the two methods, six corresponding points (points a, b and c in Figure 9) are selected for time series analysis. Points a1 and a2 fall on two plants in the emerging industrial parks. Point b1 is located on a new high-rise building, and point b2 on Jinba Railway. Points c1 and c2 are located in the old residential area in the suburbs. The velocity difference of point a is more than 10 mm/year, of point b is between 5 and 10 mm/year and of point c less than 5 mm/year. At the same time, points a and c are located in the high subsidence area, while b in the medium area. From the displacement trend during time for the corresponding points (Figure 10), it can be seen that in general, the displacement of the points exhibits a sinking trend with a nonlinear behavior over the past two years. It is worth noting that the data fluctuation, for all the six points, is larger than for SBAS data, indicating that to a certain extent, the method of PS is more sensitive than SBAS for ground subsidence detection in our study area. In the high subsidence area, from March 2017 to March 2019, the cumulative displacement of point a1 detected by PS is −102.54 mm, larger than the It is necessary to analyze the cumulative displacement of the corresponding points with different velocities for knowing the different performance of the two MT-InSAR techniques. In our research, to compare the subsidence results of the two methods, six corresponding points (points a, b and c in Figure 9) are selected for time series analysis. Points a1 and a2 fall on two plants in the emerging industrial parks. Point b1 is located on a new high-rise building, and point b2 on Jinba Railway. Points c1 and c2 are located in the old residential area in the suburbs. The velocity difference of point a is more than 10 mm/year, of point b is between 5 and 10 mm/year and of point c less than 5 mm/year. At the same time, points a and c are located in the high subsidence area, while b in the medium area. From the displacement trend during time for the corresponding points (Figure 10), it can be seen that in general, the displacement of the points exhibits a sinking trend with a nonlinear behavior over the past two years. It is worth noting that the data fluctuation, for all the six points, is larger than for SBAS data, indicating that to a certain extent, the method of PS is more sensitive than SBAS for ground subsidence detection in our study area. In the high subsidence area, from March 2017 to March 2019, the cumulative displacement of point a1 detected by PS is −102.54 mm, larger than the displacement of

Comparison of the Two MT-InSAR Techniques
In the field of land subsidence monitoring, the MT-InSAR methods used in our research have significant advantages over traditional methods or standard D-InSAR. The establishment of a traditional ground monitoring system requires a lot of manpower and material resources. In Tianjin, the leveling survey monitoring system includes 2100 benchmarks, 12 continuous GPS stations and 12 groups of multilevel layer compression monitoring wells [29]. The coverage of the monitoring area continues to increase, and the leveling survey is conducted every year. Although traditional methods have high accuracy, they are still not suitable for large-scale investigations. The stranded D-InSAR method has its inherent shortcomings. For example, Zhou et al. [39] applied a geostatistical approach of combining D-InSAR with in situ leveling measurements to quantify land subsidence in Tianjin. However, because the D-InSAR method cannot address atmospheric effects, only by ensuring a sufficiently high coherence of the associated images can it be combined with a number of in situ leveling data to reduce the measurement errors. Considering the different performance in identifying the coherent points, using only one MT-InSAR method may not necessarily provide a comprehensive

Comparison of the Two MT-InSAR Techniques
In the field of land subsidence monitoring, the MT-InSAR methods used in our research have significant advantages over traditional methods or standard D-InSAR. The establishment of a traditional ground monitoring system requires a lot of manpower and material resources. In Tianjin, the leveling survey monitoring system includes 2100 benchmarks, 12 continuous GPS stations and 12 groups of multilevel layer compression monitoring wells [29]. The coverage of the monitoring area continues to increase, and the leveling survey is conducted every year. Although traditional methods have high accuracy, they are still not suitable for large-scale investigations. The stranded D-InSAR method has its inherent shortcomings. For example, Zhou et al. [39] applied a geostatistical approach of combining D-InSAR with in situ leveling measurements to quantify land subsidence in Tianjin. However, because the D-InSAR method cannot address atmospheric effects, only by ensuring a sufficiently high coherence of the associated images can it be combined with a number of in situ leveling data to reduce the measurement errors. Considering the different performance in identifying the coherent points, using only one MT-InSAR method may not necessarily provide a comprehensive understanding of the subsidence in the study area. When studying ground deformation using the MT-InSAR technique, Huang et al. [25] found that due to the sparse distribution of the coherent points acquired from PS, it is not easy to analyze the overall regional deformation and trends. The results of different MT-InSAR methods can be mutually confirmed and complement each other. The cross-comparison analysis used in our research further suggested that the subsidence velocities at the same location obtained by PS and SBAS are likely slightly different, which may indicate which areas are more unstable. Therefore, the application of two methods for remote sensing inversion and cross-comparison of their results is helpful for an in-depth understanding of the land subsidence range and its development trend.
Similar to other studies, our research proves that the PS and SBAS methods are effective in monitoring land subsidence over a large spatial area [23,40,41]. By comparing the outcomes of the two approaches, they proved to be generally in good agreement in identifying the ground subsidence rate and the spatial distribution of the displacement in the Tianjin study area. The interpretation of the observed subsidence rates in Tianjin suggests the presence of deformation signals with different spatial extents and magnitudes, between 2017 and 2019.
Due to different algorithm mechanism, the two MT-InSAR techniques could gain different coherent points [8,11,23]. The PS technique takes advantage of evaluating stable targets in the multi-temporal interferometric SAR images [13]. Combining the remote sensing images and the distribution data of main roads in the study area, we can see that the PS points are mainly distributed in urban/ suburban buildings, railways (such as G8936, Jinba Railway, Jingshan Railway, Jinpu Railway), expressways (Jinghu Expressway, Jinbao Expressway, for example), bridges and other man-made objects (Figure 11a,b). There are more PS points in the urban area than in the suburban area because of the existence of denser man-made structures. On the contrary, PS points are relatively sparsely distributed in vegetation coverage areas such as farmland mainly characterized by the cultivated land and crops, while there are almost no PS points in rivers and lakes (Figure 11b). Unlike the PS-InSAR method requiring the targets with high coherence values, SBAS is a powerful means which can identify targets with sufficiently high coherence values [12,14]. The combination of small baselines relays on the vector involving values of the unwrapped differential interferometric phase [19]. Previous studies have shown that when comparing PS and SBAS methods, the method of SBAS can obtain considerably more coherent points than PS [23,42,43]. In our study, the annual LOS mean subsidence velocity maps obtained based on SBAS were composed of 912,385 coherent pixels, with a density of 1981 units/ km 2 , more than PS-InSAR returning 728,903 PS points with 1582 units/ km 2 . Especially in the high subsidence area (i.e., the suburbs of Tianjin), the density of coherent points acquired from SBAS is 1.32 times more than that of PS, and the extension of the subsiding area detected by means of SBAS is 43.38% larger than that of PS (Figure 4a,b). On the whole, SBAS can obtain the settlement of low coherence points in the high subsidence areas, which makes up for the shortcomings of the PS algorithm to a certain extent.
Our research showed that both the PS and SBAS methods can achieve good results in the subsidence inversion process of key zones such as emerging industrial parks, residential areas (old or new) and transportation facilities (railways, for example). Further, the subsidence value obtained by PS is generally greater than the SBAS value ( Figure 10). For example, the subsidence rate at point b1, located in a newly built high-rise building area, is −40.56 mm/year (PS) and −35.32 mm/year (SBAS); meanwhile, the subsidence rate at point b2, situated on the railway, is −41.33 mm/year (PS) and −35.82 mm/year (SBAS). In addition, in the areas with the same subsidence level, the subsidence rates gained from the two methods are slightly different due to their different sensitivity to ground objects with different structures. Both points a and c are in the high subsidence area. Point a is located in the emerging industrial parks, where the groundwater is unevenly extracted, and uneven subsidence occurs, resulting in a velocity difference between PS and SBAS greater than 10 mm/year. Point c is in the old residential area in the suburbs, where the land subsidence is more widespread and uniform with respect with the industrial area; here, the velocity difference between the two methods is less than 5 mm/year. The difference in subsidence rates obtained by PS and SBAS seems to depend on the greater capability of the PS technique to discriminate ground objects with different velocities. Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 22 Our research showed that both the PS and SBAS methods can achieve good results in the subsidence inversion process of key zones such as emerging industrial parks, residential areas (old or new) and transportation facilities (railways, for example). Further, the subsidence value obtained In our research, the standard deviations of the average subsidence rate (Figure 11c), calculated from every temporal subsidence based on ENVI SARscape, represent the precision of the mean subsidence rate for the SBAS method. It is worthy of note that the spatial distribution of PS points is in good agreement with the high-precision value of the SBAS inversion results (Figure 11b,c) and 68.45% of PS points are distributed in the area with a standard deviation value of smaller than 0.5 mm/year. In addition, in the PS map, there are a few regions with a serious loss of coherence located in the high and medium subsidence areas (I, II and III in Figure 11b,c) and correspondingly, the precision value of these areas in the SBAS results is relatively low, falling in the class 1.5-2 mm/year. Combined with remote sensing optical images, it can be seen that there is a large distribution of farmland in these areas, causing low coherence, as mentioned above. It can be concluded that in urbanized areas with artificial buildings, PS and SBAS can obtain enough coherent points, while in farmland areas, both the two methods suffer from low coherence, being the PS technique more sensitive.

Recent Ground Subsidence Characteristics and Its Influencing Factors in Tianjin
This paper proves that recently the ground subsidence velocity in most of the urban area of Tianjin is under the threshold of −10 mm/year, while there are still obvious serious settlement funnels in the suburban area (over −50 mm/year in local regions). This concurs with the findings of Guo et al. [12], who derived the land subsidence of Tianjin from 2015 to 2016 based on Sentinel-1A using the method of SBAS, showing that the subsidence rate of the central urban area varies from −6.1 to 3.9 mm/year and the subsidence rate of the high settlement area in the suburbs ranges between −126.3 and −65.7 mm/year. Due to Tianjin and the central government adopting a variety of measures to control subsidence [3,6,26] since the 1980s (such as deriving alternative surface water, restriction on groundwater extraction, manual recharge of aquifers, and establishment of a settlement monitoring network), the urban settlement in Tianjin is presently under control [27,29]. However, subsidence in the suburbs of Tianjin has not been effectively managed and the issue is becoming serious in recent years. For example, during the period from 1992 to 2000, land subsidence occurred in almost all the outskirts around the Tianjin urban center, consisting of the Wuqing, Beichen and Xiqing districts, and most subsidence rates were between −30 and −60 mm/year [44]. Since the 21st century, the subsidence in the western suburbs of Tianjin has continued to increase, and many researchers have performed inversion analysis of the area using In-SAR technology. Fan et al. [45] showed that, between 2003 and 2004, in Tianjin, the western zone is most severely affected, being Xiqing district the worst-hit area. Zhang [46] pointed out that the maximum subsidence rate of the settlement funnel in the western suburbs of Tianjin exceeded −80 mm/year from 2007 to 2010. This has been further confirmed by Zhang et al. [6], who found the south-western area of Wuqing district and western part of Beichen district displaying subsidence rates of up to −136 mm/year from 2016 to 2017. In our research, the latest land subsidence in Tianjin is basically consistent with previous studies. The land subsidence in the suburbs is still not under control, and some areas are even getting worse.
As a typical geological hazard in Tianjin, ground subsidence seriously threatens the safety of people's life and property. The subsidence has a long history and its settlement mechanism is very complicated. Based on the previous results, the causes of Tianjin subsidence can be summarized as natural and human factors. Natural factors consist of tectonic sinking and loose soil consolidation. Wang et al. [47] studied the background of tectonic subsidence in the Tianjin area, finding that the tectonic subsidence velocity is 1.7 mm/year approximately. The soil in most parts of Tianjin is under-consolidation, often with the natural water content greater than the liquid limit. Even without external load, this type of formation is still slowly consolidating under self-weight stress, causing a small amount of ground subsidence over a long period [28]. The Tianjin Land Subsidence Controlling Office showed that the natural average deformation value of the unconsolidated strata in Tianjin is about 1.0 mm/year [48]. In addition to natural factors, for many years, human activities, especially the excessive exploitation of groundwater resources, have always been the main influencing factors for Tianjin's land subsidence. Taking Wangqingtuo town (TJ10 and TJ11 in Figure 4b) as an example, its land subsidence rate has been high for many years [12,28]. As the largest bicycle production base in northern China, Wangqingtuo town has a total of 268 bicycle manufacturers with an annual output of 8 million bicycles. Many production workshops are densely distributed in the town's industrial park (Figure 5a). The production process, especially electroplating, requires a large amount of water, mostly derived from underground mining. Luo et al. [28] indicated that the shallow groundwater in Wangqingtuo town is relatively poor, and when it is excessively withdrawn, the surface formed by alluvial/diluvial sediments is prone to uneven subsidence. In our study, the town is still located in the high subsidence area, and its average subsidence rate in the past two years has exceeded −50 mm/year. Field investigations have found that long-term ground subsidence has caused some buildings to sink, and local residents have had to move away (Figure 5b,c).

Suggestions on Controlling Tianjin Suburbs Ground Subsidence
As previously described, recently, subsidence in the suburbs has been increasing. In order to control the ground subsidence, the following measures have to be taken: (1) To improve the regional ground subsidence monitoring network dynamically.
Although the Tianjin ground subsidence monitoring network is gradually improving in recent years (including 2100 benchmarks and 5900 km leveling route), the monitoring density in the suburbs is far from enough. In some newly discovered subsidence areas, there are no settlement monitoring facilities installed [49]. Level of monitoring density should be increased in high-rise buildings, public transportation and emerging industrial parks, where the ground subsidence situation is serious, and the monitoring system for layerwise marking should be established in time.
(2) To optimize the monitoring methods and improve the level of survey.
We suggest investigating the applicability of different monitoring methods in depth, choosing the appropriate monitoring means for different monitoring targets and exploring multi-technology integration methods, especially combining InSAR with traditional techniques to improve the monitoring accuracy and efficiency. A feasible "from surface to point" method is to first obtain the long-term sequence of ground subsidence development trends and spatial distribution characteristics based on InSAR technology at a large spatial scale in Tianjin, and then establish a monitoring network in key subsidence areas using traditional means, such as ground leveling, and GPS to achieve an accurate real-time measurement.
(3) To strengthen the supervision of laws and regulations, and increase publicity and education.
Though the local government has issued a number of management policies to control ground subsidence, it is still difficult to fully carry out, especially in the suburbs. Our site investigation found that the local residents know very little about the causes of ground subsidence and lack water-saving awareness. Therefore, it is necessary to strictly control groundwater exploitation, through the implementation of groundwater extraction permit systems and artificial recharge systems. Subsidence control rules should be imposed on local firms and residents, such as Management Measures for Controlling Land Subsidence in Tianjin. In addition, it is essential to establish an exhibition platform to disseminate scientific knowledge on the prevention of ground subsidence prevention to the public and enhance the awareness of disaster prevention and mitigation.

Conclusions
Our research presented the spatial-temporal distribution of land subsidence in a Tianjin study area based on 20 Sentinel-1A SLC SAR scenes collected from March 2017 to March 2019 using two MT-InSAR methods. Both PS and SBAS techniques successfully estimated the extent of land subsidence and recognized a similar subsidence pattern as a result of excessive extraction of groundwater in Tianjin. Both of the two methods suffer from low coherence, being the PS technique more sensitive. The subsidence value obtained by PS is generally greater than the SBAS value, and the magnitude of their velocity difference seems to depend on the greater capability of the PS technique to discriminate between ground objects with different velocities. The study area exhibited prominent uneven subsidence characteristics in the past two years. The subsidence velocity in the suburbs is much larger than that in the urban areas, and in particular, high subsidence areas are mainly at the junction of the Wuqing, Xiqing and Beichen districts with a subsidence velocity of over −50 mm/year. It is necessary for the local government and residents to act together to take effective measures to prevent increased subsidence in Tianjin's suburbs.