Natural Hazards and Earth System Sciences Ground Subsidence Geo-hazards Induced by Rapid Urbanization: Implications from Insar Observation and Geological Analysis

Due to the convenient transportation and construction , cities are prone to be situated in areas with flat terrain and unstable sediments, resulting in the concurrence of ground subsidence and urbanization. Here the interaction between geology, anthropogenic processes and ground subsidence geo-hazards were investigated in the Greater Pearl River Delta region of China. Geological evidences and 2006–2010 persistent scatterer data indicate that anthro-pogenic activities are dominant, although the distribution of river system and Quaternary sediments are also highly related to significant displacements (primarily at a rate of −15 to 15 mm a −1). The surface displacements derived by synthetic aperture radar interferometry suggest that the urbanization rhythm has to be routinely monitored. Considering analogous urbanization modes, particularly in developing countries, ground subsidence monitoring together with the analysis of its driving force are critical for geo-hazards early-warning, city planning as well as sustainable urbanization .


Introduction
Urbanization is a global trend, particularly in developing countries because of the productivity rise in dense agglomerations.Apart from advantages, e.g.closeness, lower carbon emission (Glaeser, 2011); urbanization is faced with challenges, including contagious disease and congestion.The acceleration of urbanization has enhanced human influence on the natural environment, requiring a detailed insight into the limitations on development that will enable us to avoid manmade disasters.For the convenient transportation and construction, cities are prone to be situated in flat areas accompanied by well-developed sediments, particularly in the vicinity of rivers and coasts.Sediments deposited by flowing water or tide dynamics, as in a riverbed, flood plain or delta, are physically unstable.Ground subsidence geo-hazards have frequently occurred in cities situated on unconsolidated alluvial deposits or reclaimed regions resulting from creeping, compaction and infrastructure construction (Meckel et al., 2006;Raucoules et al., 2008;Kim et al., 2010).
In this study, the Greater Pearl River Delta (GPRD), representative of a fast urbanization region in alluvial and lacustrine deposits of the world, is investigated.The GPRD, comprising Hong Kong, Macao and the nine municipalities (Guangzhou, Shenzhen, Zhuhai, Dongguan, Zhongshan, Foshan, Huizhou, Jiangmen, and Zhaoqing) of Guangdong Province in mainland China, is one of the most densely urbanized regions and the main hub of China's economic growth.Since the launch of China's reform program in 1978, the GPRD has set several world records, e.g. the fastest growth in gross domestic product (GDP), from US$ 8 billion in 1980 to US$ 100 billion in 2001 as well as a population explosion in GPRD cities such as the fastest one, Shenzhen, whose population increased from around 30 000 in 1980 to over 14 000 000 in 2011 -nearly 500 times.The rapid urbanization has exerted pressures (Karen et al., 2000;Zeng et al., 2008;Stone, 2009) on the sustainable development of the GPRD.Geo-hazards initiated by severe ground subsidence have been highlighted due to direct economic losses (3.47 billion RMB) and casualties (768 deaths) from 1994 to 2005 in the GRPD (Zhu et al., 2007).Interferometric synthetic aperture radar (InSAR) (Zebker et al., 1994;Jonsson et al., 2003), particularly the persistent scatterer interferomety (PSI) (Ferretti et al., 2001;Lanari et al., 2004;Tizzani et al., 2007;Lin et al., 2011), has proven to be powerful in determining a time series of displacements by exploiting stable, point-like radar targets with a large stack of SAR images (>15) (Hilley et al., 2004).Although several PSI studies were conducted (Zhao et al., 2009;Chen et al., 2010), they focused on a single city and left major gaps in terms of spatial coverage of the GPRD.Nowadays, the GPRD is integrated as a whole because of the geological environment as well as economic and policy coordination.Furthermore, the large-scale displacement trend only can be detected by investigating the entire agglomeration.Otherwise, hot-spots for potential geo-hazards can easily be left out.To address this issue and understand the exact principle of motion (local-scale and large-scale), it is essential to estimate the ground displacement and its time series for geo-hazards early-warning related to fast urbanization, via the analysis of motion spatial heterogeneity or using displacement thresholds based on engineering specifications.

Data and methodology
In total, 61 Envisat ASAR single look complex (SLC) images from two ascending tracks in image swath 2 mode (31 from Track 25 and 30 from Track 297), as listed in Table 1, were used to cover most urban regions of the GPRD (except for Zhaoqing City) of approximately 150 × 200 km 2 ; the area extends from 109 • 23 to 117 • 38 E and from 19 • 50 to 25 • 35 N (Fig. 1).All the ASAR images used in this work were acquired at the Hong Kong Remote Sensing Ground Receiving Station, which was established at the Chinese University of Hong Kong in 2006.The pixel spacing was 7.80 m in the range direction, and 4.05 m in the azimuth direction.DORIS precise orbits data provided by the ES-RIN Help desk of ESA were applied to calculate interferometric perpendicular baselines.The 3-arcsecond (∼90 m) Shuttle Radar Topography Mission (SRTM) (Farr and Kobrick, 2000) DEM data from the United States Geological Survey (USGS) were used for topographic phase estimation at the first step, and then for geocoding the InSAR products (transforming Range-Doppler coordinates into Universal Transverse Mercator map geometry system) were used.
When the temporal image co-registration was completed, interferograms were generated by conjugate multiplication of complex image pairs.The small baseline strategy (smaller than 350 m spatially and 3 months temporally) was applied in this procedure, resulting in high coherence.The corresponding interferogram formation graphs are illustrated in  Effects of flat terrain were estimated by the DORIS precise orbits data, and the topography phases were removed using the above mentioned SRTM DEM data.We used an improved PSI model developed from interferometric point target analysis (IPTA) to analyse the effects of urbanization in the GRPD, and to expose the possible geo-hazards of this runaway growth in population and industry.Generally, the PSI implementation exploits the temporal and spatial characteristics of interferometric signals collected from point-likely or high coherent targets (referred to as PS candidates hereafter for simplicity), and allows a precise estimation of deformation, topographic phase and other parameters.The PS candidates were determined by two selection criteria in this study.One was based on the spectral properties of each single SLC image; the other was derived by utilizing the average spatial coherence of all generated interferograms.In such a way, several distributed scatterers (DS) (Ferretti et al., 2011) showing similar reflectivity and coherence values could also be extracted to increase the density of PS candidates.Phase unwrapping was performed in the spatial domain using a minimum cost-flow (MCF) algorithm for sparse data (Costantini and Rosen, 1999).Then the i-th unwrapped interferometric phase φ x,i at the location x is the sum of topographic phase φ topo,x,i , ground deformation phase φ def,x,i , differential atmospheric phase φ atm,x,i and phase noise φ noise,x,i : The whole deformation rates along line of sight (LOS) from September 2006 to August 2010 were estimated by the improved PSI model, negative for subsidence, i.e. the motion of the ground away from satellite.Several key points, including residual orbits phase estimation, atmospheric phase screen mitigation, and large-scale displacement inversion with single-reference point, have been mastered: -residual orbits phase was initially estimated on the differential interferogram fringe and SLC image coregistration offsets, then a 2-D quadratic model phase function is employed for global phase trend removal; -global component of atmospheric phase screen has been partially mitigated using the above 2-D quadratic model; the residual component was further eliminated by filters in space-temporal domain during regression; -the density of measurement points increased thanks to the small baseline strategy; high density PS candidates guarantee the successful implementation of large-scale deformation inversion using a single-reference point (located in stable old towns and selected by manual) for each frame.
The main steps of the improved PSI are summarized as follows: first, a simulation model based on singular value decomposition (SVD) was used for n differential phase time series retrieval (single-reference stack).Then, a twodimensional regression based on a linear deformation model was applied on the baseline and acquisition time interval domain to estimate the residual height ε and linear displacement rate v by maximizing the following equation (Ferretti et al., 2001): The temporal ensemble coherence ( ), ranging from 0 to 1, was considered as a reliability measure in fitting the predefined displacement model.Finally, the atmospheric phase, non-linear displacement and error terms of the residual phase were further discriminated utilizing their different spatial and temporal dependencies.The derived accuracy of PSI model could reach as high as 1 to 2 mm a −1 (Teatini et al., 2005;Stramondo et al., 2008) with satisfying long time spanning, e.g.approximate 5 yr in this investigation; however, the evaluation of the error also depends on the specific data used, pre-processing, time interval, and observed scenarios.Compared with the conventional PSI processing approach, except for urban areas, the proposed PSI model is feasible for largescale natural regions monitoring thanks to the enhanced spatial density of PS candidates, which is similar to the newly developed SqueeSAR algorithm (Ferretti et al., 2011).Moreover, the improved PSI is more efficient than SqueeSAR because of its straightforward procedures at cost of resolution and accuracy loss jointly determined by the spatial averaging and phase unwrapping.
3 Results and analysis

Deformation field
Two stacks of SAR data were processed separately, than the results were integrated using the common PS candidates in the overlap region.That is, the motion values from the corresponding frame (Track 297) were calibrated based on the reference frame (Track 25).The final derived velocity field shows no large-scale subsidence bowls; however, local displacement regions of several square kilometres occurred over the study site.The corresponding values are primarily in the range of −15 mm a −1 to 15 mm a −1 .With the least square strategy, the velocity offset and standard deviation (also referred to as relative accuracy) between the two adjacent tracks were estimated, presenting 2.35 mm a −1 and 3.18 mm a −1 , respectively.It is clear that the displacements are located along the Pearl River, as marked by the blue line (interpreted by manual using the SPOT 5 image) in Fig. 3.The GPRD is formed by three major branches of the Pearl River: the Xi Jiang, Bei Jiang and Dong Jiang, and is comprised of two alluvial deltas separated by the main stream of the Pearl River.Due to the combination of fluvial outwash and tide dynamics, the depth of the silt layer is approximately 15-40 m.This kind of sediment or alluvial deposit, with a high compressibility, low permeability and bearing capacity, can be unstable when external forces are applied, easily leading to ground deformation damage, foundation sinking, as well as seismic subsidence (Zhao et al., 2009).The association of higher ground subsidence with the geology indicates that the clay deposits are unstable and alluvial outwash dynamics can play a role in the surface displacement following alternating drought-flood seasons.This physical mechanism may be the cause of building collapse along the He Jiang in August 2010, Zhaoqing City.Signatures of the uplift are jointly determined by the alluvial accumulation, phase errors from residual height and unwrapping procedures.Figure 3a and b shows the detailed surface displacements of the Guangzhou and Jiangmen urban areas.Field investigations of the past few years demonstrate that the spatial distribution of surface deformation in Guangzhou has a close relationship with human activity, including underground civil projects, construction of factories, and the expansion of development zones.Geo-hazards, including ground collapse, debris flow and building damage, occurred frequently within the subsidence zones.In 2007-2008, six ground collapse accidents fell within the subsidence region induced by the construction of metro lines (Zhao et al., 2009), marked by the white triangle in Fig. 3a; in 2010, tens of buildings in Jinshazhou, Guangzhou were seriously damaged by local ground subsidence induced by the high-speed railway underground project, resulting in an economic loss of 30 million RMB, marked by the white pentagram in Fig. 3a.The location of deformation zones in Jiangmen City corresponds with the urban expansion pattern of the last 10 yr; that is, approximately 80 % of subsidence areas are located in new development zones, implying induced consolidation of the soft soil during and after completion of the construction projects as the cause, marked by the white arrow in Fig. 3b.In addition, the elongated pattern of significantly moving PS corresponds with the location of the river network, which is partially determined by the more compressible deposits resulting from the fluvial outwash.

Deformation vs. Quaternary geology
Apart from anthropogenic activity, Quaternary geology is another prominent factor.The Pearl River basin was originally formed in the northern and central parts of the present delta as a result of faulting.During the Quaternary, a large amount of silt was deposited in the mouths of the Xi Jiang, Bei Jiang, and Dong Jiang, and eventually formed the composite GPRD (Weng, 2007).The NE-SW, NW-SE, and E-W faults intersect, arranging the geological structural units into a chessboard pattern, where the limestone corrosion fractures and the Karst caves are gestated.Taking Guangzhou City as an example, Fig. 4 shows that it is not clear whether the group of faults (Guangzhou-Conghua) in the urban area is related to the distribution of significant subsidence (extracted by a pre-defined threshold < −8 mm a −1 ), but it is highly related to the distribution of Quaternary sediments, especially to the boundary of Quaternary sediments and Tertiary rocks, see Sects. 1 and 2 of N 1 and Sect. 3 of P z1 ; many faster displacement points are located along the transition belts of different geological units.The profile sections of a-b and c-d demonstrate the correlation between the boundaries of rock distribution and PS deformation velocities.Note that the PS deformation rates in the profiles are derived from the original velocity field without the threshold filtering.Assuming identical external pressures, different sediments with different porosities can induce uneven displacement trends, which should be taken into account during city planning.

Analysis of driving forces
The driving forces of surface displacement are rather complicated, not even completely understood at present, owing to the diversity of anthropogenic activities, runoff and geology in large-scale urban agglomerations.Figure 4 shows that the displacement field and the geological structural units are weakly interrelated, implying that there was no obvious tectonic motion during the period of observation.Field investigation demonstrates that outwash dynamics and anthropogenic activities interact sometimes, although the latter plays a more significant role, including urban expansion, reclamation, infrastructure construction and underground water pumping, etc.In general, different driving forces correspond to diverse geological deformation modes, e.g.aquifer overexploitation links quadratic or linear sinking, exponential settlement links land reclamation, and sudden deformation links engineering construction.Consequently, the characteristics of historical LOS displacements from InSAR observations allow us to discriminate those external forces via geological deformation modes vs. PS motion time series (Fig. 5).New development or reclamation zones accelerate the compression and creeping of silt layers because of the downward pressure induced by buildings or other infrastructures, resulting in severe surface subsidence (−25.2 mm a −1 ) at the primary consolidation stage; this phenomenon can be explained by the Terzaghi theory (Terzaghi and Peck, 1967).However after another 2 ∼ 3 yr, secondary compression occurs and the surface displacement is nearly in balance (Fig. 5a).The surface subsidence triggered by infrastructure construction has a close relationship to the project period (Fruneau and Sarti, 2000); compared with remarkable subsidence during the working period (−31.3 mm a −1 ), the surface is relatively stable before and after the completion of a project, demonstrating cascaded motions (Fig. 5b).Industry zones need a large volume of fresh water for the manufacturing process.Groundwater overexploitation changes the local water table (Raucoules et al., 2008), which results in settlement and damage to neighbourhood buildings or sewage systems.Figure 5c

Discussion and conclusion
In the past 30 yr, starting with new cities called "Special Economic Zones" established next to Hong Kong and Macau, the spatial mode of rural urbanization in the GPRD has evolved from "micro-centralization and macro-dispersal" to "macrocentralization and micro-dispersal".The polycentric mega urban region around the estuary of the Pearl River has risen to dominate the landscape from Hong Kong to Guangzhou and Foshan in the northwest and then through Zhongshan down to Zhuhai and Macau in the southwest.Since 2005, the shortage of land in the urban centers of the GPRD has triggered urban sprawl from land reclamation and rural urbanization, which in turn shifts the infrastructure activities from the inner cities to those "micro-dispersal" suburban or development districts.The newly developed towns and displacement zones need extra attention by city planners to avoid underlying geo-hazards.The frequent concurrences of ground subsidence geohazards and fast urbanization have been the focus of modern society, and thus feasible monitoring solutions are imperative.Our study demonstrates an improved PSI method that can resolve detailed ground settlements for the urbanization process.The results in GPRD show that potential geo-hazards related to surface displacement can be precisely detected by the space-based InSAR technology.Although the driving forces of ground settlements are complicated, it is possible to discriminate some of them by using the time series analysis.We conclude that anthropogenic influences are the determinant for the ground subsidence geo-hazards in the GPRD.The distribution of river system and Quaternary sediments is highly related to significant displacements, indicating the demand for geological investigation before future city planning, owing to the well-developed, overlying silt layers.Considering the analogous geological condition (e.g.sediments or alluvial deposits) and human behaviours (e.g.urban expansion, land reclamation and underground project construction), this study describes a potential geo-hazards early-warning approach under global urbanization by jointly considering the motion spatial heterogeneity and displacement thresholds based on engineering specifications.Excessively rapid urbanization not only induces pollution and places extra stress on the area's limited resources, but also causes serious ground subsidence, hampering the economic zone's sustainable development.Utilizing the close relationship between obvious surface displacement and geological disasters (Steven et al., 2003;Chen et al., 2010), improper anthropogenic activities (e.g.resource over-exploitation and low-quality infrastructure construction) should be avoided by governments through scientific management (Glaeser, 2011) to make agglomerations be the centre of productivity and pleasure.

Fig. 1 .
Fig. 1.Location of the Greater Pearl River Delta (GPRD) region, comprising Hong Kong, Macao and the nine municipalities of Guangdong Province.Two ascending tracks of Envisat ASAR data (Tracks 25 and 297), marked by red rectangles, are analyzed for surface settlement inversion covering the whole GPRD region except for Zhaoqing City.

Fig. 2 .
Fig. 2. Baseline-time plots relevant to the ASAR images and interferograms used in this study; (a) Track 25, (b) Track 297.

Fig. 2a (
Fig. 2a (Track 25) and Fig. 2b(Track 297), respectively.Effects of flat terrain were estimated by the DORIS precise orbits data, and the topography phases were removed using the above mentioned SRTM DEM data.We used an improved PSI model developed from interferometric point target analysis (IPTA) to analyse the effects of urbanization in the GRPD, and to expose the possible geo-hazards of this runaway growth in population and industry.Generally, the PSI implementation exploits the temporal and spatial characteristics of interferometric signals collected from point-likely or high coherent targets (referred to as PS candidates hereafter for simplicity), and allows a precise estimation of deformation, topographic phase and other parameters.The PS candidates were determined by two selection criteria in this study.One was based on the spectral properties of each single SLC image; the other was derived by utilizing the average spatial coherence of all generated interferograms.In such a way, several distributed scatterers (DS)(Ferretti et al., 2011) showing similar reflectivity and coherence values could also be extracted to increase the density of PS candidates.Phase unwrapping was performed in the spatial domain using a minimum cost-flow (MCF) algorithm for sparse data(Costantini and Rosen, 1999).Then

Fig. 3 .
Fig. 3. Annual surface displacement rates over the GPRD (marked by red rectangles in Fig. 1) derived by Persistent Scatterer Interferometry (PSI) model, overlapped on a SPOT remote sensing image.The primary values are in the range of −15 mm a −1 to 15 mm a −1 in the span of 2006 ∼ 2010, and the significant displacements are located along the river network, marked by the blue line (derived by manual using the SPOT 5 image interpretation).Two detailed displacement images (a) Guangzhou (b) Jiangmen combined with field investigations demonstrate that the settlements have close relationship with anthropogenic activities.The white cross indicates the reference point for the integrated whole velocity field; symbols of pentagram, triangle and arrow indicate the locations of ground collapse accidents in 2007-2008, ground instability by the high-speed railway track in 2010, and Jiang men new development zones built in the last 10 yr, respectively.The numbers "1, 2 and 3" indicate the locations of PS time series in Fig. 5a, b and c.

Fig. 4 .
Fig. 4. Geological relationship with ground displacement in Guangzhou, covering same area as Fig. 3a.Instead of the two groups of faults, the distribution of Quaternary deposits, especially to the boundary of Quaternary sediments and late Tertiary rocks, is highly related to significant displacements (< −8 mm a −1 ) validated by the results in Guangzhou (Q d : Quaternary deposits, including gravel, sand and clay; N 1 : Neogene-Miocene conglomerate, sandstone and shale; E: Paleogene-Eocene conglomerate, sandstone and limestone; J 2 : Jurassic-middle conglomerate, sandstone and shale; J 1 : Jurassic-early conglomerate, sandstone and shale; P z1 : Permian-late shale, sandstone and siltstone; D 3 m : Devonian-late sandstone, shale and siltstone; ηγ : mica granite during Yanshanian period between Jurassic and early-Cretaceous).The profile sections of a-b and c-d are examples of rock distribution, geological and the PS deformation velocities relationships, demonstrating the correlation between the boundary of different geological units (i.e.Quaternary distribution or old rocks) and significant displacement trends.

Fig. 5 .
Fig. 5. Driving force discrimination conducted by the comparison between PS displacement time series and geological deformation modes; (a) a reclamation mode, the subsidence rate can reach as high as −25.2 mm a −1 during the primary consolidation compared with approximately zero in the secondary compression stage; (b) an infrastructure construction mode, the historical displacements are cascaded, demonstrating close relationship with the construction period; (c) a groundwater overexploitation mode, the surface subsidence is ongoing with an average rate of −18.41 mm a −1 during the whole observation span.

Table 1 .
The Envisat ASAR data used in this study.