Monitoring Land Surface Displacement over Xuzhou (China) in 2015–2018 through PCA-Based Correction Applied to SAR Interferometry

: Land surface deformation in metropolitan areas, which can cause varying degrees of hazard to both human lives and to properties, has been documented for decades in cities worldwide. Xuzhou, is one of the most important energy and industrial bases in eastern China, and has experienced signiﬁcant land subsidence due to both excessive extraction of karst underground water and exploitation of mineral resources in recent decades. Furthermore, Xuzhou has recently undergone rapid urbanization in terms of urban expansion and underground construction, which could induce additional pressure on the urban land surface. However, most previous research on land surface deformation in the Xuzhou urban areas has been conducted based on traditional ground-based deformation monitoring techniques with sparse measurements. Little is known about the regional spatiotemporal behavior of land surface displacement in Xuzhou. In this study, a detailed interferometric synthetic aperture radar (InSAR) time series analysis was performed to characterize the spatial pattern and temporal evolution of land surface deformation in central areas of Xuzhou during 2015–2018. A method based on principal component analysis was adopted to correct artifacts in the InSAR signal. Results showed the correction strategy markedly reduced the discrepancy between global navigation satellite systems and InSAR measurements. Noticeable land subsidence ( − 5 to − 41 mm / yr) was revealed widely within the Xuzhou urban areas, particularly along subway lines under construction, newly developed districts, and in old coal goafs. Remarkable consistent land uplift (up to + 25 mm / yr) was found to have signiﬁcantly a ﬀ ected two long narrow areas within the old goafs since 2015. The possible principal inﬂuencing factors contributing to the land surface displacements such as subway tunneling, building construction, mining, underground water levels and geological conditions are then discussed.


Introduction
Land surface deformation induced by endogenous forces of the earth (e.g., volcanic events, earthquakes, landslides and collapses) and/or anthropogenic activities (e.g., urban construction, mining, oil and groundwater extraction) has been observed worldwide for decades and it has become

Study Area
Xuzhou, located in the northwest of Jiangsu Province in eastern China, is bordered by Shandong, Henan and Anhui Provinces (Figure 1a,b). It lies on the North China Plain in the flood plain of the Yellow River and in the southern margin of the mountainous area of southern Shandong Province. It is not only the central city of the Huaihai economic zone (Figure 1b), but it is also one of the core cities of the three metropolitan circles planned and constructed in Jiangsu Province. It is an important transportation and railway hub, as well as one of the primary energy and industrial bases in China. (b) Location of Xuzhou (outlined in red) in the Huaihai Economic Zone (outlined in yellow). The Xuzhou major geological disaster prevention and control Area (MPCA) is outlined in blue. The black rectangle represents the outline of (c). (c) Image of mean intensity of Sentinel-1A synthetic aperture radar data of the study area. Subway Lines 1-3 are indicated by red, blue and orange lines, respectively. Pink dashed lines indicate two synclines. Filled red triangles indicate the location of the Global Navigation Satellite Systems stations used in this study. Maps are projected in geographic Lat/Lon referenced to the WGS-84 datum.

Study Area
Xuzhou, located in the northwest of Jiangsu Province in eastern China, is bordered by Shandong, Henan and Anhui Provinces (Figure 1a,b). It lies on the North China Plain in the flood plain of the Yellow River and in the southern margin of the mountainous area of southern Shandong Province.
It is not only the central city of the Huaihai economic zone (Figure 1b), but it is also one of the core cities of the three metropolitan circles planned and constructed in Jiangsu Province. It is an important transportation and railway hub, as well as one of the primary energy and industrial bases in China. Xuzhou is also an area prone to geological disasters. According to the "Geological Disaster Prevention and Control Planning of Jiangsu Province (2016-2020)" (hereafter, abbreviated to Planning 2016-2020) [38], Xuzhou is one of the seven major geological disaster prevention and control areas (that are populated and liable to geological disasters such as landslide, collapse and land subsidence) in Jiangsu Province, as indicated in Figure 1b (blue lines).
The study area, outlined by the dashed black rectangle in Figure 1b, covers the southern part of the major geological disaster prevention and control area and the urban area (34.12 • -34.45 • N, 116.93 • -117.50 • E). The landform of the study area is characterized mainly by alluvial plains, piedmont and intermontane plains with elevation of generally 30-50 m, although there are also several low mountains and hills with elevation of 100-250 m. Rivers crisscross the city. The Beijing-Hangzhou Canal passes through the north, and the Old Yellow River that flows through the west and east connects the Beijing-Hangzhou Canal to Yunlong Lake (Figure 1c). Lakes and reservoirs are scattered all over the study area composing a water reticulation system together with the rivers. The study area is covered with loose quaternary sediments with depth of 0-50 m. The regional geological structures, such as the fault zones of the Old Yellow River, are beneficial to karst development. Karst caves and underground rivers are widely formed and they contain abundant underground water. Occurrences of karst collapse in the past were spatially correlated with these fault zones [38]. Meanwhile, coal deposits are well developed within this area. Coal mining activity and underground water exploitation in Xuzhou have been reported for decades and they are considered the main factors contributing to the current geological hazards [3,34,35]. For instance, Jiuli Lake and Panan Lake in the north (Figure 1c) have both been formed from transformed coal mining subsidence areas. In recent years, there has been a succession of mine closures in Xuzhou and the exploitation of underground water has been strictly limited. Consequently, the occurrence of related geological hazards has been reduced significantly. However, with continued urban expansion and the construction of the subway network traversing the main urban area, new urban areas and the Tongshan New District (Figure 1c), the factors contributing to surface displacement in Xuzhou have become increasingly complex.

Dataset
In this study, 52 C-band Sentinel-1A SLC SAR images from the ascending orbit acquired between 27 November 2015 and 8 June 2018 were used to perform the PSI analysis. The data were acquired in the TOPS mode. The track number was 142 and the resolutions in the slant range and azimuth were approximately 2 and 14 m, respectively. The central incidence angle over the study area was 40.12 • . Detailed information regarding the SAR dataset is presented in Table 1. A 30 m shuttle radar topography mission [39] digital elevation model (DEM) from the United States Geological Survey was adopted to simulate and subtract the topographic contribution, and to georeference the interferograms. The data recorded by three GNSS stations (TSGT, CUMT, and HCXY, indicated by red triangles in Figure 1, were used to validate the InSAR results and the vector data of the subway lines used to extract the subway construction areas, and were provided by the Xuzhou Land and Resources Bureau.

Methodology
As an advanced InSAR technique, the PSI algorithm can overcome the problem of decorrelation by identifying pixels whose echo is dominated by a single stable scatterer in a series of interferograms [40]. Existing PS methods [25,26,40,41] have been used successfully in analyses of urban areas. In this study, the PSI technique was expected to be capable of detecting sufficient coherent points because of the following: (1) The dataset contained a large number of images with dense temporal spacing, and (2) the study area mostly covered an urban area affected by slow continuous surface deformation, meaning the backscattered radar signals would not vary substantially. Therefore, the PSI method was adopted to perform the time series analysis.
The data processing methodology is illustrated in Figure 2. Up to the step of the single master interferogram formation, the processing was conducted using SNAP software [42] provided by the European Space Agency. The main steps of the PSI algorithm were implemented using the StaMPS package [36]. The traditional long-wavelength artifact correction method provided by StaMPS and the artifact correction method grounded on the theory of PC analysis (PCA), recently proposed by Chen et al. [6] (Strategy 1), were both found inadequate for high-precision displacement monitoring in this study (see following subsections). The main reason is that they do not account for the spatiotemporally uncorrelated components contained in the InSAR signal. Therefore, new correction strategies (Strategy 2 and Strategy 3) by further enhancing the method addressed in Chen et al. [6] were proposed in this study. These two strategies analyze, in detail, each PC derived from the PC decomposition. The methodology is detailed in the following.

Initial InSAR Time Series Generation
First, the Sentinel-1A precise orbit files [43] were applied to refine the orbit state vectors of the SAR acquisitions. The image acquired on 13 June 2017 was selected as the master image by taking into account the overall correlation coefficient; thus, all the slave images were coregistered and resampled to the master image. The maximum temporal baseline was 564 days and the perpendicular baseline was in the range 3-84 m ( Figure 3). Debursting was applied to each image before generation of the inteferogram. Multilooking was not conducted to retain the highest spatial resolution. The topographic contribution was removed and the interferograms were then georeferenced with the aid of the 30 m DEM. The coherent pixels (referred to as PS pixels) were then selected and purified based on statistical analysis of the amplitude difference dispersion and the phase stability [40]. To reduce the amount of data for further processing, all single master interferograms were subsampled to 20 m resolution. Consequently, 430,651 PS pixels were ultimately retained. Then, the phases of the PS pixels were unwrapped using a statistical cost flow approach named the 3D unwrapping method [40,44]. The initial InSAR time series data were obtained after estimating and removing the baseline-dependent phase residuals in the interferograms induced by errors in the DEM, as well as the atmosphere and orbit error (AOE) phase of the master image present in every interferogram ( Figure 2).
Chen et al. [6] (Strategy 1), were both found inadequate for high-precision displacement monitoring in this study (see following subsections). The main reason is that they do not account for the spatiotemporally uncorrelated components contained in the InSAR signal. Therefore, new correction strategies (Strategy 2 and Strategy 3) by further enhancing the method addressed in Chen et al. [6] were proposed in this study. These two strategies analyze, in detail, each PC derived from the PC decomposition. The methodology is detailed in the following.  First, the Sentinel-1A precise orbit files [43] were applied to refine the orbit state vectors of the SAR acquisitions. The image acquired on 13 June 2017 was selected as the master image by taking into account the overall correlation coefficient; thus, all the slave images were coregistered and resampled to the master image. The maximum temporal baseline was 564 days and the perpendicular baseline was in the range 3-84 m ( Figure 3). Debursting was applied to each image before generation of the inteferogram. Multilooking was not conducted to retain the highest spatial resolution. The topographic contribution was removed and the interferograms were then georeferenced with the aid of the 30 m DEM. The coherent pixels (referred to as PS pixels) were then selected and purified based on statistical analysis of the amplitude difference dispersion and the phase stability [40]. To reduce the amount of data for further processing, all single master interferograms were subsampled to 20 m resolution. Consequently, 430,651 PS pixels were ultimately retained. Then, the phases of the PS pixels were unwrapped using a statistical cost flow approach named the 3D unwrapping method [40,44]. The initial InSAR time series data were obtained after estimating and removing the baselinedependent phase residuals in the interferograms induced by errors in the DEM, as well as the atmosphere and orbit error (AOE) phase of the master image present in every interferogram ( Figure  2). To further mitigate the impact of the remaining uninteresting contributions in the InSAR signal, we used an empirical quadratic function (Equation (1), Chen et al. [6]) to approximate the phase ramps induced by long-wavelength artifacts (possible residual orbital errors and stratified tropospheric delays):

Initial InSAR Time Series Generation
where for pixel in an interferogram, is the interferometric phase of the initial time series, and represent the pixel coordinates, is the elevation from the DEM, and , , … , ℎ are fitting parameters.
As shown in Table 2, the significant variance reduction (the root mean square (RMS) values obtained from 15.17 mm to 11.02 mm) computed in areas where no displacement was expected indicates the effectiveness of the method. Nevertheless, an uncertainty of ~11 mm is nonnegligible with respect to high-accuracy subtle displacement InSAR monitoring. It was thus necessary to further To further mitigate the impact of the remaining uninteresting contributions in the InSAR signal, we used an empirical quadratic function (Equation (1), Chen et al. [6]) to approximate the phase ramps induced by long-wavelength artifacts (possible residual orbital errors and stratified tropospheric delays): where for pixel i in an interferogram, φ i is the interferometric phase of the initial time series, x i and y i represent the pixel coordinates, z i is the elevation from the DEM, and a, b, . . . , h are fitting parameters. As shown in Table 2, the significant variance reduction (the root mean square (RMS) values obtained from 15.17 mm to 11.02 mm) computed in areas where no displacement was expected indicates the effectiveness of the method. Nevertheless, an uncertainty of~11 mm is nonnegligible with respect to high-accuracy subtle displacement InSAR monitoring. It was thus necessary to further improve the signal-to-noise ratio. The high-level phase variance in non-displacement areas (selected based on a priori knowledge) could be attributable to the following reasons. The first is that the long-wavelength artifacts were not properly approximated because of the contamination of the long-wavelength deformation signal. A similar phenomenon was found in a previous study by Chen et al. [6] over volcanic areas, where the surface displacement was correlated with elevation. Although the topography varies slightly over our study area in Xuzhou, and the effect from stratified tropospheric delays was expected to be less significant, the estimation of long-wavelength artifacts attributable to orbital errors might be biased by the long-wavelength displacement. The second reason is that turbulent tropospheric delays and/or random noise phases could play an important role. Many previous studies [31,45] used bandpass filtering to remove turbulent tropospheric delays and random noise phases that are spatiotemporally uncorrelated. However, great caution must be exercised when dealing with this type of filtering because the deformation signal is vulnerable when the parameters are not set properly. Consequently, to avoid losing any interesting signals, we preferred to explore the InSAR signal more deeply instead of applying any temporal or spatial filtering. a Correction Strategy 1 refers to the approach addressed in Chen et al. [6]. b Correction Strategy 2 refers to the enhanced approach that subtracts the PC3 and the PC4 along with the removal of the long-term displacement model (LTDM) before long-wavelength artifact correction. c Correction Strategy 3 refers to the enhanced approach that subtracts the PC3, PC4 and PC5 along with the removal of the LTDM before long-wavelength artifact correction.

InSAR Signal Analysis Based on PC Decomposition
PCA is an excellent technique for extracting information contained in data and concentrating it within a few PCs. Therefore, to further explore the InSAR signal, the PC decomposition implemented in the PCA-based inversion method package [37], was applied to the initial time series, as obtained in Section 3.2.1. The input data matrix comprised 430,651 lines (total pixel number) and 52 rows (epochs). It was decomposed into a linear combination of PCs, each of which was individually associated with a spatial function (U), significance value (S) and temporal function (V) (Equation (2)). The spatial function presents the spatial distribution of displacements, the significance value indicates the importance of the relative component, while the temporal function describes the temporal behaviors of the displacements (see Kositsky and Avouac [37] for more detailed theories, and see Lin et al. [46], Remy et al. [47] and Chen et al. [6] for examples of application).
where m is the number of pixels, n is the number of epochs, k is the number of PCs, X m×n denotes the centered data matrix, U m×k and V k×n represent the spatial function and temporal function of the components, respectively, S k×k is a diagonal matrix whose values signify the significance of each component and V T k×n denotes the transpose of V k×n . Figure 4 summarizes the temporal functions and the spatial functions weighted by the significance values of the PCs obtained from the InSAR time series. The PC1 ( Figure 4a) exhibits a strong spatially correlated long-wavelength signal that is uncorrelated in time. This long-wavelength signal could be mostly related to residual orbital errors and partially affected by topography-correlated tropospheric artifacts, according to the spatial function-topography scatter plot shown in the inset. The dominant proportion (P S1 = 0.40) of PC1 suggests the residual orbital errors could be significant in some places, even for modern satellites with precise orbits (e.g., Sentinel-1). The temporal function of PC2 (Figure 4b) behaves as a linear trend indicating high temporal correlation of the signal, which suggests that PC2 could be related to long-term displacement. To address this consideration, investigation of the initial mean velocity and the initial time series of several points was carried out for comparison with the PC2 (see Section 1 in the Supplementary Materials). Results showed the spatial pattern of PC2 is very similar to that of the initial mean velocity and most of the points reveal a similar linear displacement trend. Therefore, the signal in PC2 is related mostly to the land surface displacement affecting the study area. Conversely, both PC3 and PC4 (Figure 4c,d) present signals that are mainly uncorrelated spatiotemporally over the study area. The amplitudes in the temporal functions fluctuate randomly around zero, while the spatial pattern shows high variability with no evident correlation with topography, coinciding well with the characteristics of the turbulent component of the tropospheric artifact. The temporal behavior of PC6 (Figure 4f) shows clear seasonal variation with mainly positive values in summer and negative values in winter. The scatter plot reveals clear correlation between the spatial function values and topography from~70 to~200 m, evidencing that part of the signal presented in PC6 is attributable to the seasonality of the long-wavelength stratified component of the tropospheric artifact, although the topography (range: 16-230 m for PS pixels) varies slightly within the study area. Based on visual inspection, signals below 70 m could correspond to seasonal variation in the displacement of some local constructions related to thermal expansion and construction activity. PC5 ( Figure 4e) presents a similar spatial pattern to both PC3 and PC4 without clear correlation with topography. However, the temporal function appears to show slight seasonal behavior. This makes it difficult to classify it in any known category. It might be the combination of the turbulent tropospheric artifact, and/or seasonal displacement and/or noise. Therefore, special caution should be observed when handling this component. Detailed analysis of the initial time series based on PC decomposition helped elucidate each component of the InSAR signals; hence, the correction strategy (addressed in the next section) in terms of accurate measurement of subtle long-term displacements was proposed based on the above analysis.

Purification of Time Series Displacement
In the artifact correction approach proposed by Chen et al. [6], the displacement-related component was used to estimate the long-term displacement model (LTDM) that was removed prior to artifact estimation. Then a quadratic function (Equation (1)) was adopted to approximate the phase delays induced by the long-wavelength artifacts from the residual signals. Finally, the LTDM was

Purification of Time Series Displacement
In the artifact correction approach proposed by Chen et al. [6], the displacement-related component was used to estimate the long-term displacement model (LTDM) that was removed prior to artifact estimation. Then a quadratic function (Equation (1)) was adopted to approximate the phase delays induced by the long-wavelength artifacts from the residual signals. Finally, the LTDM was added back. Results showed that this approach was capable of obtaining high-precision displacement in their application to a volcanic area. As evidenced, PC2 is related mainly to surface displacement over the Xuzhou area; therefore, we used a linear function (dashed red line in Figure 4b) to fit the temporal function of PC2. The LTDM was estimated by reconstructing the data using U2, S2 and the best-fit linear model and then it was removed from the initial time series. Nevertheless, after having applied the same correction approach (i.e., Chen et al. [6] referred to as Correction Strategy 1 in this study), little improvement was obtained regarding the variance in non-displacement zones when compared with the result from the long-wavelength artifact correction before PC decomposition (RMS: 9.89 versus 11.02 mm; Table 2). This was largely because the long-wavelength stratified tropospheric delay was the dominant artifact in the case of Chen et al. [6], while this study found that the turbulent tropospheric delay could play an important role in the observed InSAR signal (accounting for about one quarter of the observed signal; Figure 4c,d), as analyzed in Section 3.2.3. Therefore, we added an enhancement to this approach (Figure 2). In addition to the removal of the LTDM, the spatiotemporally uncorrelated components PC3 and PC4, which are considered mostly related to turbulent tropospheric delays, were also subtracted from the initial time series. However, for the relatively critical component PC5, seasonal displacement information might have been lost if it were removed directly. Therefore, we decided to derive two independent sets of residual time series for further correction: One with PC5 (referred to as Correction Strategy 2) and one without PC5 (referred to as Correction Strategy 3). Then, Equation (1) was applied to estimate the long-wavelength artifacts that remained in each map in the two sets of residual time series. The corrected time series were finally obtained after the LTDM was added back.
As shown in Table 2, Correction Strategies 2 and 3 reduced the variance in non-displacement areas by approximately half with respect to Correction Strategy 1. To further evaluate the effectiveness of the different correction strategies, the InSAR time series observed the pixels where the three GNSS stations (TSGT, CUMT and HCXY; see Figure 1 for locations) are installed, and were compared with the measurements recorded at the GNSS stations ( Figure 5), and the RMSs of the differences between the two datasets were calculated ( Table 2). The GNSS receivers record one measurement every 30 s. The daily GNSS solutions were processed using Bernese GNSS software. The average 1σ accuracies of the daily positions were~3 and~6 mm for the horizontal and vertical components, respectively. The precise ground displacement time series data in three dimensions (up-down, U-D; east-west, E-W; north-south, N-S) at each GNSS station were then obtained after having been corrected for the plate motion using the formula in Yang et al. [48]. The three-component positional time series from the GNSS stations were projected into the radar line-of-sight (LOS) direction of the satellite before comparison with the InSAR data. Table 2 shows the artifact correction methods based on PC decomposition are more effective than the conventional long-wavelength artifact correction method that does not account for the possible impact induced by the displacement signal. In comparison with Correction Strategies 2 and 3, Strategy 1 (Chen et al. [6]), consistent with the case in non-displacement zones, makes less significant contributions to reducing the discrepancies between GNSS and InSAR observations, which is mainly because the Correction Strategy 1 is unable to handle the spatiotemporally uncorrelated InSAR signals (i.e., mostly turbulent tropospheric delays). Correction Strategies 2 and 3 both show almost equally matched performance in purifying the InSAR time series, which can be explained by the reasonable agreement between the GNSS and the corrected InSAR measurements (Table 2 and Figure 5). Although the InSAR time series corrected by Strategy 3 appear smoother in non-displacement zones when compared with those corrected by Strategy 2, it does not perform better than Strategy 2 in reducing the RMSs of the differences between the GNSS and InSAR measurements ( Table 2). This is especially evident at TSGT where the GNSS time series data show seasonal variations ( Figure 5). The InSAR time series derived by Correction Strategy 2 display consistent seasonal variations with GNSS data. However, the InSAR time series derived using Strategy 3 appear over-corrected because the seasonal variation in displacement tends to be much less apparent ( Figure 5). The analyses above confirm the assumption addressed in Section 3.2.2, that the seasonal displacement signal might be one of the compositions of PC5. Therefore, the InSAR time series data corrected by Strategy 2 are used to expand the following sections of this paper. PC5. Therefore, the InSAR time series data corrected by Strategy 2 are used to expand the following sections of this paper.

Displacement Rate Analysis
The displacement rate map of the study area derived from the corrected InSAR time series is illustrated in Figure 6. Results showed that a large part of the study area is affected by heterogeneous ground displacement in the range of approximately −41 to +25 mm/yr (in the LOS direction) during the entire study period (November 2015 to June 2018). The western part of the study area suffers surface displacement more severely than the east. The major displacement areas include the buffer zones along the subway lines and Zones 1-4. Zone 1 is located to the south of Yunlong Lake and to the west of the Tongshan New District, which is the area of most recent urban expansion (Figures 1  and 6). Zones 2-4 are coincidently located in three old goaf areas, the western goaf (WG), eastern goaf (EG) and Jiuli Hill-Ma village goaf (J-MG), respectively [49]. Of those main displacement areas, the subway buffer and Zones 1-3 are mostly characterized by surface subsidence. The most significant subsidence rates are distributed to the south of Jiuli Lake (WG) and to the east of Panan Lake (EG), i.e., maximum rates of up to −40.5 and −37.5 mm/yr, respectively. Remarkable subsidence rates up to −28.4 mm/yr are also found near the three subway stations along Line 1 to the north of Yunlong Lake (indicated by black-edged white-filled circles in Figure 6). Relatively, less significant subsidence rates ranging from a few to a dozen centimeters are found within Zone 1 and the remainder of the subway buffer is mainly concentrated within the main urban area. Conversely, unexceptional land uplift signals with deformation rates in the range 4-25 mm/yr are observed in Zone 4 and along the

Displacement Rate Analysis
The displacement rate map of the study area derived from the corrected InSAR time series is illustrated in Figure 6. Results showed that a large part of the study area is affected by heterogeneous ground displacement in the range of approximately −41 to +25 mm/yr (in the LOS direction) during the entire study period (November 2015 to June 2018). The western part of the study area suffers surface displacement more severely than the east. The major displacement areas include the buffer zones along the subway lines and Zones 1-4. Zone 1 is located to the south of Yunlong Lake and to the west of the Tongshan New District, which is the area of most recent urban expansion (Figures 1  and 6). Zones 2-4 are coincidently located in three old goaf areas, the western goaf (WG), eastern goaf (EG) and Jiuli Hill-Ma village goaf (J-MG), respectively [49]. Of those main displacement areas, the subway buffer and Zones 1-3 are mostly characterized by surface subsidence. The most significant subsidence rates are distributed to the south of Jiuli Lake (WG) and to the east of Panan Lake (EG), i.e., maximum rates of up to −40.5 and −37.5 mm/yr, respectively. Remarkable subsidence rates up to −28.4 mm/yr are also found near the three subway stations along Line 1 to the north of Yunlong Lake (indicated by black-edged white-filled circles in Figure 6). Relatively, less significant subsidence rates ranging from a few to a dozen centimeters are found within Zone 1 and the remainder of the subway buffer is mainly concentrated within the main urban area. Conversely, unexceptional land uplift signals with deformation rates in the range 4-25 mm/yr are observed in Zone 4 and along the southern and eastern border of Zone 2. Further analysis and interpretation of the surface displacement, including spatiotemporal behavior and possible origins, are presented in Section 5.

Reliability Analysis
As analyzed in Section 3.2.3, the RMS of the differences between the InSAR-derived results and the GNSS time series are in the range of approximately 3-6 mm, showing a reasonable agreement between the two datasets that quantitatively confirms the reliability of our InSAR results (Table 2 and Figure 5). The standard deviations of the displacement rates are analyzed statistically in this section, allowing us to assess the internal precision of the InSAR results for every single PS pixel. The distributions of the standard deviations of the displacement rates computed based on Bootstrap statistics from two datasets (i.e., the initial time series and corrected time series obtained from Strategy 2) are shown in Figure 7. The significant reduction in the standard deviations of the displacement rates (median from ~3 to ~0.8 mm/yr) demonstrates again, that the improved artifact correction method has apparently succeeded in purifying the InSAR displacement signal. Almost the entire population of corrected data falls into the interval 0-4 mm/yr, and about 86% of the PS pixels have a standard deviation of displacement rate of <1 mm/yr (Figure 7). Consequently, both external GNSS data and internal precision checking evidence the high reliability and precision of our InSAR results.

Reliability Analysis
As analyzed in Section 3.2.3, the RMS of the differences between the InSAR-derived results and the GNSS time series are in the range of approximately 3-6 mm, showing a reasonable agreement between the two datasets that quantitatively confirms the reliability of our InSAR results (Table 2 and Figure 5). The standard deviations of the displacement rates are analyzed statistically in this section, allowing us to assess the internal precision of the InSAR results for every single PS pixel. The distributions of the standard deviations of the displacement rates computed based on Bootstrap statistics from two datasets (i.e., the initial time series and corrected time series obtained from Strategy 2) are shown in Figure 7. The significant reduction in the standard deviations of the displacement rates (median from~3 to~0.8 mm/yr) demonstrates again, that the improved artifact correction method has apparently succeeded in purifying the InSAR displacement signal. Almost the entire population of corrected data falls into the interval 0-4 mm/yr, and about 86% of the PS pixels have a standard deviation of displacement rate of <1 mm/yr (Figure 7). Consequently, both external GNSS data and internal precision checking evidence the high reliability and precision of our InSAR results.

Surface Subsidence Associated with Urban Construction
Building construction during urban expansion and tunnel excavation during subway construction are the main causes of land subsidence in cities, as reported in previous research [14,33,50]. The results of this study show that the urban area of Xuzhou has suffered noticeable surface subsidence that could be attributed to both urban expansion and subway construction ( Figure  6). As mentioned in Section 4.1., Zone 1 is located in the western part of the Tongshan New District that has been experiencing rapid development during recent years. Many urban construction activities including the construction of commercial centers, high-rise residential buildings, hospitals and other urban infrastructure are concentrated within this area, resulting in ground subsidence over many years. The magnitude of the displacement rates in Zone 1 are mainly below -10 mm/yr. Thus, we selected only the buffer of the subway lines in addressing a detailed analysis on the relationship between surface subsidence and urban construction (Figures 8 and 9).
Construction of the three subway lines (Lines 1-3) in Xuzhou commenced in December 2015, February 2016 and August 2016, respectively, and remained ongoing at the end of our study period. Although the SAR images are incapable of acquiring information regarding the inner structures and tunnels of the subway lines, it is possible to capture the overlying surface subsidence along the paths of the subway tunnels using InSAR [12,51,52]. The displacement rates within the 1 km buffer zone of the subway lines between November 2015 and June 2018, as derived by PSI analysis, are shown in Figure 8. They indicate a spatially heterogeneous trend of subsidence development concentrated in the main urban area, mainly between Xingshanzi station and the Xuzhou railway station along subway Line 1, and from the Pengcheng Square station to Keyunbei station along Line 2. Except when traversing the main urban area, most stations of subway Line 3 were stable during the study period. This uneven behavior of displacement could be attributable to the fact that the main urban area possesses the majority of dense high-rise buildings, and that the construction of the subway lines most probably adds further pressure to the already fragile overlying strata, resulting in land surface displacement. The area of most significant displacement (approximately −10 to −30 mm/yr) was observed around the People's Square, Gongnong Road and Hanshan stations that are on the western limb of Line 1. The inset in Figure 8a shows an enlargement of the displacement rate of the area outlined by the white rectangle. To obtain a clearer idea of the displacement pattern, we constructed a profile line (P-P') shown in Figure 8b. The profile shows that the pattern of subsidence in the transverse direction to the subway line is >1 km wide and that the transverse subsidence

Surface Subsidence Associated with Urban Construction
Building construction during urban expansion and tunnel excavation during subway construction are the main causes of land subsidence in cities, as reported in previous research [14,33,50]. The results of this study show that the urban area of Xuzhou has suffered noticeable surface subsidence that could be attributed to both urban expansion and subway construction ( Figure 6). As mentioned in Section 4.1., Zone 1 is located in the western part of the Tongshan New District that has been experiencing rapid development during recent years. Many urban construction activities including the construction of commercial centers, high-rise residential buildings, hospitals and other urban infrastructure are concentrated within this area, resulting in ground subsidence over many years. The magnitude of the displacement rates in Zone 1 are mainly below -10 mm/yr. Thus, we selected only the buffer of the subway lines in addressing a detailed analysis on the relationship between surface subsidence and urban construction (Figures 8 and 9).
Construction of the three subway lines (Lines 1-3) in Xuzhou commenced in December 2015, February 2016 and August 2016, respectively, and remained ongoing at the end of our study period. Although the SAR images are incapable of acquiring information regarding the inner structures and tunnels of the subway lines, it is possible to capture the overlying surface subsidence along the paths of the subway tunnels using InSAR [12,51,52]. The displacement rates within the 1 km buffer zone of the subway lines between November 2015 and June 2018, as derived by PSI analysis, are shown in Figure 8. They indicate a spatially heterogeneous trend of subsidence development concentrated in the main urban area, mainly between Xingshanzi station and the Xuzhou railway station along subway Line 1, and from the Pengcheng Square station to Keyunbei station along Line 2. Except when traversing the main urban area, most stations of subway Line 3 were stable during the study period. This uneven behavior of displacement could be attributable to the fact that the main urban area possesses the majority of dense high-rise buildings, and that the construction of the subway lines most probably adds further pressure to the already fragile overlying strata, resulting in land surface displacement. The area of most significant displacement (approximately −10 to −30 mm/yr) was observed around the People's Square, Gongnong Road and Hanshan stations that are on the western limb of Line 1.
The inset in Figure 8a shows an enlargement of the displacement rate of the area outlined by the white rectangle. To obtain a clearer idea of the displacement pattern, we constructed a profile line (P-P') shown in Figure 8b. The profile shows that the pattern of subsidence in the transverse direction to the subway line is >1 km wide and that the transverse subsidence approximates a Gaussian distribution. The magnitude of subsidence decreases gradually with distance from the center of the subway tunnel, indicating that the pattern of subsidence is strongly associated with subway construction. Fieldwork results revealed that the pattern of subsidence has had marked effects on surrounding buildings (Figure 8c). The peak average subsidence rate of this pattern (−28.4 mm/yr) was found at point A near the People's Square station (Figure 8a). The displacement time series for points A-D located near Xuyifuyuan station in the urban area, Keyunbei station on the northern limb of Line 2, and Xinyuan Avenue station in the new urban area are illustrated in Figure 9. It can be seen that the temporal evolutions of subsidence at points A−D can be described approximately by linear trends with average subsidence rates of −28.4, −10.8, −7.1 and −17.6 mm/yr, respectively. However, two different stages of subsidence can be recognized at points B−D: (1) The ground surface was most likely stable until around mid-2016, and (2) rapid subsidence has occurred since then. At the end of the study period, the subsidence rates at the four points showed no sign of decreasing, indicating that subsidence will likely continue over the following years. The displacement time series at the four points appear to coincide with the time of subway construction work, suggesting the subsidence at these points is related to the construction work. approximates a Gaussian distribution. The magnitude of subsidence decreases gradually with distance from the center of the subway tunnel, indicating that the pattern of subsidence is strongly associated with subway construction. Fieldwork results revealed that the pattern of subsidence has had marked effects on surrounding buildings (Figure 8c). The peak average subsidence rate of this pattern (−28.4 mm/yr) was found at point A near the People's Square station (Figure 8a). The displacement time series for points A-D located near Xuyifuyuan station in the urban area, Keyunbei station on the northern limb of Line 2, and Xinyuan Avenue station in the new urban area are illustrated in Figure 9. It can be seen that the temporal evolutions of subsidence at points A−D can be described approximately by linear trends with average subsidence rates of −28.4, −10.8, −7.1 and −17.6 mm/yr, respectively. However, two different stages of subsidence can be recognized at points B−D: (1) The ground surface was most likely stable until around mid-2016, and (2) rapid subsidence has occurred since then. At the end of the study period, the subsidence rates at the four points showed no sign of decreasing, indicating that subsidence will likely continue over the following years. The displacement time series at the four points appear to coincide with the time of subway construction work, suggesting the subsidence at these points is related to the construction work.

Surface Subsidence and Uplift in Old Goafs
Xuzhou has a coal mining history extending over 130 years that has resulted in many areas of large subsidence. By 2014, the area of subsidence in the Xuzhou urban planning area related to coal mining activities was approximately 180 km 2 [49]. Most coal mines in the Xuzhou area have closed

Surface Subsidence and Uplift in Old Goafs
Xuzhou has a coal mining history extending over 130 years that has resulted in many areas of large subsidence. By 2014, the area of subsidence in the Xuzhou urban planning area related to coal mining activities was approximately 180 km 2 [49]. Most coal mines in the Xuzhou area have closed successively in recent years following the structural de-capacity policy of the Chinese coal industry (e.g., Quantai mine closed in 2010; Jiahe, Pangzhuang and Zhangxiaoloujing mines closed in 2015; Qishan mine closed in 2016; Shitun mine closed in 2017; see Figure 6 for locations). However, because of the overburden of strata, persistent subsidence of between approximately −5 and −40 mm/yr remains evident in the few years following closure. As one of the first mines to close, the area around the Quantai mine still experienced land surface subsidence in the range of −5.2 to −15.5 mm/yr during the five to eight years after closure. Approximately 85% of this area is still affected by subsidence rates of under −10 mm/yr, suggesting that the goaf of the Quantai mine might tend towards stability in the near future. The most significant subsidence rates of up to −19.0, −24.3 and −32.0 mm/yr were found for Jiahe, Pangzhuang and Zhangxiaoloujing mines, respectively. Enlarged maps of the displacement rates in the Shitun (WG) and Qishan (EG) mining areas (outlined using white rectangles in Figure 6) are illustrated in Figure 10a,b, respectively. The villages of Shixi and Cuizhuang were officially designated hidden danger sites (at extra high and high levels, respectively) of geological disaster attributed to ground subsidence, according to the Geological Disaster Verification Report of Xuzhou [49]. Similarly, they are also within the areas found most affected by land subsidence related to the Shitun and Qishan mines in our study (Figure 10a,b). The subsidence rate in Shixi increases from northwest (−8.1 mm/yr) to southeast (−40.5 mm/yr), while that in Cuizhuang increases from northeast (−13.4 mm/yr) to southwest (−37.5 mm/r). The time series for points E and F in Figure 10a,b are shown in Figure 10c,d. It can be seen that the subsidence in both Shixi and Cuizhuang follow a linear trend of temporal evolution. Severe fissures have been caused in local houses by land subsidence (see examples in Figure 10e,f). More than 1500 people in the area are threatened by the land subsidence hazard and the potential cost of the damage is 30 million RMB [49], which highlights the importance of monitoring displacement in goaf areas.
In addition to the surface subsidence that affected the old goafs in Xuzhou during 2015-2018, uplift of up to +25 mm/yr was also found, e.g., along the southern and eastern borders of Zone 2 and in Zone 4. However, it is worth mentioning that these two areas were affected by subsidence and not land uplift before 2015, according to the Geological Disaster Verification Report of Xuzhou [49], which makes 2015 a critical year. Xuzhou is the only city in Jiangsu Province that uses groundwater as its water supply. From the 1980s to the 1990s, the regional karst water level declined annually because of long-term overexploitation of karst groundwater. Regional water level drop funnels were formed in the centralized area of exploitation, which resulted in both depletion of water resources and induced geological disasters such as karst collapses. Since the beginning of the twenty-first century, following implementation of a stringent system for groundwater resource management, the regional underground water level has become controlled and the drop funnels in the urban area have gradually diminished and disappeared. Moreover, few karst collapses or land subsidence have been detected in recent years. In 2015, to further control groundwater exploitation, the "Program for Compressing Groundwater Exploitation in Jiangsu Province (2014-2020)" [53] was officially approved by the provincial government. This led to 33 karst water mining wells being permanently filled, and 66 standby karst water mining wells sealed within the Xuzhou urban planning area. We considered that the uplift phenomenon in the two old goafs since 2015 could be related to the groundwater level rise attributable to the new exploitation regime. These two uplift areas are located within old coal mining areas [54] that have experienced surface subsidence for years. Moreover, the two long narrow areas of uplift are located exactly where two northeast-southwest-oriented synclines lie (see Figures 1  and 6). Therefore, these two areas have strong affinity for water convergence when the phreatic water level rises, which can induce uplift [55,56]. Furthermore, as these areas are situated on an alluvial plain, the deposits are characterized by an unconsolidated layer of sandy clay, clay, or silty sand at a typical depth of 2-23 m below the ground surface [49]. The unconsolidated layer would inflate when experiencing rising levels of underground water, which would result in uplift. When the rate of uplift is relatively higher than the rate of subsidence induced by the persistent overburden of the strata, the uplift signal will be observed by InSAR. Equal attention should be paid to land uplift as to subsidence because both phenomena can cause geological hazards and threaten human life and property. In addition to the surface subsidence that affected the old goafs in Xuzhou during 2015-2018, uplift of up to +25 mm/yr was also found, e.g., along the southern and eastern borders of Zone 2 and in Zone 4. However, it is worth mentioning that these two areas were affected by subsidence and not land uplift before 2015, according to the Geological Disaster Verification Report of Xuzhou [49], which makes 2015 a critical year. Xuzhou is the only city in Jiangsu Province that uses groundwater as its water supply. From the 1980s to the 1990s, the regional karst water level declined annually because of long-term overexploitation of karst groundwater. Regional water level drop funnels were formed in the centralized area of exploitation, which resulted in both depletion of water resources and induced geological disasters such as karst collapses. Since the beginning of the twenty-first century, following implementation of a stringent system for groundwater resource management, the regional underground water level has become controlled and the drop funnels in the urban area have  Figure 6). The village of Shixi is outlined by dashed black lines. The map conventions are the same as in Figure 4. (b) Enlarged map of displacement rate (2015-2018) around Qishan mine (the eastern white rectangle in Figure 6). The village of Cuizhuang is outlined by dashed black lines. The map conventions are the same as in Figure 4. (c,d) Displacement time series and displacement rates of points E and F marked by black crosses in (a,b), respectively. (e,f) Damages to local houses caused by land subsidence.

Potential of the Proposed Method for High Precision InSAR Observation
This study proposed an enhanced artifact correction method based on PCA theory for high-precision InSAR time series analysis. Results showed this method could be used for comprehensive investigation of the compositions of an InSAR signal, and that it could largely improve the observation accuracy. Thus, the method was demonstrated to have considerable potential in relation to the acquisition of high-precision InSAR time series observations. Nevertheless, the proposed method does have several deficiencies. The first is that the number of PCs plays an important role in the PCA-based PSI analysis. However, there is no absolute optimal number because of the lack of physical meaning of each PC. To assist time series analysis, the PCs are abstractly linked to certain physical meanings based on their spatiotemporal characteristics. Therefore, the number of PCs is determined largely on experience as well as repetitive testing. The order of magnitude could vary depending on the research purpose. For instance, some researchers [47] have used PCA to invert displacement components and to derive information for further modeling of displacement sources. In such work, the data matrix was decomposed into dozens of PCs. First, several PCs were extracted to reconstruct a new dataset and to replace the true InSAR measurements. However, the primary purpose of this study, as well as that of Chen et al. [6] (8 PCs), was to seek a suitable PC with which to estimate the LTDM, not to filter out any uncertain components. Too many PCs tend to disperse the displacement signals, especially when the magnitudes of such signals are equivalent to or lower than the artifacts. Fewer PCs tend to mix different signals, making it difficult to obtain a PC that is spatiotemporally related to the expected ground displacement. These conclusions were drawn from repetitive testing, and then the ideal number of PCs in this study was set as six. The means of determination of the number of PCs does indeed have a certain degree of subjectivity, but it does have certain sophistication and it has been demonstrated to work well both in previous research and in this study. The second deficiency is that the method would work better when the entire study area shares similar temporal behavior of displacement (e.g., the linear trend found in this study and the exponential trend in Chen et al. [6]). The third deficiency is that the direct removal of one of the PCs could introduce the risk of losing interesting information. Therefore, external data (e.g., GNSS or leveling) are required to assist in verifying the performance of the correction strategies. The above shortcomings mean that the proposed method is not a "ready-made" approach suitable for all situations. Introducing PCA to InSAR time series analysis is an innovative idea but its functionality will vary case by case.