2.1. Time-Series InSAR Processing
On PS-InSAR processing [
41], the master image is chosen to optimize spatial-temporal de-correlation. To all the chosen Persistent Scatterer Candidates (PSCs),
, the interferometric phase of a flattened differential interferogram,
N, is given by a combination of the various contributions expressed below [
35]:
where
s and
k indicate the interferometric pairs, where
s belongs to the slave image,
k can be defined as a master image [
35]. Because of external DEM (Digital Elevation Model) inaccuracies,
,
corresponds to DEM error and is expressed as a linear rate of deformation.
denotes the atmospheric phase delay, and
expresses the temporal and geometrical decorrelation noise.
The Amplitude Dispersion Index (ADI) identifies pixels of relatively low SAR amplitude variation values to be candidates for PS-InSAR analysis, i.e., as PSCs. Hence, a higher average coherence or a lower ADI is associated with superior phase quality. The particular pixel ADI is denoted as
where
and
represent the standard deviation and the mean value of the amplitude, respectively.
The SBAS-InSAR process [
42,
43], involves four steps [
44]:
Step 1: Select a small baseline subset.
Step 2: Set the interference processing of the defined baseline sets.
Step 3: Retrieve the unwrapped deformed phases with the three-dimensional phase unfolding algorithms.
Step 4: Separate the components of error in the unfolding phase to derive the true deformation. The SBAS-InSAR formula is presented below [
45]:
Equation (
3) represents the number ranges of differential interferograms produced from
SAR images of the same region in a given time period
. Equation (
4) outlines the formation of the interferometric phase of the interferogram
j (generated from two images at
and
) in the pixel
, in which
x and
r indicate the azimuth and distance coordinates, separately. For each PS-InSAR or SBAS-InSAR, the deformation rates along the line of sight (LOS) and the deformation time-series during acquisition can be efficiently measured. When accounting for the direction of movement, the most suitable LOS for the radar sensor is along the ascending track, via which it can measure the actual displacement more reliably.
The deformations registered by PS-InSAR are distinguished by negative values. In addition, the notations
,
,
,
and
are produced by the terrain, satellite orbit error, topography, atmosphere effect, and other noises, respectively. Furthermore, the
M equations with
N unknowns can be represented in the shape of the matrix as follows.
where A denotes the
coefficient matrix.
is the
vector of unknown deformed phase value,
is the
vector of the unfolded phase value.
2.2. Proposed Method
The proposed method is divided into two parts. First, the Iterative Self-Organizing Data Analysis Techniques (ISODATA) algorithm is used to define the boundary between the slow subsidence region and the fast-subsidence region by clustering the distribution of land subsidence rates obtained by PS-InSAR. Second, when new SAR data are needed, SBAS-InSAR is used to efficiently develop small SAR data sets that include partial time-series SAR data and new SAR data. Moreover, the monitoring results are the sum of the latest results of SBAS-InSAR and the original result of the PS-InSAR for the time-series SAR datasets that are not included in the small SAR data sets. After updating with a certain amount of SAR data, the boundary and the monitoring results will be recalculated by PS-InSAR and SBAS-InSAR. The processed algorithm is partitioned into 10 steps, the procedure of which is shown in
Figure 1.
Step 1: N SAR images are acquired from the SAR images set, for which N is more than or equivalent to 20, the SAR image set is defined as depicting the identical observation scenario, and the sets are aligned.
Step 2: The SAR image set derived in Step 1 is used in PS-InSAR for land settlement monitoring and to acquire the PSCs. The SBAS-InSAR also processes the SAR image set acquired in Step 1 to perform ground deformation monitoring and acquire the Slowly Decoherent Filtering Phase (SDFP) and the corresponding pixel coordinates.
Step 3: The acquired SAR image data are filtered and put into the image data set if they meet the filtering conditions. The primary screening condition involves calculating the coherence coefficient of the acquired image. If N1 of the calculated N coherence coefficient exceeds the predetermined threshold value p, then the alignment results satisfy the preset coherence conditions, where N1 > N × 0.6; P > 0.6.
Step 4: If the number M of images in the current image datasets meets the conditions (M − N)/N > 50%, go to Step 5, Otherwise, go to Step 6.
Step 5: PS-InSAR processing is performed on the image data from Step 4, the PS surface deformation results are calculated, the PS deformation results calculated in Step 2 are updated, and the maximum comprehensive coherence coefficient derived via PS-InSAR processing 1 and the number of PS-InSAR images are updating, return to Step 3.
Step 6: N2 SAR images are selected from the image data set in Step 4 to form M1 small baseline sets and SBAS-InSAR processing is performed on the small baseline sets of the SAR image data newly acquired in Step 3 to obtain new SBAS surface deformation results, The average comprehensive coherence coefficient 2 is updated.
Step 7: The ISODATA algorithm is used for the cluster analysis of the solution results derived in Step 6 to divide the deformation area and extract the most strongly developed region.
Step 8: The new SAR images of the strongly developed land subsidence area are acquired and screened again for SBAS-InSAR processing.
Step 9: The distribution of the accelerated land subsidence area is determined through the SBAS-InSAR monitoring results derived in Step 8.
Step 10: The SBAS deformation results from Step 6 and Step 8 are fused by the method of the weighted average for the accelerated land subsidence area. The PS deformation results from Step 5 and the SBAS deformation results from Step 8 are fused by the weighted average method for the slow land subsidence area. The weighted average method’s formula is as follows:
where,
R represents the fused land subsidence monitoring results,
R1, and
R2 indicate the ground settlement monitoring results of PS and SBAS-InSAR, separately, and
w1 and
w2 are weighting factors.
2.3. Land Settlement Rate Classification and Division
The ground settlement rate classification is an important basis for designing ground settlement monitoring networks and selecting ground settlement prevention methods. The principle of classifying strong and weak areas in this paper is that if the annual average settlement rate is less than 10 mm·a
, the area is classified as a weakly developed area of ground settlement, if the annual average land subsidence rate is between 10 and 30 mm·a
, the area is defined as a moderately developed area of ground settlement, and if the annual average land subsidence rate is ≥30 mm·a
, then the area is categorized as a strongly developed area of ground settlement, as listed in
Table 1.
The ISODATA algorithm is a spatial clustering algorithm based on image segmentation and it is an unsupervised classification method. The ISODATA algorithm can use spatial information to identify and segment regions in images with highly similar target values and large differences between targets and achieve classification in a stepwise evolutionary manner. This study is used to classify land subsidence rates using three data layers simultaneously, such as SDFP, the corresponding pixel coordinates, and SABS InSAR-based ground subsidence values [
46]. The ISODATA algorithm employs the following procedure [
47]:
- (1)
Select initial parameters.
- (2)
Compute the distance index function for each of the clusters.
- (3)
Merge or split clusters in accordance with the given requirements.
- (4)
Loop iteration. Calculate the new index and define whether the result matches the clustering demands [
47].
The ISODATA algorithm proceeds as follows, as shown in
Figure 2:
Based on the unsupervised classification tool in the Arctoolbox of ArcGIS, the initial value
of the clustering center can be chosen randomly and is set to 1 here.
,
and
I are considered the default values for the ISODATA clustering algorithm in ArcGIS,
is given as 20, and the size of the neighborhood window R is a 20 × 20 grid. The input parameters are described in
Table 2.
2.4. New SAR Data Filtering, Updating, and Integrating
Due to the high temporal and spatial sampling rate, SBAS-InSAR requires few satellite images, but generally, more than eight images, to ensure monitoring accuracy [
6,
48]. The fast-rate region can be searched for by applying SBAS-InSAR to small SAR data sets, including partial time-series SAR data and new SAR data from strongly developed ground deformation areas. To avoid the event of low coherence between the new SAR data and the original SAR data, the coherence coefficient is applied to screen the new SAR data.
The maximum temporal resolution of the Sentinel-1A data is 12 days, with a critical time baseline of 4a, a critical Doppler center frequency difference of 486.486 Hz, and a critical spatial baseline of 6448 m. The screening criteria for the new SAR data are listed in Equations (
7) and (
8), and the master image is chosen according to the stacked coherence of the interferogram using Equation (
7).
where
n,
l is the serial number of the master and slave images.
represents a temporal baseline for the master and slave images.
indicates the spatial baseline of the master and slave images.
expresses the Doppler center frequency difference between the master and slave images.
indicates a critical effective baseline.
is a critical temporal baseline.
denotes the critical Doppler center frequency difference.
is the set of images meeting the filtering criteria (
).
is the master-slave image coherence coefficient (
) [
32,
49,
50].
When the SDFPs from the SBAS-InSAR data set with higher spatial coverage are applied to the subsidence regions with rapid rates, two sets of SDFPs data are superimposed according to the sampling times of the two adjacent phases of SBAS-InSAR results, due to the different time scales of the two phases. The low-quality and low-density PSCs from the rapid-subsidence regions are removed from ArcGIS, and the updated SDFPs in the rapid-subsidence regions and the PSCs in the slow-subsidence regions are overlaid to obtain the fusion results for the areas of intensely different subsidence.