Summary

The Pearl River Delta (PRD) is one of the most important economic regions with the highest population densities in China. With its dramatic increasing population and economy, hazards associated with land subsidence frequently occur here that amplify the negative effect of sea level rise. However, land subsidence has not been regularly measured in this region. Here, we use interferometric synthetic aperture radar (InSAR) to investigate the rate and extent of land subsidence in the PRD region. Assuming purely vertical displacements, multi-track interferograms from different viewing geometries are combined to estimate the linear rate map and time series at a higher resolution in time than is possible with a single track. The results show apparent subsidence along the coastal region of Shenzhen associated with rapid urban development in recent years. The average subsidence rate within 500 m of the coast is about 2.5 mm yr−1, and the maximum is up to about 6 mm yr−1 with respect to the central part of the city. Much of the land surface in the PRD is less than 2 m above mean sea level; high-precision geodetic measurements throughout the PRD region are therefore critical for conducting risk assessments, considering the rate of about 2–3 mm yr−1 of current global sea level rise.

1 Introduction

The global mean sea level is estimated to be rising at a rate of about 2–3 mm yr−1, with a significant acceleration since the early 1990s when observations from satellite altimetry are available (Cazenave & Llovel 2010; Church & White 2011; Church et al. 2011). Meanwhile, many large deltas in coastal regions undergo rapid land subsidence with the population and economic growth, which results in problems such as compaction of reclaimed land and deltaic deposits by high-rising infrastructures, over-pumping of the groundwater system (Syvitski et al. 2009). Land subsidence amplifies local relative sea level rise and puts the affected population under a greater risk of environmental problems such as flooding, shoreline erosion and saltwater intrusion. Accurate measurements of land subsidence are therefore critical for risk assessments on sea level rise in coastal regions.

The Pearl River Delta (PRD), also known as Zhujiang Delta, is the low-lying area surrounding the Pearl River estuary where the river flows into the South China Sea (Fig. 1). It is one of the most densely populated areas in the world with a population of about 56 million in 24 500 km2. About 3,720 km2 is less than 2 m above sea level. It has been the fastest growing economic region in China since the launch of China's reform program in the late 1970s. Unlike the Yangtze and Yellow river deltas that have undergone rapid regional-scale subsidence due to groundwater pumping and mining, land subsidence in the PRD was thought to be slower and mainly caused by sediment compaction and concentrated urban development. In spite of the fact that land subsidence has caused considerable damage here in recent years, the rate and extent of land subsidence are still unknown here.

Figure 1

(a) Location map of the PRD. The blue rectangles show extents of Envisat SAR imagery on one descending (T175) and two ascending (T297/T025) tracks. The red box indicates our study area. (b) Elevation from SRTM DEM superimposed on the average amplitude of 54 SAR images acquired between 2007 and 2010. The yellow lines are the city boundaries between Dongguan and Shenzhen, two megacities in the PRD.

Figure 1

(a) Location map of the PRD. The blue rectangles show extents of Envisat SAR imagery on one descending (T175) and two ascending (T297/T025) tracks. The red box indicates our study area. (b) Elevation from SRTM DEM superimposed on the average amplitude of 54 SAR images acquired between 2007 and 2010. The yellow lines are the city boundaries between Dongguan and Shenzhen, two megacities in the PRD.

Interferometric synthetic aperture radar (InSAR) measures surface displacements by differencing the phase observations between two complex radar images. Using multiple images and advanced analysis methods, InSAR is capable of remotely mapping time-varying surface deformation with mm-level accuracy and metre-level spatial resolution over large areas (see reviews in Simons & Rosen 2007; Sansosti et al. 2010; Hooper et al. 2012). It is therefore the best method to investigate regional-scale land subsidence when ground benchmarks (such as GPS and levelling) are not available. In recent years, InSAR has widely been used to measure subsidence associated with groundwater pumping (Hoffmann et al. 2001; Schmidt & Bürgmann 2003; Bell et al. 2008; Aly et al. 2009; González & Fernández 2011), tectonic activities (Lanari et al. 2004; Argus et al. 2005; Wisely & Schmidt 2010), sediment compaction (Mazzotti et al. 2009; Teatini et al. 2011) and underground mining, oil and gas extraction (Gourmelen et al. 2007; Jiang et al. 2011).

Most InSAR studies measure line-of-sight (LOS) deformation using radar images on a single track. If the target area is observed with different viewing geometries, two or three components of displacements can be inverted (Wright et al. 2004; Gourmelen et al. 2007; Wang et al. 2007). Ozawa & Ueda (2011) developed a method to estimate two components (quasi-east-west and quasi-up-down in the literature) of deformation time series by averaging multi-track data into regular time windows (46 days used in their study). This method can estimate multi-dimensional deformation time series, but it reduces temporal resolution of the original data. For most subsidence studies, the deformation is primarily in the vertical direction and the horizontal is usually negligible (e.g. Bell et al. 2008). If so, InSAR LOS observations can be projected into the vertical using local incidence angles while neglecting the horizontal components. The vertical displacements within the overlapping region can then be combined straightforwardly to estimate deformation time series at a higher resolution in time than is possible with a single track. This method can combine interferograms with different viewing geometries from multiple satellites and can significantly improve the accuracy and temporal resolution of InSAR results. The PRD region lies in the rigid South China block, where the horizontal strain rate is negligible as measured by GPS (Calais et al. 2006). Therefore, this method is suitable to study land subsidence in the PRD region.

2 Method

For a given pixel in the geocoded interferogram, the LOS displacement dlos (positive away from satellite) can be written as  

1
formula
where de, dn and du are displacements in east, north and up directions, respectively, α is the azimuth of the satellite heading (positive clockwise from the north) and ϑ is the radar incidence angle. For Envisat/ASAR data in image (IM) mode, the incidence angle varies gradually from 19° at near range to 26° at far range.

The horizontal components are often negligible for subsidence monitoring. Assuming this, the vertical component du can then be converted from dlos using the local incidence angle,  

2
formula

If the interferograms are geocoded with the same pixel spacing, du can be combined for displacement rates and time-series analysis no matter what tracks or satellites are used. Following Berardino et al. (2002) and Schmidt & Bürgmann (2003), the system of linear equations can be written as  

3
formula
where d is a vector containing the vertical displacements du.

For the time-series inversion, G and m are constructed from the time-spans and the velocities between consecutive acquisition dates. Given an interferogram Ii, j formed by two images i and j, the associated design matrix Gi,j and the model vector m can be expressed as  

4
formula
 
5
formula
where Δti =ti + 1ti, t is the acquisition date, n is the total number of the acquisitions, 0 is a zero vector indicating acquisitions which are not covered by the interferogram Ii, j, vi is the velocity of the ith time-span. Here, the acquisitions must be chronologically ordered for all the tracks (Fig. 2).

Figure 2

Schematic showing InSAR time series inversion within the overlapping region of multiple tracks. The circled numbers denote radar acquisitions chronologically ordered for all tracks.

Figure 2

Schematic showing InSAR time series inversion within the overlapping region of multiple tracks. The circled numbers denote radar acquisitions chronologically ordered for all tracks.

Because the subsets cannot be connected via cross-track interferograms, the above system of eq. (3) is rank-deficient and has infinite solutions for multi-track data sets. An optimal solution can be obtained using the singular value decomposition (SVD) method (Berardino et al. 2002) or Laplacian smoothing operator (Schmidt & Bürgmann 2003).

For the linear displacement rates inversion, G and m can be expressed as (6) and (7), respectively. The mean velocity graphic can be solved using weighted least-squares method.  

6
formula
 
7
formula

The general equation (3) can be extended to estimate DEM errors as well as linear rates or time series if the perpendicular baseline is available for each pixel. In this case, the extended matrices G′ and m′ can be written as  

8
formula
 
9
formula
 
10
formula
where Δh is the DEM error, B is the perpendicular baseline, ρ is the distance between satellite and ground scatterer and m is the number of the interferograms used in the inversion.

3 Insar Processing

Here, we use data acquired between 2007 and 2010 on three tracks by the European Space Agency's (ESA) Envisat/ASAR instrument operating at C-band with a wavelength of 5.6 cm. We select images to form interferograms according to their perpendicular and temporal baselines (Fig. 3). Envisat SAR images were only occasionally acquired in the PRD region before 2007. Most interferograms formed from these images are incoherent due to the large temporal or spatial baselines.

Figure 3

Temporal and perpendicular baselines for the interferograms used in this study on the ascending tracks T297, T025 and descending track T175. The circles denote SAR acquisitions among which 18 acquisitions are from T297, 17 from T175 and 19 from T025. The lines indicate interferometric pairs including 26 from T297, 22 from T175 and 29 from T025. The solid lines are those selected using the MST algorithm. The colour of the lines shows the fraction of successfully unwrapped pixels w.r.t. the image size.

Figure 3

Temporal and perpendicular baselines for the interferograms used in this study on the ascending tracks T297, T025 and descending track T175. The circles denote SAR acquisitions among which 18 acquisitions are from T297, 17 from T175 and 19 from T025. The lines indicate interferometric pairs including 26 from T297, 22 from T175 and 29 from T025. The solid lines are those selected using the MST algorithm. The colour of the lines shows the fraction of successfully unwrapped pixels w.r.t. the image size.

Much land surface in the PRD region is covered by dense vegetation so that many interferograms, particularly those formed with summer acquisitions, are significantly decorrelated for a perpendicular baseline more than 150 m or duration more than 1 yr. Furthermore, quite a few summer acquisitions were strongly contaminated by tropospheric delay errors due to the humid subtropical climate and lengthy monsoon season (April–September) in the PRD. We eliminate the acquisitions with variance over 100 mm2 to minimise the atmospheric delay errors, leaving 54 images for the following analysis.

We produce interferograms for each track individually using the JPL/Caltech ROI_PAC software (Rosen et al. 2004). Topographic phase is removed using a 1-arc-second DEM interpolated from the original 3-arc-second SRTM DEM (about 90 m) (Farr et al. 2007). The 1-arc-second DEM is used here to geocode radar amplitude images—a by-product of InSAR analysis—at higher resolution (Fig. 1). The coherence of the interferograms is improved by multi-looking, typically to four range looks (about 80 m), and then smoothed using a power spectrum filter (Goldstein & Werner 1998). The interferograms are unwrapped using a branch-cut algorithm (Goldstein et al. 1988). Some isolated patches have to be unwrapped manually by setting bridges due to loss of coherence. The unwrapped interferograms are then checked to exclude phase unwrapping errors using a phase closure technique: (i) We use the minimum spanning tree algorithm (MST, Kruskal 1956) to create a chain connection of the non-redundant interferograms (solid lines in Fig. 3); (ii) A closure will be generated while adding any other redundant interferogram (dashed lines in Fig. 3)—we get the closures using the Dijkstra algorithm (Dijkstra 1959) by adding these interferograms one by one; (iii) For each pixel, phase unwrapping errors are identified if the absolute sum of the closure is larger than 2π. Because of the limitations of decorrelation and strong atmospheric noise, we can only generate a few redundant interferograms (Fig. 3). Finally, we obtain 77 pairs of interferograms geocoded with the SRTM DEM coordinate system with a pixel spacing of 1-arc-second.

4 Π-Rate: Rate Map and Time-Series Analysis

Other than the deformation signal of interest, the interferometric phase contains artefacts from imperfect orbital knowledge, atmospheric propagation delays and DEM errors. Different approaches have been developed to account for these artefacts based on multi-interferogram analysis (e.g. Ferretti et al. 2001; Berardino et al. 2002; Hooper et al. 2007; Biggs et al. 2007; Hooper 2008; Gourmelen et al. 2010; Hetland et al. 2012). Π-RATE (poly-interferogram rate and time series estimator, http://homepages.see.leeds.ac.uk/~earhw/software/pi-rate/index.html) is an open-source package of Matlab codes to estimate the linear displacement rate map, time series and their associated uncertainties using a set of geocoded interferograms based on the method of Biggs et al. (2007), with further improvements by Elliott et al. (2008), Wang et al. (2009) and Wang & Wright (2012). It has been successfully used to investigate slow surface movements such as interseismic deformation with mm-level accuracy (e.g. Wang et al. 2009). Here, we improve the algorithms in Π-RATE so that it can estimate the displacement rates and time series by combination of InSAR data from multiple tracks. In the following, we describe the main processing steps in Π-RATE.

4.1 Non-redundant observations selection

For an interferogram network with N acquisitions, we can produce N(N − 1)/2 interferograms at most among which N − 1 are non-redundant. Unlike the other small baseline subset (SBAS) methods (e.g. Berardino et al. 2002; Gourmelen et al. 2010), Π-RATE estimates the time series and linear rate map using non-redundant observations. For an interferogram closure I1, 2, I2, 3 and I1, 3, for example, two interferograms have provided all the information observed by the three acquisitions I1, I2 and I3. We can find the most reliable two according to the quality of the interferograms.

Here, we adapt the MST algorithm (Kruskal 1956) to select the non-redundant observations pixel by pixel according to the unwrapped fraction of the interferograms (Fig. 4). We use all the coherent pixels in the non-redundant interferograms. For the incoherent pixels in the non-redundant interferograms, we reapply the MST algorithm to make use of the corresponding coherent data in the redundant interferograms. This strategy to select partially coherent pixels increases the coverage of the rate map and is especially useful for low coherence regions.

Figure 4

Example of selected InSAR observations (green) using MST algorithm for an interferogram closure 091209-100217-100428. The percentage is the fraction of the successfully unwrapped pixels (Fig. 3). All coherent pixels in (a) and (c) are chosen as they have higher fraction of phase unwrapping. Some coherent pixels in (b) are used if they are incoherent in (a) or (c).

Figure 4

Example of selected InSAR observations (green) using MST algorithm for an interferogram closure 091209-100217-100428. The percentage is the fraction of the successfully unwrapped pixels (Fig. 3). All coherent pixels in (a) and (c) are chosen as they have higher fraction of phase unwrapping. Some coherent pixels in (b) are used if they are incoherent in (a) or (c).

4.2 Reference point and phase

In InSAR applications, a reference point is required to fix all the interferograms to the same reference frame by setting the phase of the reference point to zero. Therefore, the reference point should be stable and persistently coherent in all the interferograms. We test a group of tie points throughout the scenes, and select one with the maximum fraction of coherent pixels and smallest phase variance within a small search window. The reference point is used when estimating the orbital and topographically correlated atmospheric delay error models.

Usually, it is unreliable to calibrate the interferograms by setting the phase of the reference point (or the average of a window around it) to zero even if the point is stable. For example, a local atmospheric pattern above the reference point could produce an offset for the entire interferogram. Following Finnegan et al. (2008), we prefer to remove the reference phase by setting the average displacement of each interferogram as zero. This method can avoid a blind choice of a single stable point and the differences of atmospheric artefacts on individual interferograms. Such reference phase will make the resultant deformation as relative to the entire scene, but the choice is reasonable here as no image-wide deformation is expected in the PRD region.

4.3 Orbital errors correction

In InSAR processing, precise orbital information is used to remove interferometric phase associated with differences in the satellites position between passes, that is, the reference ellipsoid and topography. Wang et al. (2009) show a roughly normal distribution of the parameters for orbital corrections with ERS and Envisat data, implying comparable accuracies of the precise orbits for these satellites. The best available orbits for the ERS satellites have a quoted accuracy of 5–7 cm radially and 10–15 cm cross-track (Scharroo & Visser 1998). Therefore, residual orbital phase ramps remain in the interferograms. It is relatively easy to correct orbital errors for subsidence monitoring since image-wide deformation usually is not expected. Here, we estimate orbital errors using quadratic polynomial models based on a network approach following Biggs et al. (2007) (Figs 5b–c). Apart from orbital errors, the quadratic models will remove long-wavelength components of atmospheric noise.

Figure 5

Example of InSAR artefacts correction for the interferogram 091209-100428: (a) original unwrapped interferogram; (b) quadratic orbital errors; (c) interferogram after orbital error correction; (d) topographically correlated atmospheric delay errors; (e) APS estimation; (f) interferogram corrected using atmospheric delay estimation from (d) to (e). Note that different colour scales are used for the subplots. A threshold of 35 non-redundant coherent observations is used for raw time series reconstruction and APS estimation, so there are more masked (grey) patches in (e)–(f).

Figure 5

Example of InSAR artefacts correction for the interferogram 091209-100428: (a) original unwrapped interferogram; (b) quadratic orbital errors; (c) interferogram after orbital error correction; (d) topographically correlated atmospheric delay errors; (e) APS estimation; (f) interferogram corrected using atmospheric delay estimation from (d) to (e). Note that different colour scales are used for the subplots. A threshold of 35 non-redundant coherent observations is used for raw time series reconstruction and APS estimation, so there are more masked (grey) patches in (e)–(f).

4.4 Atmospheric delay errors correction

Atmospheric artefacts in InSAR are due to the difference in propagation delay through the ionosphere and troposphere between the two SAR acquisitions. The ionospheric delay is sensitive to the wavelength of the signal. It is important for L-band data, but usually is negligible for C-band data at scales less than ∼50 km (Hanssen 2001). Tropospheric delay can be a signal of significant magnitude and is often divided into a vertically stratified component and a turbulent mixing component in InSAR studies (Hanssen 2001). The stratified delay is a function of the topographic height differences throughout the scene. We estimate the stratified delay using a linear function of the SRTM DEM based on the network approach following Elliott et al. (2008). Unlike other studies (Elliott et al. 2008; Doin et al. 2009), our stratified contribution is very small in magnitude because of the flat topography in the PRD region (Fig. 5d).

The behaviour of turbulent mixing varies in both the temporal and spatial domains. Such artefacts (often referred to as the atmospheric phase screen; APS) can be estimated using a temporal-spatial filter when tropospheric delay estimates from other techniques are unavailable (Ferretti et al. 2001), assuming that it is temporally random and spatially correlated.

A raw time series of displacements is required in order to estimate the APS. The system of eq. (3) is rank-deficient when combining multi-track interferograms for the time-series inversion. We connect the subsets by solving the velocity between consecutive SAR acquisition dates via the SVD method (Berardino et al. 2002). Basically, the SVD method fills in holes due to loss of coherence or independent subsets using a minimum-norm criteria, somewhat like a smoothing operator. It may introduce sharp patterns at the edge of the holes. To avoid such effects, we interpolate the interferograms before the raw time-series inversion and mask the holes after removing the APS.

We estimate the APS by first smoothing the raw time series in time based on a Gaussian filter with length 0.5 yr. The high-pass component is then smoothed with a spatial low-pass Butterworth filter. In stead of using a fixed size for the spatial filter, we estimate an adaptive window size for each acquisition as the e-folding wavelength α of a 1-D covariance function graphic (Parsons et al. 2006), where cjk is the covariance between pixels j and k, σ2 is the variance and djk is the distance between pixels.

Fig. 6 shows the 1-D covariance functions for all the interferograms before and after APS correction. The mean standard deviation σ and e-folding wavelength α are 4.2 mm and 2.6 km, respectively, before APS estimation (Figs 6b–c). The quadratic orbital models have removed long-wavelength atmospheric delay errors during orbital correction, so that the mean e-folding wavelength here is smaller than those estimated elsewhere (e.g. about 10 km in Hanssen 2001 and Wang et al. 2009). The σ and α are 1.5 mm and 0.6 km, respectively, after removing the APS estimates (Figs 6d–e). The APS correction has reduced the amplitude and wavelength of the covariance, on average, by about 64 and 77 per cent, respectively, comparable to those derived in Gourmelen et al. (2010).

Figure 6

(a) Exponential covariance function of the interferograms before (blue) and after (red) APS correction. (b)–(c) Histogram of standard deviation and e-folding wavelength of the exponential function before APS correction. (d)–(e) Histogram of standard deviation and e-folding wavelength of the exponential function after APS correction.

Figure 6

(a) Exponential covariance function of the interferograms before (blue) and after (red) APS correction. (b)–(c) Histogram of standard deviation and e-folding wavelength of the exponential function before APS correction. (d)–(e) Histogram of standard deviation and e-folding wavelength of the exponential function after APS correction.

4.5 Time series and rate map estimation

The effect of APS estimation depends on the window size of the spatio-temporal filters. Smaller scale or anisotropic atmospheric delay errors may remain in the interferograms even after APS correction. For the final time-series analysis, further smoothing is required to remove residual APS as well as other error sources such as decorrelation and thermal noise (Hanssen 2001). As referred to above, the system of eq. (3) is rank-deficit when combining multi-track interferograms. Following Schmidt & Bürgmann (2003), we introduce a Laplacian smoothing operator to the system to constrain the time series (eq. 11).  

11
formula
where 2 is a finite difference approximation of the Laplacian smoothing operator, and the factor κ2 determines the weight of smoothing. We select a smoothing factor for the best-fit time series using a trade-off curve between weighted misfit and solution roughness (Fig. S1).

Most radar images were not acquired at evenly spaced intervals, but it is reasonable to expect temporally smooth velocities for consecutive time-spans. Therefore, we constrain the smoothness of the velocities between the consecutive acquisitions instead of the incremental displacements used in Schmidt & Bürgmann (2003).

We solve the system of eq. (11) using a weighted least-squares method by minimising the following objective function,  

12
formula
where W is the matrix from Cholesky decomposition of the weight matrix graphic, and Cd is a full covariance matrix in time (Biggs et al. 2007). The estimated parameters and their associated errors for each pixel are given by  
13
formula
 
14
formula
where K22.

We estimate linear displacement rates based on the eqs (6)–(10) using an iterative weighted least squares method (Wang et al. 2009). We estimate the displacement rate for each pixel and evaluate the ratios between the residuals and their associated a priori standard deviations. The observed data with the largest ratio are removed to exclude outliers, then we redo the rate inversion until all the ratios are smaller than a given threshold (3 is used in this study). Such procedure is more reliable than fitting a linear trend of the cumulative time series which might be contaminated by certain epochs. The displacement rates and their uncertainties are given by (15) and (16) in a least squares sense.  

15
formula
 
16
formula

5 Results and Discussion

We first estimate the linear rate map and their uncertainties for each track independently (Fig. S2). Although all the rate maps show subsidence along the coast, the magnitudes are not consistent and some features away from the coast are unreliable. The inconsistency is mostly due to strong atmospheric artefacts since more than half of the images were acquired during monsoon season. Since atmospheric artefacts are random in time (Hanssen 2001), they can be mitigated by stacking multiple interferograms (Wright et al. 2001; Biggs et al. 2007). According to eq. (16), the accuracy of the rate map depends on the number, time-span, and accuracy of the interferograms which form matrices G′ and Cd. By combination of interferograms from three tracks, we take use of data about three times that of each track individually, resulting in cleaner and more accurate rate map (Fig. 7a). Fig. 7(b) shows 1−σ uncertainty of the resultant rate map that is lower than that of each track (Figs S2d–f).

Figure 7

(a) Mean displacement rates in vertical (positive means subsidence) and (b) the associated uncertainties. The black lines indicate two profiles in the western and central part of Shenzhen. (c) Rates along the profiles 1–7 (red) and 1′–7′ (blue). The numbers are the average velocities for different segments. (d) The gray bands show elevation along the profile 1–7. The circles denote the time when the land surface will fall to mean sea level considering the observed local subsidence rates and global sea level rise (assuming 3 mm yr−1).

Figure 7

(a) Mean displacement rates in vertical (positive means subsidence) and (b) the associated uncertainties. The black lines indicate two profiles in the western and central part of Shenzhen. (c) Rates along the profiles 1–7 (red) and 1′–7′ (blue). The numbers are the average velocities for different segments. (d) The gray bands show elevation along the profile 1–7. The circles denote the time when the land surface will fall to mean sea level considering the observed local subsidence rates and global sea level rise (assuming 3 mm yr−1).

The most conspicuous feature in the rate map is an elongated subsidence zone along the coast of Shenzhen, a new megacity developed from a small village since 1980s. The red dots in Fig. 7(c) indicate subsidence rates along the profile 1–7 within a swath width of 500 m. The average subsidence rate is about 2.5 mm yr−1, but it is non-uniformly distributed along the profile. We find the highest subsidence rates at Shajing (profile 1–2) and Fuyong (profile 3–4), both with average displacement rates larger than 3 mm yr−1. The maximum velocity is even up to 6–7 mm yr−1 here. Between these two segments (profile 2–3), the average subsidence rate is 2.9 ± 0.6 mm yr−1. The rapid subsidence is most likely related to the ongoing urbanisation in the two towns, being built as part of the Binhai (coastal) New Area of Shenzhen.

The Shenzhen Bao'an International Airport lies along the profile 4–5. It is the fifth busiest airport in China, which opened in 1991 and commenced international operations in 1993. Although built only 20 yr ago, the airport is currently subsiding at rates of 1–2 mm yr−1. Since 2008, the airport has undergone major expansion for a second runway with a length of 3.6 km and width of 60 m and a third terminal connecting the current runway. The runway was completed in June 2011 and the new terminal is expected to be completed in 2012. The new runway was built on reclaimed land extending into the PRD. Unfortunately, we were not able to measure subsidence here because of poor coherence due to rapid surface changes.

The subsidence is slower along the profile 5–6 with an average rate of 0.9 ± 1.3 mm yr−1 where there are mostly residential area, parks and schools. But the subsidence rate increases to 2.8 ± 0.7 mm yr−1 again southwards at the western side of the Binhai Avenue (profile 6–7), a newly built area near the centre of Bao'an District.

Another subsidence zone (marked as B in Fig. 7a) lies at Shatian village of Dongguan, an industrial city in the PRD region. The subsidence rates are about 2–4 mm yr−1. Here, there is the typical geomorphology of a river delta, located in an alluvial plain with braided rivers. The soft soil layer is thicker than 10 m and the maximum is about 19 m at Shatian. In the PRD region, most subsidence hazards are caused by the compaction of soft soil associated with anthropogenic activities. Therefore, extended observations should be made for the other areas, such as Guangzhou, Foshan and Dongguan, where similar conditions exist.

On the other hand, we identify that a few sites are uplifting in the study area. The most apparent two sites are denoted as D1 and D2 in Fig. 7(a), both lying near reservoirs. The uplift rates are about 3–4 mm yr−1 at both sites. We suspect that this is related to the water level changes in the reservoirs although no data are available in public yet. Another area of uplift occurs offshore of Humen marked as C in Fig. 7(a), a village of Dongguan, though the reason for this uplift is not known.

As described above, we fixed the average deformation as zero for the rate map estimation. To evaluate the effect of the reference phase, we make another profile 1′–7′ parallel to the profile 1–7 at the centre of the study area. It shows that the average subsidence rates are between −0.6 and 0.9 mm yr−1 along the profile with 1σ errors of about 0.5 mm yr−1 (blue dots in Fig. 7c). The average subsidence rate is 0.3 mm yr−1 along the full profile that is within 1σ bounds. It implies that most parts of the area are stable and the observed deformation can be regarded as absolute, assuming that there is no image-wide systematic subsidence.

Most of the area is less than 2 m above sea level along the profile 1–7 in terms of the SRTM DEM (Figs 1 and 7d). No publicly available data exist for local sea level rise in Shenzhen. Assuming a global mean sea level rise of 3 mm yr−1 (Church & White 2011), we calculate the time when the land surface will fall to sea level for 116 equally spaced points along the profile 1–7. Among them, 41 points are less than zero now, and 45 will fall to sea level in 300 yr. Our results show that many areas are subsiding at rates comparable to or even much higher than the rising sea level, thus amplifying the rates of local sea level rise. Therefore, assessments of risk due to flooding hazard based only on sea level rise have been underestimated in this coastal region.

Seasonal deformation has been found routinely for post-processing continuous GPS time series, for example, QOCA software by Dong et al. (2002). Because of low time resolution, it is difficult for InSAR to distinguish between linear and seasonal deformation although such signal has been reported before (Hoffmann et al. 2001; Schmidt & Bürgmann 2003; Lanari et al. 2004; Bell et al. 2008). Here, we combine InSAR data from three tracks to yield displacement time series (Fig. S3). The temporal resolution should be roughly three times of that from a single track (i.e. about 11 days for the Envisat data). However, our time series has an average time interval of 22 days because of the gaps in SAR acquisitions. The time-series profiles in the Fig. 8 clearly show non-linear cumulative displacements for all the sites. The blue dots in Fig. 8 indicate the residuals after removing a linear slope from the original time series. Despite different wavelengths, all the residuals show a clear cyclic deformation with amplitude of about 1–2 mm mostly associated with seasonal variations. We fit the residuals using an annual and semi-annual function (Fig. 8). Some are fitted very well (e.g. points B and C) using this model. We must be aware the trade-off for the wavelength between the temporal filter in APS estimation and the residual cyclic deformation. Periodic deformation shorter than half a year, the temporal filter cutoff used for APS estimation, would be regarded as high-frequency noise. Therefore, shorter period deformation may be found when higher temporal-resolution SAR images are available.

Figure 8

Time series of the six sites (A1–A6) along the profile 1–7 and four (B, C, D1, D2) in the other areas. The locations of these sites are denoted in Fig. 7a. The red dots indicate time series of the cumulative displacements with 1σ errors, and the blue indicate residuals after removing a linear trend from the cumulative displacements. The black curves show annual and semi-annual models to fit the residual displacements. The average velocity for each point is shown in the box.

Figure 8

Time series of the six sites (A1–A6) along the profile 1–7 and four (B, C, D1, D2) in the other areas. The locations of these sites are denoted in Fig. 7a. The red dots indicate time series of the cumulative displacements with 1σ errors, and the blue indicate residuals after removing a linear trend from the cumulative displacements. The black curves show annual and semi-annual models to fit the residual displacements. The average velocity for each point is shown in the box.

6 Conclusions

We have developed a method, implemented in the Π-RATE software, to combine multi-track interferograms for time-series analysis assuming purely vertical deformation. It can also be extended for other applications if the displacement direction is known, for example, assuming pure horizontal deformation for a given strike-slip fault. Furthermore, it can be used for measuring deformation in areas where SAR data overlap if the data from a single track are not sufficient to yield a good time series. This method enables the large archives of ERS, Envisat and other SAR satellite data to be fully exploited. The method can significantly improve temporal resolution, which is one of the most important limitations for InSAR applications in civil engineering. For the ESA Sentinel-1 satellite to be launched in 2013 with a revisit time of 12 days, the temporal resolution could be up to 3 days in the overlapping areas of neighbouring tracks observed by ascending and descending orbits (Torres et al. 2012). Because our method does not rely on a single satellite, we can combine data from multiple satellites (e.g. Sentinel-1, TerraSAR-X and COSMO-SkyMed) with different look angles and acquisition modes. By doing so, InSAR time series with daily temporal resolution will be possible in the next few years, which is enough for most engineering applications.

By combination of interferograms from two ascending and one descending tracks, we obtain the deformation rate map and time series in the PRD region. The results show an elongated zone of rapid subsidence associated with the urbanisation of Shenzhen, particularly at Shijing and Fuyong. The subsidence rates vary in space with an average of 2.5 mm yr−1. We also find subsidence at Shatian, Dongguan, due to the compaction of the thick soft soil layer, the most prevalent cause for land subsidence in the delta. Although no regional-scale subsidence was observed in our rate map, long-term routine deformation monitoring should be conducted in the PRD region where the subsiding coastal area and rising sea level are a considerable hazard.

Coastal subsidence has been a critical problem in many countries. It has even more significant impacts in China, which has a coastline length of about 180 000 km. The fast economic and population growth is mainly from the coastal cities in eastern China over the past 30 yr. Significant coastal development has been a consequence of China's fast economic growth. All the special economic zones (SEZs) and national economic development areas lie in the coastal regions. However, many megacities have suffered from land subsidence over many years due to fast development. Unlike the PRD, land subsidence hazards are widespread in the Yangtze River Delta surrounding Shanghai and Bohai Economic Zone surrounding Beijing. In both of the regions, much more SAR data have been acquired from multiple tracks. The multi-track InSAR method developed in this study can therefore play an important role for mitigating the effects of land subsidence and sea level rise.

Acknowledgments

The work was supported by Guangzhou Science and Technology Program (11G0041) and Guangzhou Urban Planning & Design Survey Research Institute. TJW and HW were supported by the Royal Society through a University Research Fellowship and an International Incoming Fellowship respectively. SAR data are copyrighted by ESA through a category-1 project 7924. We used the JPL/Caltech ROI_PAC software to produce InSAR data, and the Generic Mapping Tools (GMT) to prepare the figures. We thank Matt Garthwaite for discussions, Duncan Agnew and two anonymous reviewers for constructive comments that helped improve the manuscript. The Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET+) is part of the UK Natural Environment Research Council's National Centre for Earth Observation.

References

Aly
M.H.
Zebker
H.A.
Giardino
J.R.
Klein
A.G.
,
2009
.
Permanent Scatterer investigation of land subsidence in Greater Cairo, Egypt
,
Geophys. J. Int.
 ,
178
(
3
),
1238
1245
.
Argus
D.F.
Heflin
M.B.
Peltzer
G.
Crampe
F.
Webb
F.H.
,
2005
.
Interseismic strain accumulation and anthropogenic motion in metropolitan Los Angeles
,
J. geophys. Res.
 ,
110
,
B04401
, doi:10.1029/2003JB002934.
Bell
J.W.
Amelung
F.
Ferretti
A.
Bianchi
M.
Novali
F.
,
2008
.
Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge
,
Water Resour. Res.
 ,
44
,
W02407
, doi:10.1029/2007WR006152.
Berardino
P.
Fornaro
G.
Lanari
R.
Sansosti
E.
,
2002
.
A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms
,
IEEE Trans. Geosci. Remote
 ,
40
(
11
),
2375
2383
.
Biggs
J.
Wright
T.
Lu
Z.
Parsons
B.
,
2007
.
Multi-interferogram method for measuring interseismic deformation: Denali Fault, Alaska
,
Geophys. J. Int.
 ,
170
,
1165
1179
.
Calais
E.
Dong
L.
Wang
M.
Shen
Z.
Vergnolle
M.
,
2006
.
Continental deformation in Asia from a combined GPS solution
,
Geophys. Res. Lett.
 ,
33
,
L24319
, doi:10.1029/2006GL028433.
Cazenave
A.
Llovel
W.
,
2010
.
Contemporary sea level rise
,
Annu. Rev. Mar. Sci.
 ,
2
,
145
173
.
Church
J.A.
White
N.J.
,
2011
.
Sea-level rise from the late 19th to the early 21st century
,
Surv. Geophys.
 ,
32
,
585
602
.
Church
J.A.
, et al.,
2011
.
Revisiting the earth's sea-level and energy budgets from 1961 to 2008
,
Geophys. Res. Lett.
 ,
38
,
L18601
, doi:10.1029/2011GL048794.
Dijkstra
E.W.
,
1959
.
A note on two problems in connexion with graphs
,
Numer. Math.
 ,
1
,
269
271
.
Doin
M.-P.
Lasserre
C.
Peltzer
G.
Cavalie
O.
Doubre
C.
,
2009
.
Corrections of stratified tropospheric delays in SAR interferometry: validation with global atmospheric models
,
J. Appl. Geophys.
 ,
69
(
1
),
35
50
.
Dong
D.
Fang
P.
Bock
Y.
Cheng
M.K.
Miyazaki
S.
,
2002
.
Anatomy of apparent seasonal variations from GPS-derived site position time series
,
J. geophys. Res.
 ,
107
(
B4
),
2075
, doi:10.1029/2001JB000573.
Elliott
J.R.
Biggs
J.
Parsons
B.
Wright
T.J.
,
2008
.
InSAR slip rate determination on the Altyn Tagh Fault, northern Tibet, in the presence of topographically correlated atmospheric delays
,
Geophys. Res. Lett.
 ,
35
,
L12309
, doi:10.1029/2008GL033659.
Farr
T.G.
, et al.,
2007
.
The shuttle radar topography mission
,
Rev. Geophys.
 ,
45
(
2
),
RG2004
, doi:10.1029/2005RG000183.
Ferretti
A.
Prati
C.
Rocca
F.
,
2001
.
Permanent Scatterers in SAR interferometry
,
IEEE Trans. Geosci. Remote
 ,
39
(
1
),
8
20
.
Finnegan
N.J.
Pritchard
M.E.
Lohman
R.B.
Lundgren
P.R.
,
2008
.
Constraints on surface deformation in the Seattle, WA, urban corridor from satellite radar interferometry time-series analysis
,
Geophys. J. Int.
 ,
174
(
1
),
29
41
.
Goldstein
R.M.
Werner
C.L.
,
1998
.
Radar interferogram filtering for geophysical applications
,
Geophys. Res. Lett.
 ,
25
(
21
),
4035
4038
.
Goldstein
R.M.
Zebker
H.A.
Werner
C.L.
,
1988
.
Satellite radar interferometry: two-dimensional phase unwrapping
,
Radio Sci.
 ,
23
(
4
),
713
720
.
González
P.J.
Fernández
J.
,
2011
.
Drought-driven transient aquifer compaction imaged using multitemporal satellite radar interferometry
,
Geology
 ,
39
(
6
),
551
554
.
Gourmelen
N.
Amelung
F.
Casu
F.
Manzo
M.
Lanari
R.
,
2007
.
Mining-related ground deformation in Crescent Valley, Nevada: implications for sparse GPS networks
,
Geophys. Res. Lett.
 ,
34
,
L09309
, doi:10.1029/2007GL029427.
Gourmelen
N.
Amelung
F.
Lanari
R.
,
2010
.
Interferometric synthetic aperture radar - GPS integration: interseismic strain accumulation across the Hunter Mountain fault in the eastern California shear zone
,
J. geophys. Res.
 ,
115
,
B09408
, doi:10.1029/2009JB007064.
Hanssen
R.F.
,
2001
.
Radar Interferometry: Data Interpretation and Error Analysis
 ,
Kluwer Academic Publishers
Oxford, UK
,
Dordrecht
.
Hetland
E.A.
Musé
P.
Simons
M.
Lin
Y.N.
Agram
P.S.
DiCaprio
C.J.
,
2012
.
Multiscale InSAR time series (MInTS) analysis of surface deformation
,
J. geophys. Res.
 ,
117
,
B02404
, doi:10.1029/2011JB008731.
Hoffmann
J.
Zebker
H.A.
Galloway
D.L.
,
2001
.
Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by synthetic aperture radar interferometry
,
Water Resour. Res.
 ,
37
(
6
),
1511
1566
.
Hooper
A.
,
2008
.
A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches
,
Geophys. Res. Lett.
 ,
35
,
L16302
, doi:10.1029/2008GL03465.
Hooper
A.
Bekaert
D.
Spaans
K.
Arikan
M.
,
2012
.
Recent advances in SAR interferometry time series analysis for measuring crustal deformation
,
Tectonophysics
 ,
514-517
,
1
13
.
Hooper
A.
Segall
P.
Zebker
H.
,
2007
.
Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos
,
J. geophys. Res.
 ,
112
,
B07407
, doi:10.1029/2006JB004763.
Jiang
L.
Lin
H.
Ma
J.
Kong
B.
Wang
Y.
,
2011
.
Potential of small-baseline SAR interferometry for monitoring land subsidence related to underground coal fires: Wuda (Northern China) case study
,
Remote Sens. Environ.
 ,
115
(
2
),
257
268
.
Kruskal
J.B.
,
1956
.
On the shortest spanning subtree of a graph and the traveling salesman problem
,
Proc. Am. Math. Soc.
 ,
7
(
1
),
48
50
.
Lanari
R.
Lundgren
P.
Manzo
M.
Casu
F.
,
2004
.
Satellite radar interferometry time series analysis of surface deformation for Los Angeles, California
,
Geophys. Res. Lett.
 ,
31
,
L23613
, doi:10.1029/2004GL02294.
Mazzotti
S.
Lambert
A.
Van der Kooij
M.
Mainville
A.
,
2009
.
Impact of anthropogenic subsidence on relative sea-level rise in the Fraser river delta
,
Geology
 ,
37
(
9
),
771
774
.
Ozawa
T.
Ueda
H.
,
2011
.
Advanced interferometric synthetic aperture radar (InSAR) time series analysis using interferograms of multiple-orbit tracks: a case study on Miyake-jima
,
J. geophys. Res.
 ,
116
,
B12407
, doi:10.1029/2011JB008489.
Parsons
B.
, et al.,
2006
.
The 1994 Sefidabeh (eastern Iran) earthquake revisited: new evidence from satellite radar interferometry and carbonate dating about the growth of an active fold above a blind thrust fault
,
Geophys. J. Int.
 ,
164
,
202
217
.
Rosen
P.A.
Hensley
S.
Peltzer
G.
Simons
M.
,
2004
.
Updated repeat orbit interferometry package released
,
EOS Trans.
 ,
85
(
5
),
47
, doi:10.1029/2004EO050004.
Sansosti
E.
Casu
F.
Manzo
M.
R
L.
,
2010
.
Space-borne radar interferometry techniques for the generation of deformation time series: an advanced tool for earth's surface displacement analysis
,
Geophys. Res. Lett.
 ,
37
,
L20305
, doi:10.1029/2010GL044379.
Scharroo
R.
Visser
P.
,
1998
.
Precise orbit determination and gravity field improvement for the ERS satellite
,
J. geophys. Res.
 ,
103
(
C4
),
8113
8127
.
Schmidt
D.A.
Bürgmann
R.
,
2003
.
Time-dependent land uplift and subsidence in the Santa Clara valley, California, from a large interferometric synthetic aperture radar data set
,
J. geophys. Res.
 ,
18
(
B9
),
2416
, doi:10.1029/2002JB002267.
Simons
M.
Rosen
P.A.
,
2007
.
Interferometric synthetic aperture radar geodesy
, in
Treatise on Geophysics
 , pp.
391
446
, editor in,
chief Gerald Schubert
E.
,
Elsevier
Oxford, UK
,
Amsterdam
.
Syvitski
J.P.M.
, et al.,
2009
.
Sinking deltas due to human activities
,
Nature Geosci.
 ,
2
,
681
686
.
Teatini
P.
Tosi
L.
Strozzi
T.
,
2011
.
Quantitative evidence that compaction of Holocene sediments drives the present land subsidence of the Po Delta, Italy
,
J. geophys. Res.
 ,
116
,
B08407
, doi:10.1029/2010JB008122.
Torres
R.
, et al.,
2012
.
GMES Sentinel-1 mission
,
Remote Sens. Environ.
 ,
120
,
9
24
.
Wang
H.
Ge
L.
Xu
C.
Du
Z.
,
2007
.
3-D coseismic displacement field of the 2005 Kashmir earthquake inferred from satellite radar imagery
,
Earth Planets Space
 ,
59
(
5
),
343
349
.
Wang
H.
Wright
T.J.
,
2012
.
Satellite geodetic imaging reveals internal deformation of western Tibet
,
Geophys. Res. Lett.
 ,
39
,
L07303
, doi:10.1029/2012GL051222.
Wang
H.
Wright
T.J.
Biggs
J.
,
2009
.
Interseismic slip rate of the northwestern Xianshuihe fault from InSAR data
,
Geophys. Res. Lett.
 ,
36
,
L03302
, doi:10.1029/2008GL036560.
Wisely
B.A.
Schmidt
D.
,
2010
.
Deciphering vertical deformation and poroelastic parameters in a tectonically active fault-bound aquifer using InSAR and well level data, San Bernardino basin, California
,
Geophys. J. Int.
 ,
181
(
3
),
1185
1200
.
Wright
T.
Parsons
B.
Fielding
E.
,
2001
.
Measurement of interseismic strain accumulation across the North Anatolian Fault by satellite radar interferometry
,
Geophys. Res. Lett.
 ,
28
(
10
),
2117
2120
.
Wright
T.J.
Parsons
B.E.
Lu
Z.
,
2004
.
Toward mapping surface deformation in three dimensions using InSAR
,
Geophys. Res. Lett.
 ,
31
,
L01607
, doi:10.1029/2003GL018827.