Investigation of Ground Deformation in Taiyuan Basin, China from 2003 to 2010, with Atmosphere-Corrected Time Series InSAR

Abstract

Excessive groundwater exploitation is common through the Taiyuan basin, China, and is well known to result in ground subsidence. However, most ground subsidence studies in this region focus on a single place (Taiyuan city), and thus fail to demonstrate the regional extent of the deformation phenomena in the whole basin. In this study, we used Interferometric Synthetic Aperture Radar (InSAR) time series analysis to investigate land subsidence across the entire Taiyuan basin region. Our data set includes a total of 75 ENVISAT ASAR images from two different frames acquired from August 2003 to September 2010 and 33 TerraSAR-X scenes spanning between March 2009 and March 2010. ERA-Interim reanalysis was used to correct the stratified delay to reduce the bias expected from the systematic components of tropospheric delay. The residual delay after correction of stratified delay can be considered as a stochastic component and be mitigated through spatial-temporal filtering. A comparison with MERIS (Medium-Resolution Imaging Spectrometer) measurements indicates that our atmospheric corrections improved agreement over the conventional spatial-temporal filtering by about 20%. The displacement results from our atmosphere-corrected time series InSAR were further validated with continuous GPS data. We found eight subsiding centers in the basin and a surface uplift to the north of Taiyuan city. The causes of ground deformation are analyzed and discussed in relation to gravity data, pre-existing faults, and types of land use.

1. Introduction

Ground subsidence is the lowering of the Earth’s surface caused by natural processes as well as human activities, such as natural consolidation of soils [1], tectonic motion [1], extraction of natural resources such as groundwater [2,3], gas and coal [4], thawing permafrost [5], and construction of high-rise buildings [6]. Starting from the 1960s and continuing up to the present day, several areas in the Taiyuan basin have experienced remarkable land subsidence due to excessive groundwater pumping related to urbanization and industrialization [7,8,9]. One of the significant effects of land subsidence in this area centers on the water conveyance infrastructures of the Shanxi Wanjiazhai Yellow River Water Diversion Project (YRDP) that traverses a subsiding area [7,8]. Some parts of those structures are gravity driven through the use of very small gradients; therefore, even small changes in these gradients can cause a reduction in the water flow capacity. Moreover, structural damages to pipelines, buildings, roads, bridges, and high speed railways is an ever-present risk. Some buildings have been deemed unsafe and thus have been closed, and some of them have needed repair to keep them intact. Furthermore, some communities and villages must be relocated.

From 1956 to 2007, subsidence phenomena in the basin have been studied using spirit-leveling and GPS surveys and possible contributors to ground deformation have been documented [10,11,12]. Although the conventional geodetic tools can provide accurate measurements, they are time consuming, labor intensive and costly. Thus, the existing monitoring networks are limited in spatial and temporal resolution, lacking the capacity to fully capture the land deformation over a long period across a large area of the basin. Hence, space-borne synthetic aperture radar interferometry (InSAR) can play an important role in complementing these conventional techniques for providing maps of ground displacement over large areas (tens of thousands of square kilometer) at a spatial resolution of a few meters. Furthermore, InSAR time series techniques such as persistent scatterer InSAR (PSI) and small baseline subset (SBAS) can provide information on the slow progression of deformation over time at millimeter level accuracy by processing a multiple of SAR acquisitions [13,14].

Atmospheric phase delay is the main source of uncertainty for InSAR measurements [15]. Propagation of the microwave signal is affected by both the ionosphere and troposphere. Variations of free electrons along the travel path in ionosphere cause a phase advance or a group delay for radar signal. Several methods have been developed to correct the ionospheric effects [16,17]. Variations in pressure, temperature, and relative humidity in the lower part of the troposphere (<5 km) result in phase delay in an interferogram with a magnitude up to 15–20 cm, and can often mask or lead to false interpretation of ground displacement [15]. Correcting for tropospheric delay remains one of the major challenges in increasing the accuracy of ground deformation measurements using InSAR.

Several approaches have been recommended for mitigating the tropospheric delay in the InSAR data. One approach splits the tropospheric delay into two components: the stratified delay and non-stratified delay. Typical algorithms based on this idea are the multiscale approach [18] and the spatially variable power law method [15]. Another approach estimates tropospheric delay from GPS observations [19,20,21]. This method relies upon a dense GPS network, which is not always available in the study area. A third method is based on the atmospheric water vapor observations from satellite multispectral imagery such as MODIS (Moderate Resolution Imaging Spectroradiometer) and MERIS (Medium-Resolution Imaging Spectrometer) [22]. Multispectral imagery is limited to daylight hours and cloud-free conditions. A fourth technique exploits global atmospheric reanalysis such as ERA-Interim (ERA-I), the North American Regional Reanalysis (NARR), and the Modern Era-Retrospective Analysis for Research and Applications (MERRA) [15,23]. In spite of their coarse spatial and temporal resolution, those atmospheric models are successfully applied to mitigate the long-wavelength tropospheric stratification components in different regions.

In the absence of any auxiliary data for atmospheric mitigation, a fifth approach depends on a combination of filtering in time and space [24,25,26]. The basic idea behind the spatial-temporal filtering is that the atmospheric signal is correlated in spatial domain and uncorrelated in time domain, whereas the ground deformation is correlated in both space and time. The uncorrelated in time assumption, however, is not rigorous because seasonal variations in the moisture content of the atmosphere introduces a systematic component of tropospheric delay, which is correlated with time. The seasonality of tropospheric delay can not be effectively captured by using the spatial-temporal filtering, resulting in bias and uncertainty in the displacement time series and velocity fields. Even using more SAR acquisitions that sampled evenly in time does not necessarily reduce the bias and uncertainty [27].

This paper demonstrates the seasonal effects of atmospheric delay in our study area and develops an atmosphere-corrected time series InSAR method to reduce the delay more effectively. The stratified delay was simulated and removed by using ERA-I atmospheric model; then, the residual delays after correction of stratified delay were considered randomly in time, and were further mitigated using spatial-temporal filtering. We compared the atmospheric correction results to MERIS data and validated the InSAR displacement results using continuous GPS measurements. We discussed the causes and temporal evolution of the observed deformation, focusing on specific areas of the basin, by analyzing gravity changes, pre-existing faults, and types of land use. The availability and the detail of all of the above information has allowed us to better understand the observed subsidence and its relationships with relevant environmental variables.

2. Study Area

Taiyuan basin in central Shanxi province is situated in the north of China, in an area surrounded to the west by the Lvliang Moutains and to the east by the Taihang Mountains. Topography is high in the surrounding mountains (up to 3000 m) and is quite flat in the central basin with elevation in the range of 735 m–830 m. Such terrain features resulted from the neotectonic motion that formed the Taiyuan basin [28]. The total area of the basin is approximately 6159 km${}^2$ stretching from the northeast to southwest with a length of about 150 km in north–south direction and a width of 12–60 km in the east–west direction. This area hosts a total of eleven cities and about 5 million people, while 2 million of them are living in Taiyuan, the capital city of Shanxi province. Areas around the central basin are a broad and flat alluvial plain. The middle section of Fenhe River runs through the Taiyuan basin from north to south. Along the river bank are neatly arranged agricultural fields, primarily consisting of grasslands and various agricultural crops, making severe temporal deccorrelation for C-band SAR in this region. The topography in Figure 1a and optical image in Figure 1b show the relief feature and the distribution of main cities in the basin. The basin is characterized by a semiarid climate. The mean annual precipitation ranges from 400 to 650 mm and above 70% of the annual rainfall occurs between June and September. Winter droughts are common because this area is subjected to the full force of the dry northwest wind that blows in the winter from the Mongolian Plateau.

Taiyuan basin is a representative large-scale Cenozoic fault basin and its western and eastern sides are bounded by two boundary regional faults, the Jiaocheng Fault and Taigu Fault (see Figure 1a), respectively. The geological structures are principally controlled by the joint effect of Jiaocheng Depression of late Cenozoic Fault Zone and the western fault terraces [28]. The exposed or deeply buried stratums around it are mainly Ordovician stratums and continental stratums of the Permian system, carboniferous system, coal-bearing series and Triassic System. The interior deposits contain the stratums of Cenozoic Erathem with a depth of about 1000–2000 m. Most of the stratums locating in interior bases are the sand shale of the Permian system and Triassic system. Loose rock groups of the tertiary and Quaternary systems are distributed throughout the interior areas of the basin [28]. The area receives only a small amount of rainfall annually. Subsequently, the residents in this area have extracted a large amount of groundwater for various purposes over the last three decades, as the region experienced both severe droughts and rapid economic growth. Over time, over-exploitation caused groundwater-level declines accompanied by soil compaction in the aquifer system, resulting in land subsidence. Currently, the total area of subsidence in Taiyuan city is about 548 km${}^2$ and five subsidence centers have formed: Xizhang, Wanbailin, Xiayuan, Wujiabao, and Xiaodian (Figure 2). The land subsidence history in Taiyuan city can be divided into four stages as shown in Figure 2. From 1956 to 1981, only one subsidence center was found in Wujiabao, the maximum subsidence was up to 850 mm at the beginning of 1980s (Figure 2a). Three more subsidence centers (Xizhang, Wanbailin, and Xiayuan) were formed in the period of 1980–1989, and the accumulative subsidence was rapidly growing (over 1600 mm in Wujiabao), indicating that the subsidence rate was accelerating during this period (Figure 2b). During 1990–2000, a new subsidence center (Xiaodian) in the south of Taiyuan city began to appear and the other subsidence centers was still subsiding at an accelerating rate (Figure 2c). The subsidence center (Xiaodian) was formed during 2000–2007, which can be seen in Figure 2d, but the subsidence rates in the other subsidence centers were slowing down during this stage. In the case of the central area of Taiyuan basin, the lands are mainly used for agriculture; here, farmers highly rely upon underground water for dry land agriculture and plantation activities. Excessive pumping in such irrigated agricultural areas has caused severe land subsidence. South of the basin is one of the major coal-producing areas in Shanxi Province in places such as Xiaoyi and Jiexiu city. Mining activities affect the nearby environment, causing ground fractures and damage to the villages in this area.

3. Available Data Sets

3.1. SAR Data

Three groups of SAR images (

Table 1. SAR images used in this work, N is the number of SAR images, M is the number of interferograms.

SAR Sensor	Track, Frame	Time Span	Orbit Mode	$\mathit{\theta }$ (Degree)	N	M
ENVISAT ASAR	T75, F23	20030817-20100919	descending	23	39	181
ENVISAT ASAR	T75, F24	20040104-20100919	descending	23	36	186
TerraSAR-X	-	20090321-20100323	ascending	30	33	32

) were used to map the ground deformation over Taiyuan basin. These include 39 ENVISAT ASAR descending acquisitions (Track 75, Frame 23) collected between 17 August 2003 and 19 September 2010, 36 ENVISAT ASAR descending acquisitions (Track 75, Frame 24) collected between 4 January 2004 and 19 September 2010, and 33 TerraSAR-X ascending scenes collected between 23 March 2009 and 21 March 2010. The ASAR images were acquired at a medium resolution of 25 m and with a large ground coverage of approximately 100 km × 100 km. The two ASAR datasets from two different frames allow us to investigate the land displacement over the whole area of the basin. TerraSAR-X scenes were only about 10 km × 10 km and cover mainly Taiyuan city, but at a higher spatial resolution (3 m) and thus more suitable for infrastructures deformation monitoring [29,30]. These SAR datasets at different ground coverages, wavelengths, revisiting periods and spatial resolutions allow for a comprehensive analysis of the spatial and temporal patterns of ground deformation across the whole basin.

3.2. Continuous GPS

In 2009, a continuous Global Positioning System (CGPS) network with 67 stations was established in the Shanxi province to support a range of tasks such as geodetic measurements for crustal deformation studies and meteorological research. There are seven CGPS stations (A001, A002, K001, K002, K003, J004, J005) located in the basin. CGPS locations are indicated as the red triangles in Figure 1a. The A001 station is located in Taiyuan downtown area, while A002, K001, K002, K003, J004, J005 are located in Qingxu county, Jinzhong city, Qixian county, Jiexiu city, Fenyang city, and Jiaocheng county, respectively. Details of the processing strategy and estimation of uncertainties for CGPS measurements are described in [31]. The velocity fields of these CGPS stations in three dimensions are shown in

Table 2. CGPS site velocity (mm/year) in the ITRF (International Terrestrial Reference System) 2008 frame.

GPS Station	Longitude	Latitude	North	East	Up
A001	112.5	37.9	−11.53 ± 0.09	31.36 ± 0.09	3.42 ± 0.48
A002	112.4	37.6	−9.17 ± 0.17	28.9 ± 0.15	−35.84 ± 0.81
K001	112.7	37.7	−13.64 ± 0.10	33.31 ± 0.09	0.95 ± 0.43
K002	112.4	37.4	−10.85 ± 0.09	32.30 ± 0.09	−15.72 ± 0.52
K003	111.9	37.0	−8.34 ± 0.09	34.22 ± 0.09	−29.81 ± 0.46
J004	111.8	37.2	−9.41 ± 0.10	31.77 ± 0.09	4.91 ± 0.48
J005	112.2	37.5	−10.68 ± 0.10	37.63 ± 0.11	−5.98 ± 0.47

. We selected K001 as the reference site because it is located on stable ground with a very small vertical rate 0.95 mm/year.

3.3. Gravity Measurements

Gravity measurements at ten stations in the Taiyuan basin obtained from 2001 to 2012 were used in this study to investigate the relationship between ground subsidence and gravity variations. These ten gravity measurement points located in the basin are represented as the white circles (G1 to G10) in Figure 1a. The gravity data include relative and absolute gravity measurements. Relative gravity measuring campaigns from 2001 to 2006 were conducted three or four times annually and the average values at each point were computed as the gravity value for that year. These relative gravity data from 2001 to 2006 are shown in

Table 3. Gravity variations at ten points in Taiyuan basin [32]. Gravity values are in unit of ${10}^{-5}$ ms${}^{-2}$.

Station	Aug 2001	Aug 2002	Aug 2003	Oct 2004	Aug 2005	Aug 2006
G1	−9.9945	−9.9977	−10.0015	−10.0000	−10.0067	−10.0057
G2	−12.1628	−12.1594	−12.1641	−12.1619	−12.1587	−12.1640
G3	−42.0322	−42.0217	−42.0130	−41.9909	−41.9833	−41.9793
G4	−23.8529	−23.8363	−23.8437	−23.8571	−23.8595	−23.8382
G5	−32.5414	−32.5428	−32.5482	−32.5334	−32.5322	−32.5202
G6	−30.7725	−30.7630	−30.7628	−30.7578	−30.7608	−30.7506
G7	−44.4278	−44.4295	−44.4323	−44.4239	−44.4206	−44.4004
G8	−53.5224	−53.5077	−53.5015	−53.5002	−53.5266	−53.5017
G9	−57.0167	−57.0028	−56.9864	−56.9667	−56.9620	−56.9340
G10	−49.1414	−49.1466	−49.1387	−49.1337	−49.1309	−49.1207

. Gravity is affected by mass and distance; a change in mass or a change in surface elevation will cause a change in gravity. In Section 6.2, we use gravity changes and displacement time series from InSAR to analyze the causes of ground deformation in the Taiyuan basin.

4. Methodology

4.1. Tropospheric Delay in InSAR Measurements

The contribution from tropospheric delay to InSAR phase measurements between pixel p and q and between two interferometric acquisitions is the integral of the spatial and temporal differences of refractivity over the signal propagation path L, as [33]

$$\delta L_{p,q}^{m,s}=\left[{\int }_0^{L_p}N(L,s)dL-{\int }_0^{L_p}N(L,m)dL\right]-\left[{\int }_0^{L_q}N(L,s)dL-{\int }_0^{L_q}N(L,m)dL\right],$$

(1)

where $L_p$ and $L_q$ are the SAR signal propagation path at pixel p and q, respectively, and N is the refractivity. m and s represent the master and slave image respectively. Ignoring the ionospheric and liquid term, the refractivity N can be written as

$$N=k_1\frac{p_d}{T}+(k_2\frac{e}{T}+k_3\frac{e}{T_2}),$$

(2)

where $P_d$ is the partial pressure due to dry air in hPa, T is the atmospheric temperature in $Kelvin$, e is the partial pressure of water vapor in hPa. The atmospheric refractivity constants are $k_1=0.776$ KPa${}^{-1}$, $k_2$ = 0.716 KPa${}^{-1}$ and $k_3$ = 3.75 × 10${}^{-3}$ K${}^2$Pa${}^{-1}$. The two refractivity terms on the right-hand side of Equation (2) are referred to as hydrostatic and wet term, respectively, and substituting them in Equation (1) gives the corresponding hydrostatic and wet delay. Producing the vertical profiles of $P_d$, T, and e from the ERA-I model, the total tropospheric delay (both hydrostatic and wet) are estimated using Equations (1) and (2) and interpolated horizontally and vertically with bilinear and spline interpolation operators, respectively. The delay maps for each interferogram were produced by differentiating delay maps at each acquisition used to generate interferogram and then converted to the InSAR line-of-sight (LOS) using a simple mapping function.

4.2. Seasonal Effects of Tropospheric Delay

The tropospheric wet delay observed in interferograms can be subdivided to systematic and stochastic components in both space and time [27]. In the spatial dimension, the systematic and stochastic components refer to stratified and turbulent delay, respectively. In time, we assume the systematic component is dominated by the seasonal periodicities and the stochastic component is the non-seasonal residual. In this section, we demonstrate the systematic and stochastic components in the time domain and evaluate the performance of stratified delay correction with ERA-I model.

To investigate the seasonal effect of wet delay, we fit the time series of Zenith Wet Delay (ZWD) by weighted linear least squares with a slope rate and annual and semiannual terms as

$$ZWD\left(t_i\right)=a+b{t}_i+csin\left(2\pi t_i\right)+d\cos \left(2\pi t_i\right)+esin\left(\pi t_i\right)+f\cos \left(\pi t_i\right)+\nu \left(t_i\right)$$

(3)

with intercept and rate terms (a,b) and seasonal term coefficients (annual- c, d; semiannual- e, f). The measurement errors $\nu$ are treated as Gaussian white noise. As weather conditions remain stable over years, we expect no a trend in the time series of delay, i.e., the slope rate $b=0$. The estimated seasonal terms (c, d, e, f) are calculated using trigonometric relations, to the amplitude and phase of the annual and semiannual components.

We evaluated the seasonal periodicities using time series of zenith wet delay from radiosonde and ERA-I model. By assuming the wet delays obtained from radiosonde providing the truth, we evaluated how it can be approximated by the delay predicted from the ERA-I data. Given the zenith wet delays from radiosonde ($ZW{D}^{radio}$) and atmospheric reanalysis ($ZW{D}^{ERA-I})$), the residual zenith wet delay ($RZWD$) is computed as

$$RZWD\left(t_i\right)=ZW{D}^{radio}\left(t_i\right)-ZW{D}^{ERA-I}\left(t_i\right).$$

(4)

Our wet delay consists of daily 2004–2011 measurements calculated from the radiosonde wet delay at the ZBYN station (gray points in Figure 3a) and daily delay predicted from ERA-I (blue points in Figure 3a), both of the times are at 12:00 a.m. UTC, which is closest to the SAR acquisition times (UTC 2:43 a.m.). The wet delay between two pixels about 100 km away from each other and with an altitude difference of 730 m is dominated by periodic variations (Figure 3a). The time series of the radiosonde zenith wet delay between the two pixels is characterized by annual and semiannual amplitudes of 20 mm and 6 mm, respectively. To evaluate the component of stochastic delay, we subtracted the best fitted seasonal delay, black solid line in Figure 3a, from the zenith wet delay, transforming the skewed distribution with skewness of 1.1 (Figure 3b) to a Gaussian distribution with standard deviation of 21 mm (Figure 3c). From the frequency distributions of the residual delay, we concluded that it can be considered as a stochastic process. The difference (red points in Figure 3a) of the wet delay from radiosonde data and ERA-I model is also a Gaussian distribution (as plotted in Figure 3d). A smaller standard deviation (19 mm) in Figure 3d than in Figure 3c implies that the stratified delay correction with ERA-I not only mitigates the seasonal delay (systematic component) but also reduces some uncertainty due to the stochastic component.

The time series of zenith wet delay from radiosonde and ERA-I data in our study area show a systematic, temporally correlated component dominated by seasonal variations with annual and semiannual periodicities. The seasonal variations skew the distribution of the delay and bias the InSAR deformation time series and velocity fields. Thus, the spatial-temporal filtering strategies can produce biased displacement results. The stratified delay corrections with ERA-I mitigate the effects of the seasonal delay, largely eliminate the systematic biases, and produce the temporally uncorrelated Gaussian residual seen in Figure 3d. Therefore, the stochastic component after removing stratified delay can be properly mitigated using the spatial-temporal filtering in time series analysis.

4.3. Atmosphere-Corrected Time Series InSAR

The Stanford Method for Persistent Scatterers (StaMPS) package was chosen to process SAR data since it is more suitable to estimate ground displacement in natural areas, such as the agricultural areas in the Taiyuan basin. This package includes a PSI method that refers all interferograms to a common master image and a SBAS method that combines interferograms with multiple master images [34]. The PSI method uses PS targets containing a stable dominant scatterer within a SAR resolution cell while the SBAS method exploits distributed scattering pixels. By exploitating of multiple interferograms formed between image pairs that have small perpendicular, temporal and Doppler baselines, the SBAS method could strongly increase the number of stable pixels in natural terrains. In the StaMPS package, pixels selected from PSI and SBAS are referred to PS and slowly decorrelating filtered phase (SDFP), respectively.

After PS/SDFP selection, the wrapped phase of the selected pixels is corrected for DEM error and then unwrapped. After unwrapping, the high-pass (HP) temporal filtering and the low-pass (LP) spatial filtering are used to estimate the atmospheric disturbances assuming that atmospheric effects occur randomly in time and smooth in space. We referred the temporal-spatial filtering in StaMPS as the conventional time series method in our paper. As described in Section 4.2, the spatial-temporal filtering is less useful for removing stratification of tropospheric delays, as there is a risk of aliasing atmospheric effects into displacement results because of its seasonal characteristics, especially in cases with only a small number of SAR acquisitions available at unevenly temporal intervals. It has been found that the InSAR displacement is still retained in the residuals of the delay over volcanic areas and these effects are mainly induced by the systematic stratification of delay [35]. Therefore, we proposed a new approach that corrects the stratified tropospheric delays (systematic component) with an ERA-I model, and the residuals (stochastic component) are considered as Gaussian noise that can be effectively eliminated by time-weighted filtering. Figure 4 shows a flow chart of the proposed approach (atmospheric-corrected time series InSAR). The detailed procedures of this method are in reference [8].

5. Results

To systematically evaluate the effects of the stratified and turbulent delays, we calculated the coefficients of correlation between the uncorrected interferogram and local topography from the SRTM DEM. Larger coefficient values indicate significant stratification of the atmosphere, and therefore low levels of atmospheric turbulence. We compared the coefficient values to the reduction of standard deviation, as shown in Figure 5a. A statistical correlation was found between atmospheric stratification and the standard deviation reduction from Figure 5a. In the ASAR Frame 23 data sets (Figure 5b), we found 134 out of 181 interferograms experienced standard deviation reduction with a maximum improvement of 53% with an average of 16%, after delay correction. For the ASAR Frame 24 data sets (Figure 5c), 92 out of 186 interferograms experienced delay reduction with a maximum reduction of 35% and an average of 14%. For the TerraSAR-X interferograms (Figure 5d), 11 out of 32 interferograms showed reduced delay after correction that the maximum and the mean reduction were 26% and 10%, respectively.

As a demonstration, Figure 6 shows the delay maps of the original interferogram, the ERA-I delay predictions and their differences for three interferograms from each data sets. Both interferograms and ERA-I delays show a significant vertical stratification of atmospheric delay on the mountain regions. In the central basin area, land subsidence is clearly visible in several locations. As seen in the figure, the delays in the surrounding mountains were substantially removed in the difference maps (the third column in Figure 6). This further implies that the considerable mitigation is largely due to the removal of delays, which are strongly correlated with topography. This also can be visualized by plotting the delay versus elevation as shown in Figure 7. Nevertheless, there are some discrepancies between the ERA-I delay simulations and the original interferograms (indicated as the ellipse in the upper left of Figure 6f). Particular attention should be paid to the TerraSAR-X interferogram (Figure 6g) which only has a temporal interval of 11 days in which the deformation signal can be neglected. The southern part of this interferogram is located over the flat basin that the stratified delay is negligible and the delay variations are only caused by turbulent mixing. The magnitude of this turbulent delay is up to −60 mm, which could be mistakenly interpreted as the deformation signal. Although the spatial delay pattern in this area was partially mitigated by ERA-I (Figure 6i), the modeled lateral delay variation (Figure 6h) is much smoother than the interferogram. This implies that the magnitude of this small-scale (5–12 km) turbulence delay is underestimated by the weather model due to its limited resolution of 75 km. In general, the stratified delay is reasonably well represented by the ERA-I, whereas the turbulence delay is poorly modeled and is much smoother than the real delay in the interferograms.

The residual delays after subtracting ERA-I-derived delays were mainly the turbulence delay and were then estimated using the spatial-temporal filtering. We validated the performance of our proposed method by comparing the estimated delays with the independent delays derived from MERIS. In Figure 8, we find that the atmospheric delays estimated from conventional InSAR time series (Figure 8a,d) show disagreement with the delay obtained from MERIS data (Figure 8c,f), especially in mountain regions. These differences are mainly due to the systematic component of tropospheric delay, i.e. the stratified delay. However, the delay acquired by our proposed method (Figure 8b,e) is quite similar to that estimated from MERIS data (Figure 8c,f), even in mountain regions. For the first interferogram, the standard deviation of the residual between the conventional method and MERIS is 10.02, and it reduces to 8.27 when comparing to our proposed method (17% improvement). In the second interferogram, the standard deviation reduces from 11.28 to 9.01 (20% improvement). In addition, it is noted that the spatial variations of the delay maps from the conventional filtering methods are similar to that from our approach, except in mountain regions. This indicates that the residual delays caused by atmospheric turbulence can be effectively estimated by the spatial-temporal filtering. The main difference between the conventional InSAR time series and our proposed method comes in when dealing with the stratified delay.

After correction of the tropospheric delay, we obtained the deformation maps in the Taiyuan basin from the ASAR Frame 23, ASAR Frame 24, TerraSAR-X data sets as shown in Figure 9, Figure 10 and Figure 11, respectively. We projected the line-of-sight (LOS) measurements to the vertical direction under the assumption that there is no significant horizontal surface movement and all of the measured deformation comes from subsidence/uplift. This allows measurements from multiple SAR sensors with different imaging geometries to be compared with each other. The displacement results can also be compared with vertical surface change measured with GPS methods. The velocity map shows that the area of rocky outcrops encircling the basin did not develop subsidence. This result is consistent with previous works [36], which indicated that subsidence occurs only within the basin where a significant thickness of unconsolidated sediments exists. We observed eight land subsidence centers in the basin (black diamonds in Figure 9, Figure 10 and Figure 11). Subsidence rates in these centers are all over −20 mm/year. Seven of them occurred in or nearby the urban areas, from north to south: Wanbailin, Wujiabao, Xiaodian, Qingxu, Jiaocheng, Xiaoyi, and Jiexiu. The remaining one subsiding center is located in agricultural areas outside the urban centers, which is in the west of Taigu and Qi county and is the largest subsidence area in the basin. We also found surface uplift in one location (Xizhang, the northernmost diamond in Figure 9 and Figure 11), which is in the north of Taiyuan city, and only covered by ASAR Frame 23 and TerraSAR-X images.

6. Discussion

6.1. Comparison between InSAR and CGPS Observations

There are only about two years (2009–2010) of overlap between the InSAR and the CGPS measurements. The details of the comparison between GPS- and InSAR-derived deformation velocities are shown in

Table 4. Comparison between CGPS and InSAR time series average velocity (mm/year) of vertical deformation rate. Only overlap period between these two data sets were compared.

CGPS	GPS	ASAR Frame 23		ASAR Frame 24		TerraSAR-X
	V	(2009.084-2010.715)		(2009.084-2010.715)		(2009.216-2010.221)
		${\mathrm{V}}^a$	${\mathrm{V}}^b$	${\mathrm{V}}^a$	${\mathrm{V}}^b$	${\mathrm{V}}^a$	${\mathrm{V}}^b$
A001	1.3	5.7	4.1	-	-	1.3	1.1
2009.349-2010.998
A002	−21.3	−17.7	−22.9	−24.5	−21.7	-	-
2009.793-2010.995
K002	−16.4	-	-	−13.2	−17.0	-	-
2009.327-2010.998
K003	−30.3	-	-	−29.4	−30.2	-	-
2009.327-2010.998
J004	2.3	-	-	0.9	1.1	-	-
2009.327-2010.998
J005	−7.0	-	-	−10.6	−9.9	-	-
2009.327-2010.998

${\mathrm{V}}^a$ is the vertical velocity without atmospheric correction. $V^b$ is the vertical velocity with atmospheric correction by our proposed method.

. Site K001 was excluded from this comparison because it was used as the reference point. Small changes are seen in the InSAR displacement velocities before and after atmosphere correction because these points are located in the flat basin area and were not affected by the stratified delay. In comparison to CGPS deformation rates, we can see that the InSAR atmosphere-corrected results have a better agreement with CGPS. This implies that the significant atmosphere effects in the mountains could cause some uncertainties in the displacement results on the flat basin. Figure 12 compares the displacement time series between CGPS and InSAR. A good agreement between the two techniques is apparent in Figure 12, as the RMSE (Root Mean Square Error) is 4.6 mm, 5.3 mm, 8.5 mm, 10.5 mm, 6.2 mm, 2.0 mm for A001, A002, K002, K003, J004, J005, respectively. The comparisons between ASAR and TerraSAR-X results at the common CGPS station A001 (Table 4 and Figure 12a) indicate similar deformation time series and velocities from the two different SAR band satellites. We also compared the deformation results between ASAR Frame 23 and Frame 24 at the common CGPS stations A002 and J005 in Figure 12b,f, and very good agreements were found in both displacement rates and time series. The surface uplift in the Xizhang region was also validated by the CGPS station A001, as shown in Figure 12a, with a slight uplift rate of 1.4 mm/year.

6.2. Ground Deformation and Gravity Changes

Gravity measurements (spaceborne, airborne or ground-based) can provide useful information for understanding the mechanism of the ground movements because they reflect the underground density variations or mass movements. Observing gravity changes is one of the most powerful methods for detecting mass redistribution associated with plate tectonic motions, such as those resulting in earthquakes and volcanism [37,38]. At a given gravity station, gravity change results from the effects of free-air gravity gradient and mass redistribution. By carefully modeling the environmental gravity variations and the admittance factor between gravity and height, gravity change can be converted to height change (uplift or subsidence) [39]. The declining of the groundwater level causes land subsidence in Taiyuan basin and thus results in gravity changes. The availability of gravity values shown in Table 3 provides useful information about the mechanism of the land subsidence in our study area. For the case of the Taiyuan basin, we assumed that the gravity variations are largely due to groundwater extraction and land subsidence [33]. The relative gravity measurements were used to obtain the gravity variation rate of $(4.3\pm 0.5)\times {10}^{-8}$ ms${}^{-2}$/year in the period from 2001 to 2006 [33]. The gravity changes were then converted to the land subsidence with a subsidence rate of 21.8 mm/year assuming the gravity values does not include the groundwater change [33]. Absolute gravity measurements in Jinci shows that the gravity values increased at a rate of $1.3\times {10}^{-8}$ ms${}^{-2}$/year in the period 2009–2012, and equivalent to a subsidence rate of 6.6 mm/year [33]. The land subsidence rate during 2001–2006 inferred from gravity change is much larger than that after 2009, due mainly to the increased of the groundwater level because of the implementation Shanxi Fenhe river water re-flow project.

We analyzed the time series data of gravity changes and ground deformation for stations G1–G10 in the Taiyuan basin, as presented in Figure 13, Figure 14 and Figure 15. The displacement time series was obtained from ASAR Frame 24 and was converted into the vertical direction, considering that the horizontal component of displacement was negligible. The gravity for G1 station in Figure 13 decreased linearly in time at a rate of $-2.29\times {10}^{-8}$ ms${}^{-2}$/year, while the InSAR measurement show a land uplift rate of 11.22 mm/year. The gravity values in G2 station underwent large oscillations, most likely due to the local hydrological effects; however, the displacement time series increased in a near-perfect linear behavior. As can be seen in Figure 14, the G3 station is located close to the largest subsidence area in the Taiyuan basin and close to Pingyao fault. The gravity values of G3 increased linearly at a large rate of $11.47\times {10}^{-8}$ ms${}^{-2}$/year, while the height of ground decreased at a linear rate of 5.34 mm/year. Both the gravity values and ground deformation time series for the G4 site underwent large oscillations, most likely caused by the natural or artificial groundwater recharging in this area. The G5, G6, G7 and G8 stations are outside subsidence areas, and the gravity and height change for these stations also increased in time and with some fluctuations. In the south of the basin, two land subsidence centers (Xiaoyi and Jiexiu county) among all ten gravity stations were caused by coal mining activities, shown in Figure 15. The G9 station is located in Jiexiu county, with a maximum gravity change rate of $15.81\times {10}^{-8}$ ms${}^{-2}$/year and a maximum subsidence rate of 11.7 mm/year among these ten gravity stations. Located in Xiaoyi city but outside the subsidence area, the G10 site shows a significant increase in gravity values despite a small ground uplift. Many factors, such as tectonic movement, groundwater level change, earth’s mass variation, ground deformation, can cause ground gravity change. We cannot extract the influence of other factors from ground deformation because there is no groundwater level data for this area. In general, however, the gravity values at each station increase with time, showing traces of ground subsidence.

6.3. Correlation between Deformation and Pre-Existing Faults

One of the most significant observations from the deformation velocity maps in Figure 14 is the identification of differential subsidence across the faults. The shapes of subsidence zone near the faults clearly followed the general trends of the mapped faults. In order to investigate this phenomenon, we plotted the displacement rates along two profiles A-A’ and B-B’ (the locations shown in Figure 14) across the three pre-existing faults (Jiaocheng fault, Pingyao fault, and Taigu fault), as displayed in Figure 16. High velocity gradients of ground deformation were found across the faults, suggesting that subsidence in these locations is limited by the mapped faults. In the Qingxu and Jiaocheng county locations, we noticed a significant differential subsidence (a velocity gradient of 40 mm/year) across the faults on narrow zones (about 3 km wide) (as seen in Figure 14 and Figure 16). Traditionally, the lateral sediment faces is gradually changed and so too the ground deformation associated with their consolidation. Thus, the observed abrupt changes in ground deformation suggests the presence of mapped or unmapped faults. One problem caused by abrupt changes of land subsidence is the ground fissures that have been observed close to these two areas (Qingxu and Jiaocheng county) [40,41]. Researchers also found that fissures are spatially associated with land subsidence and pre-existing faults within the whole Taiyuan basin. The tectonic motion in the Taiyuan basin is less than 2 mm/year [11], as compared to the subsidence rate of 55 mm/year in the Qingxu and Jiaocheng subsidence center, suggesting that land subsidence is independent of tectonic movement. Because of the observed high differential deformation rates (several centimeters per year) and the small spatial range of deformation relative to the length of the fault, we believe that the movement of the fault is driven by water extraction rather than regional tectonic activity. Similar studies were made by [42] who suggested that faults act as barriers or conduits to fluid flow. In the Taiyuan basin, faults coincide with limits of the subsidence areas and thus likely serve as barriers to groundwater flow, which creates groundwater level difference across the faults. The development of rock material along fault planes with different grain size and the different characteristics of the lithologic units juxtaposed by faulting likely explain this phenomenon [43].

6.4. Migration of Subsidence Areas to Newly Urbanized Areas

Land subsidence has appeared in Taiyuan city since the 1960s, with five subsidence centers identified, including: Xizhang, Wanbailin, Xiayuan, Wujiabao, and Xiaodian, as shown in Figure 2 and Figure 9. Wanbailin, Xiayuan and Wujiabao are the downtown areas of Taiyuan city with the densest population and the highest concentration of industries. The land subsidence usually started from the downtown areas [44].

Table 5. Comparison of subsidence rate (mm/year) in different subsidence centers between different stages in Taiyuan city. The subsidence rates in period 1956–2000 are from [10,11,12].

Subsidence Centers	1956–1981	1981–1989	1989–2000	2003–2008	2009–2010
Xizhang	−3.0	−43.4	−25.0	+2.3	+2.5
Wanbailin	−2.5	−46.8	−46.7	−20.5	−12.6
Xiayuan	−3.6	−55.3	−86.0	−22.5	−15.4
Wujiabao	−33.1	−114.0	−96.2	−28.9	−30.9
Xiaodian	no subsidence	no subsidence	no subsidence	−27.9	−50.8

compares the subsidence rates in these different subsiding centers between different time periods. The ground deformation started in Wujiabao during 1956–1981 with a subsidence rate of 33 mm/year while the others (except Xiaodian) subsided slowly at only several mm per year. The subsidence became much more serious during 1981–1989 for all regions (except Xiaodian). It is believed that the rapid development of urbanization and industrialization during this period led to the accelerating subsidence rates. The land subsidence in Xizhang, Wanbailin, and Wujiabao came to be alleviated during 1989–2000, but, in Xiayuan, the situation still accelerated. After 2003, all subsidence centers were basically controlled because the local government did substantial preventive and remedial measurements to control subsidence. Even a ground rebound was detected in the Xizhang region (in the north of Taiyuan city) with an uplift rate of about 2 mm/year, possibly due to the replenishment of groundwater.

However, the Xiaodian district, in contrast, had no subsidence before 2000 and it began to subside in 2003 and the rate grew rapidly during the study period. As shown in Table 5, the subsidence velocity was only 27.9 mm/year during 2003–2008, but it increased to 50.8 mm/year during 2009–2010. In 2002, the government of Taiyuan city started to shift the economic center to the south, building a new economic development zone in the Xiaodian district. Figure 17 illustrates the land subsidence and types of land use from Google Earth optical images in the Xiaodian district during our experimental period. Before 2003, the land in Xiaodian was mainly for agricultural purposes (see Figure 17c), so the ground remained stable from January 2004 to July 2005 (see the displacement time series before the vertical yellow line). Some housing constructions and industrial facilities started to appear in 2006 (Figure 17d). This was accompanied with a large amount of groundwater pumping and thus triggered the land subsidence in July 2005. After July 2005, land subsidence developed rapidly due to more and more buildings (Figure 17e) and depletion of groundwater in this area. We fitted a linear regression line (the red dashed line in Figure 17b) to the displacement time series after July 2005 and the corresponding velocity is up to 47 mm/year. Thus, the subsidence in Taiyuan city has migrated from existing downtown areas to surrounding newly urbanized districts such as in Xiaodian. This situation is expected to continue and increase because the industrialization and urbanization will continue for decades.

7. Conclusions

Applying time series InSAR methods to ENVISAT ASAR and TerraSAR-X images, we detected ground subsidence in the whole Taiyuan basin not only in one single city as previous studies. Systematic and seasonal atmospheric effects were identified in the high mountain areas surrounding the basin, which could cause errors to InSAR displacement results. An atmosphere-corrected time series InSAR was developed to mitigate the atmospheric delay more effectively than the conventional spatial-temporal filtering. The atmospheric correction results were validated by MERIS measurements and deformation results were validated by continuous GPS observations. We found eight subsidence centers in total in the basin and a surface uplift in the north of Taiyuan city.

The analysis of the InSAR displacement time series versus the gravity changes suggested that the land subsidence was caused by soil compaction and/or groundwater level changes, which are connected with excessive groundwater extraction. Significant velocity gradients were observed between several pre-existing faults that constraint subsidence bowls. The correlation of subsidence to the locations of known active faults allows us to understand better the role of regional tectonics in land subsidence. In recent years, the implementation of groundwater withdrawal control in downtown areas (Wanbailin, Xiayuan, and Wujiabao) of Taiyuan city has significantly reduced the subsidence velocity. The most severe subsidence center in Taiyuan city has migrated from the downtown areas to a newly economic development zone—the Xiaodian district. As the rapid urbanization in those new urbanized areas continues, the prevention and control of land subsidence will continue to be a challenging task for sustainable development.