Landslide Identification and Gradation Method Based on Statistical Analysis and Spatial Cluster Analysis

Abstract

As a type of earth observation technology, interferometric synthetic aperture radar (InSAR) is increasingly widely used in the field of geological disaster detection. However, the application of InSAR in low-coherence areas, such as alpine canyon areas and vegetation coverage areas, is subject to considerable limitations. How to accurately identify landslides from InSAR measurement data in these areas remains the subject of several challenges and shortcomings. Based on statistical analysis and spatial cluster analysis, in this paper, we propose an automatic landslide identification and gradation method suitable for low-coherence areas. The proposed method combines the small baseline subset InSAR (SBAS-InSAR) method and the interferogram stacking (stacking-InSAR) method to obtain a deformation map in the study area, using statistical analysis and spatial cluster analysis to extract deformation regions and landslide polygons to propose a landslide screening model (LSM) based on multivariate features to screen landslides and reduce the interference of noise in landslide identification, in addition to proposing a landslide gradation model (LGM) based on signum function to grade the identified landslides and provide support to distinguish landslides with different deformation degrees. The method was applied to landslide identification in the upper section of the Jinsha River basin, and 47 potential landslides were identified, including 15 high-risk landslides and 13 landslides endangering villages. The experimental results show that the proposed method can identify landslides accurately and hierarchically in low-coherence areas, providing support for geological hazard investigation agencies and local departments.

1. Introduction

Due to the intensification of modern human activities, the global climate is gradually warming, which has led to a considerable increase in the frequency of landslides [1,2]. In the past 30 years, the occurrence of many large-scale landslide disasters has caused considerable casualties and property losses. The effective identification and subsequent continuous monitoring of unknown hidden landslides can protect people’s lives and property to a considerable extent.

Interferometric synthetic aperture radar (InSAR) is a new space observation technology that has been used to monitor surface deformation since the 1990s; it has also been widely used in the field of geological hazard monitoring based on its outstanding advantages, such as non-contact, all-weather, all-time, high-precision, and wide coverage operation [3,4,5]. Multitemporal InSAR (MT-InSAR) is a general term for a variety of InSAR techniques that can obtain long time-series deformation, represented by Persistent Scatterer InSAR (PS-InSAR) [6,7] and SqueeSAR [8] proposed by Ferretti, small baseline subset InSAR (SBAS-InSAR) [9,10] proposed by Berardino, the Stanford method for persistent scatterers (StaMPS) [11] proposed by Hooper, and interferogram stacking (stacking-InSAR) [12] proposed by Sandwell. By using mutitemporal SAR images, MT-InSAR can minimize the influence of spatiotemporal decorrelation [13], atmospheric delay [14,15], and satellite orbit error [16] of traditional differential InSAR (D-InSAR) [17] and can obtain high-precision and long time series of surface deformation results. The time series deformation information obtained by MT-InSAR can be used to retrieve the formation mechanism and kinematic evolution process of landslides, including landslide types, trigger factors, failure modes, depth and volume estimation, and risk assessment [18,19,20,21,22,23,24]. Scholars have carried out a considerable among research on the application of MT-InSAR for landslide monitoring. Some scholars have analyzed the influencing factors of landslides in various areas based on the results of MT-InSAR and found that rainfall is the main factor affecting the occurrence of landslides, whereas river distribution, fault zone distribution, slope, and slope direction are also important factors [25,26,27,28]. Some scholars have detected and identified landslides in large areas based on MT-InSAR results and verified the effectiveness of MT-InSAR in landslide identification combined with measured data [29,30]. Some scholars have uses InSAR to predict landslides, including early detection of unstable landslides, spatiotemporal prediction of large-scale landslides, risk assessment of landslide collapse, etc. [31,32,33].

The identification and detection of landslides in low-coherence areas, such as alpine canyon areas and vegetation coverage areas, have always been a difficult problem. On the one hand, large noise error occurs in low-coherence areas; on the other hand, the lack of data caused by decorrelation leads to landslide leakage. In the natural environment, the lack of persistent scatterers limits the application of PS-InSAR, although the SBAS-InSAR method makes up for this defect. However, the SBAS-InSAR method is also affected by spatiotemporal decorrelation. Excessive deformation or low coherence result in a lack of high-coherence points in the study area using SBAS and a consequential lack of results. Some scholars have attempted to use the stacking-InSAR method to improve the results in low-coherence areas [34,35,36]. This method suppresses the influence of random phases by weighted averaging of interferograms and does not involve the selection of coherence points so that the landslide displacement information can be better preserved in low-coherence areas. Furthermore, the technical premise of stacking-InSAR assumes that there is only linear deformation in the study area, which is consistent with the mainstream translation landslide deformation model [37]; therefore, stacking-InSAR has considerable advantages in the field of geological hazard monitoring.

In the identification of landslides using InSAR involves two important concepts: the extraction of deformation regions and the identification of suspicious landslides. Threshold segmentation is the most commonly used method to extract the deformation regions at present [38,39]. However, the geological structure, lithology, and surface cover differ depending on the region. Using the same threshold to extract the deformation from different regions leads to inaccurate results, and the use of different thresholds by different researchers results in landslide identification results that are considerably influenced by the subjectivity. Invariant thresholds are no longer used to extract deformation regions and identify landslides; instead, adaptive thresholds based on statistical characteristics of the phase are used [40,41]. Given the need for landslide identification in large areas, some automatic procedures have been proposed according to statistical indicators to automatically identify and classify landslides [42,43,44]. In terms of the identification of suspicious landslides, most of current research still employs the human visual interpretation method to identify landslides, which is not suitable for the identification and continuous monitoring of landslides in large areas. Deep learning methods have also been attempted for landslide identification and detection. For example, mapping of regional landslide susceptibility maps based on deep learning can provide a reference for landslide risk prevention in small areas [45,46]. Some morphological methods, such as hotspot analysis, are used to process deformation results, which can extract landslides with spatial aggregation characteristics [47,48]. These methods can aid in the automatic identification of landslides in a given area, but there are also some drawbacks: deep learning methods rely on the quality of prior datasets, and morphological methods cannot eliminate landslide false alarms caused by phase noise.

Based on these existing problems and previous studies, an automatic landslide identification and classification method based on statistical analysis and spatial cluster analysis is proposed in this study, which is suitable for landslide identification in low-coherence areas, such as alpine canyon areas and vegetation coverage areas. The primary contributions of this method include the following:

• SBAS-InSAR and stacking-InSAR are combined to generate deformation maps in the study area, which can obtain more information to avoid the omission of landslides in the study area;
• Statistical analysis and spatial cluster analysis are used to extract deformation regions and landslide polygons, which can provide a unified benchmark for landslide identification in different regions;
• A landslides screening model (LSM) based on multivariate features is proposed to screen landslide polygons, which can reduce the interference of noise in landslide identification; and
• A landslides gradation model (LGM) based on the signum function is proposed to divide screened landslides into four risk grades for geological disaster monitoring and risk assessment to provide effective information.

The structure of this paper is arranged as follows. In Section 2, we describe the principle and details of the proposed method. In Section 3, we introduce the experimental area. In Section 4, we demonstrate the specific implementation steps of the proposed method in the Jinsha River basin and verify the validity of the proposed model. In Section 5 and Section 6, we discuss the obtained results and present our conclusions.

2. Methods

In this section, we describe how to accurately and reliably identify potential landslides in low-coherence areas, such as alpine canyon areas and vegetation coverage areas. Figure 1 shows the workflow of the proposed method. The specific process can be divided into the following three parts.

2.1. InSAR Processing and Deformation Results Acquisition

Through D-InSAR processing, SBAS-InSAR processing, and stacking-InSAR processing, three deformation maps of the study area are obtained, including the annual velocity map obtained by SBAS, a time-series cumulative displacement map obtained by SBAS, and an annual velocity map obtained by stacking.

• D-InSAR processing.

GMTSAR software [49] is used to process Sentinel-1 images to generate interferograms. Based on the principle of the shortest total temporal baseline, the super master image is selected for registration. The spatiotemporal baseline between each image is calculated according to precise orbit determination ephemerides, which are used to select interferogram pairs according to the SBAS-InSAR algorithm strategy. The interferograms are unwrapped using the statistical-cost network-flow algorithm for phase unwrapping (SNAPHU) [50,51]. The topographic phase is removed using a digital elevation model (DEM) from the Shuttle Radar Topography Mission [52].

2.

SBAS-InSAR processing.

The basic principle of SBAS-InSAR is to minimize the error and spatiotemporal decorrelation caused by the difference in viewing angles between interferogram pairs by setting the spatiotemporal baseline threshold, and the time-series deformation results are obtained by least square solution and singular value decomposition (SVD). Through the SBAS program encapsulated by GMTSAR, four results are obtained: annual velocity, cumulative displacement, root mean square error (RMSE) during linear deformation fitting, and estimated DEM error (DEME).

3.

Stacking-InSAR processing.

Stacking-InSAR can further weaken the influence of orbit, atmosphere, and terrain errors by weighting together N unwrapped interferograms [34,53]. Stacking-InSAR assumes that the annual velocity is the ratio of cumulative displacement to time and does not consider the distribution of other noise elements in the process of phase unwrapping with different interferograms; that is, only linear deformation is considered in the study area. The mainstream landslide sliding model assumes that the landslide moves along the steepest downhill slope [37], constituting a linear deformation model. Therefore, using stacking-InSAR technology to identify landslides is a more effective method. The annual velocity in the time series is calculated using Equation (1).

$$V_{stacking}=\frac{\lambda }{4\pi }\cdot \frac{{\sum }_{i=1}^N{\phi }_i\Delta T_i}{{\sum }_{i=1}^N\Delta T_i^2}$$

(1)

where $V_{stacking}$ is the annual velocity in mm/yr, λ is the radar wavelength in millimeters, N is the number of interferogram pairs, ${\phi }_i$ is the unwrapped interferometric phase, and $\Delta T_i$ is the temporal baseline length of interferogram pairs in years.

2.2. Deformation Region Extraction and Spatial Cluster Analysis

SBAS-InSAR considerably reduces the influence of spatiotemporal decorrelation and atmospheric delay by setting a baseline threshold and spatiotemporal filtering. Stacking-InSAR also weakens the influence of orbital error, terrain error, and atmospheric error in interferograms by stacking. The residual noise phase after application of the above noise reduction methods has random characteristics.

According to the central limit theorem, the distribution of a large number of random variables is close to the normal distribution. Based on the assumption of normal distribution, deformation regions and non-deformation regions are extracted from the statistical indicators (mean and standard deviation; μ and σ, respectively) of three deformation maps. Assuming that an area without any deformation, the results of surface deformation obtained by InSAR should obey the distribution of N(0,σ). Due to the existence of a series of noise reduction methods and random noise phases, the actual obtained distribution is close to the normal distribution (N(μ,σ)). For the normal distribution, the random variables distributed in the range of μ ± 3σ account for 99.74%. After calculating μ and σ of the deformation map according to this theory, the results outside the range of μ ± 3σ are categorized as deformation regions (DR), and the results within the range of μ ± 3σ are categorized as non-deformation regions (non-DR). The deformation regions map is obtained, as well as the distribution histogram of deformation results and its fitting normal distribution curve, as shown in Figure 2.

Spatial cluster analysis is a method to quickly screen points with aggregated deformation, which can remove the interference of discrete points. A landslide is a collection of clustered points, so the spatial cluster analysis method can be used to identify suspicious landslide areas in a deformation regions map. The hotspot analysis method based on the Getis-Ord Gi* statistics [54] and the point aggregation method are used for spatial cluster analysis to obtain landslide polygons in the study area. As a hypothesis-testing method, the hotspot analysis method judges whether deformation points are significantly aggregated by checking the difference between adjacent elements within the specified critical distance range; this method can identify spatial clusters with large deformation on a deformation regions map. The hotspot analysis method uses the Getis-Ord Gi* statistic [47] to determine the spatial aggregation of the deformation points. Equation (2) can calculate a value of G_i*(d) for each deformation point; the larger the absolute value of G_i*(d), the stronger the aggregation of deformation points. In the present study, the confidence of the hypothesis test is set to 0.99, and the corresponding critical value (|G_i*(d)|) of the deformation point is 2.58.

$$G_i^{*}(d)=\frac{{\sum }_{j=1}^{n_{ij}}{x}_j+x_i-n_{ij}\overline{x}}{S^{*}\cdot \sqrt{\frac{n\cdot n_{ij}-n_{ij}^2}{n-1}}}$$

(2)

where n represents the total number of deformation points; n_ij represents the number of deformation points within the search distance (d) of the i-th deformation point (i = 1, 2, …, n; j = 1, 2, …, n_ij; j ≠ i); x represents the value of a deformation point (annual velocity and cumulative displacement); and $\overline{x}$ and $S^{*}$ represent the mean and standard deviation of all deformation point values, respectively. With a confidence of 0.99, if the |G_i*(d)| value of a deformation point is greater than 2.58, the deformation point is considered to have strong spatial aggregation.

After hotspot analysis, all feature points with high aggregation are extracted, but there are still a small number of discrete points. The point aggregation method is used to further eliminate these discrete points. The point aggregation method aggregates the points into faces by setting a certain distance threshold, which can extract landslide polygons with aggregation features on a deformation regions map, as shown in Figure 3.

Hotspot analysis and point aggregation are clustered based on the distance between deformation points. Selecting an appropriate distance threshold is the key to correctly identifying landslides. Using the smallest possible distance threshold for hotspot analysis and point aggregation can maximize the spatial differences between deformation point clusters.

2.3. Landslide Screening Based on LSM

After spatial cluster analysis of the deformation regions extracted from these three InSAR results, the landslide polygons of each deformation result are obtained, although most are phase noise, not real landslides. To remove these noises, LSM-based multivariate features is proposed to screen all landslide polygons. Figure 4 shows LSM-screened landslide polygons based on six features: time-series average coherence (γ_mean) [55], DEME, RMSE, slope, area, and optical image features.

The time-series average coherence (γ_mean) from the coherence maps of all interferogram pairs is calculated by Equation (3).

$${\gamma }_{\mathrm{mean}}=\sqrt{\frac{1}{\frac{1}{m}{\sum }_{k=1}^m\left(\frac{1-{\gamma }_k^2}{{\gamma }_k^2}+1\right)}}$$

(3)

where γ_k is the coherence map of each interferogram pair, and m is the number of interferograms. DEME and RMSE are the results generated by the SBAS program encapsulated by GMTSAR software. According to the statistics [56], landslides mostly develop on terrain slopes in the range of 25° to 45°. The slope threshold of landslide polygons is set to 25°, and the average slope value is calculated based on DEM data. The area threshold is set to 0.1 km² to screen landslide polygons larger than medium size. For each landslide polygon obtained based on SBAS-InSAR data, the μ and σ of the five indicators are calculated and used as thresholds to screen all landslide polygons. For landslide polygons obtained by stacking-InSAR, because there is no DEME and RMSE, only three indicators are used for screening: γ_mean, area, and slope. To avoid missing dangerous landslides, the landslide polygons with maximum deformation values outside the range of μ ± 6σ are also screened.

After the above processing steps, it is also necessary to observe whether the polygonal area of the landslide shows the topographic and geomorphic characteristics of landslides in optical images [57,58,59,60]. Unstable landslides usually exhibit some macroscopic features before sliding (downward gentle dip and unleveled phenomena; new gullies on the surface of landslides; loose soil and rocks in front of landslide and face the danger of river erosion; the surface of the landslide is bare and not covered by vegetation). Google Earth optical images are used to judge whether the landslide polygon has any of the above characteristics, and landslide polygons are further screened to obtain a landslide map of the study area.

2.4. Landslide Gradation Based on LGM

After the above screening process, the screened landslides are obtained based on the three deformation results and spatially fused to obtain a landslide map of the study area, as shown in Figure 5.

It is necessary to evaluate the grade of the landslide to estimate the risk of disasters after obtaining the basic landslide information. In this paper, LGM based on the signum function is proposed to grade screened landslides, as shown in Figure 6.

The signum function adds a graded label corresponding to three deformation maps to each landslide according to Equation (4).

$$classif{y}_i=\{\begin{array}{c}1,\mathrm{others}\\ 0, \mu -6\sigma <max<\mu +6\sigma \end{array}$$

(4)

$$grade={\sum }_{i=1}^{3}classif{y}_i$$

(5)

Here, μ ± 6σ is set as the threshold of deformation gradation [61]; max refers to the maximum value for the positive value of LOS direction and the minimum value for the negative value of LOS direction. The gradation result of each landslide (Equation (5)) is obtained by superimposing graded labels corresponding to three deformation maps.

According to the gradation results, the screened landslides are divided into four grades: L-I, L-II, L-III, and L-IV. The higher the grade, the greater the degree of landslide sliding and the greater the risk of landslide disaster. It is necessary to strengthen the monitoring of high-grade landslides.

3. Study Area and Datasets

3.1. Study Area

The Jinsha River constitutes the upper reaches of the Yangtze River, which flows through the four provinces of Qinghai, Sichuan, Tibet, and Yunnan. It originates from the Tanggula Mountains and joins the Yangtze River in Yibin, Sichuan. The terrain of the Jinsha River basin is undulating, and the drop of up as much as 3300 m provides considerable kinetic energy to the river water. Under the impact of the river and the uplift of the Qinghai-Tibet Plateau, a “V”-shaped structure of alpine valleys is formed. The Jinsha River basin is a high-incidence area of landslides, with many landslides, collapses, and blockages recorded throughout history.

The study area is located in the upper section of the Jinsha River basin, as shown in Figure 7. The latitude and longitude ranges are 97°52′E–100°46′E and 29°33′N–31°36′N, respectively, and the altitude span is 2300–6000 m, covering an area of 46,157 km². The study area comprises typical alpine and canyon terrain, with dense vegetation coverage, numerous snow-capped mountains, and densely developed rivers. The climatic zones of the study area are the temperate humid climate zone of the western Sichuan Plateau and the temperate semi-humid climate zone of the eastern Tibetan Plateau [36], and vertical climate change is obvious, depending on the altitude image. The study area experiences considerable rainfall in summer, with rainfall lasts for a long time due to the large area of the watershed. The flood season occurs from late June to mid-October, especially from July to September. The temperature difference between day and night in the study area is large, and the effect of freeze–thaw weathering is extremely strong. Coupled with river scouring, the stability of the rock mass in the basin is poor.

3.2. Datasets

The SAR data used in this paper comprise Sentinel-1 images. For the selected research area, we intended to used ascending and descending images for joint observation to obtain more landslide information. After actual processing, serious decorrelation of the descending images was observed, and more effective information could not be obtained. Therefore, we only used ascending images to demonstrate the proposed method.

Sentinel-1 images and the corresponding precise orbit determination ephemeris files were obtained from the European Space Agency (ESA) website. The data acquisition period was from September 2018 to February 2021, with 76 scenes in total. Details of used data and relevant processing parameters are presented in

Table 1. Data information and related processing parameters.

Sensor	Sentinel-1 A
Orbit Direction	Ascending
Path Frame	99–1280
Acquisition Time	September 2018–February 2021
Number of Images	76
Azimuth Angle	247.3°
Incidence Angle	34.34°
Ground Resolution	5 m × 20 m
Multi-looked Ratio	8:2
Temporal Baseline	25 day
Spatial Baseline	100 m
Number of Interferograms	130

, and the combinations of interferogram pairs are shown in Figure 8. The cropped DEM data of the research area were obtained from the official website of GMTSAR software, with a 1 s arc resolution corresponding to 30 m ground resolution.

4. Results and Analysis

4.1. Deformation Results

Using the SBAS-InSAR and stacking-InSAR methods, three deformation maps were obtained from the ascending data in the study area, as shown in Figure 9, Figure 10 and Figure 11. Combining the three results revealed two obvious landslide groups: G1 and G2; landslide group G1 comprises four landslides, including the Baige landslide, and landslide group G2 comprises several landslides, including the Xiongba ancient landslide. These results are consistent with previous studies [18,38,39]. The distribution histograms of the three deformation maps were generated (Figure 12), and the normal distribution curves were fitted according to the μ and σ of the three results. The fitted normal distribution curves of the three results are in agreement with the distribution histogram.

4.2. Landslides Identification and Gradation

Figure 13 shows the process of spatial cluster analysis and landslide screening on annual velocity map obtained by SBAS. Five areas presented with obvious deformation, as shown in Figure 13a, numbered as (I)–(V). The deformation regions were first extracted using statistical analysis, as shown in Figure 13b second column, with negative values in red and positive values in blue. Then, spatial cluster analysis was performed to extract landslide polygons from deformation regions, as shown in Figure 13b third column and fourth column. A total of 72,171 deformation points were identified in the deformation regions map. Via hotspot analysis, 37,674 deformation points with strong aggregation characteristics were extracted, of which 17,356 were negative deformation points (dark blue) and 20,318 were positive deformation points (dark red). Via point aggregation analysis, 595 landslide polygons with negative values were generated. Then, the proposed LSM was used to screen all landslide polygons, as shown in Figure 13b fifth column. The time-series average coherence map, DEME map, RMSE map, and slope map of the study area are shown in Figure 14. The landslide polygon statistics and model parameters are shown in

Table 2. Data statistics and threshold parameters.

Features	Research Area				Average within Polygons
	Min	Max	Mean	Std	Min	Max	Mean	Std	Threshold
γmean	0.07	0.97	0.23	0.12	0.08	0.4	0.13	0.04	>0.13
RMSE (mm/yr)	0	9.85	2.4	0.64	1.35	6.69	2.95	0.84	<2.95
DEME (m)	−48.4	90.6	11.6	10.1	−12.57	51.05	14.24	11.6	−11.6~11.6
Slope (°)	0	83	24	12.1	2.38	58.74	28.36	10.62	>25
Area (km2)	-	-	-	-	0.003	4.185	0.124	0.289	>0.1

Note: The first item represents the statistical data of each feature in the whole image range, and the second item represents the average statistical data of each feature in the landslide polygons.

. The landslide map includes 32 landslides obtained according to the steps outlined above.

Using the same method, 31 landslides were screened from 1011 landslide polygons on the cumulative displacement map obtained by SBAS, and 29 landslides were screened from 477 landslide polygons on the annual velocity map obtained by stacking. These landslides identified by three deformation maps were spatially fused to obtain the landslide map, which included 63 landslides with negative values. Then, an additional 32 landslides with positive values based on this method were included, and a complete landslide map containing 95 landslides was obtained. To reduce the repetitive workload, the optical image screening step that should have been performed before spatial fusion was performed after spatial fusion. Based on the macroscopic features of unstable landslides mentioned in Section 2.3, Google Earth optical images were used to confirm whether the screened landslides had corresponding macroscopic features. Ultimately, 47 landslides were identified, as shown in Figure 15.

These landslides were graded using the proposed LGM. After gradation, among the 47 landslides, 15 were identified as L-IV landslides, 10 as L-III landslides, 10 as L-II landslides, and 12 as L-I landslides. Among the 47 identified landslides, 13 landslides directly endangered villages. If a geological disaster occurs, the villages located below will be directly impacted. Therefore, it was also necessary to strengthen the monitoring of these 13 landslides.

4.3. Validation of Landslide Identification Results

To validate the landslide identification method proposed in this paper, two landslide groups (G1 and G2 in Figure 9) were chosen as the research objects to compare with the results obtained by other researchers. Figure 16a–d corresponds to the G2 landslide group, and Figure 16e,f corresponds to the G1 landslide group. The landslide boundaries identified in the present study were consistent with the landslide boundaries artificially identified by other researchers, and more landslide features can be observed: the pink line of landslide #19 in Figure 16 identifies the boundary reflected by the cumulative displacement map, which was significantly smaller than the blue and green boundaries, indicating that the landslide body was composed of a part of the continuous sliding body and a part of the seasonal sliding body; the pink and blue lines of landslides #17-1 and #17-2 in Figure 16d indicate that the landslide was composed of two independent sliding blocks; landslide #18 in Figure 16d and landslide #28 in Figure 16f indicate that stacking-InSAR can obtain more complete landslide boundary information than SBAS-InSAR, and the combined use of the two InSAR methods can solve the problem of omission or lack of landslide boundaries caused by decorrelation.

To validate the LGM proposed in this paper, feature points were selected among the 47 identified landslides, among which four feature points of landslide #17 were selected. Using the time-series cumulative displacement results obtained by SBAS-InSAR, the deformation on the time series of each feature point was extracted, as shown in Figure 17. The distribution of these 50 feature points is shown in Figure 18. According to the principle of least squares, the linear deformation velocity of each feature point was fitted, and the correlation coefficient (R²) between the fitted linear deformation velocity and the time series deformation was calculated. The R² under different grades was quantified, and the statistical results are shown in

Table 3. Relationship between landslide grade and correlation coefficient.

Grade	Num	Percent	R2
			Min	Max	Mean
L-I	12	25.53%	0.11	0.82	0.55
L-II	10	21.28%	0.18	0.78	0.59
L-III	10	21.28%	0.70	0.87	0.77
L-IV	15	31.91%	0.57	0.90	0.80

. With increased deformation grade, the R² of each landslide feature point increased significantly, indicating that the R² is obviously positively correlated with the landslide grade. Furthermore, the R² of the feature points can reflect the linear deformation trend of the landslide point. The larger the deformation degree and the stronger the linear deformation trend, the more urgent the need for landslide monitoring. The proposed LGM can provide a distinction basis for the monitoring of landslides in large areas.

5. Discussion

Frequent landslides in alpine canyon areas and vegetation coverage areas have attracted the research attention of scholars. Accurate identification of potential landslides and monitoring of high-risk landslides in these areas is an important task. Generally speaking, such areas are considerably affected by decorrelation (temporal decorrelation caused by vegetation coverage and spatial decorrelation caused by excessive deformation). In the present study, we combined SBAS-InSAR and stacking-InSAR to weaken the influence of decorrelation and avoid the omission of landslide identification by obtaining more data. In the process of landslide identification, accurate definition of the deformation region of landslides is very important. Based on the random characteristics of the phase noise and statistical analysis method, in the present study, we used the mean and three standard deviations of deformation results to distinguish between deformation regions and non-deformation regions to extract landslide deformation regions for subsequent landslide screening.

In the research presented in this paper, the spatial cluster analysis method was used to extract landslide polygons with clustering characteristics in deformation regions, an LSM based on multivariate features was proposed to screen the landslide polygons, and an LGM based on a signum function was proposed to grade the screened landslides. The proposed LSM uses the multivariate features of landslide polygons to evaluate the reliability of each landslide polygon as “true”, combined with optical image features to screen landslide polygons that are statistically excellent, effectively reducing false alarms of landslides. The proposed LGM uses the signum function to add results labels to each identified landslide and grade the landslides according to the degree of deformation, providing support for risk assessment of landslides in large areas and effectively avoiding waste of human and material resources. The landslide identification method proposed in this paper was automated and parameterized, which can minimize the identification error caused by human subjective factors.

Some of the landslide boundaries automatically identified in this paper were differed from real landslide boundaries. As shown in Figure 5 and Figure 15, landslide #17 comprised fusion of three landslides, the boundary of landslides #18 and #28 spanned the Jinsha River, and the boundary of landslide #33 was smaller than the real landslide boundary. A possible reason why the identified landslide boundary is larger than the real boundary range is that some noise points close to the landslide were classified into the landslide range during the spatial cluster analysis process. The possible reason that the identified landslide boundary is smaller than the real landslide boundary is the decorrelation caused by excessive deformation. To solve this problem, the next step was to introduce the idea of hydrological analysis to extract the ridge and valley lines of the study area, which were used to segment the landslide polygon, solving the problem of the landslide polygon crossing the terrain.

6. Conclusions

In this study, a method for landslide identification and gradation based on statistical analysis and spatial cluster analysis was proposed, which included a landslide screening model (LSM) based on multivariate features and a landslide gradation model (LGM) based on signum function. The proposed method was suitable for landslide hazard identification in low-coherence areas, such as alpine canyon areas and vegetation coverage areas, etc. It minimized the impact of phase noise, accurately identified landslides in large areas, and solved the problem of excessive false alarms and the considerable influence of decorrelation in the current InSAR landslide identification study. The proposed method can identify landslides with deformation characteristics and optical image characteristics in the study area and automatically extract the location, size, and boundary of potential landslides, which can provide support for geological disaster investigation departments.

Based on the landslide identification method proposed in this paper, 47 landslides were identified with obvious landslide characteristics in seven counties in the Jinsha River basin. These landslides were categorized as 15 L-IV landslides, 10 L-III landslides, 10 L-II landslides, and 12 L-I landslides using LGM. Geological disaster investigation agencies and local departments need to conduct real-time monitoring of 15 L-IV landslides and 13 landslides that endanger villages to prevent major geological disasters.

Our follow-up plan involves applying the proposed method to more areas in the southwestern mountainous region and to introduce the idea of hydrological analysis to solve the problem of landslide boundaries crossing terrain.