Land Subsidence and Uplift Surveying by Synthetic Aperture Radar Interferometry in Lianjiang Plain

In order to investigate the distribution of land subsidence in the Lianjiang Plain, the interferometric point target analysis (IPTA) method has been carried out with RADARSAT-2 data in this work. It was found that large-scale land uplift occurred in Puning City and Shantou City, and the maximum rebound velocity exceeded 100 mm/a in 2019. Then, the Sentinel-1 data collected from June 2015 to November 2019 were used to carry out a long-term sequence inversion of the ground rebound zone and further verify the RADARSAT-2 monitoring results. The results of both kinds of data show the similar distribution characteristics of the two uplift zones. After the preliminary analysis of interferometric synthetic aperture radar (InSAR) results, geological data, and local government reports, the major factor causing land subsidence and uplift was concluded to be the change in groundwater level. The thick Quaternary strata and the compressible/expandable soil layers provide conditions for rapid subsidence and uplift. Also, the south edge of each uplift zone is limited by the WNW-trending Puning-Tianxin fault, where the land deformation is accompanied by ground fissures.


Introduction
Land subsidence is a phenomenon that the elevation of the ground in the vertical direction changes due to the exploitation of underground resources or a change in stratum structure. (1)(2)(3) It is a slowly evolving geological disaster that threatens the viability and sustainability of economic development. At present, there are more than 50 cities in China suffering from land subsidence, causing direct economic losses of tens of billions of yuan each year. The total area with cumulative subsidence exceeding 200 mm is more than 79000 km 2 . (4) In these subsidence areas, the overexploitation of underground water for industry and agriculture is the main cause of regional land subsidence. In order to reduce the hazards of land subsidence, a variety of control measures have been adopted, among which the effect of prohibiting and restricting groundwater extraction is the most significant. The groundwater level has dropped more slowly and recovered in some areas where measures to control ground subsidence have been implemented. However, in some areas of southern China with groundwater exploitation management, ground uplift has been triggered. Ground uplift is also a slowly evolving geological disaster. (1,5) When ground fissures occur with ground uplift, the impact of ground uplift is more obvious. Ground fissures further aggravate damage to the land surface, which may damage the structures of buildings, reduce the bearing capacity of foundations, and cause problems such as flooding in the basements of buildings.
InSAR has been the main technology to survey and monitor regional land subsidence in China. InSAR uses the geometric relationship between sensors and ground objects to obtain large-scale, high-precision three-dimensional information and surface morphological changes. Radar differential interferometry (D-InSAR) technology is a further improvement of InSAR technology for obtaining surface deformation. It can also be considered as an extension of InSAR technology for obtaining information on surface changes, and its accuracy can reach centimeter level or higher. (6)(7)(8) After decades of technical development and engineering applications, it has been proved that InSAR technology has a wide coverage, short measurement time, and high accuracy and can overcome the traditional measurement deficiencies resulting from the weather, topography, and human factors. (9)(10)(11) In this study, 14 RADARSAT-2 images taken from December 2018 to November 2019 were used to survey the land deformation of the Lianjiang Plain in Guangdong Province based on the time series InSAR analysis method. The results show that there are two obvious ground rebound zones with areas of 56 and 13 km 2 , and the maximum annual rebound rate exceeded 10 mm/a. The land uplift was induced by the ban and restriction of groundwater use, and the rebound speed for such a large area is the highest in eastern China. The spatial characteristics of rebound are related not only to the thickness distribution of the Quaternary strata, but also to the distribution of faults. Then, the Sentinel-1 SAR data collected over more than 4 years from June 2015 to October 2019 were used to verify the ground uplift and obtain the history of ground deformation. Finally, the mechanism inducing the regional land deformation was preliminarily analyzed.

InSAR Deformation Survey Technology
In this work, the interferometric point target analysis (IPTA) method developed by GAMMA Remote Sensing, which is based on multiple master images with small spatial and temporal baselines, was used. This method integrates the advantages of the small baseline subset (SBAS) and permanent scatterer interferometry (PSI) techniques for land deformation time series analysis. Using this method, land deformation analysis can also be effectively performed with a reduced set of data, minimizing the inaccurate external digital elevation model (DEM) effect and the spatial and temporal decorrelation. Sections 2.1-2.3 describe the process flow (Fig. 1).

Interferogram stack generation
In this work, all the data pairs within a specific baseline have been combined to calculate the interferogram set. Precise orbit data are needed for SAR data co-registration, and external DEM is used for terrain phase calculation and geographic rectification. After the co-registration of SAR data and DEM, the interferometric phase is calculated by the conjugate multiplication of each pair of radar images, and the differential phase of each interferogram is obtained by the subtraction of terrain and flat phases in the image plane. The strategy of using multiple master images works with all available image pairs, in contrast to the use of one common master image in the traditional PSI technique. The use of multiple master images increases the number of interferograms used for phase regression, even with a small amount of SAR data. Moreover, the selection of a small orbital separation can minimize the effects of the DEM error and spatial decorrelation in differential interferograms, and enables more coherent targets to be preserved in the regression analysis. After the conventional D-InSAR process for each pair, the set of differential interferograms is prepared for the next step.

Point list generation
A key element of IPTA InSAR is that the interferometric analysis is only carried out for the selected targets. The intensity variability and spectral diversity criteria have been used jointly to identify the coherent targets from a stack of SAR single-look complex (SLC) images. (7) Because the distributed targets do not show speckle behavior, whereas a single coherent scatterer dominates the echo, a significantly lower temporal variability is observed for point targets than for the distributed targets. In the intensity variability selection method, this characteristic is used to identify point target candidates in large SLC data stacks. As the measure of the temporal variability, the mean/sigma ratio is used, where the mean is the temporal average of the backscattering and sigma is the standard deviation of the backscattering. As a result, these single isolated scatterers smaller than the resolution can be recognized from a larger number of images by considering the intensity standard deviation of each pixel. Point targets with backscattering above the given threshold are selected. In the case of small SLC data stacks, the temporal variability criterion used to select point target candidates becomes unreliable for statistical reasons. For the spectral diversity criterion, only pixels whose energy remains approximately the same when processing different looks with fractional azimuth and range bandwidth are expected. A list of point target candidates for each SLC is generated by this method and is merged into a one-point candidate list.
These two point target candidate lists, derived from the intensity variability and spectral diversity criteria, are merged into the final coherent point target candidate list. Then, the differential interferometric phase for each point target candidate is generated from the twodimensional pixel data obtained in step 1.

Deformation parameter estimation
IPTA phase regression is performed to estimate the deformation parameters with the point stack of differential interferograms, which supposes that the deformation is linear, with the average deformation rate map as one of the results. The point candidates are related within a Delauney triangulation, with the selected points as the nodes related by arcs. The deformation and DEM error are estimated using a linear model. Actually, the surface subsidence is nonlinear as the groundwater level changes seasonally, which can be clearly seen when the differential interferograms are reviewed for error rejection. Once the linear terms of the deformation and DEM error have been estimated, by considering the different characteristics of the atmospheric distribution and the nonlinear deformation in time and space, their respective contributions to the phase can be separated with a combination of temporal and spatial filters. A time series of unwrapped deformation phases can be generated, given the multiple temporal reference stack of unwrapped phases primarily due to deformation. Then, the deformation history is extracted by the weighted least-squares method, including the linear and nonlinear components, taking the time interval as the weight value.

Test site
The Lianjiang Plain is located in the southern part of the Chaoshan Plain in Guangdong Province, covering an area of more than 500 km 2 (Fig. 2). The Chaoshan Plain is one of the Quaternary basins along the southeast coast of China. It not only has a large distribution area, but also has a large Quaternary sediment thickness (maximum thickness of 168 m). Song et al. divided the Quaternary strata in the Lianjiang Plain into nine basic layers based on a drill histogram. (12) The bottom stratum of the Quaternary strata is a Middle Pleistocene layer of continental sedimentary rock with a thickness of about 40 m. The underpart of the basin belongs to the Medio-Pleistocene. The top part of the Medio-Pleistocene series is stable mottled clay with a thickness of about 10 cm, which reflects the fact that there is a long sedimentary interruption between the Medio-Pleistocene and the Epipleistocene. The distribution of the Medio-Pleistocene series is limited to the area between Puning and Liangying, but the Epipleistocene and Holocene series occupy a larger part of the plain, indicating that part of the fault depression starts in the Medio-Pleistocene with a large amount of sediment incorporated since the Epipleistocene. (12) The WNW-trending Puning-Tianxin fault (F2), NE-trending Puning-Chaozhou fault (F1), and Raoping-Huilai fault (F3) are located in the Lianjiang Plain. (13) Water is a very scarce resource in the Lianjiang Plain, and the per capita surface water resource is only one-fifth of that of Guangdong Province. There are a large number of printing, dyeing, textile, paper, and other plants in the area with heavy water consumption. Since the beginning of this century, ground subsidence and ground fissures have successively appeared, seriously affecting the productivity and safety of residents. According to geological experts, they have mainly been caused by the overexploitation of groundwater, (14) particularly by local private printing and dyeing enterprises, as pointed out clearly in the local government report of Chaonan District in 2011. This report proposed the strict control of the exploitation of groundwater, the strengthening of the monitoring of ground fissures, and the adoption of measures such as relocation and evasive action.

SAR data
Two kinds of data were used in this work, namely, RADARSAT-2 fine mode data of 5 m spatial resolution and 24 d revisiting time, and Sentinel-1 IW mode data of 30 m spatial resolution and 12 d revisiting time. The acquisition period of the 14 RADARSAT-2 images was from December 2018 to November 2019. Using the IPTA method, all 91 data pairs were used in the regression analysis of the coherent point target time series. The interferometric pairs of the RADARSAT-2 data combination are shown in Fig. 3(a). The subsidence rate from December 2018 to November 2019 was obtained.
The data acquisition period of the 114 Sentinel-1 images was from June 2015 to November 2019. There were 676 interferograms generated with a temporal baseline of 36 to 150 d and a spatial baseline of less than 80 m. The interferometric pairs of the Sentinel-1 data combination are shown in Fig. 3(b). After the IPTA time series analysis, the land deformation time series from June 2015 to November 2019 and four annual deformation rate maps in 2016, 2017, 2018, and 2019 were calculated, which clearly revealed the deformation history of the two rebound areas of the Lianjiang Plain.

Regional land subsidence and uplift distribution
The land deformation distribution acquired from the IPTA processing with RADARSAT-2 data in the Lianjiang Plain in 2019 is clearly shown in Fig. 4. Both land subsidence and land rebound occurred in the site. The color of the point targets in Fig. 4 indicates the land deformation rate in 2019, where a negative value indicates land subsidence and a positive value indicates ground uplift.
Large areas of slight subsidence and small centers with large subsidence can both be found in the Lianjiang Plain. The large regions of slight land subsidence are distributed along both banks of the Lianjiang River, with a subsidence rate below 30 mm/a in Heping Town, Tongyu Town, Ximapu Town, and Chendian Town. There are three small subsidence centers, which are located in Xiancheng Town and the northwest and southeast of Gurao Town with a maximum subsidence rate of more than 50 mm/a. The smaller region with an area of more than 13 km 2 is located in Liangying Town in Chaonan District of Shantou, whose uplift rate is more than 80 mm/a in Xincuo Community (23°12′12.64″ N, 116°21′57.27″ E).
In January 2020, we visited the region and a field investigation was carried out. Many ground fissures and cracks in houses were clearly observed at the southern edges of the two uplift regions where the uplift rate changed rapidly. The cracks occurred early in the period of groundwater overexploitation and developed along the WNW-trending Puning-Tianxin fault (F2 in Fig. 2). Images showing damage to the ground and a house are shown in Fig. 5.

Ground subsidence and uplift history
The IPTA processing with a long-term Sentinel-1 data sequence obtained from July 2015 to November 2019 was carried out to obtain the ground subsidence and uplift history. The annual deformation was obtained from 2016 to 2019 as shown in Fig. 6. The results verify that largescale rapid ground uplift occurred in the Lianjiang Plain in 2019.   The maximum uplift rate of Liangying Town was 60 mm/a, which is slightly lower than the result for RADARSAT-2.
The histories of the two uplift centers obtained by IPTA series analysis are drawn in Fig. 7.

Analysis of land subsidence and uplift mechanism
The validity of the results obtained from two kinds of radar data showing that Puning City had large-scale and rapid land uplift in 2019 was verified by their very similar ranges and speeds.
The surface uplift zone of Puning City, Xiajiashan Town, Junbu Town, and Zhanlong Town is the region with the fastest uplift in China. Ground subsidence and uplift are mainly caused by local groundwater exploitation. Puning City and Liangying Town are both bases of the Chinese textile industry and there are many private printing and dyeing enterprises. Groundwater pumping during the production of textiles has caused the groundwater level to drop continuously and induced land subsidence. Land subsidence occurred on the north side of the Puning-Tianxin fault, which has strongly contrasting landforms on opposite sides. The upthrown side (southwest side) of the fault is the hilly Dananshan area with an elevation of 300-600 m, while the downthrown side (northeast side) is the flat Lianjiang Plain, (15) so the fault constitutes an obvious geomorphic boundary, resulting in significant differences in land deformation on opposite sides of the fault and the appearance of a ground fissure.
With the implementation of the local government's measures to control the extraction of groundwater, the groundwater source has gradually been replenished and the groundwater level has started to rise. (5) After the recovery of the groundwater, the volume of the soil layers, especially the sand-bearing layer in the Quaternary strata, may recover after water replenishment, then the ground may rise. (16) The zones where land changed from subsidence to uplift are distributed in the region of the Lianjiang Plain where the thickness of the Quaternary strata is the largest, with a deposition thickness greater than 140 m. The drilling data in this zone shows that the main lithology of the Quaternary strata includes gravel-bearing medium-coarse sand, silty fine sand, sand clay, clay, weathered clay, and silt, and there are multiple rhythms or sedimentary cycles in the strata. (10) The multiple compressible layers and thick deposition strata in this area provide conditions for the land surface to undergo subsidence or uplift according to the groundwater level.

Conclusion
Land subsidence is a slowly evolving geological disaster caused by the overexploitation of underground resources or changes in stratum structure. Large-scale regional land subsidence reflects the bearing capacity of local water resources. Our land deformation survey using InSAR technology showed the following: (1) In 2019, there was large-scale regional land subsidence and uplift in the Lianjiang Plain.
The subsidence was mainly distributed along the two banks of Lianjiang River, and the rate was below 30 mm/a. At the same time, there were two significant uplift zones centered on Zhanlong Town in Puning and Liangying Town in Shantou, whose uplift rates exceeded 100 and 50 mm/a, respectively.  Her research interests include LiDAR data processing and application, InSAR and its application to disaster monitoring, and the inversion of biophysical parameters. She has published more than 15 papers and won one provincial award and one ministerial award.