Surface deformation evolution in the Pearl River Delta between 2006 and 2011 derived from the ALOS1/PALSAR images

This study monitors the land subsidence of the whole Pearl River Delta (PRD) (area: ~ 40,000 km2) in China using the ALOS1/PALSAR data (2006–2011) through the SBAS-InSAR method. We also analyze the relationship between the subsidence and the coastline change, river distribution, geological structure as well as the local terrain. The results show that (1) the land subsidence with the average velocity of 50 mm/year occurred in the low elevation area in the front part of the delta and the coastal area, and the area of the regions subsiding faster than 30 mm/year between 2006 and 2011 is larger than 122 km2; (2) the subsidence order and area estimated in this study are both much larger than that measured in previous studies; (3) the areas along rivers suffered from surface subsidence, due to the thick soft soil layer and frequent human interference; (4) the geological evolution is the intrinsic factor of the surface subsidence in the PRD, but human interference (reclamation, ground water extraction and urban construction) extends the subsiding area and increases the subsiding rate.


Introduction
Deltas have abundant natural resources, superior natural environment and geographical locations, so they accommodate the most intensive human activities. However, most deltas suffer from geological disasters caused by land subsidence, which damage buildings, railways, highways, subways and underground facilities (water supply pipes, gas and electrical installations), hindering the sustainable development in economy and environment (Ericson et al. 2006;Törnqvist et al. 2008;Wang et al. 2017). About 24 of the 33 major deltas in the world have surface subsidence, including the Pearl River Delta (PRD) in China (Syvitski et al. 2009). As one of the most economically developed and urbanized areas in China, the Pearl River Delta (PRD) has the largest urban agglomeration in the world and the population of 55 million. It has suffered many surface deformation disasters, including land subsidence, landslides and collapses (Syvitski et al. 2009;Chen et al. 2012;Ma et al. 2019).
The interferometric synthetic aperture radar (InSAR) technique has been widely used for monitoring surface subsidence in delta areas, including the PRD, because of its wide spatial and long temporal coverage (Higgins et al. 2014). Using C-band ENVISAT/ASAR data, Chen et al. (2012) mapped the surface deformation in the range of − 15 to 15 mm/year from 2006 to 2010. Based on the X-band and COSMO-SkyMed data of Guangzhou and Foshan, Alex et al. (2018) reported that the surface deformation was in the range of − 35 to 10 mm/year between 2011 and 2017. Ma et al. (2019) investigated the subsidence of the Guangdong-Hong Kong-Macao Greater Bay Area and found that the overall subsidence ranged from 0 to 112.3 mm/year, according to the 2015-2017 Sentinel-1 images. The results of these studies show similar subsidence patterns, but significantly different maximum subsidence orders (15 to 112.3 mm/year). Using a new kind of SAR dataset with long wavelength may address this controversy. The L-band ALOS1/PALSAR data have much longer wavelength than the C-band ENVISAT ASAR and X-band COSMO-SkyMed data, so they can detect much more reliable deformation information in densely vegetated regions, like PRD. Therefore, we use ALOS1/PALSAR data to assess the order of the subsidence in PRD during 2006PRD during -2011 The surface deformation evolution in delta areas is contributed by both natural factors and anthropogenic interference. For the PRD area, many related studies focused on the local surface deformation caused by human activities, including ground water extraction and engineering construction (Wang et al. 2012(Wang et al. , 2017Chen et al. 2012;Xu et al. 2016;Alex et al. 2018). These studies provided a lot of useful information for understanding the local subsidence pattern, but they paid little attention on the natural evolution factors, such as river distribution, coastline change and the local terrain.
In this study, we use the small baseline subset InSAR (SBAS-InSAR) method to map the surface deformation over the entire PRD (area: ~ 40,000 km 2 ) using 162 ALOS1/PALSAR images from 3 paths and 8 frames acquired between 2006 and 2011. The deformation results are verified by cross-validation of nearby path velocity rate and field investigation data. Then we analyze the temporal-spatial distribution characteristics of the deformation in the PRD, based on which we investigate the relationship between the subsidence and natural evolution factors as well as human interferences.

Geological background
The PRD is located in the south-central part of Guangdong Province, bordering the South China Sea (Fig. 1). It is a plain formed by the alluvium from the Xijiang River, Beijiang River and Dongjiang River (Huang and Zhang 2004). It covers an area of about 42, 200 km 2 from 112.5° to 114.25° E and 21.83° to 23.6° N. This delta is high in the north and low in the south, with small undulations. Influenced by the subtropical monsoon climate, the PRD has abundant rainfall and is often struck by typhoons and rainstorms. More than 300 river sections, with a total length of about 1740 km and an area of over 9750 km 2 , form a dense and complex network (Huang and Zhang 2004).
Bedrock (β), hidden karst (C), Quaternary loose soil (Q4) and soft soil form the geological structure of this area (Fig. 2). The bedrock is mainly composed of granite and partially of metamorphic sandstone, clastic rocks and basalt. The hidden karst, composed of limestone and shale, concentrates in the Guanghua Basin. The Quaternary loose soil and soft soil form alluvial layers, which are thick in the front part (such as Zhongshan city and Nansha district of Guangzhou city) and thin in the rear part of the delta. These alluvial layers were formed between the Pleistocene and the Holocene (Zong et al. 2009), with the thickness of most parts varying from 5 to 36 m, and the maximum as 55 m. These soft layers mostly distribute along the river and in the flat terrain of the delta (Li et al. 1991).

River and coastline evolution of the PRD
The rich rainfall and abundant river flow carry sediments from the upper reaches to the lower reaches. The PRD is flat, so 20-50% of the total sediment transported deposits in the river channel (Hu et al. 2011). Since the Tang dynasty, the population of the PRD started soaring. The agricultural development brought by the population growth accelerated the change of land use and led to serious soil and water loss. The increased sediment discharge contributed to the growth of the estuarine delta and outward extending of the coastline (McManus 2002;Maselli and Trincardi 2013). In the past 6000 years, the coastline of the PRD has been pushed outward 100 km (see Fig. 2) (Fyfe et al. 1999;Zong et al. 2009).
Human activities also play an important role in the morphology and geomorphology evolution of the PRD, especially in the last 100 years. The influence of human activities even exceeded that of natural factors (Ericson et al. 2006;Syvitski and Satio 2007). And some researchers, based on the remote sensing images and historical maritime map, found that the reclamation is the main factor of the PRD morphology and geomorphology changes over the past 100 years (Seto and Kaufmann 2006;Li and Damen 2010;Wang et al. 2013).

SBAS-InSAR data processing
Given a differential InSAR pair acquired at moments t A and t B , the unwrapped phase is where δφ obs is the phase of one pixel in differential interferogram, d is the line of sight (LOS) deformation of a pixel between two images, z is the topographic error, �φ atm is the atmospheric delay phase, and n j is the phase noise. Equation (1) can be written in a matrix form as follows: where v is the deformation rate. Following the SBAS-InSAR procedures (Berardino et al. 2002), we solve Eq. (2) by the least squares method, and obtain the linear deformation rate and the topographic error of each pixel. Subtracting topographic error and linear deformation phase from each differential interferogram, we get the residual phase. The residual phase consists of nonlinear deformation, atmospheric delay and noise. Taking the residual phase as the observation value, we obtain the time series phase of each SAR image by the singular value decomposition (SVD) method. We use the temporal high pass and spatial low pass filtering to obtain the atmospheric phase of each SAR image. Then, atmospheric errors are subtracted from each interferogram. Finally, the SVD method is used to calculate time series of each pixel. Figure 3 shows the flowchart of the SBAS-InSAR method.
In this study, the L-band ALOS1/PALSAR data from 3 tracks are used to obtain the subsidence that occurred in the PRD between 2006 and 2011 (see Fig. 1 and Appendix 1). The data are processed with the GAMMA software (Werner et al. 2003). The multi-look operation of 6 × 16 pixels in the range and azimuth directions is adopted to reduce phase noise and calculation burden before phase unwrapping. The result's resolution after multi-looking (1) is about 50 m × 50 m. The interferometric pairs with the temporal baseline < 700 days and perpendicular baseline < 1000 m are selected to reduce the decorrelation. The topographic phase is removed using the 1-arcsecond (~ 30 m) Shuttle Radar Topography Mission DEM (Jarvis et al. 2008). A modified Goldstein filter is used to reduce noise phases (Li et al. 2008) and a mask with a coherence less than 0.4 is adopted to weaken the effect of low coherence regions (sea, isolated islands, etc.). Finally, the phase unwrapping is performed by the minimum cost flow algorithm (Costantini 1998). The polynomial fitting is used to remove the long wavelength orbital errors in the unwrapped phase.

Surface subsidence results
Using the method described in "SBAS-InSAR data processing" section, we obtain the deformation rate map of the PRD between 2006 and 2011 (Fig. 4). Figure 4a presents the average deformation rate derived from the ALOS1/PALSAR image. The hinterland is stable and the subsidence is mainly distributed in the front part of the delta and the coastal area. The average subsidence rate ranges from 10 to 50 mm/year, with the maximum deformation rate of > 70 mm/year. Some regions have sporadic subsidence with a rate between 10 and 30 mm/year, such as the northeast part of Guangzhou city and the southeastern part of Zhaoqing city (outlined by the red dotted line in Fig. 4). Two large subsidence cones (S-1 and S-2 in Fig. 4b) are observed in the front part of the delta. The most serious land subsidence is observed in the junction area of Zhongshan city and Jiangmen city (such as Jianghai district of Jiangmen city and Henglan town of Zhongshan city), with an area of 160 km 2 . In this region, the subsidence rate of the center ranges from 30 to 50 mm/ year, even > 70 mm/year in some parts, and the average cumulative deformation is about 35 cm from 2006 to 2011 (Fig. 4d). Subsidence cone S-2 covers Jinwan district of Zhuhai city and Tanzhou town of Zhongshan city, with an area of 230 km 2 . Besides the two large subsidence cones, there are some areas with smaller subsidence rates, such as Doumen district of Zhuhai city and Xinhui district of Jiangmen city.
No large subsidence cone is found in the coastal areas. The subsidence areas scatter on the two sides of the Lingding Bay (see Fig. 4c), such as Nansha district of Guangzhou city and Nanlang town of Zhongshan city and Tangjia town of Zhuhai city. The subsidence rate is between 20 and 40 mm/year, and is over 50 mm/year in some regions. The subsidence rates in the coastal areas of Dongguan and Shenzhen are smaller. For example, the subsidence rates of Chang'an town in Dongguan and Bao'an district in Shenzhen are 10 and 25 mm/year, respectively.

Results verification
In the verification, we compared the observations of Area1 and Area2 during the same period, which are the overlapping areas of Path 461 and Path 460 (Fig. 4a). Area1 is stable but Area2 has obvious subsidence, so they can verify the results of different subsidence levels. As Fig. 4e, f show, the velocity difference in Area1 is between − 8 and 5 mm/year, and the overall mean value is − 1.3 mm/year with a variance of 2.2 mm/year. In Area2, the velocity difference is between − 10 and 6 mm/year, the mean value is − 0.4 mm/year, and the variance is 3.8 mm/year. Because of decorrelation, atmospheric delays and different incidence angles, the results of the overlapping areas of adjacent paths, like Area1 and Area2, are slightly different, so the uncertainty (2.2 mm/ year and 3.8 mm/year) between different paths is reasonable. To further verify the results, we did field investigations in some cities of the PRD, including Jiangmen city, Zhuhai city and Guangzhou city (see Fig. 5). The field observation confirmed the model results. We also reprocessed the data using the IPTA method (Werner et al. 2003). The details of data processing and result comparison are shown in Appendix 2. The result comparison also confirms that the order and patterns of surface deformation are reliable.

Discussion
Surface deformation in the PRD is usually caused by the joint force of natural geological evolution and human activities. Natural factors include sediment consolidation, tectonic activity, sea level rise and tidal load (Hu et al. 2018). Human activities are land reclamation, groundwater extraction, engineering construction and so on (Wu et al. 2016;Wang et al. 2017). Based on the subsidence results, we analyze the relationship between the subsidence and natural factors and human interference.

Surface subsidence of the PRD area between 2006 and 2011
As shown in Table 1, the average subsidence rate ranges from 10 to 50 mm/year, with the maximum rate > 70 mm/ year, and the regions with the subsidence rate > 30 mm/ year are as large as 122 km 2 .
Both the order and the area of surface subsidence we obtained are much larger than that of Chen et al. (2012) who used the ENVISAT ASAR datasets (2006)(2007)(2008)(2009)(2010) to map the surface subsidence of the PRD. In the PRD, the areas suffered rapid subsidence have Fig. 5 Photos of the field investigation in the areas denoted by the yellow stars in Fig. 4. a Xinhui port in Xinhui district, Jiangmen city (P-1 in Fig. 4a). The ground subsided about 5 cm relative to the steps; b the Kangdelai Medical city in Jinwan district, Zhuhai city (P-2 in Fig. 4a). The steps have separated from the building; c a building in Honghai village, Wanhaisha town, Nansha district, Guangzhou (P-3 in Fig. 4a). The underground pipeline has exposed because of the ground subsidence; d a house in Sanban village, Hongqi town, Jinwan district, Zhuhai city (P-4 in Fig. 4a). The windows have submerged into the ground decorrelation due to the reclamation, vegetation coverage and river distribution. ALOS1/PALSAR (L-band) have longer wavelength than ENVISAT ASAR (C-band), so their images have stronger deformation detection ability than the latter, and the detected maximum deformation is larger than that of the latter. Besides wavelength, the revisit cycle of SAR data is also very important for rapid surface deformation monitoring. The revisit cycle of the ENVISAT ASAR data is 35 days, but it is irregular in the PRD area, and sometimes is up to half a year. This may cause the high subsidence rate points in ENVISAT ASAR measurement missing or significantly less than that of ALOS1/ PALSAR data. Ma et al. (2019) used Sentinel-1A/B (2015-2017 to map the surface subsidence of the PRD, and got similar result with ours, which is much larger than that of Chen et al. (2012). The Sentinel-1A/B has a revisit cycle of 6-12 days, so their images have higher coherence and better rapid deformation detection ability than the ENVISAT ASAR dataset in this study. The consistency between our results and Ma's results indicates that the subsidence in the PRD is stable and experienced continuous and rapid surface subsidence from 2006 to 2017.

Natural factors River distribution and coastline evolution
Superimposing the river distribution map of the PRD on the deformation rate map (see Fig. 6a), we find that most subsidence areas are distributed along rivers, especially those in the front part of the PRD, where the river networks are dense (Chen et al. 2012). The coastline of the PRD continues pushing outward as silt keeps on depositing at the junction of land and sea. The geological structure formed by the highly compressible soil layers, however, is inherently unstable. The PRD is mainly composed of soft soil layers with high water content and compressibility, poor permeability, low bearing capacity and shear strength (Gebremichael et al. 2018), which easily leads to surface deformation. Most of the soft soil layers formed by sediment deposition are distributed along rivers, so subsidence usually occurs along rivers when those areas were developed by human activities and engineering construction.
The coastline formation time also has influences on the stability of the ground structure, as newly formed soil has a shorter consolidation period and may bring larger subsidence. The map of the coastline changes over the past 6000 years (see Fig. 6b) show that almost all the subsidence areas are distributed in the regions formed in the  Li et al. 1991) latest 1000 years. The land subsidence of the PRD has a high correlation with the distribution of Quaternary sediments (see Fig. 2). This has been validated in various delta regions of the world. For example, in the Fraser River Delta of Canada, the surface slowly sinks with an annual rate of 1-2 mm/year (Mazzotti et al. 2009), which is mainly caused by the consolidation of shallow Holocene sediments (soft soil). In the Mississippi Delta, the subsidence rate due to soil compaction is relatively large. The settlement speed of the soft soil area formed thousands of years ago is about 5 mm/year, and that of the area formed in the past tens to hundreds of years exceeds 10 mm/year (Törnqvist et al. 2008).

Geological structure
As shown in Fig. 7, the areas located in bedrock distribution areas (β) (see Fig. 7b) are very stable, and most surface subsidence occurred in the soft soil area. Before the modern delta was formed (6000 years ago), the bedrock areas were isolated islands, and most of them were composed of bare rocks. As the coastline moved outward, those islands gradually merged into the land. So, the geological structures of these regions are stable. As Fig. 7a, b shows, there is a strong positive correlation between the thickness of the soft soil and the deformation rate. Higgins et al. (2014) drew the similar conclusion in the study of the subsidence of Ganges-Brahmaputra Delta. Therefore, the soft soil thickness can be inverted by the correlation between surface deformation and drilling data.
Soft soil is mainly distributed along rivers and in the low-lying areas of the delta, where the erosion caused by sea level rise is serious. In Fig. 8a, the green area is lower than the average sea level, and it accounts for more than 1/3 of the entire onshore delta. The subsidence in the area below the sea level is outlined by the yellow line in Fig. 8b. The comparison between Fig. 8a, b shows that the land subsidence in the front and middle parts of the PRD concentrates in the area below the mean sea level. The annual average tidal range of the Pearl River Estuary is 1.0-1.7 m, and the load generated by tides and storm surges can cause settlements of a few millimeters to centimeters (Ye and Preiffer 1990). The low-lying areas, affected by tides, sea level rise, storms and floods, are also the areas accommodating severe subsidence disasters. According to the Intergovernmental Panel on Climate Change (IPCC), the sea level had risen 0.19 m (1.7 mm/ year) between 1901 and 2010. And the rise is accelerating (Zhang and Church 2012). These low-lying areas would suffer from both seawater erosion and surface subsidence in the future.

Human activities Reclamation
In the 1960s, the large-scale reclamation in the PRD rapidly reshaped the coastline. The reclamation rate between 1970 and 1995 reached 21 km 2 /year and then slowed down. We collect the coastline data from 1850 to 2015 with a time span of 165 years, and depict the changes in Fig. 9. The blue line is the PRD coastline in 1850, and the red lines are the coastlines of different periods. In the past 165 years, the coastline and landforms have undergone tremendous changes. A large area of land has been  reclaimed from waters and tidal flats. As the shoreline gradually moved towards the sea, more than 20 islands gradually merged into the mainland. By 2015, almost all the tidal flats in 1850 had been reclaimed into land.
As shown in Fig. 9, the reclaimed areas are also the area with the largest subsidence rate (subsidence rate: > 50 mm/year), including the coastal areas of Zhuhai, Guangzhou, Dongguan and Shenzhen. Most of these areas were reclaimed in the last 30 years (Xu et al. 2016;Huang and Zhang 2004), so the soil layers in these areas are unstable. Furthermore, the landfill materials are sediments, which causes subsidence by self-consolidation.
We select two areas, Nansha district of Guangzhou city and Jinwan district of Zhuhai city (outlined by the yellow dotted line in Fig. 9), to analyze the relationship between the subsidence and reclamation in detail. As shown in Fig. 10, these two areas were reclaimed from the sea a few decades ago. In 1984 (see Fig. 10a), Longxue and Hengmen were small islands, but they gradually expanded outwards and formed the basic shape of the current island in 2016 (see Fig. 10b-e). The area of Longxue island cofferdam grew about 150 times in 35 years. The newly reclaimed regions have serious subsidence, and the maximum subsidence rate is 75 mm/year (see Fig. 10f ). Figure 10g-k shows the variation of the reclaimed region in Jinwan district of Zhuhai city, where two small subsidence cones are obvious (see Fig. 10l). Gaolan island, which was an island before 2004, is connected to the land in 2010 (see Fig. 10j). The connection area accommodates the subsidence cones, where the maximum subsidence is 87 mm/year (see Fig. 10).
Many countries, such as the Netherlands (Hoeksema 2007), the United Kingdom (Wali 1975), Japan (Suzuki 2003) and Singapore (Douglas and Lawson 2003;Xu et al. 2016), have reclaimed land from the sea for the coastal industries and agricultural development, as well as for storm defense. Surface deformation monitoring in the Nile Delta and the Tiber Delta indicates that the reclaimed area is the most severely subsiding area (Stramondo et al. 2008;Marriner et al. 2012). So the long-term monitoring of the reclaimed area in the PRD is necessary.

Other human interventions
The PRD region is an economic center of south China, where factories, buildings, farmlands, and fisheries are densely distributed. The deposits consolidation and reclamation are the main factors for the land subsidence in deltas, but they alone cannot cause subsidence as high as 70 mm/year in small regions. We took the subsidence cone S-1 (see Fig. 4b) as an example to investigate the reason of local rapid subsidence. Through field investigation, we found that fish ponds are widely distributed in the coastal area of Jianghai (Fig. 11a), where large area land subsidence occurred. Aquaculture needs great amount of freshwater. Aggressive extraction of groundwater leads to the rapid drop of the groundwater level and causes subsidence, even settlement cones. Two profiles along S-1 (Fig. 11b), ab and cd, are used to show the settlement characteristics (Fig. 11c, d).
The subsidence caused by excessive groundwater extraction is common in delta areas. In the Mekong Delta, the average land subsidence rate has reached 16 mm/year as more than 250,000 ha of land has been converted from rice fields to shrimp ponds since 2000 (Erban et al. 2014). Oil and groundwater extraction play a key role in the land subsidence in the modern Yellow River Delta, and the land subsidence caused by aquaculture is up to 250 mm/year (Higgins et al. 2013). Therefore, we conclude that natural factors coupled with human activities cause the subsidence in the PRD. The evolution of the PRD determines the location of the subsidence areas, and human activities have become the direct cause of the settlement. Therefore, we should pay more attention to the subsidence caused by human activities.

Conclusions
We map the surface deformation of the entire PRD using the ALOS1/PALSAR images from 2006 to 2011. The results show that the areas with serious subsidence concentrate in the front part of the delta and the coastal areas. The two large settlement funnels in the Fig. 11 a Optical image of S-1 (Fig. 4). The insets are the zoom-ins of the base map, where the fish ponds are clearly visible. b Average displacement velocity map of S-1. c, d the average subsidence velocities of profiles ab and cd, respectively front part of the PRD have the average settlement rates of 30-50 mm/year, and the maximum settlement rate of larger than 70 mm/year. The surface deformation map generated from the InSAR observation shows that: (1) the subsidence area is distributed along rivers; (2) the natural evolution would be the intrinsic factor for the land subsidence of the PRD; (3) the land subsidence in the PRD is significantly related to human activities (such as land reclamation, pumping groundwater for aquaculture). Due to the natural and human factors, the settlement in the PRD region will continue. These findings have guiding contributions to disaster prevention and social construction in the PRD.

PRD: Pearl River Delta.
In order to verify the reliability of our result, we processed the ALOS1/PALSAR data with the IPTA-InSAR method and compared the results with our results in Fig. 12, which shows that the two results have good consistency. This comparison also confirms the rapid subsidence in the PRD monitored by the proposed method.