## 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 km^{2}. About 3,720 km^{2} 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.

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 *d*_{los} (positive away from satellite) can be written as

*d*

_{e},

*d*

_{n}and

*d*

_{u}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 *d*_{u} can then be converted from *d*_{los} using the local incidence angle,

If the interferograms are geocoded with the same pixel spacing, *d*_{u} 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

**d**is a vector containing the vertical displacements

*d*

_{u}.

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

*t*

_{i}=

*t*

_{i + 1}−

*t*

_{i},

*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

*I*

_{i, j},

*v*

_{i}is the velocity of the

*i*th time-span. Here, the acquisitions must be chronologically ordered for all the tracks (Fig. 2).

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 can be solved using weighted least-squares method.

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

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

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 mm^{2} 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 *I*_{1, 2}, *I*_{2, 3} and *I*_{1, 3}, for example, two interferograms have provided all the information observed by the three acquisitions *I*_{1}, *I*_{2} and *I*_{3}. 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.

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

### 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 (Parsons *et al.* 2006), where *c*_{jk} is the covariance between pixels *j* and *k*, σ^{2} is the variance and *d*_{jk} 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).

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

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,

where**W**is the matrix from Cholesky decomposition of the weight matrix , and

**C**

_{d}is a full covariance matrix in time (Biggs

*et al.*2007). The estimated parameters and their associated errors for each pixel are given by where

**K**=κ

^{2}

**∇**

^{2}.

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.

## 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 **C**_{d}. 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).

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 D_{1} and D_{2} 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.

## 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

*et al.*,

*et al.*,

*et al.*,

*et al.*,

*et al.*,