Open AccessArticle

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

Yu Chen

¹,

Kun Tan

^1,2,*,

Shiyong Yan

¹,

Kefei Zhang

^1,3

Hairong Zhang

¹,

Xiaoyang Liu

¹,

Huaizhan Li

¹ and

Yaqin Sun

Key Laboratory for Land Environment and Disaster Monitoring (Ministry of Natural Resources), China University of Mining and Technology, Xuzhou 221116, China

Key Laboratory of Geographic Information Science (Ministry of Education), East China Normal University, Shanghai 200241, China

SPACE Research Center, School of Science, RMIT University, Melbourne VIC 3001, Australia

Author to whom correspondence should be addressed.

Remote Sens. 2019, 11(12), 1494; https://doi.org/10.3390/rs11121494

Submission received: 5 May 2019 / Revised: 20 June 2019 / Accepted: 21 June 2019 / Published: 24 June 2019

(This article belongs to the Special Issue Applications of Sentinel Satellite for Geohazards Prevention)

Download

Browse Figures

Versions Notes

Abstract

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 significant 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 significantly affected two long narrow areas within the old goafs since 2015. The possible principal influencing factors contributing to the land surface displacements such as subway tunneling, building construction, mining, underground water levels and geological conditions are then discussed.

Keywords:

persistent scatterer interferometry; principal component decomposition; land surface displacement; Xuzhou; urbanization

Graphical Abstract

1. 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 one of the most prominent geological issues globally [1,2,3,4,5,6]. Land surface deformation threatens human lives and properties, particularly in metropolitan areas with a large population, because it can often result in severe damage both to buildings and to transport infrastructures.

China is among those countries that have experienced serious geological disasters because of its complex geological environment, huge population and rapidly growing economy. According to incomplete statistics, more than 96 cities in China have suffered varying degrees of land subsidence, and about 80% of them are located in the eastern and central parts of the country [7,8]. Xuzhou, as one of the important energy and industrial cities in eastern China, located in the northwest of Jiangsu Province, has undergone land surface deformation since 1949 driven predominately by both excessive extraction of karst underground water and exploitation of mineral resources. Up to 2001, the various types of geological hazards in Xuzhou had led to direct economic losses of approximately 2.04 billion RMB [9,10]. Xuzhou is one of the core cities of the Huaihai Economic Zone (Figure 1), and it was officially designated the central city in 2017 owing to the fast pace of its economic development. However, accompanying its large scale urbanization, including urban expansion, new underground construction, population migration, as well as the reclamation and exploitation of mining subsidence zones, geological hazards and regional stability problems have become increasingly serious [10]. It is worth emphasizing that Xuzhou has recently undertaken a program of subway network construction. By August 2016, the first-stage projects of the first three lines (Line 1 to Line 3) had been fully commenced, with expected completion dates in September 2019, October 2020, and December 2020, respectively. Recent land subsidence or collapse accidents caused by subway construction have been documented widely in China, including megacities such as Beijing [11], Shanghai [12,13] and Guangzhou [14,15], and many other cities. Consequently, the monitoring of land surface deformation in Xuzhou on a regular basis is an urgent requirement to prevent any associated risks and to support the integrated sustainable development of the city.

Ground-based deformation monitoring techniques such as leveling, global navigation satellite system (GNSS) measurements, and geological and geophysical investigation methods have traditionally been used to monitor land surface displacement in Xuzhou [16,17]. However, all these techniques are based on point-by-point measurements that are generally sparse in nature and involve time-consuming and labor-intensive collection, hence, they cannot meet the requirements for regional stability assessments. Compared with traditional ground-based methods, the spaceborne interferometric synthetic aperture radar (InSAR) technique can provide unparalleled resolution and large coverage of land surface displacements at a low cost [18,19]. The capability and marked advantages of the InSAR technique have been demonstrated in various applications, e.g., assessment of volcanic hazards [20], earthquakes [21], landslides [22] and land subsidence in metropolitan areas [23]. Furthermore, the advanced time series InSAR (TS-InSAR) technique, including persistent scatterer interferometry (PSI) [24,25,26], small baseline subset InSAR [27] and combined multi-temporal InSAR [28], is ideal for investigating both the spatial distribution and the temporal evolution of land surface displacement. The TS-InSAR is capable of providing displacement maps with centimeter-to-millimeter accuracy over large-scale areas and at high spatiotemporal resolution. It has been applied successfully to the study of land surface deformation in various cities such as Mexico City [5], Rome [29], Istanbul [30], Beijing [31], Shanghai [13], and Wuhan [32,33]. However, the few InSAR-based investigations [3,34,35] conducted previously in Xuzhou, focused on detecting surface subsidence in individual mining areas. Consequently, the long-term large-scale surface deformation in central parts of Xuzhou and along the subway lines under construction remains unexplored.

The main objective of this study was to characterize the spatial pattern and temporal evolution of land surface deformation in central areas of Xuzhou. The PSI algorithm, based on the Stanford method for persistent scatterers (StaMPS) package [36], was adopted to perform an initial detailed time series analysis using 52 free C-band Sentinel-1A terrain observation with progressive scans (TOPS) images acquired during 2015–2018. Then, the artifact correction approach proposed by Chen et al. [6], which is based on a principal component (PC) decomposition method [37], was further improved and used to purify the InSAR signals related to surface deformation. The land surface displacement results were illustrated and their different origins (e.g., subway construction, urban expansion, mining subsidence and geological conditions) were subsequently studied and discussed. This study represents the first investigation of land surface deformation over the Xuzhou metropolitan area using InSAR/TS-InSAR. It is hoped that this research can provide vital information for the mitigation and prevention of urban damage in the city.

2. 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.

3. Dataset and Methodology

3.1. 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.

3.2. 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.

3.2.1. 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).

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):

ϕ_{i} = a x_{i} + b y_{i} + c x_{i} y_{i} + d x_{i}^{2} + e y_{i}^{2} + f z_{i} + g z_{i}^{2} + h

(1)

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, \dots, 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.

3.2.2. 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).

X_{m \times n} = U_{m \times k} \cdot S_{k \times k} \cdot V_{k \times n}^{T}

(2)

where

m

is the number of pixels,

n

is the number of epochs,

k

is the number of PCs,

X_{m \times n}

denotes the centered data matrix,

U_{m \times k}

and

V_{k \times n}

represent the spatial function and temporal function of the components, respectively,

S_{k \times k}

is a diagonal matrix whose values signify the significance of each component and

V_{k \times n}^{T}

denotes the transpose of

V_{k \times 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_{S 1}

= 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.

3.2.3. 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.

4. Results

4.1. 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 (Figure 1 and Figure 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.

4.2. 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.

5. Interpretation and Discussion

5.1. 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 (Figure 8 and Figure 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.

5.2. 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² [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 Figure 1 and Figure 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.

5.3. 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.

6. Conclusions

This study proposed an enhanced PCA-based TS-InSAR methodology to detect ground displacement during a period of three years in Xuzhou (China) using C-band Sentinel-1A SAR images. Detailed information on the spatial pattern and temporal evolution of land subsidence and uplift within the study area was investigated for the first time, based on which the following conclusions were drawn:

The powerful potential of the enhanced PCA-based correction method for purifying InSAR displacement time series was demonstrated. The significant reduction of the variance of the interferograms and the reasonable agreement between results derived from InSAR and GNSS measurements demonstrated the method to be both efficient and effective for monitoring high-precision land surface displacement in the Xuzhou region. The success of this method suggests it might have significant potential for application in other ground displacement investigations.
Noticeable land subsidence with displacement in the range of –5 to –41 mm/yr were found widely within the urban areas of Xuzhou during the study period, particularly along the subway lines under construction, the newly developed district and in old coal goafs. This indicates that anthropogenic activities such as subway tunneling, building construction and mining could be the main factors contributing to the detected subsidence.
Remarkable long-term land uplift signals with rates of up to +25 mm/yr have begun to affect two long narrow areas within the old goafs since 2015. It is suggested that the high rate of uplift could be associated both with specific geological conditions and with rising underground water levels that could contribute to the land uplift either directly or indirectly by inflating the compositions of the unconsolidated layer.
Our results and interpretations could provide important insight into the potential instability of central areas of Xuzhou. Regular monitoring of surface displacement is needed to prevent related geohazards.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-4292/11/12/1494/s1, Figure S1: (a) Displacement rate map derived from initial InSAR time series for the entire study area during the period from November 2015 to June 2018. Positive values represent decreases in earth-satellite distance, while negative values represent increases in earth-satellite distance. Location of P1–P4 are indicated by black circles. Filled red star indicates the location of the reference area. The map conventions are the same as in Figure 4. (b) Displacement time series of points P1–P4. The red lines indicate the linear fits. Figure S2: Correlations between the PCs and the elevation.

Author Contributions

Conceptualization, Y.C. and K.T.; Data curation, Y.C.; formal analysis, Y.C.; funding acquisition, Y.C.; investigation, Y.C., S.Y., K.Z. and H.Z.; methodology, Y.C., K.T. and S.Y.; resources, Y.C. and K.T.; software, Y.C.; supervision, K.T.; validation, Y.C. and X.L.; visualization, Y.C.; writing—original draft, Y.C.; writing—review & editing, Y.C., K.T., S.Y., K.Z., H.Z., X.L., H.L. and Y.S.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 41804013), the Natural Science Foundation of Jiangsu Province (Grant No. BK20180666), and the Fundamental Research Funds for the Central Universities (Grant No. 2017XKQY019). This work was supported by a project founded by the Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions. We would like to acknowledge the support of the Jiangsu dual creative teams programme project awarded in 2017 (Grant No. CUMT07180005).

Acknowledgments

The authors would like to thank the European Space Agency (ESA) for their generosity in providing the Sentinel-1A SAR data and SNAP software for free, and Hooper and his research team at the University of Leeds for providing the StaMPS software. We also thank the Xuzhou Land and Resources Bureau for providing the vector map of the subway lines, the GNSS data and the Xuzhou urban planning area resources and environment carrying capacity evaluation report.

Conflicts of Interest

The authors declare no conflict of interest.

References

Stramondo, S.; Bozzano, F.; Marra, F.; Wegmuller, U.; Cinti, F.R.; Moro, M.; Saroli, M. Subsidence induced by urbanisation in the city of Rome detected by advanced InSAR technique and geotechnical investigations. Remote Sens. Environ. 2008, 112, 3160–3172. [Google Scholar] [CrossRef]
Remy, D.; Falvey, M.; Bonvalot, S.; Chlieh, M.; Gabalda, G.; Froger, J.L.; Legrand, D. Variability of atmospheric precipitable water in northern Chile: Impacts on interpretation of InSAR data for earthquake modeling. J. S. Am. Earth Sci. 2011, 31, 214–226. [Google Scholar] [CrossRef]
Fan, H.; Gao, X.; Yang, J.; Deng, K.; Yu, Y. Monitoring mining subsidence using a combination of phase-stacking and offset-tracking methods. Remote Sens. 2015, 7, 9166–9183. [Google Scholar] [CrossRef]
Dai, K.; Li, Z.; Tomas, R.; Liu, G.; Yu, B.; Wang, X.; Cheng, H.; Chen, J.; Stockamp, J. Monitoring activity at the Daguangbao mega-landslide (China) using Sentinel-1 TOPS time series interferometry. Remote Sens. Environ. 2016, 186, 501–513. [Google Scholar] [CrossRef] [Green Version]
Chaussard, E.; Wdowinski, S.; Cabral-Cano, E.; Amelung, F. Land subsidence in central Mexico detected by ALOS InSAR time-series. Remote Sens. Environ. 2014, 140, 94–106. [Google Scholar] [CrossRef]
Chen, Y.; Remy, D.; Froger, J.-L.; Peltier, A.; Villeneuve, N.; Darrozes, J.; Perfettini, H.; Bonvalot, S. Long-term ground displacement observations using InSAR and GNSS at Piton de la Fournaise volcano between 2009 and 2014. Remote Sens. Environ. 2017, 194, 230–247. [Google Scholar] [CrossRef]
Yin, Y.P.; Zhang, Z.C.; Zhang, K.J. Land subsidence and countermeasures for its prevention in China. Chin. J. Geol. Hazard Control 2005, 16, 1–8. [Google Scholar] [CrossRef]
He, Q.C.; Ye, X.B.; Li, Z.M.; Liu, W.B. The status and prevention strategy of land subsidence in China. Geol. J. China Univ. 2006, 12, 161–168. [Google Scholar] [CrossRef]
Hua, X.Q.; Huang, J.J.; Miao, S.X.; Zhang, L.; Wu, X. Distribution and causes of geo-hazards in Xuzhou. J. Geol. Hazards Environ. Preserv. 2015, 26, 74–80. [Google Scholar]
Hua, X.Q.; Huang, J.J.; Wei, Y.Y.; Miao, S.X.; Xu, S.Y.; Jiang, S.; Zhang, L. Demand analysis and key problems of the urban geological survey in Xuzhou. J. Geol. 2015, 39, 322–328. [Google Scholar]
Gao, M.; Gong, H.; Chen, B.; Zhou, C.; Chen, W.; Liang, Y.; Shi, M.; Si, Y. InSAR time-series investigation of long-term ground displacement at Beijing Capital International Airport, China. Tectonophysics 2016, 691, 271–281. [Google Scholar] [CrossRef]
Perissin, D.; Wang, Z.; Lin, H. Shanghai subway tunnels and highways monitoring through Cosmo-SkyMed Persistent Scatterers. ISPRS J. Photogramm. Remote Sens. 2012, 73, 58–67. [Google Scholar] [CrossRef]
Qin, X.; Yang, T.; Yang, M.; Zhang, L.; Liao, M. Health diagnosis of major transportation infrastructures in Shanghai metropolis using high-resolution Persistent Scatterer Interferometry. Sensors 2017, 17, 2770. [Google Scholar] [CrossRef] [PubMed]
Wang, H.; Feng, G.; Xu, B.; Yu, Y.; Li, Z.; Du, Y.; Zhu, J. Deriving spatio-temporal development of ground subsidence due to subway construction and operation in delta regions with PS-InSAR data: A case study in Guangzhou, China. Remote Sens. 2017, 9, 1004. [Google Scholar] [CrossRef]
Ng, A.H.-M.; Wang, H.; Dai, Y.; Pagli, C.; Chen, W.; Ge, L.; Du, Z.; Zhang, K. InSAR reveals land deformation at Guangzhou and Foshan, China between 2011 and 2017 with COSMO-SkyMed data. Remote Sens. 2018, 10, 813. [Google Scholar] [CrossRef]
Meng, L.; Feng, Q.Y.; Zhang, H.R.; Wu, Q.Q.; Lu, Y. Mechanism analysis and hazard assessment on land subsidence in Datun central district of Xuzhou city. Chin. J. Geol. Harzard Control 2008, 19, 60–63. [Google Scholar] [CrossRef]
Gao, J.X.; Hu, H. Advanced GNSS technology of mining deformation monitoring. Prodedia Earth Planet. Sci. 2009, 1, 1081–1088. [Google Scholar] [Green Version]
Bamler, R.; Hartl, P. Synthetic aperture radar interferometry. Inverse Probl. 1999, 14, 12–13. [Google Scholar] [CrossRef]
Rosen, P.A.; Hensley, S.; Joughin, I.R.; Li, F.K.; Madsen, S.N.; Rodriguez, E.; Goldstein, R.M. Synthetic aperture radar interferometry. Proc. IEEE 2002, 88, 333–382. [Google Scholar] [CrossRef]
Remy, D.; Chen, Y.; Froger, J.L.; Bonvalot, S.; Cordoba, L.; Fustos, J. Revised interpretation of recent InSAR signals observed at Llaima volcano (Chile). Geophys. Res. Lett. 2015, 42, 3870–3879. [Google Scholar] [CrossRef]
Massonnet, D.; Rossi, M.; Carmona, C.; Adragna, F.; Peltzer, G.; Feigl, K.; Rabaute, T. The displacement field of the Landers earthquake mapped by radar interferometry. Nature 1993, 364, 138–142. [Google Scholar] [CrossRef]
Cascini, L.; Fornaro, G.; Peduto, D. Advanced low- and full-resolution DInSAR map generation for slow-moving landslide analysis at different scales. Eng. Geol. 2010, 112, 29–42. [Google Scholar] [CrossRef]
Tessitore, S.; Fernandez-Merodo, J.A.; Herrera, G.; Tomas, R.; Ramondini, M.; Sanabria, M.; Duro, J.; Mulas, J.; Calcaterra, D. Comparison of water-level, extensometric, DInSAR and simulation data for quantification of subsidence in Murcia City (SE Spain). Hydrogeol. J. 2016, 24, 727–747. [Google Scholar] [CrossRef]
Hooper, A.; Segall, P.; Zebker, H. Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcan Alcedo, Galapagos. J. Geophys. Res. Solid Earth 2007, 112. [Google Scholar] [CrossRef]
Kampes, B. Displacement Parameter Estimation Using Permanent Scatterer Interferometry. Ph.D. Thesis, Technische Universiteit Delft, Delft, The Netherlands, 2005. [Google Scholar]
Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef] [Green Version]
Hooper, A. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophys. Res. Lett. 2008, 35. [Google Scholar] [CrossRef] [Green Version]
Zeni, G.; Bonano, M.; Casu, F.; Manunta, M.; Manzo, M.; Marsella, M.; Pepe, A.; Lanari, R. Long-term deformation analysis of historical buildings through the advanced SBAS-DInSAR technique: The case study of the city of Rome, Italy. J. Geophys. Eng. 2011, 8, S1–S12. [Google Scholar] [CrossRef]
Aslan, G.; Cakir, Z.; Ergintav, S.; Lasserre, C.; Renard, F. Analysis of secular ground motions in Istanbul from a long-term InSAR time-series (1992–2017). Remote Sens. 2018, 10, 408. [Google Scholar] [CrossRef]
Zhou, C.; Gong, H.; Chen, B.; Li, J.; Gao, M.; Zhu, F.; Chen, W.; Liang, Y. InSAR time-series analysis of land subsidence under different land use types in the eastern Beijing Plain, China. Remote Sens. 2017, 9, 380. [Google Scholar] [CrossRef]
Zhou, L.; Guo, J.; Hu, J.; Li, J.; Xu, Y.; Pan, Y.; Shi, M. Wuhan surface subsidence analysis in 2015–2016 based on Sentinel-1A data by SBAS-InSAR. Remote Sens. 2017, 9, 982. [Google Scholar] [CrossRef]
Bai, L.; Jiang, L.; Wang, H.; Sun, Q. Spatiotemporal characterization of land subsidence and uplift (2009–2010) over Wuhan in central China revealed by TerraSAR-X InSAR analysis. Remote Sens. 2016, 8, 350. [Google Scholar] [CrossRef]
Fan, H.; Deng, K.; Ju, C.; Zhu, C.; Xue, J.Q. Land subsidence monitoring by D-InSAR technique. Min. Sci. Technol. 2011, 21, 869–872. [Google Scholar] [CrossRef]
Liu, X.; Wang, Y.; Yan, S. Interferometric SAR time series analysis for ground subsidence of the abandoned mining area in north Peixian using Sentinel-1A TOPS data. J. Indian Soc. Remote Sens. 2018, 46, 451–461. [Google Scholar] [CrossRef]
Hooper, A.; Bekaert, D.; Spaans, K.; Arikan, M. Recent advances in SAR interferometry time series analysis for measuring crustal deformation. Tectonophysics 2012, 514, 1–13. [Google Scholar] [CrossRef]
Kositsky, A.P.; Avouac, J.P. Inverting geodetic time series with a principal component analysis-based inversion method. J. Geophys. Res. Solid Earth 2010, 115. [Google Scholar] [CrossRef] [Green Version]
Geological Disaster Prevention and Control Planning of Jiangsu Province; Department of Land and Resources of Jiangsu Province: Nanjing, China, 2017.
Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The shuttle radar topography mission. Rev. Geophys. 2001, 45, RG2004. [Google Scholar] [CrossRef]
Hooper, A.; Zebker, H.A. Phase unwrapping in three dimensions with application to InSAR time series. J. Opt. Soc. Am. A Opt. Image Sci. Vis. 2007, 24, 2737–2747. [Google Scholar] [CrossRef] [Green Version]
Lyons, S.; Sandwell, D. Fault creep along the southern San Andreas from interferometric synthetic aperture radar, permanent scatterers, and stacking. J. Geophys. Res. Solid Earth 2003, 108. [Google Scholar] [CrossRef]
Science Toolbox Exploitation Platform (STEP) SNAP Download Page. Available online: http://step.esa.int/main/download/snap-download (accessed on 15 January 2018).
Sentinel-1 POD Precise Orbit Ephemerides Products Download Page. Available online: https://qc.sentinel1.eo.esa.int/aux_poeorb (accessed on 11 February 2018).
Pepe, A.; Lanari, R. On the extension of the minimum cost flow algorithm for phase unwrapping of multitemporal differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2374–2383. [Google Scholar] [CrossRef]
Ferretti, A.; Prati, C.; Rocca, F. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2000, 38, 2202–2212. [Google Scholar] [CrossRef] [Green Version]
Lin, Y.-N.N.; Kositsky, A.P.; Avouac, J.-P. PCAIM joint inversion of InSAR and ground-based geodetic time series:application to monitoring magmatic inflation beneath the Long Valley Caldera. Geophys. Res. Lett. 2010, 37. [Google Scholar] [CrossRef]
Remy, D.; Froger, J.L.; Perfettini, H.; Bonvalot, S.; Gabalda, G.; Albino, F.; Cayol, V.; Legrand, D.; De Saint Blanquat, M. Persistent uplift of the Lazufre volcanic complex (Central Andes): New insights from PCAIM inversion of InSAR time series and GPS data. Geochem. Geophys. Geosyst. 2014, 15, 3591–3611. [Google Scholar] [CrossRef]
Yang, Y.; Zeng, A.; Wu, F. Horizontal crustal movement in China fitted by adaptive collocation with embedded Euler vector. Sci. China Earth Sci. 2011, 54, 1822–1829. [Google Scholar] [CrossRef]
Huang, J.J.; He, H.Z.; Lu, H.; Xu, S.Y.; Cui, L.Y.; Zhang, L.; Jiang, S.; Hua, M.; Wu, X.; Wei, Y.Y.; et al. Xuzhou Urban Planning Area Resources and Environment Carrying Capacity Evaluation Report; Xuzhou Land and Resources Bureau: Xuzhou, China, 2017. [Google Scholar]
Zhao, Q.; Lin, H.; Jiang, L.; Chen, F.; Cheng, S. A study of ground deformation in the Guangzhou urban area with Persistent Scatterer Interferometry. Sensors 2009, 9, 503–518. [Google Scholar] [CrossRef]
Strozzi, T.; Delaloye, R.; Poffet, D.; Hansmann, J.; Loew, S. Surface subsidence and uplift above a headrace tunnel in metamorphic basement rocks of the Swiss Alps as detected by satellite SAR interferometry. Remote Sens. Environ. 2011, 115, 1353–1360. [Google Scholar] [CrossRef] [Green Version]
Roccheggiani, M.; Piacentini, D.; Tirincanti, E.; Perissin, D.; Menichetti, M. Detection and monitoring of tunneling induced ground movements using Sentinel-1 SAR Interferometry. Remote Sens. 2019, 11, 639. [Google Scholar] [CrossRef]
Policy Archive Page of Government of Jiangsu Province. Available online: http://www.jiangsu.gov.cn/art/2015/3/16/art_65017_348358.html (accessed on 15 June 2018).
Feng, X.; Zhang, Q. The effect of backfilling materials on the deformation of coal and rock strata containing multiple goaf: A Numerical Study. Minerals 2018, 8, 224. [Google Scholar] [CrossRef]
Feng, X.; Wang, E.; Ganne, J.; Martin, R.; Jessell, M. The exhumation along the Kenyase and Ketesso shear zones in the Sefwi terrane, West African Craton: A numerical study. Geosci. J. 2018, 1–18. [Google Scholar] [CrossRef]
Feng, X.; Amponsah, P.O.; Martin, R.; Ganne, J.; Jessell, M.W. 3-D numerical modelling of the influence of pre-existing faults and boundary conditions on the distribution of deformation: Example of North-Western Ghana. Precambrian Res. 2016, 274, 161–179. [Google Scholar] [CrossRef]

Figure 1. (a) Location of Huaihai Economic Zone in China. The red rectangle indicates the outline of (b). (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.

Figure 2. Flow chart of the data processing.

Figure 3. Perpendicular and temporal baselines of interferograms. Red diamonds indicate the master image, blue diamonds represent slave images and black lines represent interferograms formed from the master image and corresponding slave images.

Figure 4. Results from the principal component (PC) decomposition applied to the initial InSAR time series, including the spatial function weighted by significance value (U*S), temporal function (V) and proportion of significance value (Ps) for each PC (sub-figures (a)–(f) for PC1–PC6). The scatter plots indicate the correlations between the PCs and the elevation (see Section 2 in the Supplementary Materials for enlarged scatter plots). Maps are overlain draped on a Google Earth image and projected in geographic Lat/Lon referenced to the WGS-84 datum.

Figure 5. Line-of-sight displacement comparison between GNSS and InSAR time series derived from different strategies at stations TSGT, CUMT and HCXY (see Figure 1 for locations). Gray lines in the background indicate GNSS measurements. Green circles, magenta squares and blue crosses indicate the measurements of the initial InSAR time series and the InSAR time series corrected by Strategies 2 and 3, respectively. Error bars represent standard deviations within 20 m neighboring areas.

Figure 6. Displacement rate map derived from InSAR time series corrected by Strategy 2, for the entire study area during November 2015 to June 2018. Positive values represent decreases in earth–satellite distance, while negative values represent increases in earth–satellite distance. The subway lines are marked by black lines. Zone 1, Zone 2 WG (the western goaf), Zone 3 EG (the eastern goaf) and Zone 4 J-MG (Jiuli-Ma village goaf) outlined by dashed black lines are the main displacement areas in addition to the subway line buffers. Location of Shitun, Jiahe, Pangzhuang, and Zhangxiaoloujing old coal mines in Zone 2 and of Qishan and Quantai old coal mines in Zone 3 are indicated by black circles. The two white rectangles represent areas which are enlarged and analyzed in detail in Section 5.2. Filled red triangles indicate the locations of the GNSS stations used in this study. Red rectangles indicate zones where no displacement is expected to occur (corresponding to “non-displacement zones” mentioned in the text). The filled red star indicates the location of the reference area. The map conventions are the same as in Figure 4.

Figure 7. Histograms of standard deviations of displacement rates derived from the initial time series and corrected time series obtained from Strategy 2.

Figure 8. (a) Displacement rate map (2015–2018) within a 1 km-wide buffer along the subway lines, extracted from Figure 6. All circles are subway station locations, among which the brown circles with station names indicate those discussed in the paper. The inset in the bottom-right corner is an enlarged map of the area marked by the white rectangle. Points A, B, C and D are discussed in the paper and are marked by black crosses. The map conventions are the same as in Figure 4. (b) The profile of P-P’ in the inset of (a). (c) Damages to a local house caused by land subsidence.

Figure 9. Displacement time series and displacement rates at points A, B, C and D shown in Figure 8a.

Figure 10. (a) Enlarged map of displacement rate (2015–2018) around Shitun mine (the western white rectangle in 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.

Table 1. Details of the processed SAR data.

Items	Description
Satellite	Sentinel-1A
Acquisition mode	TOPS
Track number	142
Orbit direction	Ascending
Polarization	VV
Wavelength (m)	0.0555
Range resolution (m)	2.33
Azimuth resolution (m)	13.94
Incidence angle (deg)	40.12
Number of images	52
Time spans	27 November 2015–8 June 2018

Table 2. Comparison of the root mean squares (RMS) in non-displacement zones and the RMS of the differences between GNSS and InSAR derived from different correction strategies (mm).

	Initial Time Series	Long-Wavelength Artifact Correction	PC Decomposition-Based Artifact Correction
	Initial Time Series	Long-Wavelength Artifact Correction	Correction Strategy 1 ^a (Chen et al. [6])	Correction Strategy 2 ^b	Correction Strategy 3 ^c
Non-displacement zones	15.17	11.02	9.89	4.75	4.67
TSGT	13.87	9.28	7.51	3.34	3.47
CUMT	15.97	12.75	9.45	3.28	3.33
HCXY	12.36	9.06	8.80	5.82	5.81

^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.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Chen, Y.; Tan, K.; Yan, S.; Zhang, K.; Zhang, H.; Liu, X.; Li, H.; Sun, Y. Monitoring Land Surface Displacement over Xuzhou (China) in 2015–2018 through PCA-Based Correction Applied to SAR Interferometry. Remote Sens. 2019, 11, 1494. https://doi.org/10.3390/rs11121494

AMA Style

Chen Y, Tan K, Yan S, Zhang K, Zhang H, Liu X, Li H, Sun Y. Monitoring Land Surface Displacement over Xuzhou (China) in 2015–2018 through PCA-Based Correction Applied to SAR Interferometry. Remote Sensing. 2019; 11(12):1494. https://doi.org/10.3390/rs11121494

Chicago/Turabian Style

Chen, Yu, Kun Tan, Shiyong Yan, Kefei Zhang, Hairong Zhang, Xiaoyang Liu, Huaizhan Li, and Yaqin Sun. 2019. "Monitoring Land Surface Displacement over Xuzhou (China) in 2015–2018 through PCA-Based Correction Applied to SAR Interferometry" Remote Sensing 11, no. 12: 1494. https://doi.org/10.3390/rs11121494

APA Style

Chen, Y., Tan, K., Yan, S., Zhang, K., Zhang, H., Liu, X., Li, H., & Sun, Y. (2019). Monitoring Land Surface Displacement over Xuzhou (China) in 2015–2018 through PCA-Based Correction Applied to SAR Interferometry. Remote Sensing, 11(12), 1494. https://doi.org/10.3390/rs11121494

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu