Understanding Land Subsidence Along the Coastal Areas of Guangdong, China, by Analyzing Multi-Track MTInSAR Data

: Coastal areas are usually densely populated, economically developed, ecologically dense, and subject to a phenomenon that is becoming increasingly serious, land subsidence. Land subsidence can accelerate the increase in relative sea level, lead to a series of potential hazards, and threaten the stability of the ecological environment and human lives. In this paper, we adopted two commonly used multi-temporal interferometric synthetic aperture radar (MTInSAR) techniques, Small baseline subset (SBAS) and Temporarily coherent point (TCP) InSAR, to monitor the land subsidence along the entire coastline of Guangdong Province. The long-wavelength L-band ALOS / PALSAR-1 dataset collected from 2007 to 2011 is used to generate the average deformation velocity and deformation time series. Linear subsidence rates over 150 mm / yr are observed in the Chaoshan Plain. The spatiotemporal characteristics are analyzed and then compared with land use and geology to infer potential causes of the land subsidence. The results show that (1) subsidence with notable rates ( > 20 mm / yr) mainly occurs in areas of aquaculture, followed by urban, agricultural, and forest areas, with percentages of 40.8%, 37.1%, 21.5%, and 0.6%, respectively; (2) subsidence is mainly concentrated in the compressible Holocene deposits, and clearly associated with the thickness of the deposits; and (3) groundwater exploitation for aquaculture and agricultural use outside city areas is probably the main cause of subsidence along these coastal areas.


Introduction
Guangdong Province is located on the South China Sea coast and has the longest coastline of any province in China, approximately 4300 km. Along this coastline, several large agglomerations (e.g., the Pearl River Delta (PRD) and Chao Shan Plain (CSP)), harbors (Zhujiang, Zhanjiang, and Shantou), and a wide range of aquaculture and agricultural areas are clustered, controlling the primary economic activity and ecological environment. Over the past few decades, this coastline has experienced rapid growth in both population and economy. The PRD has undergone a shift from an agricultural-based economy to an industrial-and technological-based economy, making the PRD Asia's fourth-largest economy [1,2]. To satisfy the needs of urban residents and export trade, an increasing amount of land is exploited through the expansion of aquaculture land use and reclamation of seashore. Among these exploitations, large-scale freshwater aquaculture has made Guangdong Province one of the

SAR Datasets and Landsat Datasets
To generate the surface motions of the whole Guangdong coastal area, we collected 253 archived SAR image datasets acquired from 2007 to 2011 by L-band ALOS1/PALSAR satellites, including 16 adjacent orbits covering an area of approximately 49,000 km 2 (represented by blue rectangles in Figure 1). Table 1 gives the imaging parameters of each orbit, including the number, time span, and systematic parameters. An average of 12 acquisitions per frame was used during the 3.5-year period.

SAR Datasets and Landsat Datasets
To generate the surface motions of the whole Guangdong coastal area, we collected 253 archived SAR image datasets acquired from 2007 to 2011 by L-band ALOS1/PALSAR satellites, including Remote Sens. 2020, 12, 299 4 of 20 16 adjacent orbits covering an area of approximately 49,000 km 2 (represented by blue rectangles in Figure 1). Table 1 gives the imaging parameters of each orbit, including the number, time span, and systematic parameters. An average of 12 acquisitions per frame was used during the 3.5-year period. In addition to the active remote sensing datasets, we also collected 8 optical Landsat images (represented by red dotted rectangles in Figure 1, see Table 2), aiming to generate land classification maps. Although a 2-year time span between SAR images and Landsat images exists, however, the changing rate of the agricultural land, forest land, and aquaculture land are only about 8.29%, 2.09%, and 2.99% respectively, according to the Statistics Bureau of Guangdong Province, see Table S1 [40]. Moreover, Landsat 8 image has a superior spectral resolution compared with Landsat 5 and 7, which is launched on 2013 [41,42].

In Situ Datasets
In addition to the remote sensing datasets, in situ datasets were collected in this study, including time-series datasets of 7 GPS stations and a dataset of Quaternary deposit thicknesses derived from boreholes (represented by black stars in Figure 1 and black circles in Figure 2, respectively). The GPS data were used to validate the InSAR-derived results, and the borehole data were adopted to generate the 2D deposit maps (Figure 2).

Methodology
MTInSAR is developed based on DInSAR, aiming to extract the deformation from mixed phases of series differential interferograms. That is, the residual topographic phase ( ), deformation phase ( ), orbit phase ( ), atmospheric phase ( ) and noise ( ) constitute the differential interferometric phase ( ).
where is a phase wrap operation. To generate high-precision deformation rate and timeseries, MTInSAR is mainly focus on the elimination of atmospheric phase and decorrelation noise based on wrapped or unwrapped phases. In this case, i.e., coastal areas with large coverage and serious temporal decorrelation, multi-master MTInSAR with a large multilooking number (3 and 8 looks in the range and azimuth directions, respectively) is adopted. Specifically, considering the complicated data conditions, such as frames with limited SAR image numbers and serious atmospheric effects in the coastal area, two kinds of MTInSAR are adopted. The first is SBAS-InSAR based on unwrapped phases, which is appropriate for most situations. The second method is an arcbased method, TCPInSAR, which is based on wrapped phases and aims to eliminate the influence of serious atmospheric effects. Considering the data processing efficiency, the latter is adopted in the case that the SAR image number is limited and the atmospheric effects are serious (tracks 453, 457, 462 and 463), see the comparison between pixel-based SBAS-InSAR and arc-based TCPInSAR in Figure S1. After generating deformation over the 16 orbits, the average velocity maps in the line of sight (LOS) direction are mosaiced after removing the relative offsets in the overlapped areas by a linear polynomial fitting [43], considering the approximate incident angle (38.7°, see Supplementary  Table S1). Moreover, the precision of the InSAR results is evaluated, including relative precision calculation by the velocity difference in the overlapped area of adjacent orbits and analysis of the absolute precision with respect to GPS measurements. Next, the land-use classification maps along the coastal areas are obtained based on the selected optical Landsat images. Finally, the deformation features are observed, and the causes of the deformation are analyzed qualitatively and quantitatively. The flowchart is shown in Figure 3.

Methodology
MTInSAR is developed based on DInSAR, aiming to extract the deformation from mixed phases of series differential interferograms. That is, the residual topographic phase (ϕ topo ), deformation phase (ϕ de f ), orbit phase (ϕ orb ), atmospheric phase (ϕ aps ) and noise (ϕ noi ) constitute the differential interferometric phase (ϕ Dint ).
where W ϕ Dint is a phase wrap operation. To generate high-precision deformation rate and time-series, MTInSAR is mainly focus on the elimination of atmospheric phase and decorrelation noise based on wrapped or unwrapped phases. In this case, i.e., coastal areas with large coverage and serious temporal decorrelation, multi-master MTInSAR with a large multilooking number (3 and 8 looks in the range and azimuth directions, respectively) is adopted. Specifically, considering the complicated data conditions, such as frames with limited SAR image numbers and serious atmospheric effects in the coastal area, two kinds of MTInSAR are adopted. The first is SBAS-InSAR based on unwrapped phases, which is appropriate for most situations. The second method is an arc-based method, TCPInSAR, which is based on wrapped phases and aims to eliminate the influence of serious atmospheric effects. Considering the data processing efficiency, the latter is adopted in the case that the SAR image number is limited and the atmospheric effects are serious (tracks 453, 457, 462 and 463), see the comparison between pixel-based SBAS-InSAR and arc-based TCPInSAR in Figure S1. After generating deformation over the 16 orbits, the average velocity maps in the line of sight (LOS) direction are mosaiced after removing the relative offsets in the overlapped areas by a linear polynomial fitting [43], considering the approximate incident angle (38.7 • , see Supplementary Table S1). Moreover, the precision of the InSAR results is evaluated, including relative precision calculation by the velocity difference in the overlapped area of adjacent orbits and analysis of the absolute precision with respect to GPS measurements. Next, the land-use classification maps along the coastal areas are obtained based on the selected optical Landsat images. Finally, the deformation features are observed, and the causes of the deformation are analyzed qualitatively and quantitatively. The flowchart is shown in Figure 3. SBAS-InSAR, developed by Berardino et al. [16] and improved by researchers in recent years, is a pixel-based approach that has been widely used for wide-ranging deformation monitoring [16,26]. In this paper, GAMMA software is used to generate the unwrapped differential interferometric phases. SAR images of each track are first coregistered with a master image, as shown in Table 1. Interferograms with a given spatial and temporal baseline threshold (850 m and 800 days respectively), aiming to reduce the influence of decorrelation and digital elevation map (DEM) error on parameter estimation. After interferograms selection, traditional DInSAR with a multilook operation (three looks and eight looks in range and azimuth direction respectively) is applied to generate the unwrapped phases, including flatten and topographic phase removing, goldstein filtering, MCF phase unwrapping. The orbit error is removed with a polynomial fitting. Then, pixels with a mean coherence above 0.6 will be selected for the following processing procedures, including residual topographic estimation, atmosphere elimination, deformation time series calculation, and average velocity calculation. Here, we used a modified method to reduce the influence of residual topography caused by the correlation of the perpendicular baseline of ALOS1 in time dimension [44]. A spatial filtering and a temporal filtering are used to eliminate the atmosphere with a window-size of 500 m and two years respectively. A bootstrap with a parameter of 500 is adopted to generate the uncertainty of the average velocity fitting, and pixels whose uncertainty exceeds three times the standard deviation (std) are deleted. We also deleted those pixels whose coherent observation number is below two thirds of the total number.
Unlike SBAS-InSAR, TCPInSAR is a method proposed based on the wrapped phase of arcs [21,45,46]. It selected temporarily-coherence points (TCPs) and the unknown parameters, i.e., deformation rate, orbit phase, residual topographic phase, can be calculated without phase unwrapping operation. TCPInSAR proposed a local Delaunay networks with a uniform sampling in the range and azimuth directions, aiming to improve the arc density of the whole image. Unknown parameters are solved with an orbit errors joint model [46] and a phase ambiguity detection model [21] based on phase difference of arcs. Here, we used a sampling spacing of 700 m in both directions. The average velocity and residual topography of arcs are solved with an ambiguity detection method, that is, arcs with residuals over a given threshold are deleted (1.5 rad is used in this case). After generating the parameters of each arc, the parameters in each pixel are calculated based on an L1 norm least-squares estimation. The nonlinear deformation is obtained after spatial-temporal filtering SBAS-InSAR, developed by Berardino et al. [16] and improved by researchers in recent years, is a pixel-based approach that has been widely used for wide-ranging deformation monitoring [16,26]. In this paper, GAMMA software is used to generate the unwrapped differential interferometric phases. SAR images of each track are first coregistered with a master image, as shown in Table 1. Interferograms with a given spatial and temporal baseline threshold (850 m and 800 days respectively), aiming to reduce the influence of decorrelation and digital elevation map (DEM) error on parameter estimation. After interferograms selection, traditional DInSAR with a multilook operation (three looks and eight looks in range and azimuth direction respectively) is applied to generate the unwrapped phases, including flatten and topographic phase removing, goldstein filtering, MCF phase unwrapping. The orbit error is removed with a polynomial fitting. Then, pixels with a mean coherence above 0.6 will be selected for the following processing procedures, including residual topographic estimation, atmosphere elimination, deformation time series calculation, and average velocity calculation. Here, we used a modified method to reduce the influence of residual topography caused by the correlation of the perpendicular baseline of ALOS1 in time dimension [44]. A spatial filtering and a temporal filtering are used to eliminate the atmosphere with a window-size of 500 m and two years respectively. A bootstrap with a parameter of 500 is adopted to generate the uncertainty of the average velocity fitting, and pixels whose uncertainty exceeds three times the standard deviation (std) are deleted. We also deleted those pixels whose coherent observation number is below two thirds of the total number.
Unlike SBAS-InSAR, TCPInSAR is a method proposed based on the wrapped phase of arcs [21,45,46]. It selected temporarily-coherence points (TCPs) and the unknown parameters, i.e., deformation rate, orbit phase, residual topographic phase, can be calculated without phase unwrapping operation. TCPInSAR proposed a local Delaunay networks with a uniform sampling in the range and azimuth directions, aiming to improve the arc density of the whole image. Unknown parameters are solved with an orbit errors joint model [46] and a phase ambiguity detection model [21] based on phase difference of arcs. Here, we used a sampling spacing of 700 m in both directions. The average velocity and residual topography of arcs are solved with an ambiguity detection method, that is, arcs with residuals over a given threshold are deleted (1.5 rad is used in this case). After generating the parameters of each arc, the parameters in each pixel are calculated based on an L1 norm least-squares Remote Sens. 2020, 12, 299 7 of 20 estimation. The nonlinear deformation is obtained after spatial-temporal filtering for the residual phase. The final deformation time series is generated by combining the linear and nonlinear deformation.
Usually, 3D deformation, that is, deformation in both the vertical and horizontal directions, is calculated by the combination of ascending and descending images, which can help us infer hazard mechanisms. However, in this case, only descending SAR images were collected. Furthermore, according to the existing research, the horizontal deformation in the coastal areas is not obvious [47]. Hence, we assume that the horizontal displacement is negligible, deformation in the vertical direction can be generated with d v = d los /cosθ, where d v and d los are the displacement in the vertical and LOS) directions, and θ is the incidence angle. The velocity values referred to in the following sections are in the LOS direction unless otherwise specified. In addition to the deformation monitoring techniques, an object-oriented nearest neighbor classification method based on multiresolution segmentation and supervised classification is also adopted for generating the land-use classification map. The details of the object classification method, which are beyond the scope of this study, can be found in previously published research [48]. As described in Section 2.2.1, because of the low change rate of land-use and the superiority of Landsat 8, the comparison between subsidence and land-use classification map is appropriate, which can help us understand the mechanism of land subsidence. The accuracy of the classification map is compared with a 30-m resolution product (with an overall accuracy of 97.29% within the range of [20 • -30 • N, 60 • -120 • E]) obtained in 2015 derived from the National Science & Technology Infrastructure of China [49,50]. The overall accuracy of the classification map in this study is about 89.5%.

Validation of InSAR Results
To quantitatively evaluate the accuracy of the InSAR-derived results, relative and absolute accuracy methods are adopted [51][52][53]. The former is also called inner precision, which is defined as the velocity difference in the overlapping areas of adjacent tracks, while the absolute accuracy is defined as the consistency with in situ measurements. Here, we used both methods to verify the deformation accuracy. Because the occurrence frequency of the large-scale terrain deformation in the coastal areas is very low, the relative accuracy is calculated by the velocity consistency in the overlapping areas between adjacent orbits after a 1 st order polynomial operation [43]. The histograms of the average velocity difference between adjacent orbits are shown in Figure 4, and the mean and std values are shown in Table 3. It is obvious that the reference precision of velocity (std value of velocity difference) is generally less than 3.1 mm/yr, with the exception of several values approximate to 4 mm/yr. The causes of these high values may be the limited number of SAR images [54].
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 20 for the residual phase. The final deformation time series is generated by combining the linear and nonlinear deformation. Usually, 3D deformation, that is, deformation in both the vertical and horizontal directions, is calculated by the combination of ascending and descending images, which can help us infer hazard mechanisms. However, in this case, only descending SAR images were collected. Furthermore, according to the existing research, the horizontal deformation in the coastal areas is not obvious [47]. Hence, we assume that the horizontal displacement is negligible, deformation in the vertical direction can be generated with = / , where and are the displacement in the vertical and LOS) directions, and is the incidence angle. The velocity values referred to in the following sections are in the LOS direction unless otherwise specified. In addition to the deformation monitoring techniques, an object-oriented nearest neighbor classification method based on multiresolution segmentation and supervised classification is also adopted for generating the land-use classification map. The details of the object classification method, which are beyond the scope of this study, can be found in previously published research [48]. As described in Section 2.2.1, because of the low change rate of land-use and the superiority of Landsat 8, the comparison between subsidence and land-use classification map is appropriate, which can help us understand the mechanism of land subsidence. The accuracy of the classification map is compared with a 30-m resolution product (with an overall accuracy of 97.29% within the range of [20°-30°N, 60°-120°E]) obtained in 2015 derived from the National Science & Technology Infrastructure of China [49,50]. The overall accuracy of the classification map in this study is about 89.5%.

Validation of InSAR Results
To quantitatively evaluate the accuracy of the InSAR-derived results, relative and absolute accuracy methods are adopted [51][52][53]. The former is also called inner precision, which is defined as the velocity difference in the overlapping areas of adjacent tracks, while the absolute accuracy is defined as the consistency with in situ measurements. Here, we used both methods to verify the deformation accuracy. Because the occurrence frequency of the large-scale terrain deformation in the coastal areas is very low, the relative accuracy is calculated by the velocity consistency in the overlapping areas between adjacent orbits after a 1 st order polynomial operation [43]. The histograms of the average velocity difference between adjacent orbits are shown in Figure 4, and the mean and std values are shown in Table 3. It is obvious that the reference precision of velocity (std value of velocity difference) is generally less than 3.1 mm/yr, with the exception of several values approximate to 4 mm/yr. The causes of these high values may be the limited number of SAR images [54].
In addition to the statistical offsets in the overlapped areas, we also compared the velocity and deformation time series with the data collected from GPS stations, as shown in Figures 4b and 5. It is noted that the velocity comparison between InSAR and GPS is in the vertical direction. The correlation of InSAR and GPS velocities is 0.71, confirming the high precision of InSAR results. Moreover, the std of deformation time series between the two measurements is also calculated, as shown in Figure 5. The GPS stations FOMO and HKSL are selected as reference stations in tracks 460 and 459 for the deformation time-series comparison, respectively, due to their different acquisition dates of SAR image series (see Table 1) and different coverage. The results show that the std of the deformation time series for stations DSMG, HKFN, HKKT, HKLT, and HKST are 5.5 mm, 8.2 mm, 6.8 mm 8.7 mm and 8.7 mm, respectively. The millimeter-level accuracy of the deformation series also proves the reliability of the InSAR-derived results.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 20 Table 3. Mean and standard deviation (std) of the average deformation rate difference in the overlapping areas of adjacent tracks.

Adjacent Tracks Mean Std Adjacent Tracks Mean Std
In addition to the statistical offsets in the overlapped areas, we also compared the velocity and deformation time series with the data collected from GPS stations, as shown in Figures 4b and 5. It is noted that the velocity comparison between InSAR and GPS is in the vertical direction. The correlation of InSAR and GPS velocities is 0.71, confirming the high precision of InSAR results. Moreover, the std of deformation time series between the two measurements is also calculated, as shown in Figure 5. The GPS stations FOMO and HKSL are selected as reference stations in tracks 460 and 459 for the deformation time-series comparison, respectively, due to their different acquisition dates of SAR image series (see Table 1) and different coverage. The results show that the std of the deformation time series for stations DSMG, HKFN, HKKT, HKLT, and HKST are 5.5 mm, 8.2 mm, 6.8 mm 8.7 mm and 8.7 mm, respectively. The millimeter-level accuracy of the deformation series also proves the reliability of the InSAR-derived results.

Results
The surface displacement map of the Guangdong coastal areas derived from MTInSAR is calculated and shown in Figure 6. Positive values (yellow colors) represent movement toward the satellite (uplift). Negative values (purple colors) indicate movement away from the satellite (subsidence). Selected areas with concentrated land subsidence rates in excess of 20 mm/yr and discrete areas with gentle deformation are analyzed to grasp the features and characteristics of surface deformation. Three obvious subsidence areas in LZP, PRD, and CSP are identified. In addition, a land use classification map of the coastal areas of Guangdong Province is also generated and shown in Figure 7. Aquaculture, agriculture, urban, forest, and water are extracted based on the collected Landsat images.

Results
The surface displacement map of the Guangdong coastal areas derived from MTInSAR is calculated and shown in Figure 6. Positive values (yellow colors) represent movement toward the satellite (uplift). Negative values (purple colors) indicate movement away from the satellite (subsidence). Selected areas with concentrated land subsidence rates in excess of 20 mm/yr and discrete areas with gentle deformation are analyzed to grasp the features and characteristics of surface deformation. Three obvious subsidence areas in LZP, PRD, and CSP are identified. In addition, a land use classification map of the coastal areas of Guangdong Province is also generated and shown in Figure 7. Aquaculture, agriculture, urban, forest, and water are extracted based on the collected Landsat images.

Results
The surface displacement map of the Guangdong coastal areas derived from MTInSAR is calculated and shown in Figure 6. Positive values (yellow colors) represent movement toward the satellite (uplift). Negative values (purple colors) indicate movement away from the satellite (subsidence). Selected areas with concentrated land subsidence rates in excess of 20 mm/yr and discrete areas with gentle deformation are analyzed to grasp the features and characteristics of surface deformation. Three obvious subsidence areas in LZP, PRD, and CSP are identified. In addition, a land use classification map of the coastal areas of Guangdong Province is also generated and shown in Figure 7. Aquaculture, agriculture, urban, forest, and water are extracted based on the collected Landsat images.

Subsidence in the Leizhou Peninsula
The LZP is located in southwest Guangdong Province. Crops are grown in more than 90% of the peninsula (see the yellow areas in Figure 7), making the LZP the granary for Guangdong Province. Following the SBAS-InSAR procedure, the result in the LZP is shown in Figure 8. The high subsidence (>20 mm/yr) in LZP is concentrated in the two bays (Zhanjiang Bay and Leizhou Bay) and the southwest of Donghai Island, as outlined by the black dotted circle in Figure 8. The subsidence area is approximately 251.5 km 2 , with an average subsidence rate of 15.1 mm/yr, and the maximum subsidence rate of 58.0 mm/yr occurs on Donghai Island. The location of these subsidence bowls is consistent with the trend of middle groundwater decline [4]. In addition to the high rates of subsidence, some other subsidence bowls with areas over 4 km 2 , located along the western coastline, are also observed. The corresponding optical images are also given, showing that obvious subsidence bowls are located in areas with aquaculture land use. Moreover, the deformation time-series of two selected points in aquaculture areas are shown in Figure 8b,c, with an approximately linear trend. Discrete subsidence with an average rate of 7.0 mm/yr is presented in inland LZP, such as the subsidence areas of Suixi and Leizhou, with a nonlinear deformation time-series trend, as shown in Figure 8e. In addition to subsidence in agricultural and aquaculture areas, small-scale displacement (with an average subsidence rate of 9.8 mm/yr) is also observed in urban areas of Zhanjiang city, especially in its industrial areas (see inset F in Figure 8a,d). In addition to the subsidence areas, a small-scale uplift (~3 mm/yr) is observed in Dengjiaolou, which is consistent with the results of traditional measurements from tectonic activity [47].

Subsidence in the Leizhou Peninsula
The LZP is located in southwest Guangdong Province. Crops are grown in more than 90% of the peninsula (see the yellow areas in Figure 7), making the LZP the granary for Guangdong Province. Following the SBAS-InSAR procedure, the result in the LZP is shown in Figure 8. The high subsidence (>20 mm/yr) in LZP is concentrated in the two bays (Zhanjiang Bay and Leizhou Bay) and the southwest of Donghai Island, as outlined by the black dotted circle in Figure 8. The subsidence area is approximately 251.5 km 2 , with an average subsidence rate of 15.1 mm/yr, and the maximum subsidence rate of 58.0 mm/yr occurs on Donghai Island. The location of these subsidence bowls is consistent with the trend of middle groundwater decline [4]. In addition to the high rates of subsidence, some other subsidence bowls with areas over 4 km 2 , located along the western coastline, are also observed. The corresponding optical images are also given, showing that obvious subsidence bowls are located in areas with aquaculture land use. Moreover, the deformation time-series of two selected points in aquaculture areas are shown in Figure 8b,c, with an approximately linear trend. Discrete subsidence with an average rate of 7.0 mm/yr is presented in inland LZP, such as the subsidence areas of Suixi and Leizhou, with a nonlinear deformation time-series trend, as shown in Figure 8e. In addition to subsidence in agricultural and aquaculture areas, small-scale displacement (with an average subsidence rate of 9.8 mm/yr) is also observed in urban areas of Zhanjiang city, especially in its industrial areas (see inset F in Figure 8a,d). In addition to the subsidence areas, a small-scale uplift (~3 mm/yr) is observed in Dengjiaolou, which is consistent with the results of traditional measurements from tectonic activity [47].

Subsidence in the Pearl River Delta
The spatial-temporal characteristic of subsidence in the PRD is much different in magnitude and range than that of the LZP, as shown in Figure 9. This result indicates that the large-scale subsidence areas with an average rate in excess of 35 mm/yr reach a total area of 263.37 km 2 , in comparison to 15.16 km 2 in LZP during the approximate period. The subsidence area in the PRD is mainly distributed along the Pearl River tributary (see Figure 9). Two large-scale subsidence bowls are observed along the Modaomen channel (outlined by B1 and B2 in Figure 9), with average subsidence rates of 32.0 mm/yr and 25.1 mm/yr, respectively, and areas of 265.43 km 2 and 387.32 km 2 , respectively. The statistical area is calculated after kriging interpolation due to the temporal decorrelation of the surrounding water. Subsidence bowls with low subsidence rates are also observed along the other three channels: the Hengmen channel, Bingqili channel, and Jiaomen channel. In addition, local subsidence areas with notable rates are distributed along the coastline, especially in reclamation areas (which are outlined between the blue and red lines during the period from 1984 to 2011 in Figure 9). For example, an artificial island located in Nansha district and a reclamation area located in Jinwan district (see dotted purple circles in Figure 9), have maximum subsidence rates of 71.9 mm/yr and 150.9 mm/yr, respectively.
Notably, some small-scale uplift is found in mountainous areas (see the black dotted circles in Figure 9). This uplift may be caused by elevation-related atmospheric effects or residual topography. Besides the average deformation map, deformation time-series of four selected points are also generated (see Figure 9b-e). The linear trend is found in both agricultural and aquaculture areas (see Figure 9f-i).

Subsidence in the Pearl River Delta
The spatial-temporal characteristic of subsidence in the PRD is much different in magnitude and range than that of the LZP, as shown in Figure 9. This result indicates that the large-scale subsidence areas with an average rate in excess of 35 mm/yr reach a total area of 263.37 km 2 , in comparison to 15.16 km 2 in LZP during the approximate period. The subsidence area in the PRD is mainly distributed along the Pearl River tributary (see Figure 9). Two large-scale subsidence bowls are observed along the Modaomen channel (outlined by B1 and B2 in Figure 9), with average subsidence rates of 32.0 mm/yr and 25.1 mm/yr, respectively, and areas of 265.43 km 2 and 387.32 km 2 , respectively. The statistical area is calculated after kriging interpolation due to the temporal decorrelation of the surrounding water. Subsidence bowls with low subsidence rates are also observed along the other three channels: the Hengmen channel, Bingqili channel, and Jiaomen channel. In addition, local subsidence areas with notable rates are distributed along the coastline, especially in reclamation areas (which are outlined between the blue and red lines during the period from 1984 to 2011 in Figure 9). For example, an artificial island located in Nansha district and a reclamation area located in Jinwan district (see dotted purple circles in Figure 9), have maximum subsidence rates of 71.9 mm/yr and 150.9 mm/yr, respectively. Notably, some small-scale uplift is found in mountainous areas (see the black dotted circles in Figure 9). This uplift may be caused by elevation-related atmospheric effects or residual topography. Besides the average deformation map, deformation time-series of four selected points are also generated (see Figure 9b-e). The linear trend is found in both agricultural and aquaculture areas (see Figure 9f-i).

Subsidence in the Chaoshan Plain
The CSP consists of three river alluvial plains, the Lianjiang River Plain (LJP), Rongjiang River Plain (RJP) and Hanjiang River Plain (HJP), represented by purple, black and red dotted circles in

Subsidence in the Chaoshan Plain
The CSP consists of three river alluvial plains, the Lianjiang River Plain (LJP), Rongjiang River Plain (RJP) and Hanjiang River Plain (HJP), represented by purple, black and red dotted circles in Figure 10, respectively. The total area of these plains reaches 2600 km 2 . The displacement results in the CSP are presented in Figure 10 and exhibit a specific spatial-temporal subsidence feature that, compared with the results of LZP and PRD, is very large and intensive. As shown in Figure 10, a large-scale subsidence bowl with an obvious boundary is observed in the LJP, outlined by B1. The area and average subsidence rate are 342.21 km 2 and 41 mm/yr, respectively. The maximum subsidence rate reaches 157.2 mm/yr. In addition to the subsidence in the LJP, some discrete small-scale subsidence areas are observed in the HJP and are marked by black circles in Figure 10, i.e., B2. The subsidence is mainly concentrated along the western side of the HJP, and the maximum subsidence rate reaches 32.9 mm/yr. In the RJP, a small subsidence area is also generated in Jieyang city, with a maximum subsidence rate of 24.2 mm/yr. Moreover, four points located in the mixed agricultural and urban areas (see Figure 10f-i) are selected to generate the deformation time-series, as shown in Figure 10b-e. The linear trend of the deformation time-series is also found.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 20 Figure 10, respectively. The total area of these plains reaches 2600 km 2 . The displacement results in the CSP are presented in Figure 10 and exhibit a specific spatial-temporal subsidence feature that, compared with the results of LZP and PRD, is very large and intensive. As shown in Figure 10, a large-scale subsidence bowl with an obvious boundary is observed in the LJP, outlined by B1. The area and average subsidence rate are 342.21 km 2 and 41 mm/yr, respectively. The maximum subsidence rate reaches 157.2 mm/yr. In addition to the subsidence in the LJP, some discrete smallscale subsidence areas are observed in the HJP and are marked by black circles in Figure 10, i.e., B2. The subsidence is mainly concentrated along the western side of the HJP, and the maximum subsidence rate reaches 32.9 mm/yr. In the RJP, a small subsidence area is also generated in Jieyang city, with a maximum subsidence rate of 24.2 mm/yr. Moreover, four points located in the mixed agricultural and urban areas (see Figure 10f-i) are selected to generate the deformation time-series, as shown in Figure 10b-e. The linear trend of the deformation time-series is also found.

Subsidence in other Coastal Areas
In addition to the three selected subsidence areas, other subsidence areas are distributed discretely along the coastline, especially in river estuaries, for example, Hailing Bay and the estuary of the Moyang River in track 463, represented by the insets A and B in Figure 6, respectively, the estuary of Guanghai Bay in track 461 (see inset D in Figure 6), and the estuary of the Huangjiang River in track 456 (inset G in Figure 6). In addition to the subsidence occurring along the river system, some local subsidence is also observed in forest areas, such as the discrete subsidence bowls in tracks T463 and T457 (see the insets C and F in Figure 6). Moreover, elongated subsidence located in the industrial petrochemical region is also observed in track 458 (see inset E in Figure 6).

Subsidence in Other Coastal Areas
In addition to the three selected subsidence areas, other subsidence areas are distributed discretely along the coastline, especially in river estuaries, for example, Hailing Bay and the estuary of the Moyang River in track 463, represented by the insets A and B in Figure 6, respectively, the estuary of Guanghai Bay in track 461 (see inset D in Figure 6), and the estuary of the Huangjiang River in track 456 (inset G in Figure 6). In addition to the subsidence occurring along the river system, some local subsidence is also observed in forest areas, such as the discrete subsidence bowls in tracks T463 and T457 (see the insets C and F in Figure 6). Moreover, elongated subsidence located in the industrial petrochemical region is also observed in track 458 (see inset E in Figure 6).

A positive Correlation between Sedimentary Thickness and Subsidence
We focused on a geological formation that is usually prone to natural hazards, namely, the Holocene series (Q4), to investigate the correlation of these two factors. As shown in Figure 1, in addition to some discrete distribution along the coastline, Q4 is mainly distributed in the LZP, PRD, and CSP. The areas of subsidence have a distribution similar to that of the Q4 deposits, especially in areas with notable subsidence rates. To clearly analyze this similarity, we calculated the histograms of the subsidence rate (>15 mm/yr) within and out of areas with Q4 deposits; the results are shown in Figure 11a. The percentage along the Y axis is calculated by the following equations: Rer Q4 = pixel number whose subsidence rate>15mm/yr within Q4 area sum total pixel number within Q4 area Rer no−Q4 = pixel number whose subsidence rate>15mm/yr out o f Q4 area sum total pixel number out o f Q4 area Subsidence is mainly located within the Q4 area, especially in areas with subsidence rates higher than 40 mm/yr (see the blue columns in Figure 11a).
To quantitatively analyze the correlation of these two factors, we also considered the thickness of Q4 deposits in the LZP, PRD, and CSP (see the Q4 depth map in Figure 2). The Q4 thickness maps of the LZP and CSP are generated by a kriging interpolation based on the depths collected from borehole measurements, which are represented by black circles [55]. We divided the thickness of Q4 into several sections, for example, 0 m to 10 m with an interval of 5 m, 10 m to 50 m with an interval of 10 m, and 50 m to 170 m with an interval of 20 m. Additionally, we calculated the average subsidence rate in these sections across the LZP, CSP and PRD. These results are represented by red diamonds in Figure 11b-d, respectively. We can see that an obvious positive correlation between Q4 thickness and subsidence rate is observed, that is, the greater the thickness is, the higher the subsidence rate, especially in the LZP and PRD. However, in the CSP, a positive correlation occurs when the thickness is larger than 50 m, and a negative correlation is generated when the thickness is low. This phenomenon may be due to the combined impact of soft soil compaction and groundwater recession.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 20 We focused on a geological formation that is usually prone to natural hazards, namely, the Holocene series (Q4), to investigate the correlation of these two factors. As shown in Figure 1, in addition to some discrete distribution along the coastline, Q4 is mainly distributed in the LZP, PRD, and CSP. The areas of subsidence have a distribution similar to that of the Q4 deposits, especially in areas with notable subsidence rates. To clearly analyze this similarity, we calculated the histograms of the subsidence rate (>15 mm/yr) within and out of areas with Q4 deposits; the results are shown in Figure 11a. The percentage along the Y axis is calculated by the following equations: Subsidence is mainly located within the Q4 area, especially in areas with subsidence rates higher than 40 mm/yr (see the blue columns in Figure 11a).
To quantitatively analyze the correlation of these two factors, we also considered the thickness of Q4 deposits in the LZP, PRD, and CSP (see the Q4 depth map in Figure 2). The Q4 thickness maps of the LZP and CSP are generated by a kriging interpolation based on the depths collected from borehole measurements, which are represented by black circles [55]. We divided the thickness of Q4 into several sections, for example, 0 m to 10 m with an interval of 5 m, 10 m to 50 m with an interval of 10 m, and 50 m to 170 m with an interval of 20 m. Additionally, we calculated the average subsidence rate in these sections across the LZP, CSP and PRD. These results are represented by red diamonds in Figure 11b-d, respectively. We can see that an obvious positive correlation between Q4 thickness and subsidence rate is observed, that is, the greater the thickness is, the higher the subsidence rate, especially in the LZP and PRD. However, in the CSP, a positive correlation occurs when the thickness is larger than 50 m, and a negative correlation is generated when the thickness is low. This phenomenon may be due to the combined impact of soft soil compaction and groundwater recession.

Quantitative Correlation between Subsidence and Land-Use Class
To analyze the correlation between subsidence and land-use class quantitatively, we generated a land-use classification map based on object-oriented segmentation, as shown in Figure 7. Four frequently used classes, agriculture, forest, urban, and water, marked by yellow, green, red, and blue shading are isolated. Notably, the urban land-use class includes manmade constructions, such as buildings and bridges, and bare soil due to their similarity in spectral information. In contrast to the previous classification results [28,33], the specific class of aquaculture is identified due to its wide range and its important role in the economy of Guangdong Province, as mentioned above. The aquaculture class is indicated by purple shading in Figure 7. The overall accuracy (OA) of the classification is about 89.5% with a land-use reference map derived from [49,50].
Comparing the subsidence map ( Figure 3) with the classification map (Figure 7), we noticed that areas with notable subsidence rate (>20 mm/yr) are mainly aquaculture and agricultural areas, as indicated by purple and yellow shading in Figure 7, respectively. The representative subsidence areas in the aquaculture areas are in the LZP ( Figure 8) and PRD (Figure 9). For example, the relevant optical images of the subsidence bowls in Leizhou Bay and Zhanjiang Bay and the other subsidence bowls along the coastline of the LZP are shown in the insets A-E of Figure 8. In the PRD, the majority of the subsidence (>20 mm/yr) is concentrated in aquaculture areas, the total area of which reaches 772.2 km 2 . However, in the CSP, the subsidence is mainly located in the mixed area with agriculture and urban classes, see the yellow and red shading in Figure 7 for the two corresponding subsidence bowls, i.e., B1 and B2 in Figure 10. In addition to these representative areas, aquaculture areas along the coastline are also undergoing subsidence but exhibit different magnitudes of subsidence (see the insets A, B, and D in Figure 6). In addition to subsidence in aquaculture and agricultural areas, notable subsidence is also observed in forest and urban areas, for example, discrete subsidence bowls in forest areas of track 463 (inset C in Figure 6) and elongated subsidence in the Nansha district of Guangzhou (inset A in Figure 9a), petrochemical area of Daya Bay (inset E in Figure 6), and Macau airport (inset B in Figure 9a).
Although the correlation between subsidence and land use is obvious after qualitative comparison, to clarify this correlation, we also give the histograms of subsidence rate in four land-use classes, aquaculture, agriculture, forest and urban, as shown in Figure 12a. Notably, pixels in water areas are ignored because of low coherence. The highest subsidence (>20 mm/yr) is mainly focused in the aquaculture areas, followed by the subsidence in urban areas, agricultural areas, and finally forest (see the red, black, blue, and green lines in the insets of Figure 12b, respectively). The total percentages of the abovementioned subsidence in the four adopted classes are 40.8%, 37.1%, 21.5%, and 0.6%, respectively. Hence, aquaculture-induced groundwater exploitation is the major cause of the notable subsidence in Guangdong Province. Moreover, we also generated the percentage of different subsidence rate sections in each land-use class (see Table 4). We can see that slow subsidence mainly occurs in agriculture areas, for example, 44.3% and 34.5% when the subsidence rate is within the range of (0, 5] and (5, 10], respectively.

Causes of Subsidence in Guangdong Province
The causes of surface subsidence are generally summarized as natural compression and artificial activities [26]. The latter is the main cause of surface subsidence in many coastal areas and deltas and includes fluid and solid resource exploitation, underground infrastructure construction and compaction of constructions overlaying the Holocene sediment. In Guangdong Province, both causes are involved and are thus qualitatively analyzed in the following points.

Subsidence Probably Caused by Groundwater Exploitation for Freshwater Aquaculture Use
The areas of subsidence rate > 20 mm/yr along the coastal area occur mainly in the aquaculture areas, as mentioned in Section 6.2. The aquaculture areas in the LZP and PRD account for 53.3% of the entire coastal area (see the purple shading in Figure 7). The corresponding subsidence areas in the LZP and PRD are 112.3 km 2 and 772.2 km 2 , respectively. The maximum subsidence rates in these two areas are 58.0 mm/yr and 150.9 mm/yr, respectively. In addition to the subsidence bowls in the LZP and PRD, local subsidence bowls are also observed in the aquaculture area, see the insets B, D, and G in Figure 6. The water sources of these aquaculture areas are mainly derived from freshwater to satisfy the water quality required by the aquatic products, such as eel and tilapia [56,57]. Moreover, the frequency of freshwater exchange is high, aiming to improve the production and quality of the products. Hence, the very high requirement for freshwater and the obvious correlation between subsidence and aquaculture land use suggest that the exploitation of groundwater for aquaculture may be the primary cause of subsidence within these areas. However, we cannot analyze the explicit relationship between subsidence and groundwater exploitation volumes due to the vast illegal exploitation. Fortunately, deformation time series analysis is commonly adopted to analyze the subsidence caused by groundwater exploitation [26,28,29,58,59]. From the results in the aquaculture area (see Figures 8b,c and 9b,c), the linear trend of the deformation time series is identified. This lack of seasonal variability is likely caused by the exploitation of groundwater from the confined aquifer instead of shallow aquifer, resulting in the phenomenon of continual subsidence.

Subsidence Probably Caused by Groundwater Exploitation for Agricultural and Residential Use
Subsidence caused by groundwater exploitation for agricultural and residential use is mainly located in the CSP and the inland area of the LZP [7,60]. Unlike the subsidence that occurs in the PRD, more concentrated subsidence at very high rates (maximum value of 157.2 mm/yr) and over large

Causes of Subsidence in Guangdong Province
The causes of surface subsidence are generally summarized as natural compression and artificial activities [26]. The latter is the main cause of surface subsidence in many coastal areas and deltas and includes fluid and solid resource exploitation, underground infrastructure construction and compaction of constructions overlaying the Holocene sediment. In Guangdong Province, both causes are involved and are thus qualitatively analyzed in the following points.

Subsidence Probably Caused by Groundwater Exploitation for Freshwater Aquaculture Use
The areas of subsidence rate > 20 mm/yr along the coastal area occur mainly in the aquaculture areas, as mentioned in Section 6.2. The aquaculture areas in the LZP and PRD account for 53.3% of the entire coastal area (see the purple shading in Figure 7). The corresponding subsidence areas in the LZP and PRD are 112.3 km 2 and 772.2 km 2 , respectively. The maximum subsidence rates in these two areas are 58.0 mm/yr and 150.9 mm/yr, respectively. In addition to the subsidence bowls in the LZP and PRD, local subsidence bowls are also observed in the aquaculture area, see the insets B, D, and G in Figure 6. The water sources of these aquaculture areas are mainly derived from freshwater to satisfy the water quality required by the aquatic products, such as eel and tilapia [56,57]. Moreover, the frequency of freshwater exchange is high, aiming to improve the production and quality of the products. Hence, the very high requirement for freshwater and the obvious correlation between subsidence and aquaculture land use suggest that the exploitation of groundwater for aquaculture may be the primary cause of subsidence within these areas. However, we cannot analyze the explicit relationship between subsidence and groundwater exploitation volumes due to the vast illegal exploitation. Fortunately, deformation time series analysis is commonly adopted to analyze the subsidence caused by groundwater exploitation [26,28,29,58,59]. From the results in the aquaculture area (see Figure 8b,c and Figure 9b,c), the linear trend of the deformation time series is identified. This lack of seasonal variability is likely caused by the exploitation of groundwater from the confined aquifer instead of shallow aquifer, resulting in the phenomenon of continual subsidence.

Subsidence Probably Caused by Groundwater Exploitation for Agricultural and Residential Use
Subsidence caused by groundwater exploitation for agricultural and residential use is mainly located in the CSP and the inland area of the LZP [7,60]. Unlike the subsidence that occurs in the PRD, more concentrated subsidence at very high rates (maximum value of 157.2 mm/yr) and over large areas (342.2 km 2 ) are observed in the mixed area of agricultural and residential land use (see B1 in Figure 10) in the CSP. Within the subsidence bowl, the percentages of manmade structures and agricultural lands are 58.6% and 30.3%, respectively, in addition to the small percentage (11.1%) of aquaculture land. Thus, groundwater exploitation for the mixture of residential and agricultural use may be the major cause of subsidence. Similarly, a linear time series instead of a seasonally variable time series is observed (see Figure 10b-e), indicating continued groundwater exploitation for residential use. In addition to the groundwater exploitation, the thick Quaternary deposits in the CSP (with a maximum thickness of 167 m) may accelerate the subsidence rate during groundwater exploitation. In addition, a wide range of slow subsidence is also observed in the agricultural land of the LZP. However, the deformation time series presents an obvious nonlinear characteristic (see Figure 8e), reflecting the seasonal irrigation in the LZP.

Subsidence Caused by Land Reclamation
In addition to the exploitation of groundwater, another dominant cause of surface subsidence is the compaction of the Q4 deposits under the overlying buildings or other infrastructures. These deposits possess a high pore ratio and low strength [61,62], resulting in subsidence due to the compaction of soft soil and aquifers. Moreover, the sea reclamation phenomenon in the PRD is widespread (see the areas between the blue and red lines in Figure 9), resulting in the accumulation and compaction of soft soil and finally leading to subsidence. For example, the deposit thickness map of the PRD in Figure 2 shows that the subsidence is mainly concentrated in the areas of soft soil, and the thicker the soft soil, the higher the subsidence rate is, as shown in Figure 11d.

Subsidence Caused by Other Reasons
In Zhanjiang and Shenzhen cities, elongated subsidence bowls with maximum subsidence rates of approximately 26.6 mm/yr and 51.4 mm/yr are observed, respectively. The bowls are located in urban areas with dense industries (see inset E in Figure 6 and the optical image in inset F of Figure 8), indicating groundwater exploitation for industrial use. Because groundwater needs for industry are high and persistent, a decline in groundwater level and continued subsidence are possible (see the corresponding deformation time series of industrial areas in Figure 8d). In addition, metro tunnel excavation is a cause of subsidence in urban areas, for example, the excavation of subway line 4 in the Nansha district of Guangzhou (see inset A in Figure).
In the LZP, discrete subsidence in a large agricultural area with a subsidence rate < 20 mm/yr is caused by not only groundwater exploitation but also natural compression, such as that due to the chemical erosion of peduncle [63].
In tracks 463 and 458, local subsidence bowls are found in forest areas and may be caused by slope instability (see insets C and F in Figure 6), potentially resulting in slow-moving landslides.

Conclusions
Through the use of two MTInSAR techniques, SBAS-InSAR and TCPInSAR, we obtained a large-scale deformation map for the coastal areas of Guangdong. We also identified three notable subsidence areas in LZP, PRD, and CSP, with maximum subsidence rates approximately 58 mm/yr, 150.9 mm/yr, and 157.2 mm/yr, respectively. The relative precision of the average subsidence velocity in the overlapped areas was generally below 3.2 mm/yr, except for several values ranging from 3.4 mm/yr to 4.3 mm/yr. The root-mean-square error (RMSE) of the velocity and deformation time series calculated between InSAR and GPS results are 1.8 mm/yr and 7.6 mm, respectively.
In addition to the subsidence map, a land-use map that identifies five classes (aquaculture, agricultural, urban, forest, and water) along the coastline was generated. The relationship between these two maps showed that a high subsidence rate (>20 mm/yr) mainly occurs in aquaculture areas, followed by urban agricultural and forest areas, accounting for 40.8%, 37.1%, 21.5%, and 0.6% of the total area with this notable subsidence rate, respectively. Subsidence due to land reclamation and manmade constructions was also monitored. Additionally, we also found that subsidence mainly occurs in the Holocene deposits, and a positive correlation between subsidence and Holocene deposit thickness is obtained after the interpolation of thickness from the selected borehole data. In addition to the average subsidence velocity, the linear deformation time series in the aquaculture and agricultural areas showed that the groundwater is probably exploited continuously, and the subsidence would continue along the coastal areas if no measures were adopted. The results highlight the aquaculture-induced subsidence, in addition to the more common agricultural-induced subsidence, along these coastal areas.