Investigating Ground Subsidence and the Causes over the Whole Jiangsu Province, China Using Sentinel-1 SAR Data

Interferometric synthetic aperture radar (InSAR) mapping of localized ground surface deformation has become an important tool to manage subsidence-related geohazards. However, monitoring land surface deformation using InSAR at high spatial resolution over a large region is still a formidable task. In this paper, we report a research on investigating ground subsidence and the causes over the entire 107, 200 km2 province of Jiangsu, China, using time-series InSAR. The Sentinel-1 Interferometric Wide-swath (IW) images of 6 frames are used to map ground subsidence over the whole province for the period 2016–2018. We present processing methodology in detail, with emphasis on the three-level co-registration scheme of S-1 data, retrieval of mean subsidence velocity (MSV) and subsidence time series, and mosaicking of multiple frames of results. The MSV and subsidence time series are generated for 9,276,214 selected coherent pixels (CPs) over the Jiangsu territory. Using 688 leveling measurements in evaluation, the derived MSV map of Jiangsu province shows an accuracy of 3.9 mm/year. Moreover, subsidence causes of the province are analyzed based on InSAR-derived subsidence characteristics, historical optical images, and field-work findings. Main factors accounting for the observed subsidence include: underground mining, groundwater withdrawal, soil consolidations of marine reclamation, and land-use transition from agricultural (paddy) to industrial land. This research demonstrates not only the capability of S-1 data in mapping ground deformation over wide areas in coastal and heavily vegetated region of China, but also the potential of inferring valuable knowledge from InSAR-derived results.


Introduction
Jiangsu province is located at the center of the eastern China coast area and in the Yangtze river delta region. Over 80% of its 107,200 km 2 area is characterized as plain, with a terrain height of less than 50 m above sea level. Given the advanced shipping network composed by the Yangtze River, the Beijing-Hangzhou Grand Canal and other numerous inland rivers/canals, and the high percentage of cultivatable land, the Jiangsu province has been the most prosperous region in China for many centuries. Since the 1980s, with fast industrialization and economic growth, ground subsidence has emerged and gradually developed into a kind of geological disaster which increases the risk of flooding and damages to buildings and infrastructure, and causes seawater intrusion and other types of losses [1]. It is reported that 5000 km 2 area has accumulative subsidence larger than 0.2 m, and the maximum accumulative subsidence reached 2 m until the year 2000 in the Suzhou-Wuxi-Changzhou region of southern Jiangsu [1].
In order to cope with the threat posed by the increasingly expanding ground subsidence, the provincial government of Jiangsu issued a 10-year plan in 2014 aiming at actions to mitigate ground subsidence, including enhancing the monitoring capability of ground subsidence. In fact, the government started to construct a ground subsidence monitoring network composed of leveling, borehole extensometers and Global Positioning System (GPS) stations in the southern part of Jiangsu as far back as the 1980s. These ground-based techniques are usually accurate in measuring subsidence, but cannot be densely deployed due to construction and maintenance costs. Because the information obtained from the ground-based network is incomplete, it has very limited value for decision-making about subsidence prevention and mitigation.
Unlike ground based methods, satellite Interferometric SAR (InSAR) has the advantage of obtaining ground deformation over an area, rather than at a discrete position. With the advent of time series InSAR (TS-InSAR) techniques, such as permanent scatterer InSAR (PS-InSAR) [2,3] and small baseline subset (SBAS) InSAR [4], InSAR has become an operationally usable tool for mapping ground deformations. Since then, TS-InSAR techniques have been used to monitor ground subsidence over many cities, such as the European cities involved in the Terrafirma project [5], Los Angeles and Houston in the US [6,7], and some Chinese cities like Beijing [8,9], Tianjing [10], Shanghai [11] and Suzhou [12]. Moreover, TS-InSAR has also been exploited for subsidence monitoring over large areas, such as the entire country of Italy [13], all of Norway (https://insar.ngu.no/), and the Beijing-Tianjing-Hebei region of China [14].
The launches of Sentinel-1A and Sentinel-1B (S-1) satellites in 2014 and 2016 respectively, further enhanced the superiority of InSAR for deformation monitoring due to the optimized features in imaging parameters, revisit time, spatial resolution, and orbital status [15,16]. The S-1 IW data, which is acquired through the novel Terrain Observation by Progressive Scans (TOPS) [17] technique, has a 250 km swath and can be regularly acquired every 6 days over most part of the world with the combination of Sentinel-1A and Sentinel-1B satellites. This superior spatial coverage, plus the free data policy, makes S-1 IW data an ideal data source for deformation monitoring over large areas. Internationally, the S-1 IW data has been used to map ground deformations over the entirety of Germany [18] and over the whole of Italy [19]. In China, most ground deformation monitoring efforts using S-1 IW data have covered relatively small areas, such as the Xiamen Airport [20] and Tianjin [21]. Therefore, the potential of S-1 IW data in monitoring ground deformation in wide areas in China has not been fully exploited.
On the other side, as a kind of geoharzards, ground subsidence may result in various damages to buildings, infrastructures and facilities. To know the spatial distribution and amplitude of ground subsidence is the first step, and many relevant InSAR researches focus on this aspect. But what is more important is to uncover some knowledge about ground subsidence, such as the causes, the tendency and so on. These knowledge can provide direct guidance for taking proper actions to prevent and mitigate subsidence hazard. In this paper, we report a research on mapping the ground subsidence and inferring the causes for the entire Jiangsu province for the period 2016-2018 with S-1 IW data.
The remainder of the paper is structured as follows. The study area and S-1 images are introduced in Section 2. Processing methodology is also given in Section 2 with emphasis on the co-registration scheme consisting of three course-to-fine registration steps and TS-InSAR retrieval over a large region. Ground subsidence results are presented in Section 3, together with the accuracy evaluation. In Section 4, the spatial-temporal characteristics and causes of subsidence over five major subsiding zones are analyzed in detail by synergetic utilization of InSAR-derived results, historical optical images and field-work findings. Finally, conclusions are drawn in Section 5.

S-1 Data and the Study Area
From the Sentinel-1 data archive accessed through the Copernicus Open Access Hub (https://scihub.copernicus.eu/), the descending S-1 acquisitions over Jiangsu province are very sparse and far less than the requirement for TS-InSAR analysis. Therefore, ascending Remote Sens. 2021, 13, 179 3 of 21 orbit S-1 IW images were chosen. A standard S-1 IW frame (or slice), covering an area of 250 * 180 km 2 in range and azimuth directions, consists of three sub-swaths with around 2 km overlap. Each sub-swath is composed of 8-10 slightly overlapping bursts. The S-1 data is provided in single-look complex (SLC) format with pixel spacing of 2.3 m and 14.0 m in slant range and azimuth directions, respectively. Six frames of S-1 IW data in 3 tracks (Path 142, 69 and 171) forms a complete coverage for the full territory of Jiangsu as depicted in Figure 1. Jiangsu province is composed of 13 prefectures with their names annotated in Figure 1. The numbers of S-1 IW images used, together with their acquisition time spans, are listed in Table 1. Due to strong atmospheric artifacts presented in some summer acquisitions, interferograms heavily contaminated by high-frequency tropospheric noise have been detected and discarded (please see the detail in Section 3.2). The 1-arcsec Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) is used in S-1 co-registration, geocoding and topographic phase removal.

S-1 Data and the Study Area
From the Sentinel-1 data archive accessed through the Copernicus Open Access Hub (https://scihub.copernicus.eu/), the descending S-1 acquisitions over Jiangsu province are very sparse and far less than the requirement for TS-InSAR analysis. Therefore, ascending orbit S-1 IW images were chosen. A standard S-1 IW frame (or slice), covering an area of 250 * 180 km 2 in range and azimuth directions, consists of three sub-swaths with around 2 km overlap. Each sub-swath is composed of 8-10 slightly overlapping bursts. The S-1 data is provided in single-look complex (SLC) format with pixel spacing of 2.3 m and 14.0 m in slant range and azimuth directions, respectively. Six frames of S-1 IW data in 3 tracks (Path 142, 69 and 171) forms a complete coverage for the full territory of Jiangsu as depicted in Figure 1. Jiangsu province is composed of 13 prefectures with their names annotated in Figure 1. The numbers of S-1 IW images used, together with their acquisition time spans, are listed in Table 1. Due to strong atmospheric artifacts presented in some summer acquisitions, interferograms heavily contaminated by high-frequency tropospheric noise have been detected and discarded (please see the detail in Section 3.2). The 1-arcsec Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) is used in S-1 co-registration, geocoding and topographic phase removal.

Data Processing
The S-1 data stacks were processed using a software named "Ground Deformation mapping with SAR Interferometry (GDEMSI)" developed by the Chinese Academy of Surveying and Mapping. It is based on a time series InSAR approach-Multiple-master Coherent Target Small-Baseline InSAR (MCTSB-InSAR), which has been used to process multi-temporal strip-mode SAR images from European Remote Sensing (ERS), Envisat and Radarsat-2 sensors over the Beijing-Tianjin-Hebei region [14].
The TOPS nature of S-1 IW data imposes stringent requirements on azimuth coregistration accuracy. According to De Zan et al. [22], at least one thousandth of one pixel co-registration accuracy is required to avoid phase discontinuity at burst edges. Therefore, a three-level course-to-fine co-registration module has been integrated into the GDEMSI software to specifically handle the S-1 co-registration.
The first-level co-registration is conducted based on a look up table (LUT) generated from the orbit information of master and slave S-1 slices and the SRTM DEM data. Given the 5 cm 3-D 1-sigma accuracy of S-1 precise orbit data [23] and the flat terrain in Jiangsu, the initial LUT can deliver around 0.5 pixel co-registration accuracy. In the second-level co-registration, offsets are measured over evenly distributed patches using cross-correlation between the master amplitude and the resampled slave amplitude. A linear model was adopted to fit these offsets measured in all bursts within the slice. Then the initial LUT was then refined accordingly. The third-level co-registration involved applying the enhanced spectral diversity (ESD) method [24] to determine the residual azimuth co-registration error to the required one thousandth of one pixel accuracy. The ESD method estimates the co-registration residual by exploiting the interferometric phase difference over the burst overlap region between master and slave bursts. The azimuth co-registration residual ∆y p at pixel p in the overlap region between the ith and (i + 1)th burst is [16]: where f az is the azimuth sampling frequency or pulse repetition frequency (PRF), ∆ f DC,p is the Doppler centroid difference of the two bursts at pixel p, ∆φ p is the unwrapped phase of a double differential interferogram like this: where m, s mean master and slave images, subscript * refers to complex conjugate. Since the azimuth co-registration residual between the two bursts can be regarded as a constant [24], it can be determined by taking into consideration all pixels with coherence larger than a threshold in the overlap region by where · denotes ensemble average. Please note, here we choose to unwrap the differential interferogram defined in Equation (2), because the interferogram is only generated for the overlap region, thus the computation cost for phase unwrapping over such a small region is trivial.
Up to now, this is the procedure to calculate the co-registration residual between one master and one slave images. When tackling time series S-1 images, we need to choose small-baseline interferograms to connect all S-1 images with redundancy and estimate the co-registration residual for every interferogram. Finally, the co-registration residual of every other S-1 image relative to a common master can be jointly estimated with the weighted least squares approach [25].
When the S-1 images have been co-registered, they are ready to be processed by GDEMSI like conventional stripmap images for time series InSAR analysis. The flowchart Remote Sens. 2021, 13, 179 5 of 21 of S-1 stack time series processing is shown in Figure 2. For more processing details, please refer to Reference [14]. When the S-1 images have been co-registered, they are ready to be processed by GDEMSI like conventional stripmap images for time series InSAR analysis. The flowchart of S-1 stack time series processing is shown in Figure 2. For more processing details, please refer to Reference [14]. As S-1 IW images can be acquired regularly every 12 days over most areas in China, there are redundant small-baseline interferogram combinations. But, we have found that atmosphere can bring in very strong artifacts on S-1 interferograms formed by two images acquired both in the summer season, that is, from June to September. This may be caused by dramatic tropospheric condition variations in summer in the mid-latitude coastal As S-1 IW images can be acquired regularly every 12 days over most areas in China, there are redundant small-baseline interferogram combinations. But, we have found that atmosphere can bring in very strong artifacts on S-1 interferograms formed by two images acquired both in the summer season, that is, from June to September. This may be caused by dramatic tropospheric condition variations in summer in the mid-latitude coastal Jiangsu province. As an illustration, three consecutive interferograms with 12-day separation are shown in Figure 3a-c. Strong and temporally uncorrelated fringes are present in the interferograms, which are obviously originated from varying atmosphere turbulence. For comparison, a good interferogram with very slight atmosphere artifacts is shown in Figure 3d, in which the master and slave images were both acquired in winter. Therefore, Remote Sens. 2021, 13, 179 6 of 21 those interferograms heavily contaminated by tropospheric anomalies were discarded, so did S-1 acquisitions not connected to the final interferograms. This is why the number of used S-1 acquisitions is smaller than the available in the archive. For comparison, a good interferogram with very slight atmosphere artifacts is shown in Figure 3d, in which the master and slave images were both acquired in winter. Therefore, those interferograms heavily contaminated by tropospheric anomalies were discarded, so did S-1 acquisitions not connected to the final interferograms. This is why the number of used S-1 acquisitions is smaller than the available in the archive. In contrast, a good quality interferogram formed by two winter images is shown in (d).
For deformation monitoring over a large area covered by multiple SAR frames, the mean deformation velocity along SAR line-of-sight (LOS) direction at selected coherent pixels (CP) over each SAR frame is estimated firstly, and then an adjustment process is implemented to harmonize the mean deformation velocity estimates from overlapping frames. Since the SAR viewing angle is different for different tracks, we need to project the LOS mean deformation velocity of overlapping tracks into a common direction. Projecting LOS mean deformation velocity into the vertical direction according to the following equation, which is actually equivalent to calculate the mean subsidence velocity, is a good solution.  (4) is only valid under the assumption that ground uplift or subsidence is the sole deformation. For this study, this assumption is reasonably true, because it has been reported that ground subsidence and uplift are the dominant ground deformation in Jiangsu province [26]. Therefore, the mean deformation velocity of the 6 S-1 frames is estimated respectively, and the mean deformation velocity is transformed as the mean subsidence velocity. Finally, the MSVs over the 6 frames are adjusted and mosaicked to generate a MSV map covering the full Jiangsu province. After knowing the adjusted mean deformation velocity and DEM error, the nonlinear deformation associated with every SAR acquisition of a CP can be estimated by temporal and spatial filtering of the residual phase. Finally, the deformation time series of a CP is obtained by combining the mean deformation velocity and the nonlinear deformations. Again we can project the LOS deformation time series into the vertical direction to obtain the subsidence time series. Figure 4 displays the resulting MSV mosaic over the full territory of Jiangsu province. The mosaic contains 9,276,214 CPs, or roughly 120 CPs in every square kilometer of nonwater area. The largest subsidence rate was 215 mm/year, located at a newly constructed dock road in the Lianyungang port within zone A of Figure 4. The majority (around 82%) For deformation monitoring over a large area covered by multiple SAR frames, the mean deformation velocity along SAR line-of-sight (LOS) direction at selected coherent pixels (CP) over each SAR frame is estimated firstly, and then an adjustment process is implemented to harmonize the mean deformation velocity estimates from overlapping frames. Since the SAR viewing angle is different for different tracks, we need to project the LOS mean deformation velocity of overlapping tracks into a common direction. Projecting LOS mean deformation velocity into the vertical direction according to the following equation, which is actually equivalent to calculate the mean subsidence velocity, is a good solution.

InSAR Results
v where v u−d (x, y) and v LOS (x, y) are the up-down mean subsidence velocity (MSV) and LOS mean deformation velocity respectively, θ (x,y) is the SAR viewing angle at position (x, y). Of course, the transformation in Equation (4) is only valid under the assumption that ground uplift or subsidence is the sole deformation. For this study, this assumption is reasonably true, because it has been reported that ground subsidence and uplift are the dominant ground deformation in Jiangsu province [26]. Therefore, the mean deformation velocity of the 6 S-1 frames is estimated respectively, and the mean deformation velocity is transformed as the mean subsidence velocity. Finally, the MSVs over the 6 frames are adjusted and mosaicked to generate a MSV map covering the full Jiangsu province. After knowing the adjusted mean deformation velocity and DEM error, the nonlinear deformation associated with every SAR acquisition of a CP can be estimated by temporal and spatial filtering of the residual phase. Finally, the deformation time series of a CP is obtained by combining the mean deformation velocity and the nonlinear deformations. Again we can project the LOS deformation time series into the vertical direction to obtain the subsidence time series.  Figure 4. Each of these major subsiding zones is located within a region with moderate subsidence, which confirms that severe subsidence is a localized phenomenon. There are 313,005 (3.37%) CPs showing uplift with MSV in the range of [5,32] mm/year. Furthermore, 80% of the uplift CPs have MSV in the range of [5,10] mm/year, which suggests that ground uplift is local and slight. The histogram of these uplift CPs is shown in Figure 5. Three uplift concentration zones are marked by dash-line ellipse in Figure 4. The causes for ground uplift is not clear. A possible reason for inland uplift is the rise of groundwater level, such as the ellipse region in Huai'an prefecture. For uplift over the other two ellipse regions along the coast of Yancheng prefecture, possible causes include tide fluctuation and estimate errors induced by atmospheric artifacts. A field campaign was conducted between April and May of 2019 to collect on-site information of ground subsidence at the five subsiding zones.

InSAR Results
of the CPs have MSV in the range of    Figure 4. Each of these major subsiding zones is located within a region with moderate subsidence, which confirms that severe subsidence is a localized phenomenon. There are 313,005 (3.37%) CPs showing uplift with MSV in the range of   5,32 mm/year. Furthermore, 80% of the uplift CPs have MSV in the range of   5,10 mm/year, which suggests that ground uplift is local and slight. The histogram of these uplift CPs is shown in Figure 5. Three uplift concentration zones are marked by dash-line ellipse in Figure 4. The causes for ground uplift is not clear. A possible reason for inland uplift is the rise of groundwater level, such as the ellipse region in Huai'an prefecture. For uplift over the other two ellipse regions along the coast of Yancheng prefecture, possible causes include tide fluctuation and estimate errors induced by atmospheric artifacts. A field campaign was conducted between April and May of 2019 to collect on-site information of ground subsidence at the five subsiding zones.  Besides the MSV mosaic covering the full province, the subsidence time series of every CP is also generated. These two kinds of product (MSV and subsidence time series) complement each other. The MSV mosaic shows the spatial pattern of ground subsidence over the full province, while the subsidence time series describes the temporal evolution of ground subsidence at an individual CP.

Accuracy Evaluation
In this research, leveling measurements are the reference data for accuracy evaluation. In Jiangsu province, the majority of leveling campaigns are sponsored by prefecturelevel governments, and not every prefecture has conducted leveling surveys. So, the point density, frequency and starting-ending date of conducted leveling campaigns vary by prefecture. We have collected 688 leveling measurements coming from four campaigns conducted in three prefectures: Changzhou, Nantong and Huai'an, to evaluate InSAR-derived MSV. Of the 688 leveling measurements, 527 are used for accuracy evaluation. The remainder are discarded because there were no CPs in the vicinity of the leveling points, that is, an 80m-radius circular neighborhood centered at the leveling point. For every used leveling point, a CP closest to the leveling point is located and its subsidence velocity is extracted from the subsidence time series according to the starting and ending dates of the leveling surveying. The statistics of the difference between leveling and InSAR subsidence velocity are listed in Table 2, which suggests a standard deviation of 3.9 mm/year over the 527 points. Besides the MSV mosaic covering the full province, the subsidence time series of every CP is also generated. These two kinds of product (MSV and subsidence time series) complement each other. The MSV mosaic shows the spatial pattern of ground subsidence over the full province, while the subsidence time series describes the temporal evolution of ground subsidence at an individual CP.

Accuracy Evaluation
In this research, leveling measurements are the reference data for accuracy evaluation. In Jiangsu province, the majority of leveling campaigns are sponsored by prefecture-level governments, and not every prefecture has conducted leveling surveys. So, the point density, frequency and starting-ending date of conducted leveling campaigns vary by prefecture. We have collected 688 leveling measurements coming from four campaigns conducted in three prefectures: Changzhou, Nantong and Huai'an, to evaluate InSARderived MSV. Of the 688 leveling measurements, 527 are used for accuracy evaluation. The remainder are discarded because there were no CPs in the vicinity of the leveling points, that is, an 80m-radius circular neighborhood centered at the leveling point. For every used leveling point, a CP closest to the leveling point is located and its subsidence velocity is extracted from the subsidence time series according to the starting and ending dates of the leveling surveying. The statistics of the difference between leveling and InSAR subsidence velocity are listed in Table 2, which suggests a standard deviation of 3.9 mm/year over the 527 points.

Discussion
In this section, detailed analyses of the spatio-temporal features and causes of the five subsiding zones marked in Figure 4 are presented. The analyses are based on CP's subsidence time series derived from the time series S-1 data, some historic optical images and the information collected during the field campaign. Understanding the factors driving the major subsiding processes is of great significance to guide the government's response to fight against ground subsidence.

Subsidence Zone A
The largest MSV of the whole province took place along a dock road within zone A, which is a part of the expansion project of the Lianyungang port. From the historic optical images shown in Figure 6a-c, the harbor dock construction did not start until 11 May 2012. The AB section of the dock road was completed through marine reclamation before October 2014, and the BC section was reclaimed more recently, that is, between October 2014 and January 2016 when the first S-1 image in the stack was acquired. Therefore, the BC section has very large MSV, that is, in the range of [−215, −100] mm/year, while the AB section has a smaller MSV varying in the range of [−80, −8] mm/year. The subsidence time series of two CPs over the span of 9 January 2016 to 12 December 2018 are presented in Figure 6f. Please note, the date format used in the X-axis of subsidence curve plot is "dd/mm/yy", and the date format in the images is "yyyymmdd," where d, m and y stand for day, month and year respectively. P1 in section BC is the CP having the largest subsidence velocity in the MSV mosaic. Therefore, the subsidence curve of P1 in the period from January 2016 to September 2018 shows a very steep gradient. The curve of P2 in section AB also shows a steady subsidence tendency but with a much smaller slope. Moreover, we can observe an obvious subsidence slow-down for point P1 and small-scale uplift for point P2 after 18 September 2018. The temporal characteristics of subsidence along the dock road is similar with the deformation detected over reclaimed land at Xiamen New Airport [20] and at the Lingang New City area of Shanghai [27]. According to Jiang and Lin [28], ground deformation associated with land reclamation can be well explained by the Terzaghi theory of consolidation [29], and typically will undergo two stages: a primary consolidation stage characterized by rapid ground subsidence with subsidence rate decreasing with the passage of time, and a second long-term compression stage with much slower subsidence rate. Therefore, it can be inferred that subsidence along the dock road in Zone A is induced by reclamation consolidation, and very likely the BC section is in late part of the primary consolidation stage, while the AB section is in the second compression stage.

Subsidence Zone B
Zone B is located in the northwest area of the Xuzhou prefecture. The enlarged MSV of zone B is shown in Figure 7a. Severe subsidence took place in two distinct areas, which are outlined by a purple-line polygon and a white-line polygon in Figure 7a, respectively. The purple-line polygon corresponds to the urban area of Feng county, which is the only urban region within Zone B. The Xuzhou prefecture is one of China's coal production centers, and the majority of active underground coal-mining fields are within the white-line polygon region. The northern part of Xuzhou prefecture has suffered freshwater shortage since the 1980s. As a result, the urban area of Feng county has many wells, as shown in Figure 7b, to support the daily water consumption of the inhabitants. The corresponding MSV is shown in Figure 7c and is in the range of [−89, −20] mm/year. We can observe a distinct spatial correlation between well distributions and ground subsidence occurrences. Moreover, the subsidence time series of two CPs: P1 and P2 marked by triangles in Figure 7c is shown in Figure 7d. Generally, the subsidence curves show steady tendency of ground subsidence, which can be explained by the steady groundwater consumption of daily life. Some smallscale oscillations are observed in both of the subsidence curves, which are very likely linked with variations of rainfall. As a summary, it can be inferred that groundwater withdrawal is the main factor accounting for the subsidence over the Feng county urban area. ence occurrences. Moreover, the subsidence time series of two CPs: P1 and P2 marked by triangles in Figure 7c is shown in Figure 7d. Generally, the subsidence curves show steady tendency of ground subsidence, which can be explained by the steady groundwater consumption of daily life. Some small-scale oscillations are observed in both of the subsidence curves, which are very likely linked with variations of rainfall. As a summary, it can be inferred that groundwater withdrawal is the main factor accounting for the subsidence over the Feng county urban area.    The underground coal mining region marked in Figure 7a is a rural landscape featured by farmlands and sparsely distributed villages. In order to avoid damages to village build-ings, the working panels (WPs) of underground mining are restricted beneath farmland patches. There are four coal mining WPs within this AOI with the parameters listed in Table 3. In this region, winter wheat is usually sowed before the end of October and reaped around the next June, then corn or rice is planted during the summer season from June to early October. An enlarged MSV image over the rectangle AOI marked in Figure 7a is shown in Figure 8a with an optical image as background. An on-site photo taken on 8 May 2019 (Figure 8b) shows the growing winter wheat on the ground above WP 74,101. Table 3. Mining parameters of the four working panels (WPs) shown in Figure 8a. The frequent land-cover changes on farmland patches makes it nearly impossible to extract CPs there, and CPs only exist on some man-made targets, like power poles and buildings. Therefore, subsidence information is missed over a large part of ground above WPs in the MSV map as illustrated in Figure 8a. Fortunately, we found that S-1 interferograms with 12-day and/or 24-day separation can keep good coherence when both images are acquired in the period between October and the next March when the ground is bare or the sowed wheat is in its early growing stage. 15 consecutive 12-day interferograms over the region of WP 92,606 in the period of 4 October 2016 to 14 April 2017 are shown in Figure 9. It can be observed from the moving deformation fringes that the underground tunneling had advanced from west to east along the strike direction of the coal seam. These 12-day and 24-day interferograms with good coherence can be reliably phase-unwrapped. Then the LOS deformation related to individual SAR acquisition can be calculated under the least squares principle by linking all of the unwrapped phases. had not yet started during the first coherent period. It can be inferred that the ground deformation highly coincide temporally with mining operations, and can reach hundreds to thousands of millimeters in several months. Ground deformation could be observed even after mining activities end, as suggested by a small-scale deformation field appeared on the upper-left side of WP 74,101 in Figure 8d. Therefore, ground deformations associated with underground mining are strongly nonlinear during the whole study period of January 2016 to December 2018, and may not be correctly retrieved by the linear phase model. As a further illustration, deformation time series of 6 points (P1, P1*, P2, P3, P4, P5) are given in Figure 10, whose locations are marked by triangles in Figure 8a,c,d respectively. From Figure 10a, although P1 and P1 * refer to the same location, the deformation calculated with the least square method over the first coherent period (P1* curve) is obviously larger than that from time series processing (P1 curve). This confirms that deformation caused by underground mining could be underestimated by time series InSAR analysis for the whole study period. On the other side, the curves at P2, P3, P4 and P5 in Figure 10b demonstrate that the magnitude of deformation triggered by active underground mining is usually very large. In summary, the deformation information in this underground mining region is missed spatially due to the lack of CPs, and underestimated in magnitude due to the temporal uneven of mining operations. is usually very large. In summary, the deformation information in this underground mining region is missed spatially due to the lack of CPs, and underestimated in magnitude due to the temporal uneven of mining operations.

Subsidence Zone C
Zone C is located in middle-south coast area of the Yancheng prefecture. The enlarged view of MSV is shown in Figure 11a. The MSV map within the white-rectangle AOI is furthered enlarged and shown in 11b with a high-resolution optical image as the background. The MSV value of Figure 11b    is usually very large. In summary, the deformation information in this underground mining region is missed spatially due to the lack of CPs, and underestimated in magnitude due to the temporal uneven of mining operations.

Subsidence Zone C
Zone C is located in middle-south coast area of the Yancheng prefecture. The enlarged view of MSV is shown in Figure 11a. The MSV map within the white-rectangle AOI is furthered enlarged and shown in 11b with a high-resolution optical image as the background. The MSV value of Figure 11b

Subsidence Zone C
Zone C is located in middle-south coast area of the Yancheng prefecture. The enlarged view of MSV is shown in Figure 11a. The MSV map within the white-rectangle AOI is furthered enlarged and shown in 11b with a high-resolution optical image as the background. The MSV value of Figure 11b is in the range of [−139, −21] mm/year. Four historical optical images corresponding to the AOI are given in Figure 11d-g. Dramatic land use changes have taken place. Paddy land and muddy beach areas in 2003, as shown in Figure 11d, had been transformed by 2017 into the industrial districts of the Fengyuan heating and power plant (FHPP) and the Yancheng steel group (YSG) respectively, as shown in Figure 11g. The subsidence time series of two CPs, one located in FHPP and another in YSG, are presented in Figure 11c. In the lower-left corner, an on-site photo taken on 21 May 2019 shows some large steel structures within FHPP. From historical images, it is certain that the heavy load of industrial equipment and buildings loaded on soft soil layers has contributed to the large-magnitude ground subsidence. However, there is no clear signal of subsidence rate slowing down with the passage of time in the subsidence evolution curves shown in Figure 11c, which is a key feature for ground subsidence merely driven by soil consolidation. Moreover, YSG and FHPP both should have great demand for freshwater supply, although we did not find any evidence of pumping out groundwater in this district during the field investigation. Therefore, it can be concluded that consolidation of underlying soft layer is certainly a cause for ground subsidence in zone C, and possibly groundwater withdrawal is another contributor. duction sites in 2018 in order to reduce groundwater extraction. In summary, ground subsidence along the coast region of Rudong county is caused by groundwater withdrawal and dominated by groundwater consumption associated with laver production.

Subsidence Zone D
Subsidence zone D is located in the coast region of Rudong county, Nantong prefecture. Three historical optical images acquired on 31 December 2008, 28 November 2013 and 7 March 2019 over a portion of zone D are shown in Figure 12a,c,d respectively. The enlarged view of the white-line AOI in (d) is shown in (b), where numerous rafts used for cultivating laver in sea water are visible. The off-shore region of Rudong county is the largest laver production base in China, accounting for 40% of the nation's annual laver yield [30]. Usually the harvest period of laver is from December to the next May. From Figure 4, we can observe a clear subsiding belt along the coast region, where the MSV value are in the range of [−31, −20] mm/year. A zoom-in view of the MSV over the red-line AOI in Figure 12d is shown in Figure 12e. To the right side of Figure 12e, two on-site photos taken within the AOI during the field campaign show a water tower and a family-use well respectively. Under a water tower there is usually a motor-pumped well which extracts much larger amounts of groundwater than a family-use well. There are many water towers and wells distributed in this region. We interviewed several local inhabitants during the field campaign, and were told that the water towers were mainly used to pump out groundwater for cleaning the laver during the harvest season and the family-use wells were used to provide fresh water for daily life of the inhabitants at all time of a year. In order to investigate the linkage between laver culture and ground subsidence, the subsidence time series of two CPs marked in Figure 12e having a MSV of −30 mm/year and −20 mm/year respectively are shown in Figure 12f. The curve of both CPs can be divided as a subsidence stage from January to September followed by an uplift stage from September to the next January during each year. This annually repeated temporal pattern can be well explained by the laver cultivation stages. The laver harvest season started in December, so large amount of groundwater was pumped out for cleaning laver starting from December and this process lasted until the next June to July, one or two months after the harvest ended in May. Therefore, we observe a subsiding period of January to early September assuming an approximate one-month time lag between groundwater withdrawal and ground subsidence [31]. Furthermore, July and August are the months usually having the largest volume of rainfall in Jiangsu province. At the same time, there is no need to pump groundwater for cleaning laver, because the harvest season already ended in May. As a consequence, the groundwater level would rebound then due to the double effect of rainfall increase and groundwater consumption decrease. Therefore, ground uplift is observed in the period of September to the next January. From the evolution curve of ground subsidence, it can also be inferred that the intensity of groundwater extraction (volume of groundwater extracted within a period) for laver cleaning is much larger than that for daily life. Moreover, there is a reduction of the rebound at the end of the monitored period in the curve of both CPs in Figure 12f. A possible explanation is that groundwater withdrawal associated with laver production was controlled in 2018. This has also been confirmed by local inhabitants during the field investigation. We were told that the scale of laver production decreased in 2018 than in 2017 because the local government noticed the damages to buildings caused by large ground subsidence and decided to shut down some laver production sites in 2018 in order to reduce groundwater extraction. In summary, ground subsidence along the coast region of Rudong county is caused by groundwater withdrawal and dominated by groundwater consumption associated with laver production. Figure 11. The enlarged view of Zone C in Figure 4 is shown in (a). The MSV within the white-line AOI in (a) is shown in (b). (d-g) are 4 historic images of the AOI which shows the land-use transition from farmland and muddy beach in 2003 to industrial districts in 2017. The yellow-rectangle region in (d) corresponds to the image area of (e-g). The subsidence time series of two CPs marked by white-triangles in (b) is shown in (c). An on-site photo of FHPP is shown to the left of (g).

Figure 12.
Three historical optical images of subsidence zone D are shown in (a,c,d). An enlarged view of the white-line AOI in (d) showing laver culture rafts in the sea is given in (b). An enlarged view of the red-line AOI in (d) with CPs' MSV overlaid is given in (e) and a further enlarged view shown in (f). The subsidence time series of two CPs marked by white-line triangle in (f) is shown in (g). In the right of (g), two on-site photos taken during the field campaign are shown, which highlight a water tower and a family-use well respectively.

Subsidence Zone E
Subsidence zone E is located in Wujiang county within the Suzhou prefecture and close to the south bank of Taihu Lake, where canals, lakes and ponds were densely distributed, and historically known as "a land of fish and rice." The enlarged view of MSV is shown in Figure 13a. The MSV within the white-rectangle AOI is furthered enlarged and shown in (f). The subsidence time series of two CPs marked by white-line triangle in (f) is shown in (g). In the right of (g), two on-site photos taken during the field campaign are shown, which highlight a water tower and a family-use well respectively.

Subsidence Zone E
Subsidence zone E is located in Wujiang county within the Suzhou prefecture and close to the south bank of Taihu Lake, where canals, lakes and ponds were densely distributed, and historically known as "a land of fish and rice." The enlarged view of MSV is shown in Figure 13a. The MSV within the white-rectangle AOI is furthered enlarged and shown in Figure 13b with a high-resolution optical image as the background. Four historical optical images corresponding to the AOI are given in Figure 13d-g. It can be seen that previous farmland, paddy and ponds as shown in Figure 13d,e were transformed to residential building region in 2017 as shown in Figure 13g. The subsidence time series of two CPs are given in Figure 13c, both of which present steady subsiding tendency with rate slightly slowing down as time passes. At the same time, the time of land use transition happed at P2 (from farmland to building) is at least 3 years earlier than that at P1(from pond to building), that is, consolidation of underlying soft soil layer at P2 started at least 3 years earlier than that at P1. Therefore, the MSV at P1 (−41 mm/year) is much larger than that at P2 (−15 mm/year). This temporal feature of the subsidence is similar with that in zone A, and also matches the Terzaghi theory of consolidation [29]. Furthermore, there is no excessive groundwater withdrawal or any other activities which may result in ground subsidence in this region. So, we can infer that ground subsidence in zone E is caused by the consolidation of underlying soft soil layers.
Remote Sens. 2021, 13, x FOR PEER REVIEW 18 of 21 shown in Figure 13b with a high-resolution optical image as the background. Four historical optical images corresponding to the AOI are given in Figure 13d-g. It can be seen that previous farmland, paddy and ponds as shown in Figure 13d,e were transformed to residential building region in 2017 as shown in Figure 13g. The subsidence time series of two CPs are given in Figure 13c, both of which present steady subsiding tendency with rate slightly slowing down as time passes. At the same time, the time of land use transition happed at P2 (from farmland to building) is at least 3 years earlier than that at P1(from pond to building), that is, consolidation of underlying soft soil layer at P2 started at least 3 years earlier than that at P1. Therefore, the MSV at P1 (−41 mm/year) is much larger than that at P2 (−15 mm/year). This temporal feature of the subsidence is similar with that in zone A, and also matches the Terzaghi theory of consolidation [29]. Furthermore, there is no excessive groundwater withdrawal or any other activities which may result in ground subsidence in this region. So, we can infer that ground subsidence in zone E is caused by the consolidation of underlying soft soil layers.

Comparative Analysis of Three Subsidence Scenarios
From the analyses on the five subsiding zones, there are three scenarios of ground subsidence in Jiangsu province in the monitoring period: subsidence linked with consolidation of underlying soft layer, groundwater withdrawal and underground mining.
Ground subsidence linked with underground mining is very characteristic of large deformation velocity. The evolution of subsidence highly coincides with mining progress. Especially, ground subsidence usually occurs immediately after underground extraction starts, and may reach a very large magnitude (usually hundreds to thousands millimeters) in months. Also, ground subsidence may still happen after mine extraction has stopped. Ground subsidence linked with soil consolidation can be well explained by the Terzaghi theory [29]. The subsidence curve of the primary consolidation stage is very distinctive, which is usually characterized by rapid ground subsidence with the rate slowing down with the passage of time. The time needed to complete the primary consolidation process varies from several years to several decades depending on the geological feature of the underlying stratum [28]. Zone A and E are typical occasions of this subsidence scenario. Ground subsidence linked with groundwater withdrawal is controlled by the intensity and period of groundwater extraction, and affected by rainfall. For example, there are two types of groundwater extraction in zone D: one for daily life of the inhabitants, the other for cleaning laver. Groundwater withdrawal for cleaning laver is more intensive than that for daily life, and presents strong periodic feature. Therefore, subsidence evolution curves in zone D show clear periodic pattern. As a contrast, groundwater withdrawal in Feng county of zone B is merely for daily life. Since freshwater demand for daily life usually keeps nearly invariant, the subsidence time series in Figure 7d presents a steady tendency, and the small-scale oscillations in the curve are very likely due to rainfall variations.

Conclusions
In this research, we report the methodology and results of monitoring ground subsidence over the whole Jiangsu province, China using S-1 IW data acquired between January 2016 and December 2018. This research is the first of its kind applications on wide-area deformation mapping with S-1 data in China.
Using time series S-1 images of 6 frames, the subsidence information for the whole Jiangsu province is retrieved. The derived MSV map has an accuracy of 3.9 mm/year based on leveling measurements. Around 82% of the territory is stable. Large subsidence mainly occurs in the north (Xuzhou prefecture), the south (Wujiang county of Suzhou prefecture) and along the coastal region of the province. Over the underground mining region of Xuzhou prefecture, ground subsidence data are largely missed due to unavailability of CPs on farmland, and InSAR-derived MSV may be underestimated.
More importantly, based on InSAR-derived subsidence time series, the causes of ground subsidence over major subsiding zones of Jiangsu province have been dig out, which include underground mining, groundwater withdrawal and consolidations of underlying soft layer associated with marine reclamation, lake/pond reclamation and land-use transition. Disclosure of the driving factors of ground subsidence has very valuable implications for the government to take proper actions. For example, it could be suggested that the efforts of mitigating ground subsidence should be focused on the management of underground coal mining and groundwater withdrawal, especially the groundwater withdrawal associated with industrial purpose, rather than daily life. It could be expected the ground subsidence related to soil consolidation would slow down and finally stabilize in some time. From this point-of-view, this research represents a significant demonstration on how the InSAR-derived deformation could be utilized for information inference and decision-making support.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.