Distributed Scatterer InSAR Reveals Surface Motion of the Ancient Chaoshan Residence Cluster in the Lianjiang Plain , China

The Lianjiang Plain in China and ancient villages distributed within the plain are under the potential threat of surface motion change, but no effective monitoring strategy currently exists. Distributed Scatterer InSAR (DSInSAR) provides a new high-resolution method for the precise detection of surface motion change. In contrast to the first-generation of time-series InSAR methodology, the distributed scatterer-based method focuses both on pointwise targets with high phase stability and distributed targets with moderate coherence, the latter of which is more suitable for the comprehensive environment of the Lianjiang Plain. In this paper, we present the first study of surface motion change detection in the Lianjiang Plain, China. Two data stacks, including 54 and 29 images from Sentinel-1A adjacent orbits, are used to retrieve time-series surface motion changes for the Lianjiang Plain from 2015 to 2018. The consistency of measurement has been cross-validated between adjacent orbit results with a statistically significant determination coefficient of 0.92. The temporal evolution of representative measuring points indicates three subzones with varied surface patterns: Eastern Puning (Zone A) in a slight elastic rebound phase with a moderate deformation rate (0–40 mm/year), Chaonan (Zone B) in a substantial subsidence phase with a strong deformation rate (−140–0 mm/year), and Chaoyang (Zone C) in a homogeneous and stable situation (−10–10 mm/year). The spatial distribution of these zones suggests a combined change dynamic and a strong concordance of factors impacting surface motion change. Human activities, especially groundwater exploitation, dominate the subsidence pattern, and natural conditions act as a supplementary inducement by providing a hazard-prone environment. The qualitative and quantitative analysis of spatial and temporal details in this study provides a basis for systematic surface motion monitoring, cultural heritage protection and groundwater resources management.


Introduction
Over the past decade, the Lianjiang Plain (LJP) has experienced rapid growth in its both population and social development.Environmental disturbances due to expanded human activity are surging [1].The increasing residential and industrial demand for water places stress on the water supply of this arid plain [2].As the primary industry of LJP, the textile and garment industry not only has a large demand for water, but it also causes pollution during the printing and dyeing process [3].Groundwater extraction in LJP has relieved the immediate problem of water supply, but it poses a potential threat to the geological environment in the future.The loss of aquifer storage results in an uneven and irrecoverable subsidence when the accumulated amount of groundwater exploitation exceeds the prescribed limit.Especially for LJP, the influenced ground targets include not only urban infrastructures, but also traditional villages (the Ancient Chaoshan Residences) which remain as the focus of cultural preservation in the Chaoshan Plain.Therefore, monitoring of surface motion changes over a wide extent is essential for probing groundwater level changes, ensuring structural health, and, most importantly, protecting heritage sites.
Surface motion can be monitored by reliable in situ measurements with high accuracy.Currently, LJP does not have such a surveying network with intensive observations because it is expensive and labor-consuming to perform.In bypassing these limitations, spaceborne synthetic aperture radar interferometry (InSAR) technology provides an alternative methodology that allows for precise deformation detections at high frequency and low cost.During the development of InSAR algorithms, Persistent scatterer interferometry methods were used to overcome the restrictions of the conventional InSAR method, including negative atmospheric effects and geometrical/temporal decorrelation [4,5].Persistent scatterers (PS), corresponding to pointwise targets with a high phase stability in all interferometric pairs, are generated by strong backscattering behaviors from artificial structures, bare rocks, and outcrops.Hence, PS methods have shown high effectiveness for usage in man-made and built-up urban scenarios.In comparison with PS, distributed scatterers (DS) correspond to distributed targets with a moderate coherence where adjacent pixels share close reflectivity values and similar statistical properties.Short vegetation, debris areas, and bare soil in nonurban environments often generate useful DS outputs.Typical algorithms based on DS, including SBAS [6,7] and SqueeSAR [8], have gained increasing attention due to their improved capabilities in spatial/temporal-decorrelated rural environments.
Many studies have investigated subsidence monitoring at a regional scale via InSAR in typical deltas, basins, and plains around the world, including the San Francisco Bay Area and Mississippi Delta of USA [9,10], Mexico City [11], Po Delta of Italy [12], Nile Delta of Egypt [13], Western Indonesia [14], Ganges-Brahmaputra Delta of Bangladesh [15], Mekong Delta in Vietnam [16], Southern Australia [17], as well as the North China Plain [18], Yangtze Plain [19,20], and Pearl River Delta Plain [21] in China.These well-developed plains have mature monitoring networks that include both ground truth and remote sensing interpretation.From a heritage-protection prospective, multi-temporal InSAR has also been proved to be an effective tool in risk assessment and damage monitoring in conserving the ancient rampart of Nanjing [22], the Summer Palace of Beijing [23], and the Historic Centre of Rome [24].Unlike these famous sites, LJP is a newly developing urban agglomeration with few existing surface motion monitoring networks, but reports of house destruction associated with surface motion change have increased recently.In the meantime, these destructions have been regarded as individual cases.Additionally, there is still no social consciousness that the potential threatened targets are not only scattered houses but also unstable ancient cultural protection sites as clusters.
To fill these research gaps, we present the first result of a DS-based InSAR (DSInSAR) prototype study in LJP that focuses on surface motion detection, as well as spatial and temporal pattern analysis.The two-tier network that combines the PS and DS approaches is applicable to complicated environments with a wide extent that contain both urban and rural areas.Furthermore, the clustered sites of the Ancient Chaoshan Residence are mainly located in suburban areas with moderate coherence where DS based methods perform well.Time-series InSAR derives deformation data within the radar line-of sight (LOS) and is projected in a vertical dimension to reveal surface subsidence and uplift in detail.A comparison between results from adjacent orbits is employed to cross-validate the consistency of the measurements.Based on our monitoring results, we explain the spatial distribution and temporal evolution by analyzing the local meteorological history, deposit characteristics, and geological structures.To quantitatively determine the relationship between observed displacement and Quaternary deposit features, a statistical logarithmic regression is built to describe their nonlinear correlation.A field investigation in selected ancient villages confirms the consequences of subsidence due to groundwater extraction.This work provides a detailed reference for systematic surface motion monitoring efforts in LJP.

Study Area
The Lianjiang alluvial plain is located in the southern part of the Chaoshan Plain in Guangdong Province, Southeast China (Figure 1).This plain evolved from a river valley, as well as an ancient coastal bay to a fault basin by sediment alluviation over tens of thousands of years [25].Presently, the plain extends 2500 square kilometers and has a thick alluvial deposit that reaches 120 m at its maximum.The terrain of the plain is tilted from northwest to southeast, with elevations ranging from 0 to 967 m.The plain is under the influence of a subtropical monsoon climate, with an annual average temperature of 21 • C and an annual average precipitation of 1750 mm.During summer and fall, tropical cyclones and typhoons frequently impact the plain and bring influent precipitation.
LJP is situated along the two sides of the middle and lower reaches of the Lianjiang River.As the main water source of the plain, the Lianjiang River has suffered severe pollution since the 1990s [26].Before the 1940s, small elevation differences and the curved channels of the river bank made it more prone to flooding during the wet season, whereas during the dry season, salt tide upwelling caused soil salinization on the plain.Later, in the 1950s, straightening and curve cut-off modifications were introduced to reform the river's course.In addition, a dam was constructed at the estuary to prevent salt tide upwelling.Thus, water pollution in natural conditions is now controlled, whereas there is still a lack of above-ground freshwater in LJP.
The administrative districts that LJP covers include Chaonan, Chaoyang, and Eastern Puning.These districts are the leading county-level districts in Guangdong Province in terms of population, with 5.34 million people in total living within the plain.The increasing population and intensive development of the textile and garment industry have caused urban agglomeration, an increased demand for water, and industry-introduced water pollution.Therefore, groundwater exploitation has been conducted, which has brought a loss of underground aquifer storage followed by a ground subsidence of unknown extent.
The Ancient Chaoshan Residences are mansions built in traditional Chinese style and have a history that can be traced back to the Ming Dynasty.Since families in the Chaoshan Plain with common ancestry and surnames usually live in groups surrounding their ancestry hall which is built in the center of their mansions, the Chaoshan Residences are always present in clusters.The aesthetic value of the Chaoshan Residences mainly rests with their specific architectural style.The overall style displays luxury and elegance with a comprehensive use of arts and crafts, like frescoes, ceramics inlay, and wood and stone carvings.Being the symbol that combines concepts of Chinese Fengshui, hierarchy, family, as well as ethics and morality, the Chaoshan Residences in the LJP are the focus of culture heritage protection due to their specific historical and artistic value.According to related research [27], there are more than a thousand Chaoshan Residences in the Chaoshan Plain.Among them, as a part of the Chaoshan Plain, LJP has 152 well-preserved ones.Figure 1 shows the location of LJP, and the distribution of the Chaoshan Residences in the plain.Figure 2 gives photos of typical mansions of the Chaoshan Residences in LJP.

Data
To reveal the surface motions of LJP, we collected 83 scenes of SAR images from Sentinel-1A, which cover two adjacent orbits.Tables 1 and 2 give the lists of the two ascending Sentinel-1A datasets.The scenes were acquired under the interferometric wide mode with the ascending mode.The two datasets have an overlap in coverage and time that is useful for cross-validation.Figure 1 (lower left) shows the coverage of the SAR images.Dataset I has 54 acquisitions from November 2015 to November 2017.We used this dataset to explore the temporal pattern with observations on a time series, although it only extends to middle and western portions of the plain.Dataset II has 29 acquisitions from September 2016 to February 2018 and has a full extent of the study area that is useful for exploring the spatial pattern of surface motion.
Figure 3 gives detailed information of the interferometric combination of the two datasets.The triangle mark and the square mark represent observations in dataset I and dataset II, respectively.The line between the marks denotes a combination of the master and slave images.The x-axis shows the spatial distribution of the multitemporal observations, and the value on the y-axis illustrates the perpendicular baseline.Except for the distribution of spatial and perpendicular baselines, the magnitude was also taken into consideration for choosing master images.A SRTM3 DEM with a sampling interval of 3 arc-seconds (nearly 90 m spatial resolution) was used to remove the topographic phase and for geocoding.

PS Approach
From N registered SLCs, interferograms were generated by master-slave pairs with an 8(range):2 (azimuth) multi-look operation where the topographic phase estimated from SRTM was removed.For the initial selection of PS candidates, the amplitude dispersion index (ADI) and averaged coherence coefficient were adopted.ADI is defined as [4]: where  and  are the standard deviation and mean value of the amplitude, respectively.Additionally, in the expectation of improving network quality, a threshold of 0.5 was set to exclude pixels with high side lobe signals.A triangulated irregular network (TIN) was established by connecting selected PS candidates.For a single arc section of the TIN, the phase difference between two endpoints can be expressed as the sum of four components: where ∆∅ denotes the phase of displacement in the direction of the satellite LOS, ∆∅ is the phase corresponding to elevation, ∆∅ is the phase term due to atmospheric turbulence, and ∆∅ is the residual phase associated with noise from different sources [28].The phase components of LOS displacement and height error can be modeled as:

PS Approach
From N registered SLCs, interferograms were generated by master-slave pairs with an 8(range):2(azimuth) multi-look operation where the topographic phase estimated from SRTM was removed.For the initial selection of PS candidates, the amplitude dispersion index (ADI) and averaged coherence coefficient were adopted.ADI is defined as [4]: where σ A and µ A are the standard deviation and mean value of the amplitude, respectively.Additionally, in the expectation of improving network quality, a threshold of 0.5 was set to exclude pixels with high side lobe signals.A triangulated irregular network (TIN) was established by connecting selected PS candidates.For a single arc section of the TIN, the phase difference between two endpoints can be expressed as the sum of four components: where ∆∅ de f denotes the phase of displacement in the direction of the satellite LOS, ∆∅ topo is the phase corresponding to elevation, ∆∅ APS is the phase term due to atmospheric turbulence, and ∆∅ noise is the residual phase associated with noise from different sources [28].The phase components of LOS displacement and height error can be modeled as: where ∆v and ∆h indicate increments of velocity and vertical height, and B t , B perp , λ, ρ, and θ represent the temporal and perpendicular baselines, the radar wavelength, the distance from satellite to target, and the incidence angle, respectively.For the atmospheric phase delay, adjacent points on an arc are assumed to have the same APS (atmospheric phase screen) where ∆∅ APS is regarded as zero.However, these long arcs connecting spatially various pixels may exhibit very different APS values.Thus, long arcs that exceed the preset threshold are excluded from the network in this approach.To estimate arc parameters and identify qualified PS pixels, a beamforming estimator was used for all the reserved arcs [29].Then, phase-unwrapping was performed in a temporal domain and the time-series of each PS was obtained.After parameter estimations, the network was adjusted to the selected reference point such that the obtained relative parameters were converted to absolute values.

DS Approach
In this study, DS measurements were exploited as an extension network to enhance the PS result.With sufficient SAR images (>20) in our study, the nonparametric Kolmogorov-Smirnov test [30][31][32] was applied to adaptively identify homogeneous neighborhoods.Each pixel was tested by KS within a designed 15 × 15 estimation window.For a hypothesis test of KS, the confidence value (level) in our case was set to 95%.After determining the SHP (statistical homogeneous pixels) family of each pixel, a coherence matrix of a reference pixel P was estimated by: in which d(P) is the column vector of a normalized complex value, and H represents Hermitian transportation.Ω is the SHP family.Suppose the optimal phase value of a DS is in the form of: where the first component of Φ is commonly set to zero.Based on the obtained coherence matrix, the maximum likelihood estimation (MLE) of Φ can be obtained by where ζ H = [1, e jϕ 1 , . . ., e jϕ N−1 ], and the mathematical operation • indicates the Hadamard product between two matrices.In this paper, we use the iterative phase linking [33] to solve y(Φ): where m/n is the row/column number of the coherence matrix, and k is the iteration step.γ PTA is then used as a "goodness of fit" index to assess the quality of phase estimation [8].After phase reconstruction by phase linking, a similar strategy to the PS approach was implemented for estimating the parameters of DS. Figure 4 shows the flowchart of the designed two-tier network combining the PS and DS approaches.

Vertical Deformation Retrieval
The above derived deformation from PS and DS approaches is in the direction of satellite LOS.According to SAR imaging geometry, a 1-D LOS deformation rate can be decomposed into a geometric sum of 3-D velocities: where v ew , v sn , and v vertical indicate two horizontal components and one vertical component, respectively.θ is the incidence angle and α is the angle between north and the projection of slant range onto the horizontal plane.In this study, we could not obtain 3-D velocities for ascending S1A datasets unless further information as an additional constraint was introduced.However, we were still able to retrieve the vertical deformation rate under the assumption that subsidence, especially in groundwater basins, is then dominated by the vertical component [34].A vertical deformation is given by: For dataset I swath 3, the incidence angle with respect to the central point was θ = 43.94• ; therefore, the vertical deformation rate was: For dataset II swath 1, the incidence angle with respect to the central point was θ = 34.03• ; therefore, the vertical deformation rate was: Remote Sens. 2019, 11 FOR PEER REVIEW 8 where v , v , and v indicate two horizontal components and one vertical component, respectively. is the incidence angle and α is the angle between north and the projection of slant range onto the horizontal plane.In this study, we could not obtain 3-D velocities for ascending S1A datasets unless further information as an additional constraint was introduced.However, we were still able to retrieve the vertical deformation rate under the assumption that subsidence, especially in groundwater basins, is then dominated by the vertical component [34].A vertical deformation is given by: For dataset I swath 3, the incidence angle with respect to the central point was θ = 43.94°;therefore, the vertical deformation rate was: For dataset II swath 1, the incidence angle with respect to the central point was  = 34.03°;therefore, the vertical deformation rate was:  Additionally, a homogeneous and stable pattern (−10-10 mm/year) were recognized in Chaoyang (Zone C).We preliminarily suggest that places with a deforming rate within the range (−80 mm/year to −140 mm/year) to be more risky in LJP due to a large number of textile factories distributing upon deposits with susceptible geology.Evidences and consequences of the destruction will be elaborated in Section 4.5.The statistical histogram of the MPs of Zone B and Zone C were plotted in Figure 6.Since we have removed the reasonable uplift MPs of Zone A from the statistics, MPs with positive values were regarded as deformations with respect to possible artifacts due to atmospheric disturbances, topographic errors, and other processing errors.As specified in Table 3, DS detected the notable number of MPs with densities of 310.39 and 346.73 MPs per square kilometer, with both densities increasing more than five times than those of the PS approach.These results show the advantage of using DS in regions with pointwise targets, as well as distributed targets.As shown in Figure 5, the DS approach exploited many low-reflectivity targets in mountainous areas in Zones A and B, as well as bare soil and submerged vegetation covered areas in Zone C.
The DS approach detected a strong reflectivity signal of the Puning Section of the Xiamen-Shenzhen Railway, as revealed in Figure 5b (lower left).The displacement pattern is given by a distinct trace along the track.With a displacement rate ranging from −10 to 12 mm/year, the track displays a relatively homogeneous deformation, and no significant heterogeneous trend could be identified along the route.With respect to spatial location, the Xiamen-Shenzhen Railway passes through the stable belt between two active deforming sections (A and B) and is not influenced by the surrounding up-and-down motion.Therefore, a reference point (0.3 mm/year) was selected along the stable route of the Xiamen-Shenzhen Railway, as highlighted by the red star in Figure 5.

Consistency Evaluation
In our study, the two datasets of similar imaging geometry have an overlap in spatial locations, which made it possible to compare the results despite a lack of ground truth infrastructure.As illustrated by Figure 5a,b, very similar surface motion patterns in the spatial dimension were detected.To verify the consistency of the InSAR-derived displacement rate, a resampling was employed where MPs were matched among results with the same geographic location.The fitting line gave the mathematical expression of linear regression by least squares estimation.As shown in Figure 7, based on the mathematical expression of linear regression, the correlation coefficient reached 0.96, with a determination R 2 of 0.92.This finding indicates a strong correlation between the independent results from the two datasets.The outliers with lower correlation results enclosed by

Consistency Evaluation
In our study, the two datasets of similar imaging geometry have an overlap in spatial locations, which made it possible to compare the results despite a lack of ground truth infrastructure.As illustrated by Figure 5a,b, very similar surface motion patterns in the spatial dimension were detected.To verify the consistency of the InSAR-derived displacement rate, a resampling was employed where MPs were matched among results with the same geographic location.The fitting line gave the mathematical expression of linear regression by least squares estimation.As shown in Figure 7, based on the mathematical expression of linear regression, the correlation coefficient reached 0.96, with a determination R 2 of 0.92.This finding indicates a strong correlation between the independent results from the two datasets.The outliers with lower correlation results enclosed by the red dashed ellipse indicate several divergent surface motion changes.The temporal mismatch of the datasets may lead to an inconsistency, because surface characteristics vary within different temporal windows in active zones.However, a cross-comparison was still performed on results from the two full datasets, rather than an estimate from two subsets within the same time.Within the overlapped temporal range, there were 36 images in subset I and only 18 images in subset II.Since insufficient observations may introduce uncertainty and, thus, decrease the robustness of the calculation, a correlation analysis with absolute accuracy cannot be achieved either based on data subset used in this study.

Spatial Pattern of Ground Subsidence
Various factors form the complex spatial pattern of surface motion in LJP.First, human activities in LJP, especially groundwater extraction, result in a high spatial variety of ground subsidence [35].As shown in Figure 8b, the interpolated surface provides detailed information, including displacement gradient, subsidence center, and the distribution of subsidence bowls.In the most active zone, Zone B, there were, in total, six subsidence bowls marked as B1, B2, B3, B4, B5, and B6.According to the density of the contour line, B1 showed the largest deforming gradient, in which the center of the fastest deforming MPs was detected at a rate of −140 mm/yr.It crossed 17 contour lines from the center to the periphery of the bowl within 0.51 kilometers (the length of EE' in Figure 8b), where the deforming gradient was the largest.For B2, the subsidence displacement was relatively high (−60 mm/yr), whereas it showed a relatively homogeneous deformation pattern from the center to the periphery.We do not consider it as a high risk zone.For B3, B5, and B6, the southern boundaries of the bowls were more compressed than other directions, which also made the subsidence center diverge from the geometric center of the bowl.B4 is a small subsidence bowl (0.67 square kilometers) with a fast deforming rate, crossing seven gradients.The administrative divisions that these six bowls are located in are Chendian, Simapu, and Gurao.These towns are the hubs of the textile and clothing industry in Chaonan, and include apparel markets, enterprises, as well as family workshops.
Second, the lithology of subsurface deposits in LJP provides an environment for the formation of such subsidence patterns.The Lianjiang Plain has a soft clay foundation with filling materials from

Spatial Pattern of Ground Subsidence
Various factors form the complex spatial pattern of surface motion in LJP.First, human activities in LJP, especially groundwater extraction, result in a high spatial variety of ground subsidence [35].As shown in Figure 8b, the interpolated surface provides detailed information, including displacement gradient, subsidence center, and the distribution of subsidence bowls.In the most active zone, Zone B, there were, in total, six subsidence bowls marked as B1, B2, B3, B4, B5, and B6.According to the density of the contour line, B1 showed the largest deforming gradient, in which the center of the fastest deforming MPs was detected at a rate of −140 mm/year.It crossed 17 contour lines from the center to the periphery of the bowl within 0.51 kilometers (the length of EE' in Figure 8b), where the deforming gradient was the largest.For B2, the subsidence displacement was relatively high (−60 mm/year), whereas it showed a relatively homogeneous deformation pattern from the center to the periphery.We do not consider it as a high risk zone.For B3, B5, and B6, the southern boundaries of the bowls were more compressed than other directions, which also made the subsidence center diverge from the geometric center of the bowl.B4 is a small subsidence bowl (0.67 square kilometers) with a fast deforming rate, crossing seven gradients.The administrative divisions that these six bowls are located in are Chendian, Simapu, and Gurao.These towns are the hubs of the textile and clothing industry in Chaonan, and include apparel markets, enterprises, as well as family workshops.
Second, the lithology of subsurface deposits in LJP provides an environment for the formation of such subsidence patterns.The Lianjiang Plain has a soft clay foundation with filling materials from land and marine sources, such as in most coastal parts of Guangdong [36].Geological explorations there report that the geology of LJP has two subdivisions in general: bedrock in mountainous and hilly environments, and alluvial deposits of the Quaternary series.The purple line shown in Figure 8a denotes the boundary between these two subdivisions given in previous studies [37].From Figure 8a we observe that most fast-deforming areas, including both uplift and downward trends, occur within the boundary of the Quaternary deposits.In contrast, the mountainous and hilly area shows a relatively stable and uniform pattern.Specifically, the Quaternary deposits can be divided into four layers with lithology compositions and thicknesses varying in different parts of LJP.The dashed black lines in Figure 8a show the division of subareas of Quaternary deposits in LJP.Table 4 summarizes detailed information of subsurface deposits in the plain.Within the full extent of LJP, the west-central subarea, where most of the surface displacements developed, has a maximum thickness of more than 140 m.
To further analyze the relationship between surface patterns and subsurface deposit characteristics, a statistical regression was implemented utilizing the borehole-derived Quaternary deposit thickness as well as InSAR-derived accumulated displacements from 2016 to 2018.The Quaternary thicknesses from sixteen boreholes in LJP (as shown in Figure 8a) were determined by several detailed geotechnical investigations [37,38].Figure 9 illustrates the plotted relationship between two variables.From Figure 9, we find that an outlier occurs at borehole no.W5.It is located in the West-Central subarea with an underlying Quaternary thickness of 97.56 m.The MP corresponding to this location experiences a severe accumulated subsidence with a rate of −56.98 mm/year.Its displacement is predominantly caused by human activities, and susceptible geology of soft deposits only partially affects its development.Except for borehole no.W5, other sites with in situ measurements have relatively moderate surface motion rates (−20-20 mm/year).Generally, there is a positive correlation between accumulated displacement and Quaternary thickness for each subarea.Excluding the outlier, a logarithmic expression was established for each subarea.Regression coefficients were estimated based on least squares estimation and evaluated by a hypothesis test.Additionally, the determination coefficient R square was employed to evaluate the goodness of fit.Equations ( 13)-( 14) display the logarithmic regression on the west, west-central, and north subareas.For the west and west-central subareas: y = 13.35ln(x) − 40.44 R 2 = 0.99.(13) For the north subarea: y = 25.45 ln(x) − 81.27 R 2 = 0.99 (14) where x represents the thickness of Quaternary deposits, and y represents the accumulated displacement.The p-values of the four coefficients were less than 0.05, meaning we can deny the null hypothesis that the logarithmic relationship between variables is false.For the south-east subarea, the fit curve tended to be horizontal, and the accumulated displacement did not change with variations in deposit thickness.We assume that the logarithmic regression is also tenable, whereas the estimated coefficients cannot pass the hypothesis test.However, the above mentioned non-linear correlations were only locally built due to the existence of several uncertainties.Firstly, it is well established that the fine-grained aquitard are mainly responsible for compaction and consequent subsidence [34,39].Whereas the unknown thickness of fine-grained aquitard may not be increase with the thickness of Quaternary series.Subareas with thick Quaternary series may contain thin aquitards.Secondly, the amount of groundwater pumping may not be spatially balanced, which would also result in biased estimation.Moreover, if we regard the Quaternary deposits in LJP as a comprehensive layer, its lithologic property determines how accurately the induced factor is reflected in the observed displacement.Additionally, the transmission of stress is not homogeneous, as it has a diminishing effect towards the bottom of the deposit.Third, geological structures control the spatial patterns of subsidence in LJP.There are strong seismic activities on the coast of Guangdong province, and activities of coastal fault basins (including LJP) have been significantly strengthened since the late Pleistocene [40].Since then, geological structures (faults) have been well developed with a series of fracture belts.In the plain, there are the NW Chaozhou-Puning (section JJ' in Figure 8a and Puning-Tianxin faults (section GG' and HH' in Figure 8a), with the latter being an NE hidden fracture of over 40 km in length.On one hand, fractures weaken the strength of the foundation.On the other hand, the development of deformations is usually consistent with the strike of the fracture, especially for those subsidences caused by human activities such as groundwater extraction [41].As shown in Figure 8a, the edge of the subsidence bowl B1 is tangential to the Puning-Tianxin fracture.The subsidence acceleration zone (with a deformation velocity of −60 mm/year to −140 mm/year) lies between sections GG' and HH' of the Puning-Tianxin fault, where the stability is the weakest.As shown in Figure 10, there is a significant acceleration of subsidence within the intersection between faults GG' and HH'.It is very possible that the expansion of the deformation is constrained by the fracture, and the fracture becomes the control boundary of the displacement.Additionally, fracture systems always serve as the partition of aquifer boundaries, thus, locally obstructing the flow of ground water.Therefore, the aquifer cannot be recharged in time through underground via natural ways.Similar observations have been explored by relevant research [39].
Groundwater extraction generally dominates the spatial distribution of land subsidence in LJP.Subsidence occurs when susceptible layers are pumped down, so water leaves their structure and the layers compact.We suggest that subsidence bowls occur where underground aquifer storage has suffered severe losses.In addition to human activity, natural conditions including alluvial deposits and faults provide an active and poorly compressive foundation, promoting the formation of such spatial patterns.These spatial patterns of subsidence have already resulted in incidents of house destruction, as reported in several news stories, which should continue to be monitored over time.In Zone A, an uplift was observed in Puning.The mechanism behind this phenomenon is not known completely, although several possible explanations can be given based on the InSAR-derived time series, as well as meteorological and geological bases.As shown in Figure 11a, the temporal evolution of P1 and P2 contains linear and seasonal trends.The linear trend reveals a generally ascending continuous uplift: P1 with a rate of 13.86 mm/yr from 2015 to 2017 and 10.76 mm/yr from 2016 to 2018; P2 with a rate of 6.84 mm/yr from 2015 to 2017 and 6.30 mm/yr from 2016 to 2018.The seasonal response appears as a nonlinear oscillation, indicating subsidence and recovery.There are in total four stages of rebound which indicate elastic properties of surface motion change.The first, third and fourth rebounds occur after or around September, which is slightly before the rainy season in Puning.Both the points surge in the fourth rebound.This large rebound may have been caused by the Super Typhoon Hato, which induced strong precipitation in August 2017.As denoted in earlier studies, sufficient rainfall recharges groundwater and can lead to uplift of surface motion [42].The  In Zone A, an uplift was observed in Puning.The mechanism behind this phenomenon is not known completely, although several possible explanations can be given based on the InSAR-derived time series, as well as meteorological and geological bases.As shown in Figure 11a, the temporal evolution of P1 and P2 contains linear and seasonal trends.The linear trend reveals a generally ascending continuous uplift: P1 with a rate of 13.86 mm/year from 2015 to 2017 and 10.76 mm/year from 2016 to 2018; P2 with a rate of 6.84 mm/year from 2015 to 2017 and 6.30 mm/year from 2016 to 2018.The seasonal response appears as a nonlinear oscillation, indicating subsidence and recovery.
There are in total four stages of rebound which indicate elastic properties of surface motion change.The first, third and fourth rebounds occur after or around September, which is slightly before the rainy season in Puning.Both the points surge in the fourth rebound.This large rebound may have been caused by the Super Typhoon Hato, which induced strong precipitation in August 2017.As denoted in earlier studies, sufficient rainfall recharges groundwater and can lead to uplift of surface motion [42].The first, third and fourth rebounds are most likely to be attributed to heavy rainfall during the rainy season.The second increase starts after February and ends before August, and the reason for this surface pattern is unknown.One potential explanation is that an earthquake hit Puning in February 2016.Although there is a general association between surface motion change and earthquakes, we still cannot determine the cause-and-effect relationship between the two.

Consequences of Subsidence
We conducted a field investigation in LJP to verify the consistency of DS result and to find out the consequences of subsidence on the spot.Four representative sites were determined by the geologic location of proposed subsidence bowls.As shown in Figure 8b, the blue stars indicate the locations of the four sites, and photos of these sites were taken, as shown in Figure 2.
According to our findings, as shown in Figure 12, these four Chaoshan Residences show evidence of house destruction, which verifies the consistency of the DS based method.The destruction is mainly manifested in cracks on the walls, which are regarded as the result of uneven surface motion changes due to land subsidence.Figure 12 gives typical examples of house destruction among the four sites.As is shown in Figure 12a-d, the degree of house destruction is diverging from slight to severe.Especially for Gurao, the damage of several houses is beyond the safe range of residence so that, a resettlement area has been provided from the government to the residents, as shown in Figure 12e.It is worth noting that, textile workshops and factories are also found very close to these ancient villages.With respect to Shanbing, the factory is less than 200 meters from the village, as shown in Figure 13a.Considering Gurao, which is famous for textile and garment industry, the whole town is surround by the textile industry, as shown in Figure 13b.

Consequences of Subsidence
We conducted a field investigation in LJP to verify the consistency of DS result and to find out the consequences of subsidence on the spot.Four representative sites were determined by the geologic location of proposed subsidence bowls.As shown in Figure 8b, the blue stars indicate the locations of the four sites, and photos of these sites were taken, as shown in Figure 2.
According to our findings, as shown in Figure 12, these four Chaoshan Residences show evidence of house destruction, which verifies the consistency of the DS based method.The destruction is mainly manifested in cracks on the walls, which are regarded as the result of uneven surface motion changes due to land subsidence.Figure 12 gives typical examples of house destruction among the four sites.As is shown in Figure 12a-d, the degree of house destruction is diverging from slight to severe.Especially for Gurao, the damage of several houses is beyond the safe range of residence so that, a resettlement area has been provided from the government to the residents, as shown in Figure 12e.It is worth noting that, textile workshops and factories are also found very close to these ancient villages.With respect to Shanbing, the factory is less than 200 m from the village, as shown in Figure 13a.Considering Gurao, which is famous for textile and garment industry, the whole town is surround by the textile industry, as shown in Figure 13b.

Conclusions
In this study, a time-series InSAR analysis using a two-tier network combining PS and DS approaches was used to detect surface motion changes in LJP.Our results demonstrate a wide distribution of land subsidence in the Chaonan District, a slight elastic rebound in the Puning District, and a relatively stable trend in the Chaoyang District.Vertical displacement data were retrieved by projecting LOS displacement to the vertical dimension, which we regard as the main component of the total deformation of groundwater exploitation-induced land subsidence.Very similar surface motion was detected from two data stacks acquired from adjacent orbits.This similarity was verified by a cross-comparison between two results by LSE-based linear regression.
The spatial distribution of surface motion change was clearly identified by interpolating timeseries InSAR results.In total, six subsidence bowls are located in towns where textile and clothing factories and workshops are clustered.This illustrates that groundwater extraction for industry use is a main factor impacting land subsidence in LJP.In addition to human-induced factors, natural

Conclusions
In this study, a time-series InSAR analysis using a two-tier network combining PS and DS approaches was used to detect surface motion changes in LJP.Our results demonstrate a wide distribution of land subsidence in the Chaonan District, a slight elastic rebound in the Puning District, and a relatively stable trend in the Chaoyang District.Vertical displacement data were retrieved by projecting LOS displacement to the vertical dimension, which we regard as the main component of the total deformation of groundwater exploitation-induced land subsidence.Very similar surface motion was detected from two data stacks acquired from adjacent orbits.This similarity was verified by a cross-comparison between two results by LSE-based linear regression.
The spatial distribution of surface motion change was clearly identified by interpolating timeseries InSAR results.In total, six subsidence bowls are located in towns where textile and clothing factories and workshops are clustered.This illustrates that groundwater extraction for industry use is a main factor impacting land subsidence in LJP.In addition to human-induced factors, natural

Conclusions
In this study, a time-series InSAR analysis using a two-tier network combining PS and DS approaches was used to detect surface motion changes in LJP.Our results demonstrate a wide distribution of land subsidence in the Chaonan District, a slight elastic rebound in the Puning District, and a relatively stable trend in the Chaoyang District.Vertical displacement data were retrieved by projecting LOS displacement to the vertical dimension, which we regard as the main component of the total deformation of groundwater exploitation-induced land subsidence.Very similar surface motion was detected from two data stacks acquired from adjacent orbits.This similarity was verified by a cross-comparison between two results by LSE-based linear regression.
The spatial distribution of surface motion change was clearly identified by interpolating time-series InSAR results.In total, six subsidence bowls are located in towns where textile and clothing factories and workshops are clustered.This illustrates that groundwater extraction for industry use is a main factor impacting land subsidence in LJP.In addition to human-induced factors, natural conditions also provide an environment for land subsidence in LJP.For geological features, the Quaternary deposits of LJP provide a soft foundation that makes LJP more prone to surface motion change.The distribution of Quaternary deposits and subsidence bowls is highly consistent.A logarithmic regression model was also built to describe the nonlinear correlation between accumulated displacement and the thickness of Quaternary deposits.For geological structures, active faults contribute to an uneven pattern of subsidence and become the boundary of deformation spreading.
The temporal analysis of represented MPs further verified the non-uniform surface patterns in Zone A, Zone B, and Zone C. The temporary uplift of Zone A indicates an elastic rebound by rainfall events recharging the groundwater.The significant subsidence of Zone B has developed from the beginning of the time-series, and accelerated after March 2017, demonstrating a potential threat to public safety over time.The mild and homogeneous time-series of Zone C were detected from the built-up area to the estuary, indicating a relatively stable trend.
The house destruction of selected sites of the Ancient Chaoshan Residence Cluster was investigated as a consequence of satellite SAR observed subsidences.The number and level of damaged mansions was shown to increase where the subsidence gradient surges.The printing and dyeing factories and workshops surround the subsidence bowls, which preliminarily verifies the association between groundwater extraction and land subsidence.For cultural heritage protection, it is highly recommended that the factories and workshops as the sources of water usage and water pollution are planned to be away from those Ancient Chaoshan Residences.
This paper demonstrates that DSInSAR is applicable for fully characterizing spatial and temporal patterns of surface motion change in alluvial plains.The study also provides a basis for future studies involving collaborations with engineering surveyors and hydrogeology experts to improve existing InSAR networks, thus improving both monitoring capabilities and providing references for cultural heritage protection.

Figure 1 .
Figure 1.Location of the Lianjiang Plain (LJP).The background is hill-shaded by the ASTER GDEM and the no data region in South China Sea is set to gray.The black dashed polygons represent the administrative division boundaries of Puning, Chaonan, and Chaoyang.The maximum elevation within the boundaries is 967 m.Each black dots represents an ancient village of the Chaoshan Residence in LJP.The lower left gives the administration boundary of China.The purple polygon in the lower inset is the coverage of the two datasets from Sentinel-1A, and the red blob is the coverage of LJP.

Figure 3 .
Figure 3. Interferometric combination of the two datasets and their perpendicular baselines.

Figure 3 .
Figure 3. Interferometric combination of the two datasets and their perpendicular baselines.

Figure 4 .
Figure 4. Flowchart of the designed two-tier network combining the persistent scatterers (PS) and distributed scatterers (DS) approaches.

Figure 5 . 10 Figure 5 .
Figure 5. Vertical subsidence rate maps of LJP.(a) Result from dataset I; (b) result from dataset II.The red star represents the reference point and the white dashed lines show the boundaries of the three zones of analysis.

Figure 6 .
Figure 6.Statistical histogram of measuring points (MPs) of Zone B and Zone C. The statistical analysis is based on the persistant scatterer interferometry using Dataset II.

Figure 6 .
Figure 6.Statistical histogram of measuring points (MPs) of Zone B and Zone C. The statistical analysis is based on the persistant scatterer interferometry using Dataset II.

11 Figure 7 .
Figure 7. Cross-validation plot between results from datasets acquired by adjacent orbits.The black line is the fitting line of the linear regression and its mathematical expression is given in the upper left of the figure.The outliers of the linear regression are circled by the red dashed ellipse, and they were not removed in the fit.

Figure 7 .
Figure 7. Cross-validation plot between results from datasets acquired by adjacent orbits.The black line is the fitting line of the linear regression and its mathematical expression is given in the upper left of the figure.The outliers of the linear regression are circled by the red dashed ellipse, and they were not removed in the fit.

Figure 8 .
Figure 8. Spatial distribution and deformation gradient of subsidence bowls in LJP.(a) The triangle marks indicate sixteen boreholes (marked as W32, W4, W6, CZK2, W5, W7, DZK1, W3, WYZK06, W1, WYZK02, ZK26, ZK15, WYZK04, ZK2, and W8) distributing in the test site.The blue lines are the sections of faults.The black dashed lines are the boundaries which divide the test site into four subareas: West, West-Central, South-East and North.(b) B1, B2, B3, B4, B5, and B6 are subsidence bowls identified by interpolation.The blue stars indicates four sites that have been visited during the field trip.

Figure 9 .Figure 9 .
Figure 9. Logarithmic regression between the accumulated displacement and thickness of the Quaternary deposit.

Figure 10 .
Figure 10.Profile of faults (section HH' and GG') and subsidence rate of section FF' in Figure 7a.

4. 4 .
Temporal Evolution of Surface Motion ChangesTo study the temporal evolution of surface motion in LJP, the time-series deformation of two MPs in each zone are shown in Figure11.Two temporal ranges (2015 to 2017; 2016 to 2018) were analyzed separately.As the accumulated displacement is a relative measure based on the first time point, the displacement value on 20 September 2016 based on the first time series is given as the initial value of the second time series.

Figure 10 .
Figure 10.Profile of faults (section HH' and GG') and subsidence rate of section FF' in Figure 7a.

4. 4 .
Temporal Evolution of Surface Motion ChangesTo study the temporal evolution of surface motion in LJP, the time-series deformation of two MPs in each zone are shown in Figure11.Two temporal ranges (2015 to 2017; 2016 to 2018) were analyzed separately.As the accumulated displacement is a relative measure based on the first time point, the displacement value on 20 September 2016 based on the first time series is given as the initial value of the second time series.

Figure 11 .
Figure 11.Temporal evolution of selected MPs.The location of the representative MPs are given in Figure 5. (a) P1 and P2 representing Zone A; and (b) P3 and P4 representing Zone B; P5 and P6 representing Zone C.

Figure 11 .
Figure 11.Temporal evolution of selected MPs.The location of the representative MPs are given in Figure 5. (a) P1 and P2 representing Zone A; and (b) P3 and P4 representing Zone B; P5 and P6 representing Zone C.

Figure 12 .
Figure 12. House destruction in (a) Shanbing; (b) Futan; (c) Simapu; and (d) Gurao.Cracks on the walls have been noticed in these Chaoshan Residences.The width of cracks varies from visible and narrow (as shown in (a)) via coin size (as shown in (b)) and phone size (as shown in (c)), to very wide (as shown in (d)).(e) Disaster resettlement area in Gurao.

Figure 12 .Figure 12 .
Figure 12. House destruction in (a) Shanbing; (b) Futan; (c) Simapu; and (d) Gurao.Cracks on the walls have been noticed in these Chaoshan Residences.The width of cracks varies from visible and narrow (as shown in (a)) via coin size (as shown in (b)) and phone size (as shown in (c)), to very wide (as shown in (d)).(e) Disaster resettlement area in Gurao.

Table 1 .
Basic information and list of ascending Sentinel-1A Dataset I in LJP: acquisition date (Date) and acquisition time in UTC (Time).The master image for PS and DS processing is 20170106.

Table 2 .
Basic information and list of ascending Sentinel-1A Dataset II in LJP.The master image for PS and DS processing is 20170910.

Table 1 .
Basic information and list of ascending Sentinel-1A Dataset I in LJP: acquisition date (Date) and acquisition time in UTC (Time).The master image for PS and DS processing is 20170106.

Table 2 .
Basic information and list of ascending Sentinel-1A Dataset II in LJP.The master image for PS and DS processing is 20170910.

Table 3 .
Statistics of MPs of the PS and DS approaches.

Table 3 .
Statistics of MPs of the PS and DS approaches.

Table 4 .
Quaternary series distribution in LJP.Layers 2 -4 are Quaternary units, and the thickness of Layer 1 (granite) is excluded from the calculation of thickness of Quaternary series.