Land deformation associated with exploitation of groundwater in Changzhou City measured by COSMO-SkyMed and Sentinel-1A SAR data

2.1 Study area

Changzhou City is the study area located on the deltaic deposits of the Yangtze River Delta on the east coast of China (Figure 1 (a) Its geographic location is 31°09’N-32°04’N and 119°08’E-120°12’E, with Taihu Lake to the south and Wuxi city to the east. The urban region of Changzhou contains Xinbei, Zhonglou, Tianning, and Wujin districts. (Figure 1(b)) The ground elevation is less than 10 m in the plain area, higher in northwest, lower in the middle and southeast. There are several bedrock hills spotted in south and west part of this city, and the groundwater elevation is about 80 to 178 m [18]. Sedimentary and igneous rocks constitute the basement of Changzhou City. This city experienced seawater intrusions four times and deposited a set of Quaternary sediments since middle Pleistocene. The thickness of Quaternary deposits strata varies, which affected by the bedrock [19]. There are four aquifers in the Quaternary alluviums, namely, one unconfined and three confined aquifers, which are denoted as U, I, II and III (Figure. 1(c)). Conventionally, groundwater with in 50 m deep is shallow, i.e., U and I, otherwise is deep, i.e., II and III [20].

Changzhou City is within a temperate climate zone of abundant rainfall and dense river network. Over the past several decades, its urban and rural areas have witnessed phenomenal economic growth, which has been accompanied by severe pollution of the surface water and general decline of the water quality of rivers. With the rapid economic development and the pollution of the surface water in Changzhou City, the groundwater usage (both civilian and industrial water) increased gradually. In our study, the second confined aquifer is the main aquifer.

Figure 1

(a) Location of the study area. (b) Red polygon denotes area under the SAR images including the major urban areas of Changzhou: Xinbei, Zhonglou, Tianning, and northeast portion of Wujin. (c) The NW-SE oriented line segment indicates the location of the geological cross section of the Changzhou City

2.2 SAR Data

All the SAR data information used in this study is show in Table 1. COSMO-SkyMed is a high-resolution radar satellite constellation that developed by the Italian space agency and defense department, and it is composed of four synthetic aperture radar satellites. COSMO-SkyMed acquires X-band images (wavelength of 3.1 cm), and is more sensitive to urban land deformation detection than the C-band (5.6 cm wavelength) and L-band (23.5 cm wavelength) radars [10]. The COSMO-SkyMed image data feature high spatial resolution, wide coverage, selectable angles of incidence and multipolarity; the image features used in conjunction with the precise orbital data allow for repeat-track interferometry in the monitoring of surface deformations. In this study, the 12 single-look complex slant (SCS) (L1A) data was acquired from the X-band SAR1 and SAR3 satellites that cover the study area. The five SAR3 images were between July and December of 2011, whereas the seven SAR1 images were between January and July of 2012. All images for this study were of HH polarization and from the ascending orbit. The look angle at the image center was about 26°, and the imaging mode was Stripmap. The Stripmap mode has two sub-modes: HIMAGE and Ping-Pong. All images for this study were of HIMAGE mode. In the HIMAGE mode, the configurations of the radar transmitting and receiving antenna do not change over time, satellites can receive the full Doppler beam width from every ground scatter, with a pattern width approximately 40 km. Currently, the X-band COSMO-SkyMed SAR data have been successfully applied to the detection of surface deformations with sub-centimeter accuracy [7], [16]. Sentinel-1A is a European C-band radar satellite launched in 3 April 2014 with a 12-day repeat cycle. The data are freely provided by the European Space Agency (ESA) (https://scihub.copernicus.eu/dhus/). It is a C-band SAR satellite. This satellite has four modes: Interferometric Wide Swath (IW), Extra Wide Swath (EW), Wave Mode (WV) and Stripmap (SM). In this study, Sentinel-1A IW Level (L1) products are used. Twenty-six scenes of Sentinel-1A IW mode (VV) data from the ascending orbit are observed and available from May 2, 2015 to April 9, 2017.

2.3 Leveling measurements

Four second-order leveling stations (labeled as L1-L4) (Figure 2) located in Tianning, Xinbei and Wujin are used to monitor land subsidence in Changhzou, with the observation once every two years. The leveling survey data is from September 2006 to September 2014. Leveling data (from 2011 to 2012) are used to validate the InSAR results.

2.4 Groundwater level

In Changzhou, the main layer from which groundwater extracted is the second confined aquifer. To explore the relationship between the variation of the water-table and the land deformation, groundwater-table data of the second confined aquifer, obtained from six water wells (labelled as W1 to W6 in Figure 2), were used in this study. These data were provided by the Jiangsu Province Hydrology and Water Resources Investigation Bureau.

2.5 SRTM DEM

The Shuttle Radar Topography Mission (SRTM) is a NASA mission conducted in 2000 to obtain elevation data for most of the world. It is the current dataset of choice for digital elevation model data (DEM) since it has a fairly high resolution (about 90 meters) and near-global coverage (from 56° S to 60° N). SRTM DEM has been publicly released in 2003, and revised many times. Data used for this study is the version 4.1 [21]. SRTM DEM was used to remove the topographic phase of interferometric phase.

Figure 2

The location of leveling stations, deformation points and groundwater monitoring stations distributed in the study area.

3.1 SBAS-InSAR principle

SBAS-InSAR is a newly developed DInSAR-based time-series analysis method first proposed by Berardino et al. [12]. Over the past decade, SBAS algorithm has been improved constantly and widely used in surface deformation studies [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,]. In the analysis, all SAR images are subject to random combinations and form interference pairs.

Supposing there are SAR images with N+1 scenes covering the same area with acquisition times of t₀, t₁ ... t_N, and an image from each scene can form an interference pair with at least one image of another scene. The SAR images whose vertical baselines are below the threshold of small baseline condition are grouped together to generate a total of L candidate images and M differential interferograms. M satisfies the following inequality (assuming N is an odd number):

If all of the differential interferograms are correctly unwrapped, they can be calibrated by a highly coherent point (x₀, r₀) that is stable or whose deformation magnitude is known. For a specific interferogram i.e., after the removal of the flat earth and topographic phases, and assume t_B> t _A (t_A, t_B are the two acquisition dates of the co-registered image pair for the interferogram i), then in a coordinate system of azimuth-range (x, r), the interferometric phase of a single-look slow-varying filtered phase (SFP) pixel (x, r) in the differential interferogram i can be expressed as follows:

where λ is the center wavelength of the imaging radar system, d(t_B, x, r)and d(t_A, x, r) represent the relative displacement in the direction of the radar line of sight (LOS) at times t_B and t _A, respectively. B⊥iΔZrsin⁡θis the phase resulted from the DEM (or height) error ΔZ and the vertical component B ^⊥_i of baseline orbit separation between the SAR image pair, and its value is directly proportional to B ^⊥_i and inversely proportional to slant range r, and it is associated with local incidence angle θ. [ϕ_atm(t_B, x, r)−ϕ_atm(t_A, x, r)] is the difference of the atmospheric delay phase of pixel (x, r) at times t_A and t_B, and Δϕ represents other noises.

After phase unwrapping, the linear deformation and DEM error phase can be used to construct a new SBAS linear equation:

where B is an M × N matrix, and C is the coefficient matrix related to the space-baseline distance. p⁰ is given by

where v represents the phase-value matrix of the average deformation rate.

The linear deformation phase and DEM errors are removed from the interferogram to generate the residual phases, which now include the atmospheric delay phase, non–linear deformation phase and decoherence noises. Subsequently, filtering is performed based on the different temporal and spatial properties of the residual phases, thereby separating the non–linear deformation phase and atmospheric delay phase. Finally, time series deformation results can be obtained by applying a least square (LS) method [22] or singular value decomposition (SVD) rule to all of the unwrapped interferograms.

3.2 SBAS-InSAR processing

Land displacement in Changzhou was mapped with the ENVI SARscape software developed by ESRI company. This software provided SBAS-InSAR processing method. The SBAS-InSAR is used to divide the generated SAR data into several groups to suppres the effect of spatial-temporal de-coherence. In the SBAS InSAR approach, thresholds of spatial-temporal baseline and Doppler centroid difference were applied to generate interferograms with SAR image pairs that maximize the InSAR coherence. In this study, only the spatial-temporal baseline was considered since the Doppler centroid difference was negligible. Spatial baseline smaller than 800 m and 300 m, the temporal interval smaller than 300 days and 500 days for images in COSMO-SkyMED and Sentinel-1A respectively. Under such design, 37 interferograms in the X-band (Figure 3(a)), and 67 interferograms in the C-band (Figure 3(b)) were generated (Figure 3(c))

Each interferometry pair that meets the temporal and spatial criteria was processed to form a DInSAR interferogram using the two-tracked method, which mainly consists of image registration and baseline estimation, interferogram generation, coherence coefficient calculation, external DEM simulation, precise baseline estimation based on precise orbits and control point data, interferogram adaptive filtering calculation [24], phase unwrapping and refinement and reflattening. To dampen the noises, boost coherence, and examine the deformation at a relatively large scale, interferograms were subject to a multi-look treatment using 11 azimuth looks and 10 range looks for COSMO-SkyMED and 1 azimuth looks and 5 range looks for Sentinel-1A. The minimum cost flow algorithm is applied for phase unwrapping [25].

Figure 3

Distribution of small-baseline pairs. (a) COSMO-SKeMED. (b) Sentinel-1A

The unwrapped phases of 16 differential interferometry pairs were processed as follows. First, the differential interferometric phase was used to estimate the linear deformation phase of the study area. Then the linear deformation phase was transformed to a LOS directional deformation rate. Subsequently, the linear deformation phase and residual topographic phase were removed from the original differential interferometric phase before the residual phase was phase unwrapped and the SVD (singular value decomposition) was used to solve the non–linear deformation phase, which culminated in the generation of a cumulative deformation time-series phase. To estimate the atmospheric delay phase, the data were subjected to high-pass (HP) filtering in the time domain followed by low-pass (LP) filtering in the space domain; subsequently the atmospheric delay phase were separated from the data.

4.1 Time evolution of land deformation

Annual land deformation maps are produced from ENVI SARscape processing. Figure 4 shows the deformation calculate with COSMO-SkyMED data (Figure 4(a)) and Sentinel-1A data (Figure 4(b)) The deformation rate from 2011 to 2012 calculated with COSMO-SkyMED data was – 55 ~ +30 mm/year. The deformation rate from 2015 to 2017 calculated with Sentinel-1A data was – 40 ~ +36 mm/year. Three subsidence zones concentrate on Xinbei, Tianning and Wujin districts marked with A, B and C respectively. There are detected in COSMO-SkyMED data. With the implementation of control measures for land subsidence by the Government, the subsidence zones of these three zones have been controlled effectively and zone C in Wujin District presents land rebound obviously (Figure 4(b))

Figure 4

Deformation rates projected onto the vertical direction for COSMO-SkyMED and Sentinel-1A data at the selected study area. A, B and C represent the regions with serious land subsidence (a) COSMO-SkyMED measurements, (b) Sentinel-1A measurements.

4.2 InSAR accuracy assessment

Land deformation calculated from SBAS-InSAR method is along LOS (line of sight) direction, but the leveling measurement is along the vertical. In order to insure the comparability of InSAR measurements and leveling measurement, the InSAR measurements (include COSMO-SkyMED and Sentinel-1A result) were projected to the vertical direction. In this study, leveling measurements from 2011 to 2012 were used to assess InSAR measurements accuracy. To compare InSAR measurements with those from leveling benchmarks, SFPs that lie within 100 mpf benchmarks are selected and then average displacement values around benchmarks. Through the above operation, leveling and InSAR measurements are in the same spatial frames. The comparison results indicated great agreement between InSAR and leveling measurements. With a minimum and maximum difference of 1.2 mm/year and 3.0 mm/year, respectively (Table 2). Measurement of land deformation based on individual interferograms can often be affected by artifacts due to temporal change in atmospheric delay and error. In satellite orbit and topographic data used, all of these factors that tend to affect the In-SAR accuracy. With SBAS-InSAR method described in section 3, atmosphere, baseline and other artifacts can be reduced effectively.

4.3 land deformation in Changzhou from 2011 to 2017

As the data acquisition time of COSMO-SkyMED is from July 21th 2011 to July 14th 2012 and Sentinel-1A data is from April 8^th 2015 to January 27^th 2017, there is a data gap between them (from July 14^th 2012 to April 8^th 2015). To fill the gap, average deformation rates during July 21^th 2011 to July 14^th 2012 are used to connect two types of land deformation of the time serious.

Figure 5

Land deformation of P1, P2, P3 and P4.

In order to better analyze the development of deformation in land subsidence areas, four land subsidence points marked with P1, P2, P3 and P4 are selected (Figure 5). The annual land subsidence rate of P1 to P4 in Figure 4 was −2.48 mm/year, −12.78 mm/year, −18.09 mm/year, and −12.69 mm/year respectively, and the accumulative land subsidence value of P1 to P4 in Figure 5 was −14.54 mm, −71.31 mm, −105.92 mm, and −74.27 mm respectively (Figure 5). The variation trends at the four points may suggest that the subsidence in these areas may continue for a while. An investigation of land subsidence in the Changzhou city by Hu indicated that over the exploitation of groundwater for industrial purposes was the main cause of land subsidence [26].

5.1 The historical land subsidence in Changzhou

According to the historical leveling survey data, land subsidence in Changzhou started in the year of 1970 and reached 50.63 mm/year between 1979 and 1983, with the maximal subsidence up to 100 mm. Between 1984 and 1991, there was a subsidence rate of 40–50 mm/year that led to a maximal cumulative subsidence of 949 mm in some of the central urban areas. The average subsidence rate was about 40 mm/year between 1993 and 1998 [26]. Since 1996, the government started limiting groundwater withdrawal [27]. Land subsidence in Changzhou showed a slowdown since 2000, at which time a subsidence rate of approximately 25 mm/year was observed. After 2004, except for the eastern part of Changzhou, land subsidence had been effectively curbed. A slight rebound occurred in some areas, and the rebound area gradually expanded. By 2007, urban land subsidence slowed down and many regions of urban areas start to rebound [26].

From July 21^th 2011 to July 14^th 2012 in this study (Figure 4(a)), land subsidence still occurred in parts of Xinbei District and Wujin District. The average subsidence rate was −7 mm/year at subsidence area of Xinbei District (e.g. the region A in the Figure 4(a)) and −15 mm/year at subsidence area of Wujin District (e.g. the region B in the Figure 4(a)). There is also a settlement phenomenon in the east part of Tianning (e.g. the region C in the Figure 4(a)), and the average subsidence rate was −5 mm/year. In addition, there was a slight rebound at an uplift rate of 4 mm/year in central urban districts (e.g., Zhonglou District and the west part of Tianning District).

From April 8^th 2015 to January 27^th 2017 results (Figure 3(b)), land subsidence in Xinbei, Tianning, Zhonglou and Wujin district has been controlled efficaciously. The average land subsidence rate of A, B and C is 0 mm/year, −1 mm/year and +6 mm/year. Although land subsidence in A, B and C region slow down, the future attention is also required to monitor land subsidence over this area.

5.2 The relationship between land deformation and groundwater exploitation

Figure 6 shows a sketch of the vertical cross-section of confined aquifers in Changzhou, which had abundant confined groundwater resources [27, 28, 29]. In the urban areas of Changzhou, the second deep confined aquifer was the main layer of groundwater extraction and the volume of groundwater exploitation has reached a considerable scale since the 1980s [27]. The mining of deep confined aquifer was at it’s minimum before 1960, but it was significantly enhanced during 1960-1970 when the social and economic development began to accelerate and led to insufficient urban water supply. Between 1970 and 1994, large water-consumption industries (e.g., printing, dyeing and textile) were developed swiftly, which caused massive exploration of groundwater resource in the surrounding towns and villages. Such exploration gradually surpassed the mining quantity of central urban districts. Thus, Changzhou experienced a drastic groundwater withdrawal during this period. After 1994, the local government realized the potential environmental damage due to excessive exploitation of groundwater, and formulated a series of rules and regulations to restrict groundwater extraction [26]. Miao et al. discovered that in Changzhou, changes in the confined water level II correlated well with the land subsidence trend. Excessive withdrawal of groundwater was the most critical factor leading to large-scale land subsidence in the city [27].

Since 2000, exploitation of groundwater has been prohibited in Changzhou city and measures were taken to recharge the depleted areas. Since 2005, a widespread increase of groundwater level has been recorded [26]. Since then, subsidence slowed down in Zhonglou and Tianning districts, alongside a slight land rebound was even observed [18, 19, 20, 21, 22, 23, 24, 25, 26, 27,]. From 2002 to 2009, Wujin had experienced a rapid industrial development that has led to a tremendous need of groundwater, lowering the groundwater table and deteriorating the land subsidence [28].

Figure 7 shows the relationship between land deformations with groundwater levels of the confined aquifer exploitation wells, (W1-W5) (Figure 7) from July 2011 to July 2012. The wells (W1-W5) were applied to exploit confined water for industrial production and the groundwater level of these wells was recorded. As the groundwater level data only covered the data acquisition time of COSMO-SkyMED, the analysis based on COSMO-SkyMED data. There is a strong correlation between the land subsidence and groundwater depth (Figure 7). The land subsidence increases with the increase of groundwater depth (i.e. decrease of groundwater level). With the groundwater rebound 2.39 m, 1.07 m, −0.51 m, 4.37 m, 4.33m and 0.9 m, the land rebound 2.95 mm, 2.95 mm, −1.79 mm, 9.14 mm, 0.92 mm and 2.60 mm in W1, W2, W3, W4, W5 and W6 respectively.

Figure 6

Cumulative displacements at locations P1-P4 (Figure 5) in the subsidence zone. The blue lines represent the water level.