An Improved Method for Automatic Identification and Assessment of Potential Geohazards Based on MT-InSAR Measurements

Abstract

Multi-temporal Interferometric Synthetic Aperture Radar (MT-InSAR) has been widely used for ground motion identification and monitoring over large-scale areas, due to its large spatial coverage and high accuracy. However, automatically locating and assessing the state of the ground motion from the massive Interferometric Synthetic Aperture Radar (InSAR) measurements is not easy. Utilizing the spatial-temporal characteristics of surface deformation on the basis of the Small Baseline Subsets InSAR (SBAS-InSAR) measurements, this study develops an improved method to locate potential unstable or dangerous regions, using the spatial velocity gradation and the temporal evolution trend of surface displacements in large-scale areas. This method is applied to identify the potential geohazard areas in a mountainous region in northwest China (Lajia Town in Qinghai province) using 73 and 71 Sentinel-1 images from the ascending and descending orbits, respectively, and an urban area (Dongguan City in Guangdong province) in south China using 32 Sentinel-1 images from the ascending orbit. In the mountainous area, 23 regions with potential landslide hazards have been identified, most of which have high to very high instability levels. In addition, the instability is the highest at the center and decreases gradually outward. In the urban area, 221 potential hazards have been identified. The moderate to high instability level areas account for the largest proportion, and they are concentrated in the farmland irrigation areas, and construction areas. The experiment results show that the improved method can quickly identify and evaluate geohazards on a large scale. It can be used for disaster prevention and mitigation.

1. Introduction

InSAR has become an important technique for large-scale ground deformation monitoring, as it has wide spatial coverage and high monitoring accuracy [1]. The fast development of the Synthetic Aperture Radar (SAR) technology provides vast amounts of SAR datasets with high spatial and temporal resolution to measure the ground deformation around the world [2]. Nowadays, the accuracy and efficiency of InSAR data processing have been greatly improved by the MT-InSAR technologies [3], such as Persistent Scatterer InSAR (PS-InSAR) [4,5], SBAS-InSAR [6,7] and SqueeSAR [8], which have overcome some limitations (temporal decorrelation, atmospheric delay, and ramp phases) inherent in Differential Interferometry SAR (DInSAR) [9,10], and have been extensively exploited to geological disaster monitoring, including urban subsidence [11,12], landslides [13,14], volcanoes [15] and earthquakes [16,17], and so on. Traditionally, the basic information (location, area, and deformation magnitude) of potential geohazards is obtained by visually interpreting the MT-InSAR measurement, which is inefficient, labor-intensive, and error-prone for large-scale area monitoring [18]. Therefore, it is necessary to develop a more efficient and comprehensive way to identify potential geohazards from the massive deformation information.

Many efforts have been made to investigate and manage geohazards, including the case study of a single geohazard, landslide identification in mountainous areas [19,20,21] and subsidence investigation in urban areas [22,23,24]. However, only a few studies investigated the methods for automatic identification and evaluation of geohazards over large-scale areas. Barra et al. (2017) proposed a method based on PS-InSAR technology to regularly update geohazard activity maps, identify active deformation areas (ADA), and update the existing geohazard database [25]. It consists of a step to detect the deformed regions and a step to evaluate the ADA quality. This method has been applied in geohazards mapping [26,27,28], geological motion assessment, and the intensity evaluation of landslides [26,29,30] based on the deformation velocity. On the basis of this method, some automatic identification and classification programs have been proposed [18,31]. In addition, some Geographic Information System (GIS) environment-based studies also combined the hazards and vulnerability of the study area to calculate the specific landslide risk [32,33].

Currently, most algorithms for geohazard detection only use the spatial information of deformation, such as the average deformation velocity, but do not consider the temporal evolution characteristics, which, however, is crucial for hazard warning. Extracting potential hazard regions from the MT-InSAR results and analyzing those temporal evolution characteristics are important for improving the interpretation efficiency and making full use of InSAR measurements. Raspini et al. [34,35] and Del Soldato et al. [36] identified anomalous changes of the time series displacements. Their methods can obtain different evolution trends in the time series displacements, such as acceleration, deceleration, and seasonal changes.

Following the previous studies of spatial analysis, this study proposes an improved method for automatically identifying the deformation regions over large-scale areas and extracting the evolution characteristics of the deformed regions by estimating the temporal development trend. This method can detect the potentially unstable areas and refine the instability boundaries of each deformation region using the spatio-temporal deformation characteristics from MT-InSAR results. We select a mountainous area (Lajia Town, Qinghai Province) and an urban area (Dongguan City, Guangdong Province) in China to validate the improved method. On the basis of the experimental results, we combine the ascending and descending data to generate the vertical deformation and the deformation in the aspect direction, and get the final results by combining the deformations in the two directions. We analyze the results of velocity gradation, temporal evolution trend estimation, instability assessment, and refinement. Finally, we discuss the advantages and limitations of the improved method.

2. The Improved Method

This section describes the improved method for automatically identifying geological disasters using the MT-InSAR measurements, which is adaptable in both mountainous and urban areas. The method has the following five modules (Figure 1):

(1)

Deformation area recognition: recognize and locate the boundary of each deformation region;

(2)

Deformation velocity gradation: classify the deformation velocity magnitude for the points in the deformed regions;

(3)

Deformation temporal evolution trend estimation: estimate the temporal evolution trend for the points in the deformed regions;

(4)

Deformation instability degree generation: combines the deformation velocity gradation in (2) and trend in (3) to obtain the instability degree for the points in the regions by on the quantitative instability matrix;

(5)

Deformation instability refinement: refines the deformation boundary of instability degree in (4).

2.1. Deformation Area Recognition

This module uses the deformation velocity obtained from the MT-InSAR to automatically identify the deformation regions (Figure 2a). It includes the following four steps:

(1)

Data preprocessing: A noise removal method is used to remove the points that are sparsely distributed or significantly different from the neighboring points in displacement [37] to obtain spatially consistent deformation results.

(2)

Cluster analysis on deformation point attributes: Set a deformation velocity threshold and the expansion radius parameters. If the point has a deformation velocity larger than the threshold, it is defined as an active point, otherwise it is a stable point. The active points are buffered according to the expansion radius, and attribute clustering is performed following the principle of spatial proximity relationship [29]. Up to this step, the deformation regions formed by clustering of active points can be obtained (Figure 2b). In order to refine the boundary shape, a smoothing process is performed to the expansion area.

• Velocity threshold: Standard deviation, $\mathsf{\sigma }$, calculated by the deformation velocity values from MT-InSAR. It reflects the sensitivity and noise level of the results. It is used to reclassify the points by activity [25,28]. In this study, the value of 3σ is adopted to determine the active point and stable point.
• Buffer distance: Deformation points within the buffer distance are considered belonging to the same deformation area. The buffer distance can be adjusted according to the spatial resolution of the MT-InSAR result and geological conditions of the study area. In this study, we set the buffer distance as 30 m.

(3)

Preliminary identification of the deformation regions: Identify multiple deformed regions according to the cluster analysis result in step (2). Some deformation regions may be noise or error. Therefore, we eliminate the deformation regions smaller than the given minimum size to improve the distribution of the results.

• Minimum size: set the minimum area of deformation regions (km²) set according to the local geological background and remove regions smaller than the threshold, as they may be caused by data processing errors, atmospheric effects, and so on. In this study, we set the minimum area as 0.05 km² and 0.1 km² for mountainous and urban areas, respectively.

(4)

Deformation result refinement: Using the preliminary results obtained in step (3) and the Delaunay triangulation algorithm [38] to connect all the boundary points of a single area to generate the boundary of the smallest convex polygon. Then, generate the boundary of each deformation region, considering the topology relationship of boundaries.

2.2. Deformation Velocity Gradation

In order to identify the deformed regions and evaluate their deformation states, we should grade the deformation velocity gradation thresholds and classify points by activity. Some studies have shown that the deformation velocity threshold should be determined according to the geological conditions of the study area [39,40,41].

In this study, the points are grades into three levels by the deformation velocity: stable point (low), active point (moderate), and highly-active point (high) (

Table 1. The thresholds for deformation velocity gradation.

Gradation	Threshold (mm/yr)
Stable point (low)	$\left|V\right|$ < 3σ
Active point (moderate)	3σ <$\left| V\right|$ < 6σ
Highly-active point (high)	$\left| V\right|$ > 6σ

Note: $\left|V\right|$ represents the absolute value of the deformation velocity obtained from the MT-InSAR results. The thresholds can be adjusted according to the deformation characteristics (density and uncertainty of point targets).

). In Section 2.1, we identify the deformation regions based on point buffered distance, so some stable points are included in the deformation regions (Figure 2c). In order to distinguish the stable and active points, we choose 3σ and 6σ as the velocity gradation threshold [18,25]. When the absolute value of velocity $\left|V\right|$ of a point is less than 3σ, the deformation velocity gradation is low, and the point is stable; when $\left|V\right|$ is greater than 3σ and less than 6σ, the velocity gradation is moderate, and this point is active; when $\left|V\right|$ is greater than 6σ, the velocity gradation is high, and the point is highly-active. The result of velocity gradation is shown in Figure 2d.

2.3. Deformation Temporal Evolution Trend Estimation

The improved method not only uses the deformation velocity to evaluate the deformation intensity of the area, but also uses the time series displacement to estimate the temporal evolution trend. According to the fitting curve of the time series displacement, the development trend can be divided into four categories: the exponential acceleration model (EAM), which shows an acceleration trend and has a high unstable level; the unary linear model (ULM), which shows a deformation with uniform rate and has a moderate unstable level; the exponential mitigation model (EMM), showing a deceleration trend, and low unstable level, gradually tending to a smooth deformation; the category of no specific deformation characteristics (NSC), which has a low unstable level, cannot be fitted by specific function expressions and without regular displacement pattern. The first three models are shown in Figure 3.

The specific procedure for estimating the deformation development trend is shown in Figure 4a. First, the first three trend models are all fitted to the target point time series, and the Root-Mean-Square Errors (RMSEs) and correlation coefficient (Rs) are calculated, which are denoted as RMSE_EAM, RMSE_ULM, RMSE_EMM, R_EAM, R_ULM, and R_EMM. The absolute value of Rs is between 0 and 1. Generally speaking, the closer $\left|\mathrm{Rs}\right|$ to 1, the stronger the correlation between time series deformation and fitting curve. In this study, we tend to choose 0.8 as the correlated threshold and 0.01 as the model similarity threshold, which are the most appropriate value for the experimental area. Therefore, if all the Rs are less than 0.8, the point trend is considered as NSC; if the RMSE_EAM is the smallest, the point trend is considered as EAM; if the RMSE_ULM is the smallest, and the difference between R_ULM and the R_EAM exceeds 0.01, the point trend is considered as ULM; if the RMSE_EMM is the smallest, and the difference between R_EMM and R_ULM exceeds 0.01, the point trend is considered as EMM. Figure 4b–d are the examples of ULM, EAM, and EMM. The threshold mentioned above can be adjusted according to the local situation.

When the trend of the fitting result of the target point is similar, that is RMSE_ULM or RMSE_EMM is the smallest and R_ULM or R_EAM is infinitely close to the R_EAM or R_ULM, respectively, the point is considered having the higher-level change trend, so as to achieve a warning effect. For example, if the value of RMSE_ULM is the smallest, and the difference between R_ULM and R_EAM exceeds 0.01, the point trend is considered as ULM, otherwise, it is EAM. After evaluating the trend, we summarize and analyze the points with the same level to obtain the distribution of different trend levels in each region. The result of trend estimation is shown in Figure 2e.

2.4. Deformation Instability Degree Generation

Using the deformation velocity gradation obtained in Section 2.2 and the development trend obtained in Section 2.3, we establish a 3 × 4 qualitative instability matrix, which uses four levels to assess the instability of the points in the deformation regions (Figure 5). The results of instability degree will be used in Section 2.5.

2.5. Deformation Instability Refinement

Each point in the deformation regions has three attributes, which are the deformation velocity graduation, trend, and instability. It is heterogeneous in the space domain. To provide useful information for disaster prevention and mitigation, we should refine the deformed regions. Large-scale monitoring areas usually have a large number of deformation points, so we propose to refine the deformation area based on the results of spatiotemporal instability. This process is similar to the recognition algorithm mentioned in Section 2.1 theoretically, but its recognition area is much smaller. In order to locate the unstable regions in each small area and improve the recognition accuracy, the instability boundary refinement method is targeted at the cluster analysis of the points in each deformation area.

The main steps include: clustering and filtering of the deformation points with the same instability degree, and refining the deformation region boundaries The refinement results are the instability classification boundary of each deformed area (Figure 2f), which is a part of the final products. This method can greatly improve the resolution of locating highly unstable areas and evaluating potential hazards.

3. Experiment and Data Processing

3.1. Study Area and Datasets

The specific location and Sentinel-1 image coverage of our study areas are shown in Figure 6. The mountainous area, north of Lajia Town, is located at the junction between Golog County and Hainan County in Qinghai Province, in the upper reaches of the Yellow River basin, with an elevation of 3000–4000 m. This area is mountainous and has plateaus, featured by rugged and steep terrain. The slopes here are generally above 30 degrees, and even greater than 50 degrees in some areas. The valleys are mostly narrow and V-shaped, with only small patches of bottomland or no bottomland. This area is especially developed with loess and other disaster-prone rocks [42]. The study area is close to the government of Lajia Town, which has dense residential areas, so serious geological disasters will cause huge damages.

The urban experimental area is Dongguan City, located at the center of the east coast of the Pearl River Delta, with a total area of 2465 km² and a permanent population of 8,464,500. The landform is dominated by hills and alluvial plains. The southeast is mountainous with large undulations, and the elevation is mostly 200–600 m. The northeast close to the riverside of East River has gentle terrain, which is easy to cause seeper. The northwest and southwest are alluvial plains with low and flat terrains, so the strata in this area have poor stability, high water content, and high compressibility. Dongguan has mild temperatures and rich rainfall throughout the year. The average annual temperature in 2019 is about 23.9 °C and the rainfall is 1900 mm (January 2019–December 2019 from http://www.dg.gov.cn/). Moreover, it suffers frequent extreme weather (such as typhoons and heavy rainfall), which, coupled with a large number of engineering projects, cause a lot of land subsidence occurring with wide distribution in Dongguan City [43,44].

Using the SBAS-InSAR technology (Figure 7), we obtain 71 descending (18 January 2017–13 July 2019) and 73 ascending Sentinel-1 images (24 January 2017–18 August 2019) covering Lajia Town, and 32 Sentinel-1 images covering Dongguan (8 December 2018–27 December 2019) to extract the time series deformation. The acquisition time and basic parameters of the used Sentinel-1 data are shown in

Table 2. Parameters and acquisition dates of the images used in this study.

Study Areas	Direction	Orbit	Heading	Incidence	Pixel Spacing (Rg × Az)	Number of Images
Lajia Town, Qinghai Province (mountainous area)	Descending	T33	−169.73°	33.69°	2.33 × 13.96 m	71
	Ascending	T26	−10.05°	33.81°	2.33 × 13.98 m	73
Dongguan, Guangdong Province (urban area)	Ascending	T11	−10.55°	34.04°	2.33 × 14 m	32

and

Table A1. Parameters and acquisition dates of the images used in this study.

Study Areas	Parameters		Acquisition Date
Lajia Town, Qinghai Province (Mountainous area)	Direction	Descending	2017/01/18, 2017/02/11, 2017/03/25, 2017/04/06, 2017/04/18, 2017/04/30, 2017/05/12, 2017/06/05, 2017/06/17, 2017/06/29, 2017/07/11, 2017/07/23, 2017/08/04, 2017/08/16, 2017/08/28, 2017/09/09, 2017/09/21, 2017/10/03, 2017/10/15, 2017/10/27, 2017/11/08, 2017/11/20, 2017/12/02, 2017/12/14, 2017/12/26, 2018/01/07, 2018/01/19, 2018/01/31, 2018/02/12, 2018/02/24, 2018/03/08, 2018/03/20, 2018/04/01, 2018/04/13, 2018/04/25, 2018/05/07, 2018/05/19, 2018/05/31, 2018/06/12, 2018/06/24, 2018/07/06, 2018/07/18, 2018/08/11, 2018/08/23, 2018/09/04, 2018/09/16, 2018/09/28,2018/10/10, 2018/10/22, 2018/11/03, 2018/11/15, 2018/11/27, 2018/12/21, 2019/01/02, 2019/01/14, 2019/01/26, 2019/02/07, 2019/02/19, 2019/03/03, 2019/03/15, 2019/03/27, 2019/04/08, 2019/04/20, 2019/05/02, 2019/05/14, 2019/05/26, 2019/06/07, 2019/06/19, 2019/07/01, 2019/07/13
	Orbit	T33
	Heading	−169.73°
	Incidence	33.69°
	Pixel Spacing (Rg × Az)	2.33 × 13.96 m
	Number of images	71
	Direction	Ascending	2017/01/24, 2017/02/05, 2017/02/17, 2017/03/25, 2017/04/06, 2017/04/18, 2017/04/30, 2017/05/12, 2017/05/24, 2017/06/05, 2017/06/17, 2017/06/29, 2017/07/11, 2017/07/23, 2017/08/04, 2017/08/16, 2017/08/28, 2017/09/21, 2017/10/03, 2017/10/15, 2017/10/27, 2017/11/08, 2017/11/20, 2017/12/02, 2017/12/14, 2018/12/26, 2018/01/07, 2018/01/19, 2018/01/31, 2018/02/12, 2018/02/24, 2018/03/08, 2018/03/20, 2018/04/01, 2018/04/13, 2018/04/25, 2018/05/07, 2018/05/19, 2018/05/31, 2018/06/12, 2018/06/24, 2018/07/18, 2018/07/30, 2018/08/11, 2018/08/23, 2018/09/04, 2018/09/16, 2018/09/28, 2018/10/10, 2018/10/22, 2018/11/03, 2018/11/15, 2018/11/27, 2018/12/21, 2019/01/02, 2019/01/14, 2019/01/26, 2019/02/07, 2019/02/19, 2019/03/03, 2019/03/15, 2019/03/27, 2019/04/08, 2019/04/20, 2019/05/02, 2019/05/14, 2019/05/26, 2019/06/07, 2019/06/19, 2019/08/18
	Orbit	T26
	Heading	−10.05°
	Incidence	33.81°
	Pixel Spacing (Rg × Az)	2.33 × 13.98 m
	Number of images	73
Dongguan, Guangdong Province (Urban area)	Direction	Ascending	2018/12/08, 2018/12/20, 2019/01/01, 2019/01/13, 2019/01/25, 2019/02/06, 2019/02/18, 2019/03/02, 2019/03/14, 2019/03/26, 2019/04/07, 2019/04/19, 2019/05/01, 2019/05/13, 2019/06/06, 2019/06/18, 2019/06/30, 2019/07/12, 2019/07/24, 2019/08/05, 2019/08/17, 2019/08/29, 2019/09/10, 2019/09/22, 2019/10/04, 2019/10/16, 2019/10/28, 2019/11/09, 2019/11/21, 2019/12/03, 2019/12/15, 2019/12/27
	Orbit	T11
	Heading	−10.55°
	Incidence	34.04°
	Pixel Spacing (Rg × Az)	2.33 × 14.00 m

.

3.2. MT-InSAR Data Processing

Considering the small perpendicular baseline distribution of Sentinel-1 satellites, we only need to consider the temporal baseline for choosing interferometric pairs. Each image is combined with its following two images to form interferometric pairs (Figure A1). To reduce the noise and average the resolution in the azimuth and range directions, a multi-look operation (range:azimuth = 5:1) is applied. Adaptive Goldstein filter is performed to reduce the phase noise and the Minimum Cost Flow (MCF) method is used to conduct phase unwrapping. We use polynomial fitting to remove the phase ramp in each interferogram. The mean coherence, intensity, and amplitude dispersion of each multi-looked pixel are used to choose high-quality pixels for the subsequent time-series displacement calculation. Multiple linear regressions are performed to estimate the topographic error of each pixel. In addition, the topographic error of each pixel is subtracted from the unwrapped phase. Then, singular value decomposition (SVD) is used to obtain the time series phases of each pixel [6]. We use the least square to estimate the linear deformation rate of each pixel, and subtract the linear deformation from the unwrapped phase. The components of the remaining phases are nonlinear deformation, atmospheric delay, and noise. Atmospheric delay phases are highly correlated in space while lowly correlated in time. Noise phases are lowly correlated in space and time. Thus, we use the temporal low-pass filter to remove atmospheric delay and noise and obtain the nonlinear deformation. At last, the final time series deformation is obtained by combining the linear deformation and nonlinear deformation.

The average deformation rate of each pixel can be estimated by the least square method. The deformation results, including the experimental area deformation map, the deformation rate of each point, and the deformation sequence information, are used for the automatic identification of active deformation regions.

3.3. Decomposition of the Slope Deformation

In urban areas, the surface deformation is dominated by subsidence, so using a single observation geometry can identify potential hazards well. In the case of Dongguan City, we test the improved method by a single orbit SAR observation. However, in the mountainous areas with complex terrains, using a single observation geometry can hardly obtain the complete deformation data, especially for the landslide monitoring, because the deformation of landslide is very sensitive to slope, aspect, and incidence angle [45]. In this study, we jointly use the ascending and descending data to retrieve two-dimensional (2D) deformation in the aspect (horizontal) and vertical directions.

The deformation in the line of sight (LOS) direction, ${\mathrm{d}}_{\mathrm{Los}}$, can be obtained from MT-InSAR, which is then decomposed into the deformation in the vertical direction ${\mathrm{d}}_{\mathrm{U}}$, the north–south direction ${\mathrm{d}}_{\mathrm{N}}$, and the east–west direction ${\mathrm{d}}_{\mathrm{E}}$. The relationship between ${\mathrm{d}}_{\mathrm{Los}}$ and the three deformation components can be expressed as follows [46]:

$${\mathrm{d}}_{\mathrm{Los}}={\mathrm{d}}_{\mathrm{U}}\cos {\mathsf{\theta }}_{\mathrm{inc}}-{\mathrm{d}}_{\mathrm{N}}{\sin \mathsf{\theta }}_{\mathrm{inc}}\cos ({\mathsf{\alpha }}_{\mathrm{azi}}-\frac{3\mathsf{\pi }}{2})-{\mathrm{d}}_{\mathrm{E}}\sin {\mathsf{\theta }}_{\mathrm{inc}}\sin \left({\mathsf{\alpha }}_{\mathrm{azi}}-\frac{3\mathsf{\pi }}{2}\right),$$

(1)

where ${\mathsf{\theta }}_{\mathrm{inc}}$ and ${\mathsf{\alpha }}_{\mathrm{azi}}$ represent the radar incidence angle and satellite heading angle, respectively, so ${\mathsf{\alpha }}_{\mathrm{azi}}-3\mathsf{\pi }/2$ is the angle between the north direction and range direction.

In this study, we use the monitoring results from ascending and descending data to decompose the 2D deformation. The horizontal deformation of the landslide usually occurs on the slope direction, so we assume that the landslide deformation can be decomposed into the horizontal deformation ${\mathrm{d}}_{\mathrm{hor}}$ and ${\mathrm{d}}_{\mathrm{U}}$. The decomposition relationship can be expressed by Formula (2):

$$\left(\begin{array}{c}{\mathrm{d}}_{\mathrm{U}}\\ {\mathrm{d}}_{\mathrm{hor}}\end{array}\right)=\left(\begin{array}{ll}\cos {\mathsf{\theta }}_{\mathrm{inc}}^{\mathrm{A}}& -{\sin \mathsf{\theta }}_{\mathrm{inc}}^{\mathrm{A}}\cos \left(\mathsf{\delta }-{\mathsf{\alpha }}_{\mathrm{azi}}^{\mathrm{A}}+\frac{3\mathsf{\pi }}{2}\right)\\ \cos {\mathsf{\theta }}_{\mathrm{inc}}^{\mathrm{D}}& -{\sin \mathsf{\theta }}_{\mathrm{inc}}^{\mathrm{D}}\cos (\mathsf{\delta }-{\mathsf{\alpha }}_{\mathrm{azi}}^{\mathrm{D}}+\frac{3\mathsf{\pi }}{2})\end{array}\right)\left(\begin{array}{c}{\mathrm{d}}_{\mathrm{Los}}^{\mathrm{A}}\\ {\mathrm{d}}_{\mathrm{Los}}^{\mathrm{D}}\end{array}\right),$$

(2)

where ${\mathsf{\theta }}_{\mathrm{inc}}^{\mathrm{A}}$ and ${\mathsf{\theta }}_{\mathrm{inc}}^{\mathrm{D}}$ are the radar incidence angles from the ascending and descending orbits, respectively; ${\mathsf{\alpha }}_{\mathrm{azi}}^{\mathrm{A}}$ and ${\mathsf{\alpha }}_{\mathrm{azi}}^{\mathrm{D}}$ are the satellite heading angles from the ascending and descending orbits, respectively; $\mathsf{\delta }$ represents the aspect of the point. Based on the monitoring results of the ascending and descending data, the 2D deformation in the vertical direction and along the slope can be obtained. The final deformation is the vector sum of the horizontal deformation and the vertical deformation. The relationship between the three deformations can be expressed by Formula (3):

$${\left|{\mathrm{d}}_{\mathrm{fin}}\right|}^2={\left|{\mathrm{d}}_{\mathrm{hor}}\right|}^2+{\left|{\mathrm{d}}_{\mathrm{U}}\right|}^2,$$

(3)

where ${\mathrm{d}}_{\mathrm{fin}}$ is the final deformation calculated.

4. Results and Analysis

4.1. Deformation Extraction Results of Lajia Town

Using the SBAS-InSAR technology, we obtained the ground deformation in Lajia Town from ascending and descending data, shown in Figure 8a,b. We selected points by three thresholds, which are the coherence threshold 0.5, the intensity threshold 0.5 and the amplitude dispersion threshold 1.0. We obtained 409,072 and 464,444 target points from the ascending and descending data, respectively, with a density of 2678/km² and 3041/km². The automatic extraction method was employed to investigate the potential hazards based on the MT-InSAR results, where ${\mathsf{\sigma }}^A$ of the ascending results was 3.6 mm/yr and ${\mathsf{\sigma }}^D$ of the descending results was 5.0 mm/yr. Finally, we identified 14 and 18 deformation regions in the results of ascending and descending data, respectively. Combining the results of these two orbital data, 23 deformation regions were determined, with an area of 7.7 km² (each area 0.05~2.30 km²). As Figure 8c,d shows, the extraction results of different orbital data are quite different.

Substituting the basic parameters, such as the incidence angle, heading angle, and landslide aspect angle into formula (2), the deformation velocities in the horizontal and vertical direction can be calculated (Figure 8c,d). The total number of the identified target points is 26, 982 (388~8644 points in each region). The maximum annual average deformation along the horizontal direction has reached 100 mm/yr, and the maximum annual average deformation velocity along the vertical direction is 30 mm/yr. As the terrain is rugged, we only decompose the 2D deformation of the points in the deformation regions. This calculation greatly depends on topographic factors, so some deformation points obtained by those two different geometries can be quite different.

The final deformation velocity of the detected areas is shown in Figure 9a, which is also the basis for the velocity gradation. After identifying the deformation regions, the method will automatically complete the velocity gradation, trend estimation, instability degree generation, and deformation boundary refinement. The point deformation trend is estimated on the basis of the time series from the ascending and descending data, respectively, and the method selects a higher-level trend as the trend estimation result at the points. The results are shown in Figure 9b–d.

Table 3. Basic information and assessment results of the deformation regions in Lajia Town.

Index	Level	Low	Moderate	High	Very High
Velocity gradation	Number of points	3049	13,093	10,863
	Proportion (%)	11.3	48.5	40.3
Trend	Number of points	1379	554	25,072
	Proportion (%)	5.1	2.1	92.9
Instability	Number of points	1192	3274	12,111	10,428
	Proportion (%)	4.4	12.1	44.9	38.7

Note: The text in italics and bold indicates the level of the largest proportion of deformed points.

shows the assessment results based on deformation points. Most potential hazards maintain the deformation accelerating trend and have a high or very high instability level. In Lajia Town, the high instability level deformation points have the largest proportion (44.89%), followed by the very-high level (38.65%). On the whole, the test area is in an active state, which needs timely monitoring and attention.

In Lajia Town, there are three regions with very high instability levels, which could be potential landslides, as shown in Figure 9a. We performed a detailed analysis on them and the results are shown in Figure 10 (Landslide 1, 2, and 3 in Figure 9). Landslide 1 in Figure 10a is close to the residential areas, and it spans from the flanks of the valley to Lajia Town at the bottom of the valley, posing great threats to local people. For Landslide 1, the deformation along the horizontal direction is much larger than that in the vertical direction, as the landslide is creeping continuously towards the bottom of the slope. It has a deformation velocity exceeding 70 mm/yr. We select the time series of three points from a factory (P1), a residential area (P2), and the road (P3), to estimate their deformation trends. Precipitation is one of the important factors for most landslides. The monthly precipitation data (January 2017−July 2019 from https://gpm.nasa.gov/) in the study area [47] show that most rain falls from May to November, and less from December to April. The largest precipitation is from September to October. Figure 10d shows that the deformation trends of the three points fit the EAM model well with the maximum cumulative deformation exceeding 180 mm, indicating that the deformation is continuously accelerating during the monitoring period. The landslide has no obvious acceleration during the rainy season, and the correlation between its deformation and precipitation is weak, indicating that precipitation is not the dominant factor in landslide deformation. The deformation of the far river section of the landslide was greater than that in the near river section. P3, which is located in a very high instability area, has a cumulative deformation > 180 mm. The deformation center of the region with a very high instability in Figure 10 (outlined by purple dashed line) is about 2.17 km length and 0.82 km width, with an area of more than 1.2 km², accounting for 95% of the entire landslide. A road traverses the most unstable area in the middle of the landslide and leads directly to the residential area near the river. The time-series results also show that the landslide deformation process is in an accelerated deformation stage, which brings dangers to the residents and traffic inside the landslide. The other two landslide areas (Figure 10b,c) are distributed along the river. The highest instability areas are located in the center of the landslides, and landslide 2 is also located near the residential area. We plot the time series displacements of two points in landslide 2 and 3 (P4, P5), which shows that the annual deformation of the two points exceeds 40 mm, and the acceleration does not coincide with the increase in precipitation. These areas have high or very high instability. Therefore, continuous monitoring is necessary for those landslides.

4.2. Deformation Extraction Results of Dongguan City

In the Dongguan experiment, using the same thresholds for point selection in Lajia Town, we obtained a total of 4,101,295 target points, with a density of 1662/km². The $\mathsf{\sigma }$ of the results was 3.2 mm/yr. Therefore, the values of 3$\mathsf{\sigma }$ (9.6 mm/yr) and 6$\mathsf{\sigma }$ (19.2 mm/yr) were adopted to classify the deformation velocity as low, moderate, and high. Using the parameters obtained in Section 2.1, 221 deformation regions were automatically identified, with a total area of 40.663 km² (each area 0.098~3.834 km²). The total number of the identified target points is 67,204 (each area 85~7908 points).

The obtained deformation map is shown in Figure 11a, and its settlement shows a scattered distribution pattern, with greater deformation in the suburb and urban fringe. More than 80% of the deformation regions have a moderate intensity, but a few deformation regions have a deformation velocity exceeding 30 mm/yr (Figure 11b,c).

Table 4. Basic information and assessment results of the deformation regions in Dongguan City.

Index	Level	Low	Moderate	High	Very High
Velocity gradation	Number of points	47,837	76,816	8957
	Proportion (%)	35.9	57.6	6.7
Trend	Number of points	47,401	8932	77,277
	Proportion (%)	35.5	6.7	57.9
Instability	Number of points	50,849	19,261	54,975	8525
	Proportion (%)	38.1	14.4	41.2	6.4

Note: The text in italics and bold indicates the level of the largest proportion of deformed points.

lists the results calculated from the points in each area. Specifically, more deformation areas (111 regions, 50.2% of the deformation regions) are around the parks or reservoir slopes on the southern boundary of the urban area, which are covered by scrubs and vegetation. The velocity level of the points in these regions is mostly low or moderate. A few deformation areas are distributed along the river in construction areas and residential areas (42, 19.0%), farmland irrigation areas (49, 22.2%), and mixed areas (including crubs, paddy field areas, and construction areas) (19, 8.6%), most of which have moderate and high intensity (Figure 11b,c and Table 4). Although most of the deformation points have moderate intensity, they have an accelerating trend. Therefore, the high instability points account for the largest proportion.

Since most of the detected areas in Dongguan have high instability, we select two for the field survey to validate the improved method (Area 1 and 2 in Figure 12). The two detected areas are mainly composed of paddy fields, factory buildings, and construction sites. Several points show high velocity, accelerating trend, and very high instability on the factory building and construction site. Area 1 (Figure 12a,d) is close to the Pearl River estuary with three very high unstable blocks (outlined by the purple line in Figure 12a) and a long-time displacement. The ground subsidence caused a slight separation between the nearby road surface and the building foundation, which in turn led to the cracks and deformation of buildings (such as walls and pillars). Area 2 (Figure 12b,e) is located near a welding workshop, whose north and east sides are large open-air parking lots, and the west side is warehouses and roads. This area is close to the river, and it has thick, soft soil and an unstable foundation. Therefore, the separation between the road and the workshop building foundation is serious. P6 and P7 are the points in the very-high instability areas, showing a tendency of deformation acceleration, whose settlement exceeds 60 mm and 80 mm, respectively. In addition, some deformed regions are distributed in the farmland irrigation or paddy field with high instability, in which the groundwater extraction may cause consolidation of the soil layer.

In general, the surface subsidence in Dongguan City is due to the unstable strata and human activities. Dongguan City is a coastal city and has soft soil layers, so there are many regions that have deformation due to external influences. Additionally, Dongguan City is an important economic city in China, so human activities are frequent, such as engineering construction, factory operations, and groundwater extraction. These factors, together, contribute to the continuous subsidence [12,48], which deserve timely monitoring.

5. Discussion

The improved method can provide robust deformation results. DInSAR technology is prone to temporal decorrelation and atmospheric delay [18]. MT-InSAR results have higher accuracy and spatial consistency. The first ADA extraction algorithm proposed by Barra et al. [25] is based on PSI technology, which could generate and update the list of geological hazards. Furthermore, some researchers used satellite images of single orbits, so they can only obtain the LOS deformation, which is insufficient for landslide research in mountain areas, because of the rough terrain. Jointly using both ascending and descending data [28] to automatically detect potential deformation areas in mountain areas can improve the accuracy, so we employed both ascending and descending data in this research. This research also eliminates the points with poor spatial consistency to avoid misclassifying deformation areas, which is useful for the MT-InSAR result caused by the temporal decorrelation in vegetation and water area.

The improved method provides trend estimation in the time domain. The deformation information in the time dimension is very important for predicting the deformation development trend and providing disaster warnings. Some methods have identified anomalous changes to obtain different trends based on the time series deformation [34,36]. In this study, we used four trend evolution models to fit time series deformation, which can effectively avoid the influence of anomalous deformation caused by atmospheric delay. Furthermore, the improved method can show the development trend of the deformation points in the time domain, and obtain the deformation development trend for disaster prevention.

The improved method also provides instability evaluation for the deformation region in the spatiotemporal domain for the deformation area and accurate locations of highly unstable areas. To assess the instability of geological hazards, we should consider both the temporal and spatial displacement characteristics of the region. Some methods focus on the vulnerability and exposure assessment, including the quantification of the potential loss by a hidden landslide [29], the GIS-based procedure for analyzing specific risk [32]. In addition to using external data or models, some researches may only utilize the deformation characteristics in the spatial dimension. In this study, we combine the velocity gradation and the time evolution characteristics to establish a qualitative instability map, which can automatically obtain the spatiotemporal characteristics of deformation points. The instability refinement experiments show that the improved method can classify the deformation regions by instability and locate the highly unstable regions. In addition, the improved identification algorithm is based on clustering of active points, so it greatly improves the identification efficiency, computing speed, and automation degree, and is suitable for recognizing and assessing the deformation area obtained by MT-InSAR technology. Therefore, the improved method can be widely used in the identification of potential hazards in large-scale areas.

The improved method is also impacted by the inherent limitations of MT-InSAR measurements, such as shadow overlay [20,28] and insensitivity to the north–south deformation [9]. Actually, besides the four evolutionary categories (EAM, ULM, EMM, NSC), there are displacements caused by sudden changes or periodic deformation trends, which are excluded or classified into one of the categories by our method. Furthermore, the estimation of deformation development trend estimation is strongly dependent on the number of available SAR data. Less image samples in the time domain would lead to low accuracy of the development trend estimation. Additionally, the improved method only identifies the potential hazards in large-scale monitoring, and the subsequent monitoring of a single hazard requires external data for assistance. Using external auxiliary data, such as regional DEM and geological data, may improve the deformation area identification and evaluation, which also can improve the result interpretation of deformation regions.

6. Conclusions

In this study, we introduced an improved method for automatic identification and assessment of potential geohazards, which can show the spatial and temporal information of the deformation, and quantitatively evaluate the instability of the deformed area. Its final outputs include a deformation velocity map, a deformation regions distribution map, the results of velocity gradation and trend evaluation, an instability result that consists of the instability level and statistical information, and the results of boundary refinement.

The method has been tested in the mountainous area of Lajia Town in Qinghai Province, China, and the urban area of Dongguan City in Guangdong Province, China. In Lajia Town, it identified 23 landslides and found that most potential hazards have very high-level instability. In the point-based instability evaluation, the high-level unstable deformation points take the largest proportion (44.9%). In Dongguan City, 221 potential deformation regions were identified, distributing less in the city center and more in the suburb. Most deformation points in the detected regions (41.2%) have high-level instability, which is mostly located in construction sites, residential areas, and paddy fields. The deformation in Dongguan City is mainly caused by soil consolidation from frequent engineering construction and human activities. The deformation areas of Lajia Town have the highest instability in the center and the instability becomes lower gradually outward. However, the instability of the deformation areas in Dongguan City shows no specific rule.

Using the InSAR monitoring results, this method improves the large-scale deformation monitoring capabilities of InSAR technology. It can extract the deformation area from different terrains. It has advantages in urban areas, where the deformation points are scattered and difficult for manual interpretation. The instability refinement method can identify the highly unstable areas caused by engineering construction. It also helps users focus on the regions deserving more attention, so as to provide a reference for disaster prevention, reduction, and relief.