Next Article in Journal
Microwave Brightness Temperature (MBT) Background in Bayan Har Block, Qinghai-Tibet Plateau and Its Importance in Searching for Seismic MBT Anomalies
Next Article in Special Issue
InSAR Modeling and Deformation Prediction for Salt Solution Mining Using a Novel CT-PIM Function
Previous Article in Journal
Semantic Segmentation of Very-High-Resolution Remote Sensing Images via Deep Multi-Feature Learning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Monitoring and Stability Analysis of the Deformation in the Woda Landslide Area in Tibet, China by the DS-InSAR Method

1
School of Land Science and Technology, China University of Geosciences, Beijing 100083, China
2
Zhejiang Huadong Mapping and Engineering Safety Technology Co., Ltd., Hangzhou 310005, China
3
The Department of Systems Design Engineering, University of Waterloo, Waterloo, ON N2L 3G1, Canada
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(3), 532; https://doi.org/10.3390/rs14030532
Submission received: 29 December 2021 / Revised: 20 January 2022 / Accepted: 21 January 2022 / Published: 23 January 2022

Abstract

:
The Woda area in the upper Jinsha River has steep terrain and broken structures, causing landslide disasters frequently. Here, we used the distributed scatterer interferometric SAR (DS-InSAR) method to monitor and analyze the Woda landslide area. With the DS-InSAR method, we derived the deformation of the Woda landslide area from 106 Sentinel-1A ascending images acquired between 5 November 2014 and 4 September 2019 and 102 Sentinel-1A descending images acquired between 31 October 2014 and 11 September 2019. The obvious advantage of the DS-InSAR method compared to the persistent scatterer (PS) InSAR (PS-InSAR) method is that the densities of the monitoring points were increased by 25.1% and 22.9% in the ascending and descending images, respectively. The two-dimensional deformation of the landslide area shows that the maximum surface deformation rate in the normal direction was −80 mm/yr, and in the east–west direction, 118 mm/yr. According to the rescaled range (R/S) analysis, the Hurst index values of the deformation trends were all greater than 0.5, which means the deformation trend will continue for some time. In addition, we analyzed the influencing factors and the deformation mechanism of the Woda landslide area and found that the surface deformation is closely related to the geological structure and precipitation, among which precipitation is the main factor triggering the deformation. Our monitoring results will help the local government to conduct regular inspections and strengthen landslide disaster prevention in low-coherence mountainous areas.

1. Introduction

Vigorous tectonic movements and natural erosions shaped the complex landform of Southwest China, where mountains, plateaus, and vertical and horizontal surface gullies are widely distributed. The Qinghai–Tibet Plateau, known as the ‘roof of the world’, has obvious crustal uplift and strong and rapid rivers, resulting in frequent geological disasters on both sides of the rivers [1,2]. Landslide disasters are easily triggered in areas with lush natural or planted forests or in steep areas, especially after extreme events. In eastern Tibet, continuous landslide disasters have seriously affected social and economic development and people’s lives [3]. The region along the Jinsha River and its tributaries has steep terrains and broken structures, so geological disasters such as ground fissures and landslides have occurred frequently. Since the 1980s, 61 landslides have been recorded [4]. The Woda landslide area in Yanbi Township, Jomda County, eastern Tibet, is located in the upper reaches of the Jinsha River. It mainly includes two landslide bodies, which started slowly deforming in 1985. Remote sensing interpretation and ground survey both show risks of further deformation and sliding that may cause river closure [5]. On 11 October 2018, a major landslide occurred in Baige village in the lower reaches of the Jinsha River, about 80 km from the Woda village [6]. The slope blocked the Jinsha River and damaged many roads, bridges, and buildings. Such major geological disasters may impact the areas hundreds of kilometers upstream and downstream. Therefore, monitoring the long-term surface deformation of the Woda landslide area is of great significance for preventing and controlling the occurrence of geological disasters.
In the past few decades, leveling [7], three-dimensional laser scanning [8], global navigation satellite system (GNSS) [9], and unmanned aerial vehicle (UAV) photogrammetry [10] have been widely used in deformation monitoring. However, they cannot provide satisfactory results in the identification and monitoring of landslides in large areas due to their long measurement cycles, small spatial coverage, large demand of human and material resources, and high-risk factors for close-range work. Synthetic aperture radar (SAR) remote sensing has the advantages of wide coverage as well as all-day and all-weather observation, so it has been widely used in geological disaster monitoring [11,12]. With the deepening of interferometric SAR (InSAR) research, differential InSAR (D-InSAR) uses SAR images before and after the deformation to monitor the surface deformation [13,14], but it has low monitoring accuracy due to time and geometric decoherence and atmospheric delay, meaning it cannot meet the long-term monitoring requirements. In 2000, Ferretti et al. proposed the theory and method of persistent scatterer InSAR (PS-InSAR) [15], suitable for monitoring targets with stable radar scattering characteristics, such as buildings and bridges in urban areas, but it is difficult to effectively apply it in areas such as farmland and woodland. For this reason, subsequently, Berardino et al. contributed many theoretical innovations and algorithm improvements to the PS-InSAR method and proposed some new time series InSAR methods, which greatly enriched the time series InSAR technical theory system [16,17,18,19]. In recent years, the time series InSAR method has been widely used in landslide investigation and identification [20], landslide sliding evaluation [21], the establishment of early warning systems for landslides [22], and landslide surface deformation monitoring [23,24], which is powerful technical support for landslide disaster research.
In 2011, Ferretti et al. proposed the second-generation permanent scatterer technology (SqueeSAR) and related scholars around the world began to shift their research focus to distributed targets with relatively weak backscatter [25]. In recent years, the distributed scatterer InSAR (DS-InSAR) method based on SqueeSAR technology was introduced to the research in low-coherent areas with weak backscattering. Goel et al. proposed the Anderson–Darling (A–D) test method to select statistical homogeneous pixels (SHPs), which effectively solved the high error rate of the K–S test in small sample data [26]. Gianfranco et al. decomposed the covariance matrix of DS points and selected the eigenvector with the largest value as the optimal phase of the DS point [27]. Jiang et al. proposed fast SHPs, modified fast SHPs, and hypothesis test of confidence interval (HTCI) algorithms, which effectively improved the SHP selection efficiency and phase optimization quality and greatly promoted the progress and development of the DS-InSAR method [28,29]. This method combines PS points and DS points to significantly increase the distribution density of monitoring points. Compared to other time series InSAR methods, it has obvious advantages and is more suitable for the surface deformation monitoring of landslides in mountainous areas [30,31].
In addition, time series InSAR technology can only extract the deformation rate and the deformation time series along the line of sight (LOS) direction [32]. In fact, real space deformation is three-dimensional, namely in the vertical, east–west, and north–south directions. The vertical direction is usually longer than the other two directions. Since the flight direction of the Sentinel-1A satellite is close to the north–south direction, the satellite heading angle is relatively small, and it is insensitive when calculating the north–south deformation [33]. However, it is not enough to only know the deformation information in the LOS direction, and it is difficult to meet the needs of surface deformation research. Generally, for research on surface deformation, it is best to perform a comprehensive analysis with ascending and descending SAR data and calculate the vertical and east–west two-dimensional deformation, further enriching and expanding the research work of the surface deformation [34,35]. In landslide areas where data sources are scarce and it is difficult to obtain leveling, GNSS, and other data, the two-dimensional monitoring results combining ascending and descending SAR data can help distinguish different landslide deformation directions, reflect different deformation rate modes, and reduce the influence of geometric distortion with high reliability and practicability [36,37].
In view of the fact that the time series and two-dimensional studies of the Woda landslide area are relatively few, in this study, the DS-InSAR method was applied to monitor and analyze the deformation in the Woda landslide area using the ascending and descending Sentinel-1A images. The HTCI algorithm with the highest computational efficiency was used to select SHPs, and the adaptive spatial nonlocal filtering method was constructed based on the results of the SHP selection to perform phase optimization. We analyzed the deformation characteristics and inducing factors of the landslide area, which can provide reference and guidance for the early investigation, monitoring, and timely warning of densely vegetated landslides.

2. Study Area and Datasets

2.1. Study Area

The study area is the Woda village and its surrounding areas in the upper reaches of the Jinsha River. The Woda village is located in Yanbi Town, Changdu City, Tibet, China, with a longitude of 98 4 9 38.6 3 E and latitude of 31 2 5 56.3 1 N (Figure 1a). There are two closely distributed landslides in the area, named landslide bodies 1 and 2, along the Jinsha River (Figure 1b). The front edges of the landslides are steep and reach the Jinsha River valley at the toe, but the middle and rear parts are relatively flat (Figure 1c). The elevation of the Jinsha River surface is about 2950 m, and that of the trailing edge apex of the landslide area is about 4000 m. Such a great elevation difference may contribute to geological disasters, such as large-scale landslides, for example, the Baige landslide, which blocked the main stream of the Jinsha River to form a dammed lake, causing huge property losses.

2.2. Datasets

We used the Sentinel-1A data in the ascending and descending orbits. Sentinel-1A is a C-band synthetic aperture radar (SAR) satellite with a new TOPSAR imaging mode. It has a 12-day return cycle and a large-scale spatial coverage of 250 km × 250 km at medium resolution. At present, the latest global real-time observation data are free of charge. A total of 208 scenes of Sentinel-1A interferometric wide (IW)-swath mode data were adopted in this study. Detailed information is shown in Table 1.
In addition, the SRTM1 30 m DEM with a 1 arc-second spatial resolution from National Aeronautics and Space Administration (NASA) was used to remove topographic phase and geocoding. A Google optical image was also used to show the final results and the analysis of the landslide area.

3. Methods

In this study, we used the DS-InSAR method to monitor the Woda landslide area, and the specific data processing flowchart is shown in Figure 2. Firstly, SAR images were registered and clipped by using precise orbit and DEM data, and then differential interferogram pair sequences were generated. When setting the temporal and space baseline thresholds to 72 days and 200 m, the 106-scene ascending SAR images generated 209 interference pairs; when setting the temporal and space baseline thresholds to 96 days and 200 m, the 102-scene descending SAR images generated 201 interference pairs. The generated interference pair information is shown in Figure 3. Then, we selected the PS and DS points and extracted their phase and elevation information to calculate the deformation. Finally, we combined the precipitation data and geological map to analyze the deformation of the landslide area in detail. The precipitation data in this paper can be downloaded from the official Soil and Water Assessment Tool (SWAT) website at https://swat.tamu.edu/ (11 October 2021), and the geological map is from the hydrogeological map of China (1:200,000).

3.1. PS Selection

The amplitude threshold and the amplitude dispersion index threshold were jointly used to select the PS point. The amplitude threshold was used to select pixels with high amplitude values as candidate points, from which those with smaller amplitude deviation values were determined as the final PS points. The pixels with high amplitude values indicate better coherence, and a smaller amplitude deviation index indicates better stability. Therefore, the selected PS points have stable radar scattering characteristics, which can improve the accuracy of the results.

3.2. DS Selection

DS point selection includes SHP selection and phase optimization. SHPs were selected according to the SAR image intensity information, which was used to optimize the interference phase of the SAR image. The phase optimization was evaluated by calculating the temporal coherence and setting the temporal coherence threshold to 0.4 to select high-quality DS points [12,38].

3.2.1. SHP Selection

There are dozens of SHP selection algorithms with different accuracies. Assuming that the SAR intensity dataset obeys the exponential distribution and the amplitude dataset obeys the Rayleigh distribution, the algorithm can be divided into two categories—non-parametric hypothesis testing methods, such as the KS test [39] and the BWS test [40], and the parameter test algorithm, such as HTCI. The HTCI algorithm has higher selection efficiency and accuracy than the KS and BWS test methods. In this study, we used the HTCI algorithm to select SHPs.
Assuming that N SAR images obey complex Gaussian distribution in the time dimension and the variance of random variables is σ 2 / 2 , the SAR image intensity sequence obeys Rayleigh distribution.
f ( I ) = 1 σ 2 e I σ 2 , I 0
According to the central limit theorem, the sample mean obeys a normal distribution of N ( σ 2 , σ 4 / ( N L ) ) . The interval estimation analytical equation is as follows:
P { σ 2 z α / 2 σ 2 / N L < σ ¯ < σ 2 + z α / 2 σ 2 / N L } = 1 α
where Z α / 2 is the percentage of the standard normal probability density function α / 2 . Defining the center pixel in the window as the reference pixel and the remaining pixels as the neighboring pixels, if the sample mean of the neighboring pixels is within the given confidence interval, the neighboring pixels are considered the SHPs of the reference pixel. However, in actual data processing, the available SAR images are very few, and the Gaussian hypothesis is often difficult to establish. The Gamma distribution of the accumulated N scene intensity images is N σ G ( N a , b ) , where a = L , b = σ 2 / L (where L denotes the multi-look number). Set Y = N σ / b , so Y obeys the standard Gamma distribution with the shape parameter N a , which is Y G ( N a , 1 ) [41]. Let g α ; N a be the α percentage of the standard Gamma distribution of the known parameter N a , and we get the following:
P { g α / 2 ; N a < Y < g 1 α / 2 ; N a } = 1 α
From which, we obtain the fixed interval analytic equation as
P { g α / 2 ; N L σ 2 / ( N L ) < σ ¯ < g 1 α / 2 ; N L σ 2 / ( N L ) } = 1 α
If there is less available time series SAR image data, Equation (4) is better than
Equation (2), as the interval value is asymmetric and the tail feature is emphasized. In summary, the algorithm defines a fixed window for each pixel of the image. Then, the mean temporal intensity of the reference pixel and that of the neighborhood pixel are taken as the values of σ 2 and σ ¯ in Equation (4), respectively. Finally, the neighborhood pixels located in the confidence interval are selected as the SHPs of the reference pixels. The algorithm selects SHPs according to the intensity mean of the local spatial pixels and the local minimum mean square error, which effectively improves its robustness and selection accuracy.

3.2.2. Optimal Phase Estimation

In the selection window, the interference phases of SHPs are affected similarly by noise, are highly correlated, and have distributed characteristics. In this study, we used the adaptive spatial nonlocal filtering method based on the results of SHP selection to perform phase optimization. In a certain window, the reference pixels were filtered by adaptive spatial nonlocal filtering using the selected SHPs to reduce the phase noise of the interferogram.
At any pixel p in the interferogram, the estimated value of the adaptive spatial nonlocal filtering is as follows [42]:
N L ( p ) = m Ω w ( p , m ) u ( m )
N L ( p ) is the phase value of the filtered pixel p , and is the phase value of the homogeneous pixel m in the search range Ω . The weight value w ( p , m ) depends on the similarity between pixels p and m , and satisfies the following conditions: 0 w ( p , m ) 1 , m w ( p , m ) = 1 .
d ( N p , N m ) = u ( N p ) u ( N m ) 2 , a 2 1 + u ( N p ) u ( N m ) 2 , a 2
Z ( p ) = m exp ( d ( N p , N m ) h 2 )
w ( p , m ) = exp ( d ( N p , N m ) h 2 ) / Z ( p )
where d ( N p , N m ) is the Gaussian weighted Euclidean distance, which is defined as the similarity between the distance weighted windows; Z ( p ) is the normalization parameter; h is the attenuation factor of the control exponential function; a is the standard deviation of the Gaussian kernel function. The larger the d ( N p , N m ) value is, the greater the weight. The value of h is usually determined according to the noise level. A value that is too large may filter out some useful features and textures, but a value too small makes it difficult to filter out noises. Therefore, selecting the value of h becomes a key factor that affects the nonlocal filtering. In this paper, it is defined as the following:
h = 1 1 + D A
where D A is the amplitude deviation index of the estimation window. In areas with large terrain undulations, a small h will be selected to retain as much texture information as possible; in areas with relatively flat terrain, a large h will be used to remove noise.
After phase optimization, the quality of the optimized phase was evaluated by Equation (10). The difference between the optimized interference phase and the original interference phase was fitted to obtain a fitting metric value of γ D S [25].
γ D S = 2 N 2 N Re n = 1 N k = n + 1 N e i ϕ n k e i ( ϑ n ϑ k )
ϑ n and ϑ k are the optimized interference phases of a single master image, and ϕ n k is the original interference phase. γ D S can also be regarded as the temporal coherence of the interference pair sequence and was used to describe the quality of the optimized phase. By setting the temporal coherence threshold, the points with high temporal coherence were selected as the DS points in the subsequent analysis.

3.3. Combined Data Processing

Combining the selected PS and DS points can effectively improve the insufficient point selection density of using the PS-InSAR method alone. The phases of the combined coherent points are differentiated and unwrapped to construct triangulation networks. Through iterative calculation and regression analysis, we separated the linear deformation, the elevation residual, and the residual phase. According to the different characteristics of the residual phase in the time and space domains, multiple filtering in the frequency domain can separate nonlinear deformation, atmospheric delay phase, and noise phase. The singular value decomposition was used to obtain the deformation rate and accumulative deformation during the study period.
Then, using the deformation obtained from the ascending and descending orbit data, we calculated the vertical and east–west deformation characteristics in the landslide area. Finally, the stability analysis of the study area was carried out according to the two-dimensional deformation results. The details are as follows.

3.3.1. Two-Dimensional Deformation Calculation

According to the geometric relationship of radar imaging, the surface deformation monitored by the ascending and descending SAR images were both in the LOS direction ( D L O S ). D L O S can be decomposed into the deformation in the vertical ( D v ), east–west ( D e ), and north–south ( D n ) directions (Figure 4).
D L O S can be expressed as [43,44] the following:
D L O S = [ cos θ sin θ cos α sin θ sin α ] [ D v D e D n ]
where θ represents the radar incident angle, and α represents the satellite heading angle. Since the flight direction of the Sentinel-1A satellite is close to the north–south direction, the deformation monitoring results in this direction are not sensitive. Therefore, we calculated the deformation in the vertical and east–west directions.
[ D L O S 1 D L O S 2 ] = [ cos θ 1 sin θ 1 cos α 1 cos θ 2 sin θ 2 cos α 2 ] [ D v D e ]
A v , e = [ cos θ 1 sin θ 1 cos α 1 cos θ 2 sin θ 2 cos α 2 ]
where A v , e is the coefficient matrix. The radar incidence angle information of each pixel of the ascending and descending SAR images is shown in Figure 5. In areas with large terrain undulations, the incidence angles of neighboring pixels are greatly different, and even the incidence angles of the same pixel in the ascending and descending orbits are significantly different. In this case, the deformation results are invalid in the overlapped and shadow areas, and the reliability in the perspective shrinkage area is reduced [45]. Considering the actual situation, we used the incidence angle of each pixel to calculate the two-dimensional deformation by the least square method as follows.
D v , e = ( A v , e T P A v , e ) 1 A v , e T P D L O S 1 , 2
where P is the weight matrix, the reciprocal of the residual standard deviation of the monitoring result is used as the weight, and D v , e = [ D v   D e ] T ; D L O S 1 , 2 = [ D L O S 1   D L O S 2 ] T . Finally, we obtained the deformation in the vertical direction and east–west direction of the monitoring point of its tangent plane. However, it should be noted that since the radar incident angle information of each pixel was considered, the vertical direction deformation obtained was actually the deformation of the surface normal direction.

3.3.2. Stability Calculation

We analyzed the stability of the Woda landslide area by the rescaled range (R/S) analysis method [46,47]. The R/S analysis method is a time series fractal statistical method based on fractal theory. It estimates the deformation trend over time and predicts the continuity of future deformation by calculating the Hurst index value H of a non-linear time series. A long time series { x N } consisting of N scenes of SAR images can be divided into m sub-intervals with the length L . For each sub-interval, set
X t , m = u = 1 L ( x u x ¯ u )
where x ¯ u is the average value of the n-th ( n m ) sub-interval x u , and X t , n is the cumulative deviation of the n-th sub-interval. Then,
R = max ( X t , n ) min ( X t , n )
S = 1 L u = 1 L ( x u x ¯ u ) 2
where R is the range of the interval, and S is the standard deviation of the interval. R/S is
R / S = k ( n ) H
where k is a constant with the empirical value of 0.5; H is the Hurst exponent. Taking the logarithm of both sides of the above equation, we get
log ( R / S ) n = H log ( n ) + log ( k )
From the above equation, the Hurst index value H can be estimated by the least square regression. The calculated Hurst index value ranges between 0 and 1. Only when the value is greater than 0.5, the monitoring time series deformation has a long-term correlation, and the current deformation trend will be maintained in the future. The closer the Hurst index value is to 1, the more reliable and stable the correlation and consistency are.

4. Results and Analysis

4.1. Monitoring Results and Analysis

From the 106 ascending and 102 descending SAR images, we derived the surface deformation of the Woda landslide area using the above method. On the basis of the SHPs selected by the HTCI algorithm, we optimized the interferograms by the adaptive nonlocal filtering method. The time coherence threshold was defined as 0.6 to obtain the DS point; the amplitude dispersion index threshold was set as 0.4 to obtain the PS point. Both the PS point and the DS point were adopted for the phase unwrapping and residual phase removal, and finally, the surface deformation during the study period was obtained.
Figure 6 shows the LOS deformation rate of the study area obtained by different methods. The Woda landslide area is mainly mountainous and lacks strong scatterer point targets, such as buildings and bare rock. So, the point targets obtained by the traditional PS-InSAR method are sparse (Figure 6b,d) and cannot be used for analysis. The results of DS-InSAR greatly improve the spatial distribution density of point targets than that of PS-InSAR, especially in the areas with weaker coherence. Specifically, DS-InSAR selected a total of 162,992 point targets from the ascending data and 149,518 point targets from the descending data, with point densities of 25.5% and 23.4%, respectively. The PS-InSAR method selected a total of 2293 point targets from the ascending data and 3088 point targets from the descending data, with densities of 0.4% and 0.5%, respectively. In addition, as Figure 6a,c shows, the results from the ascending and descending data are significantly different in spatial distribution and magnitude. The maximum deformation rate of the ascending data is −156 mm/yr, but for the descending data, it is −129 mm/yr. The reasons for this could be (1) due to the differences in the observation direction, the incidence angle and heading angle of the ascending and descending radar satellites, the projection of the deformation on the LOS of different attitudes at the same location must be different; (2) mountainous areas with great terrain undulations result in geometric distortion in SAR images, which will directly affect the deformation results; (3) due to the changes in slope and aspect, the results of the deformation in the ascending and descending LOS directions of the same location may be opposite.
In this study area, the PS-InSAR method selected no point target at the landslide bodies and could not perform subsequent experimental analysis, indicating that the method has great defects in low-coherence research areas, such as for landslide monitoring. The DS-InSAR method selected dense point targets and provided rich deformation details in the area, which was very helpful for subsequent deformation analysis and discussion.
Due to the absence of leveling and GNSS data, we verified the reliability of the DS-InSAR monitoring results by its correlation with the results of PS-InSAR [48,49]. Figure 7 shows the correlation between the PS-InSAR and DS-InSAR deformation rate results. Since both the PS-InSAR and DS-InSAR results are dependent variables, the correlation coefficient (Corr) proposed by Karl Pearson [50] was used to verify the results of the DS technology. The Corr is a statistical indicator for reflecting the correlation between variables. The closer the Corr value is to 1, the closer the correlation is between the two results. In the ascending orbit, the Corr is 0.83 and the root mean square error (RMSE) is 4 mm/yr; in the descending orbit, the Corr is 0.77 and the RMSE is 5 mm/yr. The results of these two methods are highly correlated, which proves the reliability of DS-InSAR in monitoring the surface deformation of the landslide area.
After analyzing and verifying the deformation rate results, we analyzed the time series deformation in the landslide area to understand the surface deformation characteristics more comprehensively. Figure 8 is a three-dimensional optical image feature map of the landslide area. The yellow curves outline the identification boundaries of the two landslide bodies 1 and 2. Considering the geomorphology and slopes, we selected three points from each landslide body and analyzed their time series deformation trends to understand the law of surface change in the landslide area.
Figure 9 shows the time series deformation trends of points P1–P6 derived from ascending and descending data between 5 November 2014 and 4 September 2019. The six feature points all have varying degrees of deformation and the cumulative deformation statistics of feature points are shown in Table 2. P1 is located in the west-central part of landslide body 1. The ascending result shows great deformation, with a cumulative deformation of −342 mm. The descending result shows a slight uplift (33 mm), which may be a result of the electromagnetic wave emission direction and the radar incidence angle. P2 is located in the south-central part of landslide body 1. The ascending result shows large deformation, with a cumulative deformation of −375 mm. The deformation in the descending result is smaller, with a cumulative deformation of −85 mm. P3 is located in the east-central part of landslide body 1. It has smaller deformation in the ascending result (−163 mm) than P1 and P2, but the change trend is obvious. The descending orbit result is smaller and fluctuates greatly, indicating that the monitoring results were greatly affected by noise. P4 is located in the upper part of landslide body 2. The deformation of the ascending result is large (−209 mm), but the deformation trend is not smooth. The deformation in the descending result is −73 mm. P5 is located in the middle of landslide body 2. The deformations in the ascending and descending results are −223 mm and −220 mm, respectively. P6 is located in the lower part of landslide body 2. The deformation in the ascending orbit result is −284 mm. The deformation in the descending result has a similar trend, but the value is smaller, at −214 mm.
As Figure 9 shows, the six feature points all experienced deformation acceleration in August 2017, but the deformations in the ascending and descending results of each feature point are quite different. Thus, it was difficult to get a unified summary for the surface deformation law according to the time series deformation trend. Therefore, we established a two-dimensional deformation model using the results of the ascending and descending orbits to further analyze the landslide area.

4.2. Two-Dimensional Deformation Results and Analysis

As the geographic locations of the monitoring points in the ascending orbit SAR image and those in the descending orbit SAR image may not be exactly the same, the homonymy points had to be selected before the two-dimensional deformation modeling. In this study, if the monitoring points on the two orbits had a distance of less than 10 m, they were selected as the homonymy points, and a total of 76,943 monitoring points were obtained (Figure 10).
According to the distribution results of the homonymy points, we obtained the normal deformation and the east–west deformation by Equations (12) and (14). As shown in Figure 11a,b, the two landslide bodies have obvious deformations, but the other regions are stable. The maximum surface deformation rate in the landslide bodies was −80 mm/yr in the normal direction and 118 mm/yr in the east–west direction. To analyze the surface deformation characteristics more clearly and intuitively, the region in the red dashed box was displayed in a three-dimensional map, as shown in Figure 11c,d. In the normal direction, the deformation rate was between −70mm/yr and −40mm/yr, and higher in the east–west direction. The regions with deformation rates greater than −70mm/yr are concentrated in the middle of landslide body 2. In the east–west direction, the surface deformation rates in the center area of the two landslide bodies were high (60 mm/yr to 90 mm/yr), especially in the middle part of landslide body 1 (>90 mm/yr). The deformation rate of the rest area is about 30 mm/yr and not significant.
Through the two-dimensional calculation of the landslide area from the normal and east–west directions, and combined with the actual surface model, the surface deformation analysis was more comprehensive, specific, and vivid. Therefore, compared to the single LOS direction deformation analysis, the two-dimensional deformation is beneficial for comprehensively analyzing the surface deformation characteristics of the landslide area, summarizing the surface deformation laws, as well as providing help for the subsequent discussion of the influence of precipitation and the deformation mechanism.

4.3. Stability Analysis

In order to further study and analyze the surface stability of landslide bodies 1 and 2, we analyzed the stability of the P1–P6 feature points. For the R/S analysis of feature points, we acquired 106 images from 5 November 2014 to 4 September 2019, which imposed a huge calculation burden. We divided the study period into six segments and calculated the Hurst index of each segment separately. The results in the normal direction and east–west direction are shown in Table 3 and Table 4, respectively. The distribution of the calculation results is shown in Figure 12.
From Table 3 and Figure 12a we found that, in the normal direction, the Hurst index of each feature point is greater than 0.5 at each segment, indicating that the time series deformation results have a long-term correlation. From 10 February 2017 to 26 September 2017, the Hurst index values of each feature point are low, and the Hurst index values of P2, P4, and P5 are all lower than 0.9, which is a little lower than the results of the other periods and also confirms the acceleration trend of deformation after the precipitation peak in 2017 obtained in the previous study. From 7 January 2019 to 4 September 2019, the Hurst index value of each feature point is greater than 0.9, indicating that the normal direction deformation trend will continue for some time. From Table 4 and Figure 12b, we found that in the east–west direction, the Hurst index calculated at each feature point is also greater than 0.5 at each segment, indicating that the time series deformation results have a long-term correlation. From 10 February 2017 to 26 September 2017, the Hurst index values of P5 and P6 are obviously relatively low, at 0.736 and 0.748, respectively, while the other points are greater than 0.9. From 7 January 2019 to 4 September 2019, the Hurst index value of P4 is 0.858, and those of the other points are all above 0.9, indicating that the east–west deformation trend will continue for some time.
As Figure 12 shows, the Hurst index values of P4–P6 have larger fluctuations than those of P1–P3, especially between 7 January 2019 and 4 September 2019. This means that after the peak of precipitation in 2017, the influence of precipitation on landslide body 2 was greater than on landslide body 1. In addition, the Hurst index values of the P1–P6 feature points are higher than 0.5 in both the normal and east–west directions, indicating that the current deformation trend will continue. It is safe to infer that the surface deformation of the two landslide bodies may develop along the current trend. In this study, we adopted the R/S analysis method, which not only analyzes the stability and reliability of the surface deformation trend during the study period but also predicts the deformation trend of the subsequent landslide surface and provides effective theoretical support for local disaster warning and reduction.

5. Discussion

5.1. Influence of Precipitation

In order to explore the inducing factors of the surface deformation of the Woda landslide area, we analyzed the relationships between the precipitation and normal deformation and the east–west deformation time series. Since the time span of the ascending orbit SAR images is not completely consistent with that of the descending orbit SAR images, based on the time span of the ascending SAR images (5 November 2014–4 September 2019), we performed time linear interpolation on the time series deformation results of the descending SAR images. We calculated the two-dimensional time series deformation of the ascending and descending data. The precipitation data were processed according to the temporal baseline length between the SAR images of the ascending orbit, which makes the analysis more reasonable.
Figure 13 shows the normal direction time series deformation and precipitation distribution map of P1–P6. Figure 14 shows the east–west direction time series deformation and precipitation distribution map of P1–P6. As the two figures show, the precipitation was mostly concentrated in May to September each year. Therefore, in this study, we focused on the deformation trend after the precipitation peak period from 2016 to 2018, as shown by the green dashed boxes in Figure 13 and Figure 14. As Figure 13a and Figure 14a show, after the precipitation peak period in 2016, the normal direction deformations of P1 and P3 had accelerated trends. After the precipitation peak period in 2017, the normal direction and east–west direction deformations of P1 and P2 had obvious acceleration. The east–west direction deformation of P3 accelerated in May, in the rainy season. After the precipitation peak period in 2018, the east–west deformation of P1 accelerated, and the deformation of the other points did not change significantly. In Figure 13b and Figure 14b, after the precipitation peak period in 2016, only the normal deformation of P6 had an acceleration trend. After the precipitation peak period in 2017, the normal direction deformation and east–west direction deformations of P4–P6 show an obvious acceleration. P6 shows the most obvious acceleration and maintained this trend for some time. After the precipitation peak period in 2018, the deformations of the P4–P6 feature points are not obvious. In summary, the P1–P6 feature points were most affected by the precipitation peak period in 2017.
To further analyze the relationship between precipitation and the time series deformation, more comprehensive statistics of precipitation data are listed in Table 5. In May–September from 2015 to 2019, the precipitation in May–September 2017 was the highest, reaching 676.17 mm, which is 28.5% higher than in 2016. Table 5 also lists the maximum daily precipitation rates throughout the whole study period. Three of the five days with the highest precipitation are in 2017. The daily precipitation of 13 May 2017 reached 113.69 mm, which is a heavy rainstorm; the daily precipitation for the two days of 7 September 2017 and 4 July 2017 also reached the heavy rain level. Combining the deformation trends of P1–P6 in the LOS direction (Figure 9), normal direction (Figure 13), and east–west direction (Figure 14), from Table 5, we found that the deformation of each point accelerated during August to November 2017, indicating that precipitation was an important factor leading to the deformation acceleration; the precipitation from May to September in 2015 and 2016 was less than that in 2017 and 2018, so the effect of precipitation on the deformation acceleration is not significant; extremely rainy weather and high precipitation are some of the main reasons for accelerated surface deformation in this landslide area. This provides a useful reference for analyzing and judging the trend of surface deformation in the landslide area.

5.2. Discussion of Deformation Mechanism

The Woda landslide area is located in a typical alpine and canyon landform, through which the Jinsha River runs. In this area, the Jinsha River is about 115 m wide, and the terrain of the river is gentle. The landslide bodies have circle chair shapes, with multi-level terraces [5]. A large amount of landslide material is accumulated on the slope. The average slope of the landslide bodies is 20°–25°, and the maximum slope at the front edge of the landslide bodies can reach 70°. In addition, the area has a plateau cold temperate semi-humid climate, with uneven precipitation throughout the year. Heavy rains occur occasionally in the area, which has become a significant cause of geological disasters.
As Figure 15 shows, the lithology of this area is mainly composed of schist, slate, and shale of the Upper Triassic Lanashan Formation and of crystalline limestone and feldspar quartz sandstone of the Upper Triassic Tumgou Formation. Under the influence of long-term tectonic movement, weathering, and rainfall, rock mass joints and fissures are abundant and the structure is broken. Therefore, the deformation of the Woda landslide area is controlled by geographical factors, such as stratum lithology and geological structure. In addition, external factors, such as steep terrain, climate environment, and rainfall, promote deformation.
Once a large-scale slide occurs, the materials may block the Jinsha River, forming a dammed lake and inducing indirect disasters, such as river overtopping, dam breaking, and flooding. Therefore, long-term monitoring of the surface deformation with the InSAR method is necessary, which can be used for analyzing and judging the landslide scale, deformation trends, and influencing factors. It plays an important role in improving emergency management and the prevention and mitigation of landslide disasters.

5.3. Limitations and Prospects

In this study, surface deformation monitoring and stability analysis of the Woda landslide area were successfully carried out. It is undeniable that this study has certain limitations. (1) Due to the influence of geographical location and topography, there is a lack of leveling and GNSS data for support and verification; (2) only SAR data was used, and there was a lack of optical satellite and remote sensing images for collaborative monitoring and analysis; (3) restricted by relevant management regulations, the collection of geological data was not rich enough, and there was no perfect in-depth analysis of the geological structure and lithology.
In future research, with the development of new observation platforms, such as the BeiDou Navigation Satellite System (BDS) and UAV, it is possible to combine navigation and positioning data, UAV photography data, SAR data of different polarization modes, and multi-platform spaceborne SAR data to carry out surface deformation monitoring fieldwork. In the field of landslide surface deformation monitoring, a scientific and reliable monitoring and evaluation system will be formed.

6. Conclusions

This study adopted the state-of-the-art of DS-InSAR method to extract the surface deformation of the Woda landslide area. The HTCI algorithm was used to select SHPs, and the adaptive spatial nonlocal filtering method was combined with the SHP selection result to optimize the phases. Both DS and PS points were used to increase the density of monitoring points in the study area.
From the 106 Sentinel-1A ascending images and 102 Sentinel-1A descending images, we derived the time series deformation of the Woda landslide area for nearly 5 years by the DS-InSAR method. The two landslide bodies have obvious deformations, but the other areas are stable. Compared to the traditional PS-InSAR method, the DS-InSAR method increased the monitoring points density significantly. The results of the monitoring points were verified by the Karl Pearson correlation coefficients, which are 0.83 for the ascending result and 0.77 for the descending result. The obtained two-dimensional deformation model shows that the maximum deformation rate of the normal direction is −80 mm/yr, and for the east–west direction, it is 118 mm/yr. The stability analysis results show that the current deformation trend is sustainable. The time linear interpolation analysis of the six feature points suggests that the deformation is controlled by geographical factors, such as stratum lithology and geological structure, and it is affected by external factors, such as precipitation. The research results and analysis method in this paper provide reference and guidance for the early investigation, monitoring, and timely warning of dense forest landslides.

Author Contributions

Y.L., H.Y., and S.W. conceived of the original design of this paper. Y.L. and H.Y. improve the structure of the paper. L.X. and J.P. provided comments on this paper. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (42174026), the National Key Research and Development Program of China (2021YFE011004), and the Open Fund of State Key Laboratory of Coal Resources and Safe Mining (Grant No. SKLCRSM20KFA12).

Data Availability Statement

Data sharing is not applicable.

Acknowledgments

The authors would like to thank the European Space Agency for providing Sentinel-1A data and the National Aeronautics and Space Administration for providing SRTM1 DEM data. In addition, thanks to Google for the optical imagery.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

  1. Liu, C.; Li, W.; Wu, H.; Lu, P.; Sang, K.; Sun, W.; Chen, W.; Hong, Y.; Li, R. Susceptibility evaluation and mapping of China’s landslides based on multi-source data. Nat. Hazards 2013, 69, 1477–1495. [Google Scholar] [CrossRef]
  2. Delaney, K.B.; Evans, S.G. The 2000 Yigong landslide (Tibetan Plateau), rockslide-dammed lake and outburst flood: Review, remote sensing analysis, and process modelling. Geomorphology 2015, 246, 377–393. [Google Scholar] [CrossRef]
  3. Fan, X.; Xu, Q.; Alonso-Rodriguez, A.; Subramanian, S.S.; Li, W.; Zheng, G.; Dong, X.; Huang, R. Successive landsliding and damming of the Jinsha River in eastern Tibet, China: Prime investigation, early warning, and emergency response. Landslides 2019, 16, 1003–1020. [Google Scholar] [CrossRef]
  4. Sijing, W.; Guohe, L.; Qiang, Z.; Chaoli, L. Engineering Geological Study of the Active Tectonic Region for Hydropower Development on the Jinsha River, Upstream of the Yangtze River. Acta Geol. Sin. 2000, 74, 353–361. [Google Scholar] [CrossRef]
  5. Feng, W.; Dun, J.; Yi, X.; Zhang, G. Deformation analysis of Woda village old landslide in Jinsha river basin using SBAS-InSAR technology. J. Eng. Geol. 2020, 28, 384–393. (In Chinese) [Google Scholar]
  6. Ouyang, C.; An, H.; Zhou, S.; Wang, Z.; Su, P.; Wang, D.; Cheng, D.; She, J. Insights from the failure and dynamic characteristics of two sequential landslides at Baige village along the Jinsha River, China. Landslides 2019, 16, 1397–1414. [Google Scholar] [CrossRef]
  7. Stiros, S.C.; Vichas, C.; Skourtis, C. Landslide Monitoring Based On Geodetically Derived Distance Changes. J. Surv. Eng. 2004, 130, 156–162. [Google Scholar] [CrossRef]
  8. Du, J.-C.; Teng, H.-C. 3D laser scanning and GPS technology for landslide earthwork volume estimation. Autom. Constr. 2007, 16, 657–663. [Google Scholar] [CrossRef]
  9. Muntean, A.; Mocanu, V.; Ambrosius, B. A GPS study of land subsidence in the Petrosani (Romania) coal mining area. Nat. Hazards 2016, 80, 797–810. [Google Scholar] [CrossRef]
  10. Casagli, N.; Frodella, W.; Morelli, S.; Tofani, V.; Ciampalini, A.; Intrieri, E.; Raspini, F.; Rossi, G.; Tanteri, L.; Lu, P. Spaceborne, UAV and ground-based remote sensing techniques for landslide mapping, monitoring and early warning. Geoenvironmental Disasters 2017, 4, 9. [Google Scholar] [CrossRef]
  11. Bayer, B.; Simoni, A.; Mulas, M.; Corsini, A.; Schmidt, D. Deformation responses of slow moving landslides to seasonal rainfall in the Northern Apennines, measured by InSAR. Geomorphology 2018, 308, 293–306. [Google Scholar] [CrossRef]
  12. Liu, Y.; Fan, H.; Wang, L.; Zhuang, H. Monitoring of surface deformation in a low coherence area using distributed scatterers InSAR: Case study in the Xiaolangdi Basin of the Yellow River, China. Bull. Eng. Geol. Environ. 2021, 80, 25–39. [Google Scholar] [CrossRef]
  13. Yin, Y.; Zheng, W.; Liu, Y.; Zhang, J.; Li, X. Integration of GPS with InSAR to monitoring of the Jiaju landslide in Sichuan, China. Landslides 2010, 7, 359–365. [Google Scholar] [CrossRef]
  14. Chen, Y.; Yu, S.; Tao, Q.; Liu, G.; Wang, L.; Wang, F. Accuracy verification and correction of D-InSAR and SBAS-InSAR in monitoring mining surface subsidence. Remote Sens. 2021, 13, 4365. [Google Scholar] [CrossRef]
  15. Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
  16. Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef] [Green Version]
  17. Werner, C.; Wegmuller, U.; Strozzi, T.; Wiesmann, A. Interferometric point target analysis for deformation mapping. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, Toulouse, France, 21–25 July 2003; Institute of Electrical and Electronics Engineers: Muri, Switzerland; pp. 4362–4364. [Google Scholar]
  18. Zhang, L.; Lu, Z.; Ding, X.; Jung, H.-S.; Feng, G.; Lee, C.-W. Mapping ground surface deformation using temporarily coherent point SAR interferometry: Application to Los Angeles Basin. Remote Sens. Environ. 2012, 117, 429–439. [Google Scholar] [CrossRef]
  19. Osmanoğlu, B.; Sunar, F.; Wdowinski, S.; Cabral-Cano, E. Time series analysis of InSAR data: Methods and trends. ISPRS J. Photogramm. Remote Sens. 2016, 115, 90–102. [Google Scholar] [CrossRef]
  20. Del Soldato, M.; Riquelme, A.; Bianchini, S.; Tomàs, R.; Di Martire, D.; De Vita, P.; Moretti, S.; Calcaterra, D. Multisource data integration to investigate one century of evolution for the Agnone landslide (Molise, southern Italy). Landslides 2018, 15, 2113–2128. [Google Scholar] [CrossRef] [Green Version]
  21. Rosi, A.; Tofani, V.; Tanteri, L.; Tacconi Stefanelli, C.; Agostini, A.; Catani, F.; Casagli, N. The new landslide inventory of Tuscany (Italy) updated with PS-InSAR: Geomorphological features and landslide distribution. Landslides 2018, 15, 5–19. [Google Scholar] [CrossRef] [Green Version]
  22. Intrieri, E.; Raspini, F.; Fumagalli, A.; Lu, P.; Del Conte, S.; Farina, P.; Allievi, J.; Ferretti, A.; Casagli, N. The Maoxian landslide as seen from space: Detecting precursors of failure with Sentinel-1 data. Landslides 2018, 15, 123–133. [Google Scholar] [CrossRef] [Green Version]
  23. Komac, M.; Holley, R.; Mahapatra, P.; van der Marel, H.; Bavec, M. Coupling of GPS/GNSS and radar interferometric data for a 3D surface displacement monitoring of landslides. Landslides 2015, 12, 241–257. [Google Scholar] [CrossRef]
  24. Fobert, M.-A.; Singhroy, V.; Spray, J.G. InSAR monitoring of landslide activity in Dominica. Remote Sens. 2021, 13, 815. [Google Scholar] [CrossRef]
  25. Ferretti, A.; Fumagalli, A.; Novali, F.; Prati, C.; Rocca, F.; Rucci, A. A New Algorithm for Processing Interferometric Data-Stacks: SqueeSAR. IEEE Trans. Geosci. Remote Sens. 2011, 49, 3460–3470. [Google Scholar] [CrossRef]
  26. Goel, K.; Adam, N. A Distributed Scatterer Interferometry Approach for Precision Monitoring of Known Surface Deformation Phenomena. IEEE Trans. Geosci. Remote Sens. 2014, 52, 5454–5468. [Google Scholar] [CrossRef]
  27. Fornaro, G.; Verde, S.; Reale, D.; Pauciullo, A. CAESAR: An Approach Based on Covariance Matrix Decomposition to Improve Multibaseline–Multitemporal Interferometric SAR Processing. IEEE Trans. Geosci. Remote Sens. 2014, 53, 2050–2065. [Google Scholar] [CrossRef]
  28. Jiang, M.; Ding, X.; Hanssen, R.F.; Malhotra, R.; Chang, L. Fast Statistically Homogeneous Pixel Selection for Covariance Matrix Estimation for Multitemporal InSAR. IEEE Trans. Geosci. Remote Sens. 2014, 53, 1213–1224. [Google Scholar] [CrossRef]
  29. Jiang, M.; Guarnieri, A.M. Distributed Scatterer Interferometry With the Refinement of Spatiotemporal Coherence. IEEE Trans. Geosci. Remote Sens. 2020, 58, 3977–3987. [Google Scholar] [CrossRef]
  30. Dong, J.; Zhang, L.; Tang, M.; Liao, M.; Xu, Q.; Gong, J.; Ao, M. Mapping landslide surface displacements with time series SAR interferometry by combining persistent and distributed scatterers: A case study of Jiaju landslide in Danba, China. Remote Sens. Environ. 2018, 205, 180–198. [Google Scholar] [CrossRef]
  31. Wang, J.; Wang, C.; Xie, C.; Zhang, H.; Tang, Y.; Zhang, Z.; Shen, C. Monitoring of large-scale landslides in Zongling, Guizhou, China, with improved distributed scatterer interferometric SAR time series methods. Landslides 2020, 17, 1777–1795. [Google Scholar] [CrossRef]
  32. Hanssen, R.F. Radar Interferometry; Springer: Dordrecht, The Netherlands, 2001. [Google Scholar] [CrossRef] [Green Version]
  33. Fan, H.; Wang, L.; Wen, B.; Du, S. A New Model for three-dimensional Deformation Extraction with Single-track InSAR Based on Mining Subsidence Characteristics. Int. J. Appl. Earth Obs. Geoinf. 2021, 94, 102223. [Google Scholar] [CrossRef]
  34. Béjar-Pizarro, M.; Notti, D.; Mateos, R.M.; Ezquerro, P.; Centolanza, G.; Herrera, G.; Bru, G.; Sanabria, M.; Solari, L.; Duro, J.; et al. Mapping Vulnerable Urban Areas Affected by Slow-Moving Landslides Using Sentinel-1 InSAR Data. Remote Sens. 2017, 9, 876. [Google Scholar] [CrossRef] [Green Version]
  35. Confuorto, P.; Di Martire, D.; Centolanza, G.; Iglesias, R.; Mallorqui, J.J.; Novellino, A.; Plank, S.; Ramondini, M.; Thuro, K.; Calcaterra, D. Post-failure evolution analysis of a rainfall-triggered landslide by multi-temporal interferometry SAR approaches integrated with geotechnical analysis. Remote Sens. Environ. 2017, 188, 51–72. [Google Scholar] [CrossRef]
  36. Sun, Q.; Hu, J.; Zhang, L.; Ding, X. Towards Slow-Moving Landslide Monitoring by Integrating Multi-Sensor InSAR Time Series Datasets: The Zhouqu Case Study, China. Remote Sens. 2016, 8, 908. [Google Scholar] [CrossRef] [Green Version]
  37. Herrera, G.; Gutiérrez, F.; García-Davalillo, J.C.; Guerrero, J.; Notti, D.; Galve, J.P.; Fernández-Merodo, J.A.; Cooksley, G. Multi-sensor advanced DInSAR monitoring of very slow landslides: The Tena Valley case study (Central Spanish Pyrenees). Remote Sens. Environ. 2013, 128, 31–43. [Google Scholar] [CrossRef]
  38. Li, T.; Zhang, H.; Fan, H.; Zheng, C.; Liu, J. Position inversion of goafs in deep coal seams based on DS-InSAR data and the probability integral methods. Remote Sens. 2021, 13, 2898. [Google Scholar] [CrossRef]
  39. Abrahams, J. Probability, Random Variables, and Stochastic Processes by Athanasios Papoulis. J. Am. Stat. Assoc. 1984, 79, 957. [Google Scholar] [CrossRef]
  40. Neuhäuser, M. Exact tests based on the Baumgartner-Wei-Schindler statistic—A survey. Stat. Pap. 2005, 46, 1–29. [Google Scholar] [CrossRef]
  41. Jiang, M.; Miao, Z.; Gamba, P.; Yong, B. Application of Multitemporal InSAR Covariance and Information Fusion to Robust Road Extraction. IEEE Trans. Geosci. Remote Sens. 2017, 55, 3611–3622. [Google Scholar] [CrossRef]
  42. Deledalle, C.A.; Denis, L.; Tupin, F. NL-InSAR: Nonlocal interferogram estimation. IEEE Trans. Geosci. Remote Sens. 2011, 49, 1441–1452. [Google Scholar] [CrossRef]
  43. Hu, J.; Li, Z.-W.; Ding, X.; Zhu, J.; Zhang, L.; Sun, Q. Resolving three-dimensional surface displacements from InSAR measurements: A review. Earth Sci. Rev. 2014, 133, 1–17. [Google Scholar] [CrossRef]
  44. Fuhrmann, T.; Garthwaite, M.C. Resolving three-dimensional surface motion with InSAR: Constraints from multi-geometry data fusion. Remote Sens. 2019, 11, 241. [Google Scholar] [CrossRef] [Green Version]
  45. Colesanti, C.; Wasowski, J. Investigating landslides with space-borne Synthetic Aperture Radar (SAR) interferometry. Eng. Geol. 2006, 88, 173–199. [Google Scholar] [CrossRef]
  46. Hurst, H.E. Long-Term Storage Capacity of Reservoirs. Trans. Am. Soc. Civ. Eng. 1951, 116, 770–799. [Google Scholar] [CrossRef]
  47. Pei, L.; Chen, J.; Zhou, J.; Huang, H.; Zhou, Z.; Chen, C.; Yao, F. A Fractal Prediction Method for Safety Monitoring Deformation of Core Rockfill Dams. Math. Probl. Eng. 2021, 2021, 6655657. [Google Scholar] [CrossRef]
  48. Lv, X.; Yazici, B.; Zeghal, M.; Bennett, V.; Abdoun, T. Joint-Scatterer Processing for Time-Series InSAR. IEEE Trans. Geosci. Remote Sens. 2014, 52, 7205–7221. [Google Scholar]
  49. Samiei-Esfahany, S.; Martins, J.E.; Van Leijen, F.; Hanssen, R.F. Phase estimation for distributed scatterers in InSAR stacks using integer least squares estimation. IEEE Trans. Geosci. Remote Sens. 2016, 54, 5671–5687. [Google Scholar] [CrossRef] [Green Version]
  50. Osmani, S.A.; Banik, B.K.; Ali, H. Integrating fuzzy logic with Pearson correlation to optimize contaminant detection in water distribution system with uncertainty analyses. Environ. Monit. Assess. 2019, 191, 441. [Google Scholar] [CrossRef]
Figure 1. (a) Location of the study area in Tibet Autonomous Region, China, and SRTM-derived topographic map of the Woda landslide area (Jomda County). The red and blue boxes indicate the coverage of Sentinel-1A ascending and descending data, respectively. (b) Two-dimensional Google map imagery of the study area. (c) Three-dimensional Google map imagery of the study area.
Figure 1. (a) Location of the study area in Tibet Autonomous Region, China, and SRTM-derived topographic map of the Woda landslide area (Jomda County). The red and blue boxes indicate the coverage of Sentinel-1A ascending and descending data, respectively. (b) Two-dimensional Google map imagery of the study area. (c) Three-dimensional Google map imagery of the study area.
Remotesensing 14 00532 g001
Figure 2. The flowchart of landslide monitoring by the DS-InSAR method.
Figure 2. The flowchart of landslide monitoring by the DS-InSAR method.
Remotesensing 14 00532 g002
Figure 3. Spatial and temporal baselines of the interferograms obtained by (a) ascending data; (b) descending data.
Figure 3. Spatial and temporal baselines of the interferograms obtained by (a) ascending data; (b) descending data.
Remotesensing 14 00532 g003
Figure 4. Three-dimensional spatial decomposition model of the LOS deformation. D L O S 1 and D L O S 2 are the LOS deformations of the ascending and descending orbits, respectively; θ 1 and θ 2 are the incidence angles of the ascending and descending SAR images, respectively; α 1 and α 2 are the heading angles of the ascending and descending SAR satellites.
Figure 4. Three-dimensional spatial decomposition model of the LOS deformation. D L O S 1 and D L O S 2 are the LOS deformations of the ascending and descending orbits, respectively; θ 1 and θ 2 are the incidence angles of the ascending and descending SAR images, respectively; α 1 and α 2 are the heading angles of the ascending and descending SAR satellites.
Remotesensing 14 00532 g004
Figure 5. Radar incidence angle distribution map: (a) ascending orbit SAR image; (b) descending orbit SAR image.
Figure 5. Radar incidence angle distribution map: (a) ascending orbit SAR image; (b) descending orbit SAR image.
Remotesensing 14 00532 g005
Figure 6. The LOS deformation rate maps of the Woda landslide area (a) derived from ascending data using the DS-InSAR method; (b) derived from ascending data using the PS-InSAR method; (c) derived from descending data using the DS-InSAR method; (d) derived from descending data using the PS-InSAR method.
Figure 6. The LOS deformation rate maps of the Woda landslide area (a) derived from ascending data using the DS-InSAR method; (b) derived from ascending data using the PS-InSAR method; (c) derived from descending data using the DS-InSAR method; (d) derived from descending data using the PS-InSAR method.
Remotesensing 14 00532 g006
Figure 7. Correlation between PS-InSAR and DS-InSAR deformation rates in the (a) ascending orbit and the (b) descending orbit.
Figure 7. Correlation between PS-InSAR and DS-InSAR deformation rates in the (a) ascending orbit and the (b) descending orbit.
Remotesensing 14 00532 g007
Figure 8. Three-dimensional optical image feature map of the Woda landslide area.
Figure 8. Three-dimensional optical image feature map of the Woda landslide area.
Remotesensing 14 00532 g008
Figure 9. Time series deformation of points P1–P6 (af) based on ascending and descending results.
Figure 9. Time series deformation of points P1–P6 (af) based on ascending and descending results.
Remotesensing 14 00532 g009
Figure 10. The distribution of the homonymy points: (a) ascending orbit results; (b) descending orbit results.
Figure 10. The distribution of the homonymy points: (a) ascending orbit results; (b) descending orbit results.
Remotesensing 14 00532 g010
Figure 11. Two-dimensional deformation results (a) in the normal direction and (b) in the east–west direction. Three-dimensional model of the deformation (c) in the normal direction and (d) in the east–west deformation.
Figure 11. Two-dimensional deformation results (a) in the normal direction and (b) in the east–west direction. Three-dimensional model of the deformation (c) in the normal direction and (d) in the east–west deformation.
Remotesensing 14 00532 g011
Figure 12. Hurst index values of P1–P6 feature points (a) in the normal direction and (b) in the east–west direction.
Figure 12. Hurst index values of P1–P6 feature points (a) in the normal direction and (b) in the east–west direction.
Remotesensing 14 00532 g012
Figure 13. Time series deformation in the normal direction and precipitation distribution map of the P1–P6 feature points. (a) includes P1–P3 feature points; (b) includes P4–P6 feature points. The green dashed boxes are the deformation trend of the feature points after the period of concentrated precipitation.
Figure 13. Time series deformation in the normal direction and precipitation distribution map of the P1–P6 feature points. (a) includes P1–P3 feature points; (b) includes P4–P6 feature points. The green dashed boxes are the deformation trend of the feature points after the period of concentrated precipitation.
Remotesensing 14 00532 g013
Figure 14. Time series deformation in the east–west direction and precipitation distribution map of the P1–P6 feature points. (a) includes P1–P3 feature points; (b) includes P4–P6 feature points. The green dashed boxes are the deformation trend of the feature points after the period of concentrated precipitation.
Figure 14. Time series deformation in the east–west direction and precipitation distribution map of the P1–P6 feature points. (a) includes P1–P3 feature points; (b) includes P4–P6 feature points. The green dashed boxes are the deformation trend of the feature points after the period of concentrated precipitation.
Remotesensing 14 00532 g014
Figure 15. Geological map of the study area. (1) Upper Triassic; (2) Upper Triassic diorite; (3) Upper Triassic ultrabasic rock blocks; (4) diorite vein; (5) stratigraphic occurrence; (6) actual blatt flaws; (7) supposed blatt flaws; (8) stratigraphic boundary; (9) identified landslides.
Figure 15. Geological map of the study area. (1) Upper Triassic; (2) Upper Triassic diorite; (3) Upper Triassic ultrabasic rock blocks; (4) diorite vein; (5) stratigraphic occurrence; (6) actual blatt flaws; (7) supposed blatt flaws; (8) stratigraphic boundary; (9) identified landslides.
Remotesensing 14 00532 g015
Table 1. Detailed information of the dataset.
Table 1. Detailed information of the dataset.
SensorTemporal CoverageImage NumberOrbitTrackPolarization
Sentinel-1A5 November 2014–4 September 2019106Ascending99VV
Sentinel-1A31 October 2014–11 September 2019102Descending33VV
Table 2. Cumulative deformation statistics of feature points.
Table 2. Cumulative deformation statistics of feature points.
Feature PointsCumulative Deformation (mm)
AscendingDescending
P1−34233
P2−375−85
P3−163−42
P4−209−73
P5−223−220
P6−284−214
Table 3. Hurst index calculation results in the normal direction.
Table 3. Hurst index calculation results in the normal direction.
No.Period of TimeHurst Index
P1P2P3P4P5P6
15 November 2014–18 December 20150.9340.9190.8420.9380.9530.847
218 December 2015–10 February 20170.9560.9490.9070.9420.9320.939
310 February 2017–26 September 20170.9480.8180.9430.7740.7260.900
426 September 2017–12 May 20180.9320.9390.8170.9320.9050.935
512 May 2018–7 January 20190.8170.9490.9560.8650.7480.924
67 January 2019–4 September 20190.9340.9450.9350.9040.9370.925
Table 4. Hurst index calculation results in the east–west direction.
Table 4. Hurst index calculation results in the east–west direction.
No.Period of TimeHurst Index
P1P2P3P4P5P6
15 November 2014–18 December 20150.9320.9380.8670.9630.8400.931
218 December 2015–10 February 20170.9340.9380.8280.9450.8760.936
310 February 2017–26 September 20170.9090.9110.9560.9420.7360.748
426 September 2017–12 May 20180.9250.9290.9270.9080.9590.922
512 May 2018–7 January 20190.9470.9040.9270.8810.8490.915
67 January 2019–4 September 20190.9460.9080.9350.8580.9280.905
Table 5. Detailed statistics on precipitation data.
Table 5. Detailed statistics on precipitation data.
Annual Precipitation from May to SeptemberPrecipitation per Day Ranking
DatePrecipitation (mm)DatePrecipitation (mm)
2015597.5913 May 2017113.69
2016526.277 September 201791.71
2017676.172 July 201986.79
2018642.354 July 201765.02
2019559.486 August 201545.61
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, Y.; Yang, H.; Wang, S.; Xu, L.; Peng, J. Monitoring and Stability Analysis of the Deformation in the Woda Landslide Area in Tibet, China by the DS-InSAR Method. Remote Sens. 2022, 14, 532. https://doi.org/10.3390/rs14030532

AMA Style

Liu Y, Yang H, Wang S, Xu L, Peng J. Monitoring and Stability Analysis of the Deformation in the Woda Landslide Area in Tibet, China by the DS-InSAR Method. Remote Sensing. 2022; 14(3):532. https://doi.org/10.3390/rs14030532

Chicago/Turabian Style

Liu, Youfeng, Honglei Yang, Shizheng Wang, Linlin Xu, and Junhuan Peng. 2022. "Monitoring and Stability Analysis of the Deformation in the Woda Landslide Area in Tibet, China by the DS-InSAR Method" Remote Sensing 14, no. 3: 532. https://doi.org/10.3390/rs14030532

APA Style

Liu, Y., Yang, H., Wang, S., Xu, L., & Peng, J. (2022). Monitoring and Stability Analysis of the Deformation in the Woda Landslide Area in Tibet, China by the DS-InSAR Method. Remote Sensing, 14(3), 532. https://doi.org/10.3390/rs14030532

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop