Large-Area Landslides Monitoring Using Advanced Multi-Temporal InSAR Technique over the Giant Panda Habitat, Sichuan, China

Abstract

The region near Dujiangyan City and Wenchuan County, Sichuan China, including significant giant panda habitats, was severely impacted by the Wenchuan earthquake. Large-area landslides occurred and seriously threatened the lives of people and giant pandas. In this paper, we report the development of an enhanced multi-temporal interferometric synthetic aperture radar (MTInSAR) methodology to monitor potential post-seismic landslides by analyzing coherent scatterers (CS) and distributed scatterers (DS) points extracted from multi-temporal l-band ALOS/PALSAR data in an integrated manner. Through the integration of phase optimization and mitigation of the orbit and topography-related phase errors, surface deformations in the study area were derived: the rates in the line of sight (LOS) direction ranged from −7 to 1.5 cm/a. Dozens of potential landslides, distributed mainly along the Minjiang River, Longmenshan Fault, and in other the high-altitude areas were detected. These findings matched the distribution of previous landslides. InSAR-derived results demonstrated that some previous landslides were still active; many unstable slopes have developed, and there are significant probabilities of future massive failures. The impact of landslides on the giant panda habitat, however ranged from low to moderate, would continue to be a concern for conservationists for some time in the future.

1. Introduction

Landslides are defined as the movement of a mass of rock, earth, or debris on hillslopes due to gravity and other triggering factors [1]. As one of the major natural disasters in the world, they cause serious damage to the natural environment and infrastructure, and considerable loss to the economy and of human and animal lives. In the Sichuan province China, home to China’s Giant Panda Habitat World Heritage site, widespread active faults, frequent earthquakes, cracked lithology, steep and rough topography, and localized torrential downpours render the province particularly vulnerable to landslides. To-date, more than 100 thousand landslides affecting over 100 cities have been recorded [2]; almost half of the total number of county-towns in Sichuan and habitats utilized by the giant panda has been affected. The destructive Wenchuan earthquake (Ms 8.0) that occurred on 12 May 2008 resulted in thousands of casualties and rendered a million people homeless [3]; it induced more than 5600 landslides covering a total area of 41,750 km² [4]. The strong tremors loosened and shattered bedrock, produced numerous cracked slopes that became unstable after the quake, and released a large number of falling stones. The frequency and scale of landslide occurrence increased significantly in the post-seismic period, compared to the pre-seismic period [5].

To reduce the effect of landslides, risk analysis and susceptibility assessment are crucial [6]. Deformation of unstable slopes is often used as the indicator for landslide characterization, mapping and prediction [7]. Resolving the kinematics before or after an event is a priority for better understanding potential landslides. However, the use of traditional investigation methods, such as leveling, GPS and geophysical prospecting [8], are severely limited in mountainous areas where there are no access roads or paths. Due to their all-weather, day-night, broad coverage, and high spatiotemporal resolution capabilities [9], differential synthetic aperture radar interferometry (DInSAR) has great potential for landslide deformation monitoring [10]. The first application of DInSAR dates back to the mid-1990s [11]; since then it has become widely used and, more recently, multi-temporal DInSAR (MTInSAR) techniques (represented by Permanent Scatterers—PSInSAR [12], Small Baseline Subset—SBAS [13] and SqueeSAR [14]), capable of high-precision measurements (cm to mm range) and time series analysis, have become popular. MTInSAR techniques rely on processing of time series, repeat-pass SAR images, and identification of radar targets providing regularly measurable phase signals. They have been successfully applied in many areas of landslides [7,15,16,17,18,19,20,21,22,23,24].

Classical techniques have limitations for complex regions where few SAR images are available. PSInSAR is not applicable to outside urban regions due to the scarcity of permanent scatterers (PS). SBAS focuses on coherent scatterers (CS) by utilizing high-quality interferograms with small spatial and temporal baselines. Its use for sites with rough terrain and vegetation cover is quite challenging. Analyzing distributed scatterers (DS), allow the extraction of plenty of point targets in vegetated regions using the SqueeSAR approach; but a large number of SAR images (e.g., more than 25 scenes) is needed for a robust inversion. When using DInSAR for landslide applications several challenges, such as serious decorrelation caused by the vegetation cover, rough terrains, and large deformation gradients [25,26,27], phase propagation anomalies due to atmospheric variability [28,29], and geometric distortions in mountainous areas caused by the side-looking imaging must be addressed and resolved. Moreover, the complexity of landslide dynamics and various error sources could greatly increase the probabilities of deformation estimations and interpretations being incorrect [30].

In this study, in order to monitor slope instability phenomena in the vegetated, steep terrain of Sichuan giant panda habitats, we developed an enhanced landslide-oriented MTInSAR methodology for areas where availability of SAR images are limited by taking advantage of existing SBAS and SqueeSAR technologies. We derived a surface deformation field in the area using our methodology. Up to now, in our study area remote sensing techniques have been used to investigate landslide distribution and their destructive effects mainly for past landslides [31,32,33,34]. In contrast, our methodology enables the detection of slope instabilities for potential, future landslides by monitoring displacement anomalies, thereby significantly improving hazard prediction and risk mitigation capabilities.

Due to its coherence preservation capability (as shown by many successful landslide monitoring studies [35,36,37]), we used l-band ALOS PALSAR images. In this paper, we first introduce the study area and the data used; then we explain the entire processing chain of our approach and the improvements and advantages it brings to landslide detection and monitoring in comparison to classical techniques. We have analyzed the application of MTInSAR to landslide monitoring in combination with other multi-source data (e.g., Google earth images, ICESat/GLAS altitude data, survey control points, and field surveying) and discussed the deformation characteristics of landslides and their impacts on the giant panda habitat. Based on our analyses we have drawn conclusions that may be investigated in future studies.

2. Study Area and Data

2.1. Study Site

The study area is situated in the Dujiangyan City and Wenchuan County, Sichuan province, China (see inset in Figure 1). It is near the earthquake’s epicenter and belongs to the maximum meizoseismal area, especially close to the Yingxiu Town, where 80% of the houses collapsed due to the quake [38]. In a total area of about 4200 km², the terrain is very rugged and the elevation ranges from 500 m to more than 5000 m, with slope gradients reaching 72° in some places (Figure 1). The Minjiang River flows north-to-south through the deep valley and, after the Zipingpu reservoir, it drains into the Chengdu plain. Longmenshan Fault, trending northeast and consisting of the Mao-Wen, Yingxiu-Beichuan, and Guanxian-Anxian Faults, passes through the study area and runs parallel to the ridges and valleys. The rocks mainly comprise shale, sandstone, limestone, slate, granite, and basalt, from Precambrian to Cretaceous times [31]. The area is subtropical and receives humid monsoonal rains; annual mean air temperature is about 15 °C and yearly precipitation 1200 mm. 50%–70% of the rainfall occurs during the period from June to September [38], which is a significant factor triggering instability of the slopes and frequent landslides [39].

Figure 1. (a) SRTM-derived topographic map of the study area. The black triangles indicate locations of landslides; (b) Distribution of giant panda habitats in China and the location of our study area.

Figure 1. (a) SRTM-derived topographic map of the study area. The black triangles indicate locations of landslides; (b) Distribution of giant panda habitats in China and the location of our study area.

The study area is well known as the Giant Panda Habitat World Heritage site. The Wenchuan earthquake is estimated to have destroyed about 20% of the forest in the Wolong protection zone which is a key component of the World Heritage site [40]. Loss of forests fragmented the giant panda habitat and increased the number of habitat patches to more than six times than that before the earthquake [41]. The post-seismic disasters also threaten the survival of giant panda, and it’s reported that at least one of them was killed by the debris flow. The photos in Figure 2 illustrate the scale of frequent large-area landslides in our study area.

Figure 2. Photos of landslides from locations along the Minjiang River indicated by black triangles in Figure 1: (a) Yangdianjie; (b) Luoquanwan; (c) Yinxing; (d) Yingxiu.

Figure 2. Photos of landslides from locations along the Minjiang River indicated by black triangles in Figure 1: (a) Yangdianjie; (b) Luoquanwan; (c) Yinxing; (d) Yingxiu.

2.2. Data

In order to monitor potential landslides and evaluate their influences on human and giant panda habitats and lives, this study employed 11 scenes of l-band (radar wavelength of 23.6 cm) post-seismic ALOS PALSAR SLC images for SAR interferometry analysis (see Table 1). They were acquired in the ascending mode with a center incidence angle of about 34°, during the period from 22 June 2008 to 13 February 2011. Among them, six were acquired in the fine beam single polarization (FBS, range bandwidth of 28 MHz) (HH) mode and the others were acquired in fine beam dual polarization (FBD, 14 MHz) (HH + HV) mode. In order to carry out interferometric processing with two different modes, FBD HH images were doubly oversampled to FBS geometry with pixel spacing being 3.17 m in the azimuth direction and 4.68 m in the slant range direction, respectively. The three arc-second (spatial posting of 90 m) shuttle radar topography mission (SRTM) digital elevation model (DEM) [42] was collected for the topographic phase removal and InSAR-derived productions’ geocoding. Google earth images and ICESat/GLAS data were used for the underlying graphs of exhibition, analysis, and validation of the results.

Table 1. ALOS PALSAR images used in this study and their spatial/temporal baselines in comparison to the master image (20091110).

No.	Acquisition Time (yyyymmdd)	Polarization	Perpendicular Baseline (m)	Temporal Baseline (days)
1	20080622	HH	244.8	−506
2	20090207	HH	−1102.5	−276
3	20090625	HH/HV	−439.5	−138
4	20090810	HH/HV	−662.9	−92
5	20090925	HH/HV	−184.1	−46
6	20091110	HH	0	0
7	20091226	HH	316.2	46
8	20100210	HH	807.6	92
9	20100628	HH/HV	1157.7	230
10	20101229	HH/HV	1954.1	414
11	20110213	HH	2573.3	460

Table 1. ALOS PALSAR images used in this study and their spatial/temporal baselines in comparison to the master image (20091110).

No.	Acquisition Time (yyyymmdd)	Polarization	Perpendicular Baseline (m)	Temporal Baseline (days)
1	20080622	HH	244.8	−506
2	20090207	HH	−1102.5	−276
3	20090625	HH/HV	−439.5	−138
4	20090810	HH/HV	−662.9	−92
5	20090925	HH/HV	−184.1	−46
6	20091110	HH	0	0
7	20091226	HH	316.2	46
8	20100210	HH	807.6	92
9	20100628	HH/HV	1157.7	230
10	20101229	HH/HV	1954.1	414
11	20110213	HH	2573.3	460

3. Methodology

To monitor landslides over a large region of rough terrain and dense vegetation, an enhanced MTInSAR methodology was developed by taking advantage of existing approaches to overcome certain limitations. The entire data processing chain is shown in Figure 3. Basic interferometric processing (including SLC SAR image calibration, SAR image co-registration, multi-looking, interferogram generation, removal of topography and flat earth phase, as well as 3D phase unwrapping) was carried out using GAMMA software. For other processing requirements, including extraction of CS/DS, phase optimization, reforming interferograms, modeling and removal of orbit/topography related errors and deformation inversion, we developed appropriate programs implemented in MATLAB and C++.

Figure 3. Data processing flowchart of the enhanced MTInSAR methodology combining coherent scatterers (CS) and distributed scatterers (DS).

Figure 3. Data processing flowchart of the enhanced MTInSAR methodology combining coherent scatterers (CS) and distributed scatterers (DS).

3.1. Combination of CS and DS

Measurable targets are few and unevenly distributed over vegetated mountainous regions. Their scarcity hampers the application of prevalent MTInSAR techniques in slope instability monitoring. Generally, at least a density of five points per km² is needed to ensure reliable estimation of atmospheric delay phase and deformation parameters. PS points are few in the natural scenario of the study site; and the occurrence of CS points is also limited due to severe temporal decorrelation caused by vegetation cover. DSs reflect surfaces extending uniformly over a large area and where neighboring pixels of their images have similar backscattering characteristics. They are widespread particularly on slopes covered by low vegetation or debris.

To increase the point density per area, we synergistically analyzed CS and DS. A total of more than 1.8 million measurable points (CS plus DS) covering the study area was extracted except for unfavorable radar looking regions, e.g., layovers and shadows. Several studies had proposed identification and denoising algorithms for such unfavorable areas [14,43,44]. In this study, analogous approaches of SBAS and SqueeSAR were applied for the extraction of CS and DS points. To identify CS, a complex multi-look operation (two in range and six in azimuth) was first performed producing pixels with ground range dimensions of about 20 m × 20 m. Then we generated a total of 36 interferometric pairs with small orbital and temporal separations (spatial baseline threshold 1500 m and temporal threshold 365 days) (Figure 4). Coherence maps were estimated for each interferogram and pixels with mean values higher than the assumed threshold of 0.3 were picked out as CS. Results showed that their distribution is extremely uneven; density differences ranged from hundreds of points in the flat urban area to a single point in mountainous regions.

Figure 4. Spatial and temporal baselines of the interferograms. Circles represent SAR acquisitions.

Figure 4. Spatial and temporal baselines of the interferograms. Circles represent SAR acquisitions.

To extract DS, homogeneous pixels should be identified first, for which a Kolmogorov-Smirnov test was applied to the time series amplitude data. In the 20 × 20 (pixels) moving windows the pixels with more than 180 homogeneous neighbors were taken as DS candidates, and the number of points extracted was more than that for CS. Unlike PS and CS, no dominant point targets exist in the dimensions of DS pixels, and they usually maintain low coherence and signal-to-noise-ratio (SNR) in the interferograms. To improve the qualities of DS-derived parameter estimations, a sophisticated phase filtering processing is needed [43,45,46]. Ferretti et al. [14] used the maximum likelihood (ML) method to estimate the likeliest phases of each image to replace the original values as follows:

$$\widehat{\lambda }=\arg \underset{\lambda }{\min }\left\{{\text{Λ}}^H({\left|\widehat{C}\right|}^{-1}\circ \widehat{C})\text{Λ}\right\}$$

(1)

where $\widehat{\lambda }={\left[0,{\vartheta }_2,\dots {\vartheta }_N\right]}^T$ represent the new phase values of N images to be solved when the phases of first image are assumed to be zeros. $\text{Λ}=\text{exp}\left(i\text{λ}\right)$, λ is the original value of $\widehat{\lambda }$, $\widehat{C}$ is the complex coherence matrix based on the homogeneous pixels. To solve this equation, LBFGS (Limited memory Broyden-Fletcher-Goldfarb-Shannon) algorithm was adopted, which is a popular optimization algorithm in familiar quasi-Newton methods and has great advantages in dealing with large nonlinear problems due to its efficiency and lower consumption of computer memory [47,48]. The effect of phase optimization is shown in Figure 5.

Figure 5. Comparison of original and optimized phases. Results covering part of the study area are shown, superimposed on the averaged intensity imageries. (a)–(c) are the original phases of selected regions in our study area, and (d)–(f) are their corresponding optimized results. The interferometric pairs are indicated by the black fonts.

Figure 5. Comparison of original and optimized phases. Results covering part of the study area are shown, superimposed on the averaged intensity imageries. (a)–(c) are the original phases of selected regions in our study area, and (d)–(f) are their corresponding optimized results. The interferometric pairs are indicated by the black fonts.

3.2. Parameter Estimation by SBAS Solution

After synergistically analyzing CS and DS, unknown parameters (e.g., motion and residual heights) were estimated by the SBAS solution. Differential interferograms were available after the removal of topography and earth flat phases. Estimating the integer ambiguities of the wrapped phase, so-called phase unwrapping, is the most crucial step in deriving surface deformation estimates. Decorrelation and phase discontinuity in the order of 2π caused by the combination of low density of points, large topography and high displacement rates are mainly responsible for unwrapping errors. 3D unwrapping algorithm (including in temporal and spatial domains) [49] was applied owing to its tested performance for large amount of small-baseline interferograms. A reference point assumed to be stable (in the Dujiangyan city and away from landslides) was required to transform the differential phase to absolute value.

After the processing mentioned above, the deformation phase was still mixed up with orbit errors, atmospheric artifacts and residual topography errors. For the long-wavelength ALOS satellite, the orbit errors could reach 30 cm [50], and would seriously bias estimates of changes in topography. In the mountainous areas such errors could be amplified several folds. Atmospheric path delay can be divided into two parts: one arising from the turbulent mixing process and the other due to elevation-dependent stratified component [25]. The latter is usually misinterpreted as either topography or deformation [29,51]. Topography errors, related to perpendicular baselines, were firstly removed using the algorithm of SBAS. Then, we applied a combined model [52] comprising a biquadratic model for orbital phase errors [53,54] and a linear model for the elevation dependent errors (stratified atmospheric delay and topography errors) [55] to further reduce the artifacts from orbit, residual topography and atmospheric disturbance. The combined model is described as follow:

$$\text{ø}\left(x, y\right)\text{ = }{a}_0\text{ + }{a}_{1}x + a_{2}y \text{+ }{a}_{3}xy + a_4{x}^2\text{ + }{a}_5{y}^2\text{ + }{a}_{6}h\left(x, y\right)\text{ + ε}\left(x, y\right)$$

(2)

where $x$ and $y$ indicate the range and azimuth coordinates; $h$ is the elevation; $\text{ε}$ represents the random phase error; $a_i$ are coefficients to be estimated. After correctly unwrapping and removing phase errors mentioned above, a temporal and spatial filter was applied again. Finally, deformation rates and time series displacements were derived by the least squares solution.

3.3. Improvements and Advantages

The improvements and/or innovations of the proposed MTInSAR processing chain can be summarized as follows. First, taking advantage of SBAS and SqueeSAR, a landslide oriented MTInSAR chain was developed for monitoring slope instability in a vegetated, steep terrain including important habitats of the giant panda. A ground deformation field in the area was derived due to the MTInSAR data processing chain that we have developed. A large number of measurable points (including CS and DS) were extracted. The point density is an important factor that restricts the application of MTInSAR in suburban vegetated regions. By synergistically analyzing CS and DS, we extracted 1.8 million points (that is four times more than SBAS). The LBFGS algorithm was introduced to improve the efficiency of the procedure of phase optimization. The merit of SBAS was introduced to generate multi-master interferograms, overcoming the limitation of SqueeSAR applications for SAR images based on more than 20–25 scenes. We also introduced an enhanced model to remove elevation-dependent phase errors that makes our methodology more robust than classical MTInSAR approaches, particularly for studying mountainous terrains.

4. Results, Validation and Discussion

4.1. Application to Landslide Monitoring

There are many factors that need to be taken into account when using SAR interferometry to monitor landslides. Besides point selection, density, removal of the various errors and incomplete unwrapping mentioned above, the relationship between SAR imaging geometry and topography requires careful interpretation. DInSAR techniques only detect the change of the sensor-target distance, corresponding to the component of actual displacement projected onto LOS direction. This single-direction measuring mode limits the displacement sensibility to the plane approximately orthogonal to orbit direction [30]. ALOS satellite operates on a nearly polar orbit with an inclination angle of 98.2°, and the geometry of its right-looking ascending track enables easy monitoring of gentle slopes facing east; but deformations along the north-south direction become very difficult to detect. Based on the geometric relations, the regions of low sensibility can be extracted using the following equations (see Figure 6):

$$\text{Facing north:}\left(\text{β > θ/2}\right)\text{ \& [336.8° < δ < 360° || δ < 6.8°]}$$

(3)

$$\text{Facing south:}\left(\text{β > θ/2}\right)\text{ \& [156.8° < δ < 186.8°]}$$

(4)

where $\text{β}$ is slope gradient, $\text{θ}$ is LOS direction and $\text{δ}$ represents the slope aspect clockwise from the north. The slopes facing west and steep slopes facing east are severely affected by layover and shadow [56], and we extracted them with GEO module in the Gamma software (see Figure 6). In these regions SAR signal is useless. It is clear that low sensibility and layover/shadow are widespread, and account for about 12% and 9% of the area with slope gradients bigger than 10°, respectively. This means that assuming an even distribution of landslides at least one in five landslides would go undetected.

Figure 6. Layover/shadow and low sensitivity regions in our study area. Red rectangles, indexed from I to V indicate the regions where deformation results were illustrated in detail in Section 4.3. Routes of field campaigns after the Wenchuan earthquake are shown by white lines.

Figure 6. Layover/shadow and low sensitivity regions in our study area. Red rectangles, indexed from I to V indicate the regions where deformation results were illustrated in detail in Section 4.3. Routes of field campaigns after the Wenchuan earthquake are shown by white lines.

Additionally, the slope gradient and aspect affect the presence of measurable points, and about 70% of them are on slopes of gradients 10°–50° and aspects 30°–180° (northeast to south). Very few points are located on aspects of 180°–300° due to the effect of layover; and this influence becomes greater as the slope gradients increase (see Figure 7a).The incidence angle (about 34°) of PALSAR sensor determines that radar is more sensitive to the vertical than horizontal deformations [25]. More than three independent datasets (including ascending and descending) are needed to acquire the 3D displacement. For landslide monitoring, the LOS displacement is usually transferred to the down-slope direction assuming that the movement is translational along the slope, and in accordance with the equations used in [37]. In this paper, LOS displacement was preserved for landslide analysis.

Figure 7. (a) Distribution of detected points related to slope gradient and aspect; the right legend indicates the occurrence percentage of the total number of points; (b) point height differences between SRTM DEM and InSAR-derived elevations based on the ICESat/GLAS data; (c) InSAR-derived deformation rates on survey control points.

Figure 7. (a) Distribution of detected points related to slope gradient and aspect; the right legend indicates the occurrence percentage of the total number of points; (b) point height differences between SRTM DEM and InSAR-derived elevations based on the ICESat/GLAS data; (c) InSAR-derived deformation rates on survey control points.

4.2. Validation

Several methods were used to validate the InSAR-derived results. First, we used Google earth images as the base maps for deformation illustration and landslide identification, on which landslide regions are clearly visible with much brighter tones than the original surfaces. Our study demonstrates that the occurrence of unstable points is consistent with the distribution of large-scale landslides from the past. At a local scale, we found that slope instabilities often occur close to, or in areas with similar contours to previous landslides. Through these analyses we could make sure that the InSAR-derived results are validated. Several field surveys were undertaken after the Wenchuan earthquake with the specific routes marked by white lines in Figure 6. Though away from the main disaster regions, we still could see landslides everywhere and some of them even blocked the highways.

Additionally, ICESat/GLAS elevation data was used to check the results of height correction that is not only the auxiliary product of MTInSAR but also an important index of accuracy assessment. Especially in mountainous regions, like our study area, wrong height is the primary error source of deformations. The GLAS footprints on the ground are approximately 72 m in diameter and 172 m in spatial interval. Its accuracy of elevation is probably 14 cm [57], but the value becomes extremely unreliable in mountainous regions. Considering this problem and the distance between laser footprints and measureable InSAR points (smaller than 100 m), only 50 points were preserved. Figure 7b shows the height differences of SRTM DEM and InSAR-corrected heights based on ICESat/GLAS data. It is clear that the InSAR-corrected values are more consistent with the reference values (that is, 16.7 m averaged absolute difference with 9.4 m standard deviation of InSAR VS. 20.5 m and 14.2 m of SRTM DEM), indicating the validity of height correction obtained. There were seven survey control points (black triangles in Figure 6) in the study area, which were assumed to be absolutely stable. To compare it with the MTInSAR results, we calculated the mean deformation values of five adjacent InSAR measurable points. The results show that the mean absolute deformation rates of control points is 0.46 cm/a (see Figure 7c), which confirms the reliability of MTInSAR deformation results.

4.3. Deformation and Landslides Analysis

Among the InSAR products, the backscattering characteristics, coherence map and DEM are usually used to distinguish sudden and massive landslides [32,37], while the deformation map was good for identifying the potential, gradual, and longer-term slope failures. In this experiment, using MTInSAR-derived deformation we have detected dozens of landslides in an area of approximate 4200 km², mainly distributed along the Minjiang River and the Longmenshan Fault, as well as some high-altitude regions. The LOS deformation on important regions “I, II, III, IV and V” (see Figure 6) rather than the whole study area are shown in this paper for convenience, with velocity rates primarily ranging from −7 to 1.5 cm/a. The minus values mean that targets are moving away from the radar sensor, and considering the radar geometry and slope aspect, they indicate potential landslide events. A further analysis of regional deformation and landslides is shown below.

4.3.1. Minjiang River

Minjiang is an important left-bank tributary of the upper Yangtze River. Its 1279-km-long main stream flows in a southern direction. The watershed is larger than 130,000 km². Large volume and an over-3-km drop of the water flow has made possible tens of hydropower stations, including the ancient Dujiangyan irrigation system (another World Heritage site) and the modern Zipingpu reservoir. In the upper reaches of Minjiang River, just before the Yingxiu Town, earthquakes and landslides are frequent. Lava, clasolite, and carbonatite dominate the lithology in this region. Translational landslides, i.e., mass movements with little rotation or backward tilting [1], occurred on almost every slope in the river valley, resulting in considerable surface fragmentation and serious deterioration of vegetation. Due to the geometric distortion (layover on the east bank, shadow on the west bank) and decorrelation caused by drastic displacements, the measurable points are sparse, and mainly concentrated on the unbroken surface and/or the landslide regions where relatively stable conditions had been maintained during the period for which SAR images were obtained. Through the MTInSAR-derived deformation map (see Figure 8a), we could identify dozens of potential landslides along the 40-km-long valley, such as those shown in Figure 8b–d. A large number of landslides of various magnitudes develop at altitudes from 2000 down to 1200 m. Small landslides tend to converge, generating large landslides. In the 5-km-long valley shown in Figure 8b, most serious landslides occurred at low altitudes close to the river, covering more than half the surface, with the largest size up to 400 m wide and 700 m long. Mean LOS deformation rate with approximately −3 cm/a pointed to the instability of this region and the recurrent possibility of landslides. The whole surface of a slope (southeast aspect with approximately 40° average gradient, Figure 8c), from the mountain-top to the foot, had been cut into pieces by landslides, with maximum lengths up to 1 km. Deformation results reveal that the points lying just below the landslides are more unstable than those close to the river valley. The slope in Figure 8d is located near Yinxing Town, facing east and characterized by the mean gradient of approximately 50°. It can be divided into two parts, both of which have been severely damaged by falling stones. The lower boundaries (black lines shown in Figure 8d) of previous landslides are exactly parallel to the boundary of InSAR-derived unstable slopes. Mean deformation rate of approximately −5 cm/a indicated that previous landslides here were still active, and massive movement could occur continuously, because large numbers of stones and lots of mud had been accumulated after the previous landslide events, reactivating or triggering new translational landslides because of the downward loading of loose overlying mantles. We chose three profiles in representative regions for the deformation analysis (EE’, FF’ and GG’, see Figure 8d), as illustrated in Figure 9. Potential landslides can be clearly distinguished through deformation differences, such as sections marked by dotted red-lines along profiles EE’, FF’ and GG’, respectively. In addition, a deformation zone with values ranging from −2 to −4 cm/a in the lower-left of Figure 8d was also observed and corresponds to areas impacted by engineering and construction activities such as roads and houses. Figure 8 shows the displacement time series of three typical points from unstable slopes with velocity rates up to −4.3, −1.7 and −3.4 cm/a respectively. It is clear that there is a deceleration trend due to the degraded effect of Wenchuan earthquake after 2009. Hence, Minjiang River valley is the region most affected by the landslides, which can be attributed to erosion due to the river, fragile lithology and steep topography. Fortunately, most residential areas are located at the gentle slopes away from the possible landslides. But infrastructure in the valley floor is vulnerable to landslides.

Figure 8. Deformation maps along Minjiang River (region I); (a) whole map superimposed on the averaged SAR intensity imagery; (b), (c) and (d) are images of three typical landslide regions superimposed on Google earth images, with their locations indicated by 1, 2, and 3 in (a) respectively; temporal displacement evolution of the points A, B, and C (indicated in Figure b, c and d respectively); The blue line in (a) is the Minjiang River; In (d) black solid lines represent the lower boundaries of previous landslides, and dotted blue lines (EE’, FF’ and GG’) indicate the profiles of deformation analysis (see also Figure 9).

Figure 8. Deformation maps along Minjiang River (region I); (a) whole map superimposed on the averaged SAR intensity imagery; (b), (c) and (d) are images of three typical landslide regions superimposed on Google earth images, with their locations indicated by 1, 2, and 3 in (a) respectively; temporal displacement evolution of the points A, B, and C (indicated in Figure b, c and d respectively); The blue line in (a) is the Minjiang River; In (d) black solid lines represent the lower boundaries of previous landslides, and dotted blue lines (EE’, FF’ and GG’) indicate the profiles of deformation analysis (see also Figure 9).

Figure 9. (a)–(c) illustrate the variation in point LOS deformation rates (in cm/a) along profiles EE’, FF’ and GG’ (in Figure 8d), respectively at their corresponding altitudes (in m). The sections marked by dotted red lines indicate the boundaries of the potential landslide regions.

Figure 9. (a)–(c) illustrate the variation in point LOS deformation rates (in cm/a) along profiles EE’, FF’ and GG’ (in Figure 8d), respectively at their corresponding altitudes (in m). The sections marked by dotted red lines indicate the boundaries of the potential landslide regions.

4.3.2. Longmenshan Fault

The Longmenshan Fault comprises three parts; where the Wenchuan earthquake caused no ruptures (Mao-Wen), major ruptures (Yingxiu-Beichuan) and minor ruptures (Guanxian-Anxian), respectively [31]. The dominant lithology is demarcated by the faults, mainly comprising lava, clasolite, limestone, malmstone, and mudstone. In the segment from Wenchuan to Caopo Town, the Mao-Wen Fault is undercut by the Minjiang River, rendering rocks in the Fault exposed and prone to erosion. Figure 10 shows the region southwest of Caopo, where few landslides occurred. The distribution characteristics of previous landslides here is discontinuous, small scale and low density, probably due to the occurrence of shallow valley landscapes. On the east side, the points are sparse due to the radar layover effect, while on the west side the gentle slopes are facing east and good results were obtained. The deformation result reveals that there are two extremely unstable valleys (see Figure 10a, region (1) and (2)), each a few-kilometers long and where average rates exceed −5 cm/a. More than ten villages are located at the valley in Figure 10b, with the maximum slope gradient exceeding 60°. In the 6 km-long valley, dozens of landslides occurred with lengths of individual slides varying from several hundred meters to over 1 km.

The valley presented in Figure 10c is about 8 km long with the altitudes ranging from 1600 to 2400 m and slope gradients from 20° to 50°. Villages named Dengcaoping, Changheba, and Zhangyagang are located here. This place used to be safe and suitable for human habitation, and very few landslide imprints were found in Google earth images. Because of this, a large number of houses have been built in the flat valley floor after 2008 for the homeless people displaced by the Wenchuan earthquake. However, deformation results showed that the slopes were becoming unstable, especially in the region with low altitude and close to places where engineering activities were concentrated. The appearance of the instability can be attributed to the consequences of human activities. Figure 10d shows a slope near Sanshenghao Village, facing south, with a gradient of about 40°, and it is divided by a deep groove. On the west part, the points are situated above the previous landslides. The mean deformation rate is approximately −2.3 cm/a, and the absolute values increase at lower altitudes indicating that regions closer to landslides were more unstable. The same phenomenon can be seen in the east; the unstable points are mostly concentrated around landslides rather than on sections above them. It is also worth noting that the deformation on this slope is significantly underestimated because radar is not sensitive to the component along the north-south direction.

Figure 10. Deformation maps along Mao-Wen Fault (region II). (a) The whole map superimposed on the averaged SAR intensity imagery; (b)–(d) show three typical landslide regions superimposed on Google earth images, with their locations indicated by 1, 2, and 3 in (a) respectively; Blue lines in (d) indicate the regions where previous landslides were concentrated.

Figure 10. Deformation maps along Mao-Wen Fault (region II). (a) The whole map superimposed on the averaged SAR intensity imagery; (b)–(d) show three typical landslide regions superimposed on Google earth images, with their locations indicated by 1, 2, and 3 in (a) respectively; Blue lines in (d) indicate the regions where previous landslides were concentrated.

The deformation map of region III, including the town of Yingxiu and the Zipingpu reservoir, was shown in Figure 11. Yingxiu-Beichuan Fault and Zipingpu reservoir are linked to the epicenter of the Wenchuan earthquake, and this region is characterized by seismically-induced giant landslides, often more than 2 km in length and approximately 800 m in width. The landslide distribution varies greatly in space; much less occurred on the southeast of Yingxiu due to the gentle gradient. Generally, the frequencies of landslides decrease as altitudes reduce. From the deformation results it is clear that some regions are still unstable. The most significant region is Yingxiu (see Figure 11b), which was the most vulnerable area in 2008. Few landslide imprints as well as large deformations can be found on adjacent slopes. Numerous points of negative deformations are observed at the valley floor, where a large number of disaster relief-houses were constructed. Considering that this region is gentle and close to the river, and that the houses were built in a hurry, the negative impacts could probably be attributed to engineering induced subsidence rather than landslides. Figure 11e shows the deformation variation along a profile through Yingxiu with subsiding regions highlighted by the red ellipses. Translational landslides were prevalent in areas surrounding the Zipingpu reservoir, particularly on the steep north shore where massive seismically-induced landslides were clearly visible. Rockfalls, i.e., abrupt movements of rocks and boulders due to the presence of interstitial water from precipitation, were also widespread. A comparison between the Zipingpu Dam, a stable engineering construct (11c) and a typical landslide body on the south of the reservoir (Figure 11d) was also undertaken.

Figure 11. Deformation maps along Yingxiu-Beichuan and Guanxian-Anxian Fault (region III). Figure (a) the whole map superimposed on the averaged SAR intensity imagery; Figures (b)–(d) show the typical regions superimposed on Google earth images, with their locations indicated by 1, 2, and 3 in (a) respectively; In (d) the red ellipse indicates a landslide area; Figure (e) shows the variation in point LOS deformation rates (in cm/a) along profile AA’, and red ellipses indicate the subsiding regions.

Figure 11. Deformation maps along Yingxiu-Beichuan and Guanxian-Anxian Fault (region III). Figure (a) the whole map superimposed on the averaged SAR intensity imagery; Figures (b)–(d) show the typical regions superimposed on Google earth images, with their locations indicated by 1, 2, and 3 in (a) respectively; In (d) the red ellipse indicates a landslide area; Figure (e) shows the variation in point LOS deformation rates (in cm/a) along profile AA’, and red ellipses indicate the subsiding regions.

4.3.3. High-Altitude Mountains

Some potential landslides were also found in high-altitude mountain region, such as those shown in Figure 12. The region in Figure 12a is located about 9 km southeast of Wenchuan County, and the terrain is rather rough with the altitudes up to 4500 m and perennial snow cover in shady areas. Many slopes have been broken up due to landslides of various sizes, especially at the mountain-tops where most landslides originated. From the deformation results we can clearly discriminate the unstable regions, which were mainly located at the upper section of the slopes. This is in accordance with the distribution of the previous landslides, and it means that potential landslides are originating from mountain-tops. Figure 12c shows a mountain-top with an elevation of 4000 m, and giant landslides were evident on all the slopes with various orientations. The deformation results showed that this large region was characterized by a rate of about −3 cm/a but included areas with values exceeding −5 cm/a. The occurrence of unstable points implies that the previous landslides were still active. Figure 12d shows another mountain-top with imprints of previous landslides clearly visible on Google earth images. Unstable points were detected on the east side and no points laying on the west side due to the layover.

Figure 12. Deformation maps of region IV (a) and region V (b) superimposed on the averaged SAR intensity imageries; (c–f) showing the typical landslide regions superimposed on Google earth images; with their locations indicated by 1 and 2 in (a); 3 and 4 in (b); respectively. The deformation time series of a point marked by the red circle is illustrated in (e); indicating an acceleration motion trend since 2010. Blue lines in (d) and (e) indicate the regions where previous landslides were concentrated.

Figure 12. Deformation maps of region IV (a) and region V (b) superimposed on the averaged SAR intensity imageries; (c–f) showing the typical landslide regions superimposed on Google earth images; with their locations indicated by 1 and 2 in (a); 3 and 4 in (b); respectively. The deformation time series of a point marked by the red circle is illustrated in (e); indicating an acceleration motion trend since 2010. Blue lines in (d) and (e) indicate the regions where previous landslides were concentrated.

The region in Figure 12b has a similar topography and landscape condition as that in Figure 12a, but no massive landslides were detected there. Additionally, the deformation results demonstrate that this region is still stable except for a few small slopes, such as those present in Figure 12e,f. The unstable slope in Figure 12e faces southeast and would at times be covered by a thin layer of snow during winters. A few small landslides measuring dozens of meters occurred under those unstable points. Note that the thin snow cover (e.g., less than 10 cm) would not affect the deformation monitoring for the slope considering the penetration capability of l-band SAR waves. However, the melt-water from snow could accelerate the erosion of steep slopes, increasing slope instabilities which show an acceleration motion trend since 2010 as measured on a representative point (marked by the red circle) in Figure 12e. In contrast, high elevation, north facing slopes in Figure 12f were covered by thick snow all year-long, and several hundred-meters-long landslide imprints appeared on the snow surface, indicating the potential of future landslides in the region.

4.3.4. Giant Panda Habitat

In Sichuan province the giant panda habitat is affected by human activities like logging, grazing, construction of roads, poaching as well as natural disasters such as earthquakes and landslides. Past landslides affected 1.06% of the total giant panda habitat up to 2006 [58]. They were mainly small-to-medium landslides, resulting in the reduction of habitat quality and an occasional death of individual animals. The Wenchuan earthquake resulted in the disappearance of 5.95% of the total area of the giant panda habitat (65,584 hm²) [59]. The habitat was also more fragmented because movement corridors and access to water-sources were cut-off by landslides [60]. Note that these negative impacts would persist over the short-to medium-term.

Figure 13 shows the circumstances of the giant panda habitat and landslide distribution in the study area. Two major giant panda habitats, the Minshan Mountain, and Qionglai Mountain, account for 41.66% and 26.47% of total habitat areas, and their densities of giant panda population are 0.074 and 0.072 per km², respectively. Based on the deformation results, we extracted approximately 50 potential landslides, among which 20 were located within the Giant Panda Habitat World Heritage site, and all of them were confirmed by Google earth images. Most of them were distributed along the river and fault valley, bringing in little effect to the giant panda because pandas normally prefer higher altitudes. Compared with the landslide occurrence in the Minjiang River valley (defined as the reference layer with extremely-dense landslides), the density of landslides in giant panda habitats can be qualitatively divided into five grades (from rare to extremely dense) by analyzing the spatial distribution of past landslide via Google earth images and InSAR-derived results. Statistical analysis of the distribution of the five grades of landslides showed that: (1) the landslides in Minshan Mountain were much denser than that in Qionglai Mountain, indicating that the giant pandas in Minshan were more vulnerable to impacts of landslides; (2) there was a correlation between the landslide density and the altitude; density of landslides decline as altitudes decrease. However, it is fortunate that most giant pandas frequented areas where sparse or medium-density landslides occurred.

In brief, landslide analysis by MTInSAR in this study could facilitate the effective management of the giant panda habitat, particularly in identifying priority areas for surveys following an earthquake and for designing and implementing measures to mitigate excessive habitat fragmentation and facilitate panda movements during post-seismic periods.

Figure 13. Giant panda habitats in Minshan/Qionglai Mountains and landslide distribution derived by analyzing Google earth images and InSAR-derived results.

Figure 13. Giant panda habitats in Minshan/Qionglai Mountains and landslide distribution derived by analyzing Google earth images and InSAR-derived results.

5. Conclusions

In this paper we report findings of studies on landslides over a large region including in the Giant Panda Habitat World Heritage site in Sichuan China. MTInSAR techniques used in this study are powerful in monitoring subtle deformations as well as the landslides. Other approaches have limitations when the monitored site is large and complex as the one in our study and the SAR dataset available is small. An enhanced landslide-oriented MTInSAR approach described in this paper could be effective for landslide monitoring in Sichuan including in the giant panda habitats by an integrated use of CS and DS points and the analyses of ensuing results. Dozens of potential landslides were detected in detail in a large area of approximately 4200 km², distributed mainly along the Minjiang River and Longmenshan Fault, as well as in some high-altitude mountain regions at and above 4500 m. Considering their capability to detect and monitor potential landslides, MTInSAR approaches have value for hazard prediction and risk mitigation. Our results confirm that places of past landslides would continue to be geohazards in the short-to-medium term. New slides could easily be triggered in such locations. Even where landslides were small there is the risk of several small landslides aggregating to create a massive one. The impact of landslides on the giant panda habitat however, was low-to-moderate; nevertheless impacts of future landslides on giant-panda habitats may be of a different scale, and hence continuous monitoring and regular assessments are encouraged.