Landslide Detection Based on Multi-Direction Phase Gradient Stacking, with Application to Zhouqu, China

Abstract

Landslides are a common geological disaster, which cause many economic losses and casualties in the world each year. Drawing up a landslide list and monitoring their deformations is crucial to prevent landslide disasters. Interferometric synthetic aperture radar (InSAR) can obtain millimeter-level surface deformations and provide data support for landslide deformation monitoring. However, some landslides are difficult to detect due to the low-coherence caused by vegetation cover in mountainous areas and the difficulty of phase unwrapping caused by large landslide deformations. In this paper, a method based on multi-direction phase gradient stacking is proposed. It employs the differential interferograms of small baseline sets to directly obtain the abnormal region, thereby avoiding the problem where part of landslide cannot be detected due to a phase unwrapping error. In this study, the Sentinel-1 satellite ascending and descending data from 2018 to 2020 are used to detect landslides around Zhouqu County, China. A total of 26 active landslides were detected in ascending data and 32 active landslides in the descending data using the method in this paper, while the SBAS-InSAR detected 19 active landslides in the ascending data and 25 active landslides in the descending data. The method in this paper can successfully detect landslides in areas that are difficult for the SBAS-InSAR to detect. In addition, the proposed method does not require phase unwrapping, so a significant amount of data processing time can be saved.

1. Introduction

Landslides are a common geological disaster, usually occurring in steep mountainous areas due to the inclination of the slope and the soil’s own gravity, which may lead to downward sliding. Additionally, external factors such as earthquakes, rainfall, and man-made activities may lead to accelerated sliding of landslides [1,2,3]. Therefore, in order to avoid casualties and the economic losses caused by landslides, it is important to understand the locations of landslides and whether they will cause disasters. The most important aspect is the detection of landslide locations.

Remote sensing is widely used in landslide disaster detection because of its non-contact and large-area characteristics [4]. After the occurrence of landslides, the surface soil or rock is often exposed, resulting in spectral differences between the landslide area and the surrounding environment. Therefore, the landslide can be identified from optical images by the differences in spectral characteristics [3]. However, because some landslides do not expose soil or rock before the disaster, such landslide sites cannot be detected by optical images. SAR data, on the other hand, has certain advantages in this regard, and is widely used in the monitoring of various types of geologic hazards, such as earthquakes [5], volcanoes [6], floods [7,8], and so on. As a high-precision earth observation technique, InSAR has a unique advantage in the detection of this kind of slow-moving landslide because of its ability to obtain millimeter-level deformation values [9,10,11]. InSAR mainly obtains surface deformation using phase information to find abnormal areas on the surface. The main steps are as follows: Long-term SAR data and DEM are collected, and the average deformation rate of the surface is acquired by InSAR time series analysis techniques [12,13,14,15]. Then, the landslide in the region is comprehensively judged by combining related data (DEM, optical images). Finally, the temporal deformation characteristics of the single landslide are obtained and the concrete parameters of the landslide are calculated using multi-track data [16,17,18,19].

There have been numerous studies on obtaining large-scale surface information through InSAR to identify landslide disasters [20,21], but there are still some difficulties in interpreting the obtained deformation results in a large range. Therefore, some methods that can automatically extract the deformation region have been applied to the interpretation of deformation results. In 2017, Barra [22] proposed a method to automatically update the deformation region by using Sentinel-1 data and the PS-InSAR method, which has been applied to various deformation monitoring scenarios. This method can solve the difficulty in interpreting the final results. However, since the specific location of the landslide cannot be known during the landslide survey, it is necessary to process a large range of data simultaneously during data processing, which not only has low computational efficiency but also leads to the redundancy of results in most areas without deformation. With the development of SAR satellites, SAR data has increased explosively, and therefore this disadvantage will become more obvious. To address this, some scholars proposed that the deformation results of a large range be obtained by D-InSAR or Stacking methods first, and that key areas be processed by refined timing methods [23]. However, this method still has some defects. First, D-InSAR is the deformation information obtained between two time periods, and atmospheric errors cannot be avoided. Secondly, the D-InSAR and Stacking methods also require phase unwrapping, which cannot avoid the difficulty of unwrapping caused by steep mountain terrain. Some methods that do not need phase unwrapping have been gradually used to identify landslide hazards [24,25], but it is difficult to obtain enough monitoring points in low-coherence areas. Therefore, how to avoid phase unwrapping errors and how to detect landslides in both low coherence regions and large deformation regions are the main difficulties in landslide detection using InSAR.

The phase gradient method is gradually being used to locate deformation regions. It has been used to detect the deformation information associated with seismic fault zones [26,27] and faults [28]. The deformation range caused by these cases is relatively large, often extending for thousands or tens of kilometers in space, and these regions can be identified by the phase gradient method. Compared with large-scale deformation, these landslides are generally small, and most are located in a slope unit [29]. Some scholars have also made attempts at using the phase gradient for landslide recognition. First, for the identification of reservoir bank landslides, which tend to have obvious deformation before and after impounding, obvious signals are found in the differential interferogram, so the phase gradient can quickly locate them [30]. Secondly, earthquakes often lead to slope instability, and some slopes may move slowly or accelerate in earthquake areas, therefore such landslides are often dense [31]. Similarly, some scholars have used the NASA Uninhabited Aerial Vehicle SAR (UAVSAR) data and the phase gradient stacking method to explore the deformation boundary of monomer landslides [32]. Additionally, some scholars have applied the deep learning method to the phase gradient stacking results for the automatic identification of landslides [33,34].

The current use of interferogram phase gradient stacking for the detection of deformation regions is more likely to use the results of phase gradient stacking in the range and azimuth directions. However, this approach does not make full use of the neighboring information, as a central pixel has at least eight neighbouring pixels. Thus, some landslides cannot be detected by this method due to the fact that the phase gradient stacking results in the range and azimuth directions are not obvious because of the direction of the slopes and the imaging parameters of the SAR satellites. Consequently, to ascertain landslides based on phase gradient information while leveraging the entirety of temporal signal sequences, a novel approach is proposed: a multi-direction phase gradient stacking technique tailored for the identification of slow-moving landslides. This method combines short baseline sets of interferograms over the entire time series, attenuating atmospheric effects by means of averaging over a large number of interferograms, with fast access to slow-moving landslides through the results. Subsequently, this method is applied to Zhouqu County, Gansu Province, China, to validate its efficacy. Furthermore, we employed the Small Baseline Subset (SBAS) method to determine the average deformation rate within this area, enabling a comparative analysis between the proposed method and the SBAS method in landslide detection, elucidating their disparities and respective advantages.

2. Methods

2.1. Main Components of the Method

Initially, SAR data encompassing the designated study area along with an external Digital Elevation Model (DEM) are gathered, and the data set is processed through registration and difference. Subsequently, by employing the phase gradient calculation formula, the phase gradient is computed and a stacking process is executed to yield the phase gradient stacking outcome. This outcome integrates the phase gradient in each direction, culminating in the fusion of these directional phase gradients. The landslide detection methodology, delineated in three primary stages (as depicted in Figure 1), unfolds as follows:

(1)

Data Collection and Pre-processing: SAR imagery and the DEM are acquired, with subsequent geocoding of the primary image to simulate the topographic phase under SAR coordinates using the DEM. Interferometric pairs are selected based on predefined spatial and temporal baseline thresholds, following a construction approach akin to the Small Baseline Subset (SBAS), resulting in acquisition of differential interferograms composed of small baseline sets. Processing was performed in the GAMMA202202 software.

(2)

Multi-direction phase Gradient Stacking and Fusion: the phase gradient calculation method is applied to calculate the phase gradients in the four directions of the interferograms, and then the phase gradients in the four directions of the interferograms are stacked with the temporal baseline as the weighting factor. The noise is removed by means of spatial filtering, and the stacked results of the phase gradients are fused in the four directions, and the fused stacked results of the phase gradients are obtained.

(3)

Landslide detection: leveraging the fused multi-direction phase gradient stacking results, regions exhibiting outliers are identified. These outliers predominantly arise from deformation-induced anomalies, where the interference pattern stacking across the time series mitigates atmospheric influences. We consider these anomalous regions as target areas, especially regions that indicate landslides.

2.2. Principles Involved in the Method

2.2.1. Phase Gradient Stacking and Fusion

The phase gradient can be obtained without phase unwrapping, allowing us to acquire the true gradient value, as demonstrated by Cao [31], who used the phase gradient to calculate the deformation gradient. Moreover, since the phase gradients of many interferograms do not influence each other, they can compensate for the incoherence present in some regions caused by temporal and spatial baselines. By stacking the phase gradients, we consistently obtain effective results for the entire region, enabling us to observe regions in low coherence.

Here, the central difference method is used to calculate the phase gradient. The center pixel has eight neighboring pixels, so using the central difference method only needs to calculate the phase gradient in four directions (phase gradients in opposite directions are negative to each other). Taking the center pixel (i,j) as an example, the phase gradients in the four directions are calculated as follows:

$$\left\{\begin{array}{c}\begin{array}{c}∆{\mathsf{\phi }}_0(\mathrm{i},\mathrm{j})={\mathsf{\phi }}_{\mathrm{i},\mathrm{j}+1}-{\mathsf{\phi }}_{\mathrm{i},\mathrm{j}-1}\\ ∆{\mathsf{\phi }}_{45}(\mathrm{i},\mathrm{j})={\mathsf{\phi }}_{\mathrm{i}-1,\mathrm{j}+1}-{\mathsf{\phi }}_{\mathrm{i}+1,\mathrm{j}-1}\\ ∆{\mathsf{\phi }}_{90}(\mathrm{i},\mathrm{j})={\mathsf{\phi }}_{\mathrm{i}-1,\mathrm{j}}-{\mathsf{\phi }}_{\mathrm{i}+1,\mathrm{j}}\end{array}\\ ∆{\mathsf{\phi }}_{135}(\mathrm{i},\mathrm{j})={\mathsf{\phi }}_{\mathrm{i}-1,\mathrm{j}-1}-{\mathsf{\phi }}_{\mathrm{i}+1,\mathrm{j}+1}\end{array}\right.$$

(1)

∆φ_k represents the phase gradient in a particular direction of the central pixel. We define the range phase gradient as ∆φ₀, the upper right corner phase gradient as ∆φ₄₅, the orientation direction as ∆φ₉₀, and the upper left corner phase gradient as ∆φ₁₃₅. Since the wrapped differential interferogram has values in the range (−π, π], the phase gradient calculation results lie within the range (−2π, 2π]. However, results in the range (−2π, −π] and (π, 2π] occur due to phase jumps. To obtain the true gradient value within (−π, π], we wrapped the phase gradient as follows:

$$∆{\mathsf{\phi }}^{\prime }=\left\{\begin{array}{cc}∆\mathsf{\phi }+2\mathsf{\pi },& \hspace{1em}∆\mathsf{\phi }<-\mathsf{\pi }\\ ∆\mathsf{\phi },& \hspace{1em}-\mathsf{\pi }\le ∆\mathsf{\phi }\le \mathsf{\pi }\\ ∆\mathsf{\phi }-2\mathsf{\pi },& \hspace{1em}∆\mathsf{\phi }>\mathsf{\pi }\end{array}\right.$$

(2)

∆φ’ is the phase gradient of the difference interferogram. At each pixel point (i,j), phase gradients in four directions can be obtained: ∆φ0 (i,j), ∆φ45 (i,j), ∆φ90 (i,j), and ∆φ135 (i,j). Here, the phase gradient is calculated using differential interferograms in the SAR coordinate system, with the 0° direction indicating the ranging direction and the 90° direction indicating the azimuthal direction.

In the final phase gradient result, there are inevitable random phase gradients caused by random noise points and local micro gradients caused by the atmosphere, etc. Here, we adopt two strategies to solve this problem. First, we mask the low coherence region based on coherence, because the low coherence region is often noisy so it is difficult to obtain effective information; secondly, we stack and sum the phase gradients of interferograms constructed by short baseline sets, which can weaken some atmospheric effects on large spatial scales. However, due to the differences in temporal and spatial baselines, the gradient results obtained by the differential interferogram of different temporal baselines are also different. Therefore, the method for different temporal baselines is used here, and the temporal baselines are used as a weighting factor to stack the phase gradients of different temporal baselines. The final Stacking method is shown in Equation (3):

$${\mathrm{G}}_{\mathrm{k}}=\sum _{\mathrm{n}=1}^{\mathrm{M}}{\mathrm{T}}_{\mathrm{n}}∆{\mathsf{\phi }}_0^{\mathrm{n}}/\sum _{\mathrm{n}=1}^{\mathrm{M}}{\mathrm{T}}_{\mathrm{n}}^2$$

(3)

where G_k represents the phase gradient results stacked in each direction, k equals 0, 45, 90, 135, representing the phase gradient stacking results in four directions, M represents the number of differential interferograms we used for phase gradient stacking, and T represents the temporal baseline of the two-scene complex image of the differential interferogram.

Although some of the noise and the atmosphere are removed by masking the low coherence region and stacked summation, there are still some phase gradient anomaly results caused by random noise. Therefore, in order to better identify the region of interest, we use spatial filtering to remove these noises. Of course, due to the side-view imaging of SAR and the relief of the terrain, there will inevitably be some shadows and overlying areas. DEM and SAR imaging parameters are used to calculate areas of geometric distortion to help identify landslides and reduce landslide detection errors.

By stacking the phase gradients in the four directions, the phase gradient stacking results in the four directions can be obtained respectively. Due to the relationship between the morphology of landslide and satellite imaging, the landslide signals in the four results show different forms, and it may be impossible to distinguish the landslide signals from other signals in one result. Thus, the results of the four directions are fused through Equation (4):

$$\mathrm{G}=\left(\right|{\mathrm{G}}_0|+{|\mathrm{G}}_{45}|+{|\mathrm{G}}_{90}|+{|\mathrm{G}}_{135}|)/4$$

(4)

where G is the combination of the phase gradient stacking results in four directions, and the phase gradient stacking results were averaged in four directions to obtain a more complete landslide boundary signal, which is an improvement over the results obtained in a single direction.

2.2.2. Landslide Detection

Following the fusion of gradients to yield the final result, G, variations emerge due to disparities in the number of differential interferograms and temporal baselines. This discrepancy hinders the stabilization of G, prompting the necessity for normalization. Peristaltic landslides typically exhibit prolonged, subtle deformations, fostering the persistence of stable deformation signals across numerous differential interferograms. In contrast, noise remains stochastic and can be effectively eliminated through coherence analysis. Consequently, within the result G, regions depicting landslides manifest significantly heightened gradient values compared to other areas. In addition, there is a spatial aggregation tendency in these regions. These areas are defined as landslides.

3. Study Area and Used Data

3.1. Study Area

Zhouqu is located in the southern part of Gansu Province, China (104.371° E, 33.785° N), in the transition zone between the first and second steps, on the edge of the Tibetan Plateau. The Minshan mountain system runs from southeast to northwest, and the Bailong River is the main river that passes through Zhouqu County. The topography of this area is quite undulating, with an altitude of 1000–4000 m, a relative elevation difference of more than 2500 m, and a gradient concentrated between 40 and 70°. Figure 2b is the elevation map of the study area. Because the study area is located in the seismic zone, the regional tectonic movement is active, which makes the slope prone to instability. The climate of this region is characterized by a concentrated rainy season and frequent downpours. The monthly precipitation distribution in the study area from March 2018–September 2020, as shown in Figure 3, is at its peak from May to September each year. The unique geographical location and the influence of climate have led to a large number of landslides in this region, and landslide disasters have occurred in recent decades. The huge debris flow on 8 August 2010 caused huge casualties and property losses [35]. Meanwhile, in recent years, some ancient landslides have intensified and some small landslide disasters have occurred [36,37].

3.2. Used Dataset

Sentinel I imagery from orbit 55 (frame 107) and orbit 62 (frame 479) was collected from ESA between 22 March 2018 and 19 September 2020. A total of 71 images were collected in the ascending orbit and 74 images in the descending orbit.

Table 1. SAR data parameters.

Senor	Band	Track	Pass Direction	Pixel Spacing (m) in Rg × Az	Number of Images	Temporal Coverage
Sentinel-1A	C	55	Ascending	2.33 × 13.95	71	22 March 2018–19 September 2020
Sentinel-1A	C	62	Descending	2.33 × 13.95	74	22 March 2018–19 September 2020

shows the relevant parameters, and Figure 2a shows their spatial coverage. The majority of the images were obtained within 12 days, although some time points are still missing. To remove topographic phase from the interferograms, we used a DEM with a spatial resolution of 30 m provided by the Shuttle Radar Topography Mission (SRTM).

3.3. Data Processing

SAR images were processed to select the primary image from the ascending and descending data. The remaining images were then aligned with the primary image. To improve the image quality, the images were multi-viewed in both the range and azimuth directions (multi-view ratio of 8:2), with a ground resolution of about 30 m, which is close to the SRTM acquired. Secondly, to suppress the phenomenon of incoherence, we built a short baseline set, using 36 days and 100 m as the temporal and spatial baseline thresholds to select interference pairs for ascending and descending. The temporal and spatial baseline combination of the constructed interference pairs is shown in Figure 4. Based on the interference pair selected from the temporal and spatial baseline, we used the DEM with a resolution of 30 m to generate the simulated topographic phase, and used the difference to eliminate the influence of terrain to obtain the differential interferograms. Meanwhile, the Goldstein filtering method [39] was used to filter the differential interferograms to improve the quality of the interferograms. Finally, the filtered interferograms were used for the final phase gradient calculation and stack fusion. Here, 0.15 was used as the coherence threshold mask to remove the phase gradient calculated by the corresponding interference pair.

4. Result

4.1. Landslide Detection Based on Multi-Direction Phase Gradient Stacking

Utilizing the differential interferograms from the ascending and descending tracks, our method derived phase gradient stacking results for both ascending and descending data in the study area, illustrated in Figure 5. The null regions within these results denote areas affected by geometric distortion, varying between the ascending and descending tracks due to their distinct monitoring angles. While most regions exhibit low values, indicating no deformation, conspicuous red aggregation areas highlight anomalies starkly contrasting with their surroundings. A total of 26 landslides, labeled A1–A26, were identified in the results of the ascending data (Figure 5a), and a total of 32 landslides, labeled D1–D32, were identified in the results of the descending data (Figure 5b). Notably, A1–A19 and D1–D19 represent the same landslides detectable via both ascending and descending data. These landslides span significant areas and exhibit minimal geometric distortion in both sets of results, facilitating clear detection. Conversely, A20–A26 are landslides detectable solely by ascending data, while D20–D32 are detectable only by descending data. These differences may be due to smaller landslide areas or different geometric distortions that result in some landslides only being detectable through a single track. Of course, it is also possible that the data processing process is responsible for this discrepancy.

Most of the landslides detected on the ascending and descending tracks were located on both sides of the Bailong River. Within the orbit ascent and descent results, the most prominent regions displaying robust landslide detection signals include the Suoertou landslide, the Nanyu landslide group, and the Zhongpai landslide. Locations A1 and D1 in Figure 5 represent the Suoertou landslide, exhibiting consistent deformation regions in both ascending and descending track results. The Nanyu landslide group comprises numerous smaller landslides, such as A7 and D7 along the Xieliupo landslide, A9 and D9 at the Mentou Ping landslide, and D27 at the Jiangdingya landslide, all positioned along the Bailong River. These landslides exhibit evident signals in our results due to slope instability induced by persistent erosion at their bases from river water. Notably, significant disparities exist between the monitoring outcomes of the ascending and descending data in this region. While landslides 7–9 exhibit detectable signals in both the ascending and descending tracks, additional landslides, like D25–D29, are visible only in the descending track results. Geometric distortion primarily accounts for the absence of D28 and D29 in the lifting track data. Conversely, landslides D25–D27 evade detection due to compression within a few pixels in the orbital ascent data, a situation mirrored in the SBAS results [16,40].

A total of 39 landslides with deformation can be detected in this region by applying the ascending and descending data through the method in this paper, and their specific information is shown in

Table 2. Landslide Detection Information Statistics for the Ascending and Descending Tracks.

Landslide	Longitude	Latitude	Area (km2)	Elevation (m)	Type	Landslide	Longitude	Latitude	Area (km2)	Elevation (m)	Type
A/D1	104.330	33.801	3.144	1863	Earth-flow	A21	104.226	33.811	0.209	1625	Slide
A/D2	104.328	33.790	1.676	1937	Slide	A22	104.294	33.769	0.276	1742	Slide
A/D3	104.392	33.772	0.288	1597	Slide	A23	104.370	33.804	0.192	1647	Slide
A/D4	104.414	33.747	1.672	1768	Slide	A24	104.276	33.805	0.951	1908	Slide
A/D5	104.400	33.739	0.258	1522	Slide	A25	104.285	33.810	0.182	1930	Slide
A/D6	104.418	33.730	0.645	1743	Slide	A26	104.367	33.773	0.456	1626	Slide
A/D7	104.401	33.639	1.222	2179	Earth-flow	D20	104.390	33.781	0.342	1774	Slide
A/D8	104.523	33.725	0.391	1566	Slide	D21	104.408	33.740	0.104	1357	Slide
A/D9	104.516	33.742	0.204	1665	Slide	D22	104.450	33.717	0.0797	3108	Slide
A/D10	104.340	33.620	0.204	2043	Slide	D23	104.447	33.726	0.348	1676	Slide
A/D11	104.500	33.736	3.690	1682	Slide	D24	104.463	33.729	0.409	2136	Slide
A/D12	104.385	33.647	0.440	2196	Slide	D25	104.436	33.791	0.221	1564	Slide
A/D13	104.384	33.641	0.412	2164	Slide	D26	104.245	33.859	0.172	2356	Slide
A/D14	104.510	33.623	0.180	1608	Slide	D27	104.259	33.747	0.252	2096	Slide
A/D15	104.543	33.644	0.424	1286	Slide	D28	104.436	33.643	0.913	1876	Slide
A/D16	104.537	33.783	0.435	1425	Slide	D29	104.450	33.697	0.098	1720	Slide
A/D17	104.502	33.791	0.539	2096	Slide	D30	104.421	33.737	0.081	1984	Slide
A/D18	104.416	33.759	0.470	1789	Slide	D31	104.416	33.759	0.470	1838	Slide
A/D19	104.507	33.690	0.159	1718	Slide	D32	104.320	33.822	0.4523	2294	Slide
A20	104.240	33.795	0.965	1593	Slide

The criteria for classifying landslides are documented by Hungr [41].

. The two largest landslides, the A/D1 Suoertou landslide and A/D11 Zhongpai landslide, both have an area of more than 3 km², and it can be seen that these landslides are mainly distributed at an elevation of about 2000 m. The highest landslide is D22, which is located at 3000 m above sea level. The detected landslides are mainly shallow landslides, with some debris flow landslides but no landslides with sudden deformations such as rock falls, mainly because the InSAR technology detects creeping deformations, and landslides with deformations that are too large cannot be detected by the phase information [10].

4.2. Landslide Detection by SBAS

Deformation results for the survey area were derived using the traditional temporal InSAR method (SBAS). Figure 6 illustrates the SBAS results from both the ascending and descending tracks. Negative blue values represent movement away from the satellite, while positive red values indicate deformation towards the satellite. Notably, most of the outlined landslides in both the ascending and descending SBAS results exhibit deformation, while areas lacking outlines generally show no deformation. This demonstrates the efficacy of our method in accurately detecting regions experiencing deformation. In the ascending results, two regions stand out with substantial deformation: the Suoertou landslide (A1) and the Zhongpai landslide (A11), as depicted in Figure 5a. The region with the most significant deformation showcases a rate of approximately −90.97 mm/year, while the other region registers a rate of about −83.88 mm/year. Similarly, in the orbit descent results shown in Figure 5b, these same landslides exhibit noticeable deformation. Additionally, the D28 landslide in the orbit descent results also demonstrates considerable deformation, measuring approximately −57.14 mm/year.

To compare the results of the proposed method and SBAS, we conducted a comparative analysis of the detected landslides and summarized them in

Table 3. The phase gradient stacking method and the SBAS method landslide quantity detection statistics.

Method	Number of Landslides Detection
	Ascending	Descending
Multi-direction phase gradient stacking	26	32
SBAS	19	25

. The multi-direction phase gradient stacking method identified 26 landslides in the ascending data, which is more than the SBAS, which only identified 19 landslides in the ascending data. Similarly, for the descending data, the method detected 32 landslides, while the SBAS detected 25. These results show that the multi-direction phase gradient stacking method detects more landslides on both the ascending and descending orbits than the SBAS. The difference in detection results is mainly due to the limitations of SBAS in effectively capturing deformations in low-coherence regions that have been represented as real landslides in previous studies [42]. However, the method in this paper is excellent at detecting landslides in these challenging areas, and helpsto identify more landslides.

5. Discussion

5.1. Comparison of the Difference and Efficiency of Multi-Direction Phase Gradient Stacking and SBAS for Typical Landslide Detection

The Suoertou landslide is an earth-flow landslide which has been active since its reactivation in the 1970s, and the entire length of the landslide body is about 3000 m [43]. A1 in Figure 7a,c and D1 in Figure 7b,d represent the Suoertou landslide, which shows a wide upper part and a narrow middle and lower part. There are two main deformation zones in the whole slope, with a significant deformation in the upper and lower parts, while the middle part is more stable. In the SBAS results, the deformation rate of the slope is about −81.9 mm/yr, and according to the distribution of deformation points, the deformation is greatest in the upper part of the slope bordering with the middle part of the slope. This leads to the fact that not all of the deformation points can be detected in this area, but using this paper’s method, the area is detected. This limitation is most evident in the A2/D2 landslide, where a clear signal exists throughout the landslide, but effective deformation points are difficult to detect in the SBAS results. This situation cannot be avoided during the phase deconvolution of conventional time-series InSAR methods, where deformation is too large to cause incoherence [10,44], and therefore the only way to remove low-coherence regions is to mask them during phase deconvolution. However, the method in this paper applies differential interferograms for phase gradient calculation and stacking, so these low-coherence regions may retain signals in some interferograms, and all are able to be detected in a sufficiently large number of regions.

Figure 8 depicts another significant landslide accumulation area, namely the Nanyu landslide group, known for its historical presence. While the overall slope has maintained relative stability in recent decades, localized small landslides, such as the Mentou Ping landslide and the Jiangdingya landslide [37], have emerged within certain sections. It is important to note that the Xieliupo landslide is a large landslide located north of the Nanyu landslide on an active fault line [40]. Continuous erosion at the base of the Bailong River sustains its active state. Studies indicate significant deformation within this landslide slope, posing challenges in effectively monitoring adequate deformation points [45,46], evident in the SBAS results shown in Figure 8c,d. The SBAS method primarily captures signals at the base of the slope, but struggles to detect effective deformation points in the upper regions. Conversely, the multi-direction phase gradient stacking method, showcased in the ascending and descending results in Figure 8a,b, reveals activity across the entire slope. This contrast highlights the challenge faced by the SBAS in capturing comprehensive deformation signals along the entire slope.

In other small landslides, such as A8 in Figure 8a,c, the landslide boundaries detected by this paper’s method are in high agreement with the areas with deformations obtained by the SBAS. This is also true of D27, D28, and D29 in Figure 8b,d. Comparing the results of the method in this paper with the SBAS method, it is easy to find that the results of phase gradient stacking are positively correlated with the magnitude of deformation in the detectable region. Landslides with large deformations, such as D28, also show relatively high values in the results of phase gradient stacking. In contrast, landslides with smaller deformations, such as D27, also show relatively low values in the results of phase gradient stacking. This is mainly due to the fact that landslides with large morphology variables also have more drastic phase changes in the interferograms, so the phase gradients are relatively large. This is why the A7/D7 Xieliupo landslides have the highest values of stacked phase gradients.

This method quickly extracts the landslides from the differential interferograms, providing a reference for quickly locating the landslide with deformation and improving the efficiency of large-scale landslide surveys. Since the method is based on wrapped differential interferograms, phase unwrapping is not required during data processing. The elapsed time of the multi-direction phase gradient stacking method and the SBAS method are listed in

Table 4. Comparison of efficiency between the multi-direction phase gradient stacking method and the SBAS method.

Method	Total Time (Seconds)
	Ascending	Descending
Multi-direction phase gradient stacking	404	375
SBAS	3158	3427

. There are a total of 109 interferograms in the ascending data and 107 interferograms in the descending data. The size of the study area is 1507 × 1045 pixels, about 927.19 km², and the geometric distortion area is removed. The effective gradient points of the ascending and descending tracks are about 1,325,600 and 1,388,900, respectively. Approximately 297,400 effective deformation points were calculated for the ascending track and 429,500 effective deformation points for the descending track. As shown in Table 4, the time applied by the multi-direction phase gradient stacking method in the two data sets is much lower than the time of the SBAS, which is consistent with our expectations and can save a lot of time on the basis of obtaining effective landslide results, thereby improving landslide identification efficiency.

5.2. Comparison of Range and Azimuth Phase Gradient Stacking

The application of the range and azimuth direction to seismic fault zones and active faults has been investigated [26,27,28], but the detection of landslides using only the range or azimuthal direction is subject to partial erroneous detection of landslides [34]. This is mainly due to the fact that gradients mainly detect boundary information, so landslides with deformation appear as two neighboring regions in the detection results of a single direction phase gradient stacking [34], which leads to a situation where some landslides are not signaled in one direction but are signaled in other directions. For example, the two small landslides in Box 1 of Figure 9a–e cannot be detected in G₀, while the landslide boundaries can be clearly detected in other directions. As for the landslide in Box 2 in Figure 9, the landslide boundary information is detected in the results of all four directions (Figure 9a–d), while in the fused results, the results of all four directions are combined. Since there is always a difference between the landslide direction and the irradiation direction of the radar imaging [10], the landslide morphology detected by applying phase gradient stacking in different directions is obviously different, such as the landslide in Box 3 in Figure 9. Although the landslide information can be detected in the results of all four directions, the landslide is distributed along the range towards an angle of about 45°, so it shows more prominence in G₄₅. Secondly, the landslide boundary information detected in G₁₃₅ can better complement the results of G₄₅, so the fusion of all the results can obtain more complete landslide boundary information. Overall, the fused results can combine information from all directions, and the fusion of the four methods can attenuate some of the noise signals and better highlight the landslide signals.

5.3. The Effect of Interference on the Selection of the Result

This method involves using short baseline sets, where different temporal baseline thresholds can affect the results obtained. We evaluated the temporal baseline thresholds at 12, 24, 36, 48, and 60 days, leveraging the stable 12-day revisit period of Sentinel-1A. Figure 9 shows that as the temporal baseline threshold increases, subtle signals indicating small landslides gradually become more apparent. In Figure 10, two small landslides within the dotted circle were not detected using a 12-day temporal baseline threshold. However, faint signals appeared at the 24-day threshold. At 36 days, distinct signals emerged, clearly distinguishable from the non-landslide areas. These signals became more pronounced at the 48- and 60-day thresholds. The capability of this method to detect landslides relies on deformation signals within the differential interferogram. Some smaller landslides exhibit weak deformations, which are imperceptible in the 12-day differential interferogram, leading to undetected landslides. Longer time differentials showcase more evident deformation signals, enabling this method to detect smaller landslides. Interestingly, as the temporal baseline threshold increases, the signal strength of larger landslide points gradually diminishes. Hence, in selecting a temporal baseline, using a 36-day threshold for constructing short baseline sets proves effective. This threshold allows the detection of smaller landslide signals while retaining most signals from larger landslides.

It is evident that the phase gradient stacking method is closely related to whether there are deformation signals in the temporal baseline and interferogram. According to this feature, we can select data from a fixed period of time, such as the Sentinel-1A data of 6 months or one year, and initially screen out the regions where deformations may exist over a large range. However, this approach cannot directly obtain temporal deformation information to understand the motion state. Ultimately, time series InSAR technology is needed to obtain temporal deformation of specific landslide points. This approach focuses on regions with deformation, eliminating a large number of regions without deformation, and thereby improving the efficiency with which we process and analyze data.

6. Conclusions

In this paper, a multi-direction phase gradient stacking method based on differential interferograms is proposed for landslide detection, and landslides in Zhouqu County and its surrounding area are successfully detected. The method is based on the differential interferogram of a small baseline set, and the abnormal region is obtained by multi-direction phase gradient stacking. The landslides detected by this method are compared with that detected by SBAS. It is found that the landslide detected by this method is highly consistent with that detected by SBAS, indicating that landslides can be successfully detected by this method. At the same time, landslides can also be detected by this method in the incoherent region. More importantly, this method has great advantages in dealing with large deformation landslides, and can obtain complete information of the entire landslide. Additionally, the data processing efficiency of this method is much higher than that of the SBAS method, and the landslide location can be obtained efficiently in large-scale landslide detection. However, it should be noted that the method in this paper does not obtain the deformation variable of the landslide, and can only indicate whether the landslide is slippery or not, so the deformation characteristics of the landslide cannot be analyzed. The research results also show that the selection of the temporal baseline threshold will greatly affect the results of landslide detection. If the temporal baseline threshold is too small, some landslides with weak deformations cannot be detected, whereas if the temporal baseline threshold is too large, some landslide signals with large deformations will be weakened. For example, the Sentinel-1 data, with a 36-day temporal baseline threshold, proved to be the optimal temporal baseline threshold. However, whether a 36-day temporal baseline threshold is optimal for all regions is not clear and will be addressed in subsequent studies.

The continuous development of SAR satellites will lead to the continuous surge of SAR data. Traditional temporal InSAR technology needs to re-process a large range of SAR data for newly acquired data in order to determine whether there are new geological hazard points based on deformation. Compared with the traditional temporal InSAR which needs to be untangled, this method can greatly improve the efficiency of identification. Therefore, in a large-scale landslide survey, this method can be used to quickly extract the landslide aggregation area, and the scale and aggregation degree of the landslide can be known, providing a guideline for landslide hazard evaluators to assess the hazards of landslides. Sentinel-1 data has a stable and short revisit period, which provides great convenience for our method, which can collect a large amount of data in a short time for the acquisition of results. However, the use of long-period data in this method is unknown. Follow-up research should utilize deep learning to realize large-scale automatic landslide detection.