Multi-Temporal Loess Landslide Inventory Mapping with C-, X- and L-Band SAR Datasets—A Case Study of Heifangtai Loess Landslides, China

Abstract

The Interferometric Synthetic Aperture Radar (InSAR) technique is a well-developed remote sensing tool which has been widely used in the investigation of landslides. Average deformation rates are calculated by weighted averaging (stacking) of the interferograms to detect small-scale loess landslides. Heifangtai loess terrace, Gansu province China, is taken as a test area. Aiming to generate multi-temporal landslide inventory maps and to analyze the landslide evolution features from December 2006 to November 2017, a large number of Synthetic Aperture Radar (SAR) datasets acquired by L-band ascending ALOS/PALSAR, L-band ascending and descending ALOS/PALSAR-2, X-band ascending and descending TerraSAR-X and C-band descending Sentinel-1A/B images covering different evolution stages of Heifangtai terrace are fully exploited. Firstly, the surface deformation of Heifangtai terrace is calculated for independent SAR data using the InSAR technique. Subsequently, InSAR-derived deformation maps, SAR intensity images and a DEM gradient map are jointly used to detect potential loess landslides by setting the appropriate thresholds. More than 40 active loess landslides are identified and mapped. The accuracy of the landslide identification results is verified by comparison with published literatures, the results of geological field surveys and remote sensing images. Furthermore, the spatiotemporal evolution characteristics of the landslides during the last 11 years are revealed for the first time. Finally, strengths and limitations of different wavelength SAR data, and the effects of track direction, geometric distortions of SAR images and the differences in local incidence angle between two adjacent satellite tracks in terms of small-scale loess landslides identification, are analyzed and summarized, and some suggestions are given to guide the future identification of small-scale loess landslides with the InSAR technique.

1. Introduction

Loess occupies an area of about 630,000 km², equivalent to approximately 6.6% of the total area of China [1]. This includes the loess plateau, an important base of energy sources and chemical industries, which amounts to approximately 0.4 million km² in area [2]. Since loess has distinctive physical and mechanical properties—including loose texture, vertical joints, macro-pores, high water sensitivity, and collapsibility—it is prone to geological hazards triggered by earthquakes, rainfall and irrigation [3,4,5]. One third of the geological hazards in China occur in loess plateau, e.g., landslides, collapse, ground fissures and erosion and collapsible settlement [1,6]. Over the past few decades, these geological hazards have caused enormous numbers of casualties and economic losses, such as the destruction of civil houses and factories, the destruction of major transport trunk lines, damage to oil and gas routes and the reduction of cultivated lands [7,8]. Loess landslides, one of the severe geological hazards in loess plateau, have the special characteristics of high speed, a small spatial scale, long runout distance, group occurrence, and recurrence [2,9].

Landslide inventory maps record the location and (where known) the types of mass movements and the date of occurrence that have left discernable traces in an area. It is of great significance to investigate the distribution, type, recurrence rate and statistics of landslides, and to determine their susceptibility, deadliness, vulnerability and risk [10]. Traditional methods for landslide inventory mapping mainly depend on field geological surveys, visual interpretation of stereoscopic aerial images and unmanned aerial vehicle (UAV) photogrammetry. For instance, Peng et al. [11] obtained the distribution and type of the landslides in the Heitai area of China using UAV photogrammetry and field investigations. Furthermore, Zhuang et al. [5] mapped the locations of landslides at the scale of 1:50,000 in the Shaanxi province of China via field investigations and aerial photographs. Due to the labor-intensiveness and low efficiency of field investigation method, and low accuracy of UAV DEM measurement, the inventory mapping of large area, small-scale loess landslides is still challenging. Interferometric Synthetic Aperture Radar (InSAR) can overcome the aforementioned limitations as it can remotely sense landslides over as large as over tens to ten thousands of square kilometers area with good measurement accuracy (e.g., [12,13,14]).

Over the past decades, a series of advanced multi-temporal InSAR (MT-InSAR) algorithms, including Persistent Scatterers Interferometry (PSI) [15,16,17,18,19,20,21,22] and the Small Baseline Subset (SBAS) [23,24,25,26], have been developed to overcome the major limitations of the traditional InSAR method and hence improve the applicability of the InSAR technique. These advanced MT-InSAR techniques have been widely used in the identification of landslides [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]. However, due to the special characteristics of loess landslides, which include small spatial scale and abrupt occurrence and recurrence [44], especially for multi-temporal landslide inventory mapping, these advanced MT-InSAR techniques still have some difficulties, including the decorrelation caused by surface vegetation and soil moisture changes, unwrapping errors and DEM errors caused by abrupt landslides, and (often) atmospheric artifacts (including stratified atmosphere delay and turbulence).

Owing to the increasing availability of Synthetic Aperture Radar (SAR) datasets, both archived and operational, the combination of multi-source, multi-resolution and multi-temporal SAR data for landslide identification should constitute the present strategy. Herrera et al. [45] used multi-sensor and multi-temporal SAR data to detect and monitor slow-moving landslides in the Upper Tena Valley of the Central Spanish Pyrenees. Mulas et al. [46] integrated multi-sensor SAR data to characterize the slow-moving landslide in Badia village, up-stream of Corvara in the Dolomites region of Italy. Sun et al. [47] integrated multi-sensor SAR data to detect and characterize the slow-moving landslides in Zhouqu, China. Multi-sensor InSAR datasets including L-band ALOS/PALSAR and ALOS/PALSAR-2 data, X-band TerraSAR-X data and C-band Sentinel-1A/B data with various spatial-resolutions, revisiting periods and wavelengths, make it possible to detect small-scale loess landslides and generate multi-temporal inventory maps.

2. The Study Area

2.1. Geological Setting

The Heifangtai loess terrace is located in Yongjing county, Gansu province, northwestern China (Figure 1a), and consists of Heitai and Fangtai separated by Hulang gully, occupying a surface area of about 13.7 km² [9,11,48,49]. The Yellow River and the Huangshui River surround this terrace. The elevation difference from the bottom to the top of the terrace is about 199.1 m, and the slope gradient at the edge of the terrace varies from 20 to 70° (Figure 1b). Based on borehole drilling and geophysical exploration data, the terrace has been determined to be mainly composed of Aeolian loess [50]. Figure 2 shows the engineering geological profile of the terrace along the line A–A’ in Figure 1a in the azimuth of 60° direction, which can be divided into four layers from up to the bottom of Malan loess ($Q_3^{eol}$), alluvial clay ($Q_3^{al}$), fluvial gravel ($Q_3^{al}$) and bedrock ($K_{1}hk$) [11,50]. The top layer is Malan loess with a thickness of 30–50 m, which is mainly composed of silt-sized loess with high porosity and open-paced fabric, resulting in a high susceptibility to collapse upon wetting [11].

The average annual evaporation and precipitation in the study area are about 1593 mm and 287.6 mm, respectively [1]. According to Xu et al. [48], the annual irrigation volume was about 0.7 million cubic meters in the 1980s, and increased to 0.88 million cubic meters in 2009. As a result, the groundwater level has risen by approximately 20 m at an average rate of 0.18 m per year since the 1980s [11]. The basal zone of loess is softened by saturation and prone to deformation under overburden stress [49], which has led to about 200 loess landslides and a number of dense cracks distributed at the edge of the terrace [11]. Additionally, the terrace cliff retreated by about 126 m due to the frequent occurrence of landslides between 1977 and 2010, and is still retreating at an estimated rate of 0.024 km² per year [1,49]. As such, these loess landslides seriously restrict the development of the regional economy and threaten the safety of the local residents, factors which have been extensively investigated by numerous researchers in China and abroad [52,53,54,55,56].

2.2. Landslides

The distribution of landslides in the whole Heifangtai terrace are investigated based on digital elevation model (DEM), aerial images and the results of field work [11,51], as shown in Figure 1a. A total of 69 landslides are detected in the Heitai area. These landslides are classified into two groups [11]: Loess-bedrock landslide (failure surface generates in the loess layer and then extends downwards to the bedrock layer) and loess landslide (failure surface is entirely developed in the whole loess layer). The landslides can be further divided into nine landslide groups in the spatial domain [11]: G1 (Fangtai landslide group), G2 (Xinyuan landslide group), G3 (Dangchuan landslide group), G4 (Huangci landslide group), G5 (Yehugou landslide group), G6 (Jiaojiaya landslide group), G7 (Jiaojia landslide group), G8 (Chenjia landslide group) and G9 (Moshigou landslide group). Loess landslides are developed on the Fangtai landslide group, Dangchuan landslide group, Yehugou landslide group, Jiaojiaya landslide group, Jiaojia landslide group, Chenjia landslide group and Moshigou landslide group, while loess-bedrock landslides are distributed on the Fangtai landslide group, Xinyuan landslide group, Huangci landslide group and Jiaojia landslide group. Detailed information on landslide distribution, name, length, width and volume is provided in Ref. [11]. Over the past three decades, more than 120 failures have occurred in the study area [1]. Recently, a large landslide occurred on the Dangchuan landslide group on 29 April 2015, which released about 35 × 10⁴ m³ of debris with a maximum sliding distance of 182 m, and destroyed three factories and fourteen houses. Three large landslides subsequently occurred on the Dangchuan landslide group on 1 October 2017 with a total volume of about 33.9 × 10⁴ m³ [11,57].

3. Datasets and Methodology

3.1. Datasets

In order to determine the distribution of landslides in the Heifangtai loess terrace and reveal the temporal evolution of landslides over the last 11 years, a total of 111 SAR images from four different satellites and seven different tracks were used in the study (Figure 3). The main parameters of the SAR images used in the study are shown in

Table 1. Main parameters of the SAR datasets used for Heifangtai landslide inventory mapping.

Satellite	Band	Track	Geometry	Incidence Angle (°)	Resolution in Azimuth (m) and Range (m)	Number of SAR Images
ALOS/PALSAR	L	473	Ascending	38.7385	3.16 × 4.7	22
ALOS/PALSAR	L	474	Ascending	38.7263	3.14 × 4.7	12
ALOS/PALSAR-2	L	146	Ascending	40.5539	3.25 × 4.29	6
ALOS/PALSAR-2	L	39	Descending	40.5551	3.25 × 4.29	5
TerraSAR-X	X	21	Ascending	41.1669	1.26 × 0.91	23
TerraSAR-X	X	165	Descending	41.8010	1.26 × 0.91	19
Sentinel-1A/B	C	135	Descending	33.7927	9.32 × 13.97	24

. The periods of each dataset are presented in Figure 4. Note that the periods of ascending and descending TerraSAR-X data are less than one year, while those of ALOS/PALSAR and ALOS/PALSAR-2 images are around four and three years, respectively. Unfortunately, no achieved abovementioned SAR images are available during the years of 2011 to 2014.

In order to obtain high-accuracy InSAR measurement results, for ALOS/PALSAR-2, TerraSAR-X and Sentinel-1A/B datasets, a three-meter spatial resolution digital elevation model (DEM) (hereinafter UAV DEM) generated by unmanned aerial vehicle (UAV) photogrammetry with an optical camera with the type of Sony ILCE-7R_FE35mmF2.8ZA_35.0_7360x4912 (RGB) is used to differentiate the topographic phase and to perform the result analysis. For ALOS/PALSAR datasets, one arc-second advanced land observing satellites (ALOS) global digital surface model (AW3D30 DSM) data is adopted for topographic phase removal and result analysis. A multilooking factor of two in azimuth is used in the processing of ALOS/PALSAR and ALOS/PALSAR-2 data, while the multilooking factors of 2 × 2 and 4 × 1 are used in the processing of TerraSAR-X data and Sentinel-1A/B data, respectively. The spatial resolution of multi-looked ALOS/PALSAR images is about 7.5 m in both range and azimuth directions, that of TerraSAR-X images is about 3 m in both directions, that of ALOS/PALSAR-2 images is about 6.5 m in both directions, and that of Sentinel-1A/B images is about 16 m in both directions. SAR images with these spatial resolutions have the best ability to detect small-scale loess landslides and to monitor landslide deformation with large gradients.

The SBAS method [23] is adopted independently for different SAR data to generate all possible interferograms firstly by setting the appropriate temporal and spatial baseline thresholds. Then, the coherence of interferograms are checked, especially in landslide-prone regions. The interferograms with heavy decorrelation noises caused by agricultural irrigation, heavy rainfall in summer and dense vegetation are removed. Finally, the following interferograms are selected: A total of 38 interferograms for ALOS/PALSAR images from track 473, 19 interferograms for ALOS/PALSAR images from track 474, 30 interferograms for ascending TerraSAR-X images, 25 interferograms for descending TerraSAR-X images, 8 interferograms for ascending ALOS/PALSAR-2 images, 7 interferograms for descending ALOS/PALSAR-2 images and 44 interferograms for descending Sentinel-1A/B images. Spatial–temporal baselines and high-quality InSAR combinations are shown in Figure 5. Figure 5a shows the configuration of the temporal and spatial baseline of the interferometric pairs produced by ascending ALOS/PALSAR datasets from tracks 473 (black line) and 474 (red line). Figure 5b shows the configuration of the temporal and spatial baseline of the interferometric pairs produced by ascending and descending ALOS/PALSAR-2 datasets from tracks 146 (black line) and 39 (red line). Figure 5c shows the configuration of the temporal and spatial baseline of the interferometric pairs produced by ascending and descending TerraSAR-X datasets from tracks 21 (black line) and 165 (red line). Figure 5d shows the configuration of the temporal and spatial baseline of the interferometric pairs produced by descending Sentinel-1A/B datasets from track 135. The X-band SAR data has weaker penetrability than the C- and L-band data, especially in vegetation-covered regions [58], and accordingly the interferograms generated from both ascending and descending TerraSAR-X datasets are divided into two independent subsets. The detailed data processing parameters used in the study are listed in

Table 2. The detailed data processing parameters used in the study.

Datasets	Track	Multilooking Factor	Temporal Threshold (day)	Baseline Threshold (m)	Interferograms Used
ALOS/PALSAR	473	1 × 2	500	2000	38
ALOS/PALSAR	474	1 × 2	500	2000	19
ALOS/PALSAR-2	146	1 × 2	600	500	8
ALOS/PALSAR-2	39	1 × 2	600	500	7
TerraSAR-X	21	2 × 2	60	250	30
TerraSAR-X	165	2 × 2	60	250	25
Sentinel-1A/B	135	4 × 1	48	200	44

.

3.2. Flowchart of Landslide Identification with the InSAR Method

InSAR images and derived products have proven useful for mapping and characterizing landscape changes and ground deformation [58]. In this study, a topographic slope map, InSAR deformation maps and SAR intensity images are jointly employed to detect active small-scale loess landslides.

The detailed flowchart of small-scale loess landslide detection based on multi-source and multi-temporal SAR datasets is depicted in Figure 6. For one specific SAR dataset, a suitable master image is first selected considering both the temporal and spatial baseline and the variation of the Doppler central frequency [59], and all slave images are subsequently co-registered to it. Secondly, all possible interferograms are generated by setting appropriate temporal and spatial baseline thresholds, the topographic phase is removed with respect to external DEM, interferograms are filtered [60] and unwrapped with the minimum cost flow (MCF) method [61], the component of stratified atmosphere is removed by estimating a phase delay elevation profile for each interferogram and the residual orbital ramp is removed using refined baselines. Lastly, the average deformation rate for each pixel is calculated by weighted averaging (stacking) of the high-quality interferograms [62]. It should be noted that three key steps—the selection of high-quality interferograms, DEM error correction and phase unwrapping error detection and correction—must be carefully considered in order to detect small-scale loess landslide.

The derived average deformation rate map is used to detect active landslides. It is initially regarded as an active landslide if the annual average deformation rate is greater than 2 cm or less than −1 cm in the line-of-sight direction. In order to exclude the influence of layover for landslide identification, DEM gradient map and SAR intensity images are used to further validate the identified active landslides [27]. The value of the DEM gradient threshold is set as 10° based on in situ geological investigation. It is identified to be an active landslide if the DEM gradient value is greater than 10°. The coordinates and shapes of active landslides are ascertained. The same procedure is conducted for SAR data from other satellite tracks. The independent InSAR results from different satellite tracks are cross-validated with respect to their geometries. Eventually, the identified landslides from different SAR datasets are combined into a mosaic to generate the final multi-temporal landslide inventory maps.

3.3. DEM Error Correction

A newly generated, high-resolution UAV DEM is used to eliminate the topographic phase for landslide identification. However, due to the surface change caused by repeated landslides, no DEM is available during the SAR acquisition periods. DEM errors often exist in differential interferograms, and will inevitably result in errors in the final deformation results. According to the height-to-phase formula [63], the topographic phase is proportional to the perpendicular baseline $B_{\perp }$ between two SAR images, and the DEM error $\Delta Z$ can be estimated with following equation:

$$\delta \varphi =\left[A,B\right]\left[\begin{array}{l}V\\ \Delta Z\end{array}\right]+\delta {\varphi }_n$$

(1)

where $\delta \varphi$ is the unwrapped phase, $A$ and $B$ are coefficient matrices (for the detailed form of $A$ and $B$, see [27]), $V$ is the deformation rate and $\delta {\varphi }_n$ is the residual phase including nonlinear deformation phase, turbulent atmospheric noise and phase noise. As the phase of DEM error is sensitive to the long perpendicular baseline, interferograms with small temporal baselines and long spatial baselines are selected to calculate DEM error under the norm of least squares [2]. Moreover, high quality interferograms without unwrapping errors and atmospheric perturbations are needed. DEM errors are then added back to the original DEM to update the DEM as follows:

$$Z_{topo}=Z_{DEM}+\Delta Z$$

(2)

where $Z_{DEM}$ is the height of the DEM and $Z_{topo}$ is the updated DEM. The value of interferogram coherence is used to determine whether the DEM error is reasonably corrected. The abovementioned procedure can be operated iteratively.

3.4. Phase Unwrapping Error Detection and Correction

The principal observations of InSAR are wrapped interferometric phases (modulo 2π) and the 2π ambiguity is removed in practice by phase unwrapping step to generate a continuous map of phase values. However, in the case of loess landslide identification, phase unwrapping errors often occur due to the large fringe numbers presented on the interferograms caused by large deformation gradient and coherence loss.

Some methods have been developed to evaluate phase unwrapping quality and correct phase unwrapping error. For instance, López-Quiroz et al. [64] used the root mean square (RMS) of interferometric system disclosure to detect the phase unwrapping error, while Yang et al. [65] used the region-grow method to detect and correct phase unwrapping error. In this study, unwrapping errors are firstly detected and located by the vision and closed loop methods [66]. Two different methods are then applied to correct the detected unwrapping errors: (1) When conducting the phase unwrapping, the correct unwrapping phase can be obtained by testing different coherence thresholds; (2) As the phase contributions (deformation, atmospheric and orbital) behave in a conservative manner, e.g., ${\varphi }_{AB}+{\varphi }_{BC}={\varphi }_{AC}$, where ${\varphi }_{AB}$ is the phase contribution of interferogram $AB$ generated by acquisitions $A$ and $B$, once there exits disclosure, unwrapping errors occur. The 2π integer times phase offset is then added back to the wrongly unwrapped phase to achieve the correct unwrapped phase.

4. Results

4.1. Error Detection and Correction

4.1.1. DEM Error Correction

The UAV DEM was acquired over the study region in 2015. Some landslides occurred in the following years, which caused topography change. Hence, DEM errors occur in differential interferograms generated from ALOS/PALSAR-2, TerraSAR-X and Sentinel-1A/B datasets. As an example, two typical interferograms with short perpendicular baseline and long perpendicular baseline are shown in Figure 7a,c. Figure 7a shows the original interferogram with a temporal baseline of 22 days and a perpendicular baseline of −23.217 m. Figure 7c shows the original interferogram with a temporal baseline of 11 days and a perpendicular baseline of 133.176 m. It is clear that the region indicated by the white dotted rectangle in Figure 7c is affected by DEM errors. It is worth noting that the region indicated by the white dotted rectangle in Figure 7a is dominated by unwrapping errors rather than DEM errors, which will be discussed in the next section. DEM errors are estimated with high-quality interferograms and removed from all interferograms. The results of the DEM error correction are given in Figure 7d. It can be seen that the DEM error is well corrected while the deformation signal is well retained.

4.1.2. Phase Unwrapping Error Correction

Serious unwrapping errors occur in this region for two main reasons. The first is the losses of coherence caused by dense vegetation coverage, rainfall and agricultural irrigation. The second is the dense fringes presented in the interferograms caused by DEM errors and large deformation gradient. After the correction of the DEM errors, the interferometric fringes caused by DEM errors were successfully removed and the number of fringes was significantly reduced, which makes the phase unwrapping much easier. The remaining unwrapping errors caused by low coherence and large deformation gradient were then corrected by two different methods mentioned above. The unwrapped interferograms before and after phase unwrapping error correction are shown in Figure 8. In Figure 8a,c, it is evident that the deformation signals of the region indicated by the white dotted rectangle are obscured by unwrapping errors. The unwrapping errors in Figure 8a are corrected by the method of adding 2π integer times phase offset. However, for the unwrapping errors in Figure 8c, after several tests, an optimal coherence threshold is set to obtain the corrected unwrapping phase. The corrected results are given in Figure 8b,d, where it can be seen that the phase unwrapping errors are well corrected and the deformation signals are retained.

4.2. Identification of Potential Landslides

Average line-of-sight deformation rate maps for the Heifangtai loess terrace from December 2006 to November 2017 are obtained with different SAR datasets. Note that negative values (red color) represent movement away from the satellite and positive values (blue color) indicate movement towards the satellite. The blue polygons indicate the distribution of landslide groups. Figure 9 shows the results calculated with different ascending datasets with different time spans. Figure 9a shows the deformation rate map calculated with ALOS/PALSAR data from track 473 from December 2006 to March 2011. It can be seen that the whole landslide groups including Fangtai, Yehugou and Moshigou all suffered deformation. Additionally, some large deformation occurred in the local area of other landslide groups, such as Dangchuan, Huangci, Jiaojia and Chenjia. Figure 9b shows the deformation rate map calculated with ALOS/PALSAR data from track 474 from March 2007 to October 2009. It can be seen that the active deformation areas are basically consistent with those obtained from track 473 datasets. A large deformed funnel was observed in the Huangci landslide group (i.e., Huangci No.3 landslide), the main reason for which is a landslide which occurred in May 2008. Figure 9c shows the deformation rate map calculated with TerraSAR-X data from track 21 from February 2016 to November 2016. It can be seen that some new deformed regions appeared in the study area compared with Figure 9a,b, such as in the Xinyuan, Dangchuan, Jiaojia and Chenjia landslide groups. The whole Dangchuan landslide group was deformed and two new deformation regions appeared in the Xinyuan landslide group. Additionally, it can be seen that the deformed funnel in the Huangci landslide group is not presented in Figure 9c, which suggests that there was no new deformation in this area after the occurrence of the Huangci No. 3 landslide. Figure 9d is the deformation rate map calculated with ALOS/PALSAR-2 data from track 146 from November 2014 to November 2017. It can be seen that the obtained active deformation regions are basically consistent with the results obtained from TerraSAR-X data in most regions of the study area (i.e., the Fangtai and Moshigou landslide group). Due to the low spatial resolution of ALOS/PALSAR-2 data (compared with TerraSAR-X data) and the influence of temporal decorrelation for one-year interferogram, it is evident that some potential landslides were missed in the Xinyuan, Dangchuan and Chenjia landslide groups. Figure 10 shows the results calculated with different descending datasets with different time spans. Figure 10a is the average deformation rate map calculated with TerraSAR-X data from track 165 from January 2016 to November 2016. Similarly, it can be seen that there are large deformations in the whole of the Dangchuan and Moshigou landslide groups. Other large deformations occurred in the local area of the Xinyuan, Yehugou, Jiaojia and Chenjia landside groups. Furthermore, it can be seen that some omissions for landslides appeared in the Fangtai and Chenjia landslide groups due to geometric distortions of SAR images and temporal decorrelation. Figure 10b shows the average deformation rate map calculated with ALOS/PALSAR-2 data from track 39 from May 2015 to July 2017. It can be seen that the results are in good agreement with those obtained from descending TerraSAR-X data. Some differences were observed in the Xinyuan, Yehugou and Jiaojia landslide groups, which may have been caused by the low spatial resolution of ALOS/PALSAR-2 images and the different temporal spans of the two datasets. Figure 10c is the deformation rate map calculated with Sentinel-1A/B data from track 135 from September 2016 to October 2017. It can be seen that the whole of the Dangchuan and Moshigou landslide groups, and the local areas of the Yehugou and Chenjia landslide groups, deformed continuously. The absence of deformation in some regions (i.e., the Xinyuan and Jiaojia landslide groups) may be due to the low spatial resolution of Sentinel-1A/B data resulting in the omission of some small-scale landslides.

Table 3. The changes of area in different landslide groups identified by different SAR datasets during different periods in the Heifangtai terrace.

Landslide Group	Detected SAR Data	Track	Period	Area (km2)
Fangtai	ALOS/PALSAR	473	December 2006–March 2011	0.0795
	ALOS/PALSAR	474	March 2007–October 2009	0.0708
	ALOS/PALSAR-2	146	November 2014–November 2017	0.0678
	ALOS/PALSAR-2	39	May 2015–July 2017	0.0042
	TerraSAR-X	21	February 2016–November 2016	0.0582
	TerraSAR-X	165	January 2016–November 2016	0.0076
Xinyuan	ALOS/PALSAR	473	December 2006–March 2011	0.0059
	ALOS/PALSAR	474	March 2007–October 2009	0.0069
	ALOS/PALSAR-2	146	November 2014–November 2017	0.0064
	ALOS/PALSAR-2	39	May 2015–July 2017	0.0108
	TerraSAR-X	21	February 2016–November 2016	0.0075
	TerraSAR-X	165	January 2016–November 2016	0.0165
Dangchuan	ALOS/PALSAR	473	December 2006–March 2011	0.1024
	ALOS/PALSAR	474	March 2007–October 2009	0.0951
	ALOS/PALSAR-2	146	November 2014–November 2017	0.1016
	ALOS/PALSAR-2	39	May 2015–July 2017	0.1672
	TerraSAR-X	21	February 2016–November 2016	0.0933
	TerraSAR-X	165	January 2016–November 2016	0.1391
	Sentinel-1A/B	135	September 2016–October 2017	0.1610
Huangci	ALOS/PALSAR	473	December 2006–March 2011	0.1168
	ALOS/PALSAR	474	March 2007–October 2009	0.1862
Yehugou	ALOS/PALSAR	473	December 2006–March 2011	0.0703
	ALOS/PALSAR	474	March 2007–October 2009	0.0634
	ALOS/PALSAR-2	146	November 2014–November 2017	0.0313
	ALOS/PALSAR-2	39	May 2015–July 2017	0.0184
	TerraSAR-X	21	February 2016–November 2016	0.0311
	TerraSAR-X	165	January 2016–November 2016	0.0323
	Sentinel-1A/B	135	September 2016–October 2017	0.0420
Jiaojiaya	ALOS/PALSAR	473	December 2006–March 2011	0.0018
	ALOS/PALSAR	474	March 2007–October 2009	0.0007
	TerraSAR-X	21	February 2016–November 2016	0.0023
Jiaojia	ALOS/PALSAR	473	December 2006–March 2011	0.0425
	ALOS/PALSAR	474	March 2007–October 2009	0.0544
	ALOS/PALSAR-2	146	November 2014–November 2017	0.0202
	ALOS/PALSAR-2	39	May 2015–July 2017	0.0043
	TerraSAR-X	21	February 2016–November 2016	0.0309
	TerraSAR-X	165	January 2016–November 2016	0.0413
	Sentinel-1A/B	135	September 2016–October 2017	0.0035
Chenjia	ALOS/PALSAR	473	December 2006–March 2011	0.0361
	ALOS/PALSAR	474	March 2007–October 2009	0.0622
	ALOS/PALSAR-2	146	November 2014–November 2017	0.0130
	ALOS/PALSAR-2	39	May 2015–July 2017	0.0275
	TerraSAR-X	21	February 2016–November 2016	0.0345
	TerraSAR-X	165	January 2016–November 2016	0.0128
	Sentinel-1A/B	135	September 2016–October 2017	0.0413
Moshigou	ALOS/PALSAR	473	December 2006–March 2011	0.0300
	ALOS/PALSAR	474	March 2007–October 2009	0.0402
	ALOS/PALSAR-2	146	November 2014–November 2017	0.0691
	ALOS/PALSAR-2	39	May 2015–July 2017	0.0675
	TerraSAR-X	21	February 2016–November 2016	0.0509
	TerraSAR-X	165	January 2016–November 2016	0.0569
	Sentinel-1A/B	135	September 2016–October 2017	0.0814

shows the changes of area in different landslide groups in the Heifangtai terrace identified by different SAR datasets during different periods. The following conclusions can be drawn:

(i)

It can be seen from Figure 9 and Figure 10 that the deformation regions calculated with different SAR datasets during a similar period are in good overall agreement, meaning that slow-rate surface deformation can be measured by different SAR data once it can be mapped in the line-of-sight direction. Furthermore, the existence of differences among different SAR datasets are mainly caused by different SAR acquisition parameters, including SAR acquisition date, wavelength, spatial resolution, local incidence angle and satellite tracking direction, which can be used to obtain more detailed information on each individual landslide.

(ii)

The deformation region of the Heifangtai terrace before 2011 was significantly different than that after 2014. From December 2006 to March 2011, the deformation regions of the Heifangtai terrace were mainly concentrated in the Fangtai, Yehugou and Moshigou landslide groups. However, after 2014, the main deformation region extended to Dangchuan landslide group, where large deformation occurred. Furthermore, some new deformation regions developed in the Xinyuan, Jiaojia and Chenjia landslide groups, which suggests an increasing trend of landslide distribution from December 2006 to November 2017.

(iii)

Loess landslides and loess-bedrock landslides have different formation processes. In loess landslides, the deformation of the landslide does not terminate even if it has already slid, but the deformation continues to occur on the back edge of the landslide to form a new landslide (as seen in the Dangchuan landslide group, where three large landslides occurred in 2015, and some deformation can still be monitored subsequently), which is closely related to the retrogressive failure mode of loess landslides [49]. However, in loess-bedrock landslides, no obvious deformation can be monitored in a short period of time once a landslide has already occurred.

(iv)

It can be seen from Table 3 that the area of active landslides identified by ALOS/PALSAR datasets from track 473 and track 474 are highly different in the Huangci, Jiaojia and Chenjia landslide groups. The area of active landslides identified on track 474 images is larger than that identified on track 473 images. This difference may be due to the different acquisition times of the two datasets. According to Peng et al. [11] the largest number of landslides occurred in 2007 and 2008, and the acquisition time of track 474 images is mainly concentrated in this period.

Three InSAR-derived products, including the deformation maps, SAR intensity images and the DEM gradient map are combined to detect active landslides. We set a stability threshold range of −1 cm to 2 cm, and set the threshold of DEM slope at 10°. First, active landslides are provisionally identified by independent InSAR measurements from the single satellite orbit datasets. Afterwards, landslide candidates from different satellite tracks are combined into a mosaic. Some landslide candidates are further cross-validated by independent InSAR observations from different satellite tracks.

A total of 32 potential landslide areas are identified from December 2006 to March 2011 (the first stage), and a total of 48 active landslides are identified from January 2016 to November 2016 (the second stage), which are shown in Figure 11. Note that the deformation rate calculated with L-band ALOS/PALSAR data is merged together and the single small landslides cannot be distinguished due to their low spatial resolution with respect to X-band TerraSAR-X data. Therefore, the potential landslide areas in Figure 11a may contain several small active landslides. It can be seen that the Fangtai landslide group contains one large landslide area for the first stage and three small landslides for the second stage; the Xinyuan landslide group contains one landslide for the first stage and three small landslides for the second stage; the Dangchuan landslide group contains one large landslide area for the first stage and eleven small landslides for the second stage; the Huangci landslide group contains one landslide for the first stage, and no detectable active landslides for the second stage; the Yehugou landslide group contains four landslides for the first stage and six landslides for the second stage; the Jiaojiaya landslide group contains two landslides for the first stage and one landslide for the second stage; the Jiaojia landslide group contains five landslides for the first stage and three landslides for the second stage; the Chenjia landslide group contains ten landslides for the first stage and nine landslides for the second stage; and the Moshigou landslide group contains six landslides for the first stage and nine landslides for the second stage. The distribution of landslides identified by the InSAR technique is in good agreement with that identified by field geological survey [11].

4.3. Landslide Identification Based on SAR Intensity Changes and DEM Errors

On 1 October 2017, three large landslides, namely Dangchuan No.4 (DC#4), No.5 (DC#5) and No.9 (DC#9) landslides, occurred simultaneously and released a volume of about 33.9 × 10⁴ m³ of debris. The photo of the scene after the landslides is shown in Figure 12a. Figure 12b shows an unwrapped interferogram formed from two SAR images of 28 September 2017 and 10 October 2017, which covered the date of the landslides. No useful information can be obtained due to the serious decorrelation over the landslide area caused by the mass transformation. In this case, SAR intensity change method is applied to locate the landslide area. The two closest SAR intensity images from the Sentinel-1A/B data covering the landslide event are selected, shown in Figure 12c,d. It can be seen from these figures that the intensity information of the landslide area clearly changed after the landslide occurred. The three landslides were clearly recorded in the intensity map in Figure 12d. However, due to the low spatial resolution (16 m) of Sentinel-1A/B images with respect to the landslide extent (about 140 m in width and 270 m in length), the whole landslide area cannot be precisely delineated by the SAR intensity change map. In order to quantitatively estimate the volume change caused by the landslides, DEM error is estimated with the interferograms generated after these events. Six high-quality interferograms from Sentinel-1A/B data with short temporal baseline and long spatial baseline are selected to calculate DEM error based on Equation (1), and are shown in Figure 13. Note, negative values (red color) represent a decrease in elevation (the source area) and positive values (blue color) indicate an increase in elevation (the deposit area). It is evident that serious DEM errors appeared in both the source area and deposit area. Two areas are completely and precisely mapped based on the DEM error estimation.

5. Discussion

In order to assess and generate the multi-temporal inventory maps, some concerns are discussed with respect to the effects of wavelength, resolution, track direction and local incidence angle.

5.1. The Effects of Different SAR Sensors and Bands

The main features of SAR sensors include the resolution, wavelength and data duration. The Jiaojia landslide group is taken as an example to illustrate the effects of different SAR sensors and bands. The line-of-sight deformation rates of the landslide group (Figure 14) are independently calculated by three different satellite datasets with different spatial resolutions and wavelengths, including L-band ALOS/PALSAR-2 datasets (Figure 14a), X-band TerraSAR-X datasets (Figure 14b) and C-band Sentinel-1A/B datasets (Figure 14c). It can be seen from Figure 14 that the deformation results for SAR sensors with different spatial resolutions and wavelengths are significantly different in the landslide area. Three active landslides, namely Jiaojia No.5 (JJ#5), Jiaojia No.6 (JJ#6) and Jiaojia No.9 (JJ#9) (as shown in Figure 14d) were successfully identified by X-band TerraSAR-X datasets, but failed to be identified from both L-band ALOS/PALSAR-2 images and C-band Sentinel-1A/B images. In 2017, Peng et al. [11] identified three active landslides based on aerial images, DEM and field investigations. The landslides identified by X-band TerraSAR-X data in the present study are in good agreement with the results obtained by Peng et al. We conducted a field investigation in August 2018 to further validate the reliability of TerraSAR-X data identified landslides. The scene photos of the JJ#5 and JJ#6 landslides are shown in Figure 15, where dense cracks were distributed. A crack with a maximum width of about 55 cm was found at the back edge of the JJ#6 landslide. Three landslides have lengths of about 118 m, 326 m and 100 m, and widths of about 49 m, 108 m and 66 m, respectively, and the multi-looked pixel sizes of ALOS/PALSAR-2 images, TerraSAR-X images and Sentinel-1A/B images are about 6.5 m, 3 m and 16 m, respectively. Additionally, L-band and C-band data have lower sensitivity to the surface deformation than X-band data. Therefore, it can be inferred that the L-band and C-band data are not suitable for the detection of small spatial scale and small deformation loess landslides. Nevertheless, L-band ALOS/PALSAR-2 and C-band Sentinel-1A/B data show higher coherence than X-band data, which can be applied for long term surface deformation monitoring even with a large deformation gradient.

5.2. The Effects of Satellite Track Direction

The Xinyuan No.2 landslide is taken as an example. Peng et al. [11] identified this landslide on the basis of aerial images, DEM and field investigations in 2017. Figure 16 shows the line-of-sight deformation rates of the Xinyuan No.2 landslide obtained from ascending (Figure 16a) and descending (Figure 16b) TerraSAR-X datasets from January 2016 to November 2016. It can be seen that the landslide is identified independently by both ascending and descending TerraSAR-X data, but that there are great differences in the range and magnitude of the deformation. A scene photo of the landslide was taken in August 2018, as shown in Figure 16d. Dense cracks are visible on the body of the landslide, implying that the rupture surface had run through. The imaging geometries of ascending and descending data and the main landslide direction are shown in Figure 16c. The azimuth of the landslide ranges from 80 to 220° (the main azimuth of the landslide is about 130°), and the slope angle ranges from 10 to 40° above the horizontal surface. The flight directions of ascending and descending of TerraSAR-X satellite are −9.7° and 189.6°, respectively, and the incidence angles are 41.2° and 41.8°, respectively. It can be seen from Figure 16c that the landslide moving in a downslope direction shows movement away from the ascending line-of-sight direction (in red), while movement towards the descending line-of-sight direction (in blue). The different shapes of landslide in geodetic coordinate between ascending and descending SAR data indicate the different sensitivity of two imaging geometries to different local deformation directions. For example, the small deformation region of descending data indicates that the edge of the landslide was approximately parallel to the descending satellite flight direction, and hence no deformation can be sensed. In this case, the ascending and descending SAR images have been combined to achieve the complete deformation fields and to identify the whole landslide regions.

5.3. The Effects of the SAR Geometric Distortions

Due to the side-looking imaging geometry of SAR satellites, geometric distortions including layover, shadow and foreshortening are inevitable in mountainous regions, which will cause some blind areas in single-orbit observations and seriously decrease the capability of landslide identification. The degree of geometrical distortion is determined by the local slope of the topography and the SAR satellite geometry. For landslide identification by ascending and descending SAR datasets, geometric distortions will occur in different areas, which can lead to inconsistency in the results from ascending and descending tracks [47]. Figure 17a,b shows the deformation rate maps for the Dangchuan landslide group obtained from ascending and descending TerraSAR-X datasets, respectively. Figure 18a,b shows the deformation rate maps for the Fangtai landslide group. It can be seen from Figure 17a,b and Figure 18 that the deformation rates calculated with ascending and descending data are inconsistent. Geometric distortions are mapped for the ascending and descending TerraSAR-X images by exploiting local terrain information (i.e., local terrain slope angle and local terrain azimuth) and SAR imaging geometries (i.e., the azimuth and local incidence angle of the satellite), which are shown in Figure 19a,b, respectively. It can be observed in Figure 19 that some regions in the Dangchuan landslide group are affected by severe geometric distortions (i.e., layover) in the ascending image while some regions in the Fangtai landslide group are affected by severe geometric distortions (i.e., layover) in the descending image. The geometric distortions for the ascending and descending data are quite complementary. Therefore, we can infer that the deformation rates of the Dangchuan landslide group calculated with ascending data are contradictory to the deformation rates calculated with descending data because of the effects of the geometric distortions. The geometric distortions in the descending images caused some blind areas in the Fangtai landslide group, which results in some active landslides being missed when using only descending data. Therefore, it is demonstrated that the combination of ascending and descending InSAR measurements can effectively reduce geometric distortion effects in rugged terrain. Thus, the identification of loess landslides can be effectively achieved.

For the Dangchuan landslide group, previous landslide identification studies mainly focused on aerial images, DEM, optical remote sensing and field investigations [11,67,68,69]. In this study, the InSAR technique is used for the first time to detect potential landslides. The results of landslide identification using the InSAR technique are in good agreement with those in the literatures. A field investigation was conducted in August 2018 showed that dense cracks and sliding scarps were found in the whole of the Dangchuan landslide group. The scene photos after the occurrence of the Dangchuan No.3 (DC#3) and Dangchuan No.4 (DC#4) landslides are shown in Figure 17c,d. Two long cracks (as shown in Figure 17d) with widths of about 30 cm and 40 cm were found at the back edge of DC#3 and DC#4, respectively.

5.4. The Effects of Differences in Local Incidence Angle between Two Adjacent SAR Satellite Tracks

The line-of-sight deformation rate maps for the Dangchuan landslide group between 2007 and 2009 are calculated using ALOS/PALSAR datasets from the adjacent tracks 473 and 474, as shown in Figure 20a,b, respectively, where it can also be seen that a good agreement in the spatial distribution of the active deformation areas can be obtained from both adjacent track datasets. The maximum deformation rates calculated using track 473 and track 474 are −90 mm/year and −100 mm/year, respectively. Scatter plots of the deformation rates obtained from tracks 473 and 474 are shown in Figure 20c. The root mean square error (RMSE) of the difference and the correlation coefficient are 7.3 mm/year and 0.81, respectively, which shows a good consistency of deformation rates. However, it can be observed in Figure 20c that a systematic offset still exists between the two results. In order to explain this difference, the differences in the flight direction and local incidence angle between the two adjacent satellite track data are considered. The flight directions of the ALOS/PALSAR satellite in tracks 473 and 474 are −10.3° and 9.9°, respectively, and the local incidence angles of tracks 473 and 474 at the Dangchuan landslide group are 37.2° and 39.8°, respectively. The deformation rates in the line-of-sight direction and in the downslope direction along line BB’ are shown in Figure 20d, where it can be seen that two large deformed funnels correspond to the Dangchuan No.2 and Dangchuan No.3 landslides, which occurred in April 2015 and August 2015, respectively. Furthermore, the discrepancy between the two profiles in the downslope direction is much less than that in the line-of-sight direction, which indicates the high-accuracy of InSAR measurements.

6. Conclusions

The InSAR technique is firstly explored to detect small-scale loess landslides reliably and accurately over a large area. The flowchart of multi-temporal loess landslide inventory mapping is proposed considering its special characteristics. Some key data processing strategies are summarized as follows: (1) It is very important to select SAR images and external DEM data with meter-level resolution for small-scale loess landslide detection. Additionally, dense SAR acquisition time and small multi-look number are two key factors for improving the capability of landslide detection; (2) Two special data procedures including DEM error correction and phase unwrapping error detection and correction must be conducted in the data processing. Furthermore, high-quality interferograms should be selected to obtain surface deformation; (3) For single-orbit InSAR observations, it is difficult to detect landslides with deformation parallel to the direction of the satellite flight or located in regions of SAR geometric distortion.

The combination of multi-sensor and multi-temporal SAR datasets is used to investigate small-scale loess landslides in the Heifangtai loess terrace from December 2006 to November 2017. The InSAR-derived products including deformation map, SAR intensity image and external DEM are combined to detect and map potential loess landslides and already-occurred landslides. Three different wavelengths of SAR datasets are compared in detail in term of landslide inventory mapping. The phase information and intensity information of SAR images are quite complementary for the identification of potential and already-occurred small-scale loess landslides. As for the potential for loess landslide detection, all SAR datasets can measure most landslides. Some big differences mainly result from the different satellite tracks, i.e., ascending and descending tracks, and different SAR wavelengths and spatial pixel sizes. Some slight differences are caused by the slight variation in local incidence angle between two adjacent tracks.

Two main stages of surface deformation of the whole Heifangtai loess terrace from December 2006 to November 2017 are obtained by the InSAR technique. A total of 32 active landslide areas are detected from December 2006 to March 2011, and 48 active landslides are detected from January 2016 to November 2016. These potential landslides are all distributed on the edge of the Heifangtai terrace. The spatiotemporal evolution of landslides is analyzed in details. Additionally, it is observed that loess landslides and loess-bedrock landslides show different formation processes. This research is of great significance to the early-warning and prevention of loess disasters in the study area.