**<1 km**is obtained, which gives a sound precursor for small-scale loess landslide detection. Moreover, the longer and continuous land subsidence has been monitored while deformation starting point for the landslide is successfully inverted, which is key to monitoring the similar loess landslide. In addition, the accelerated landslide deformation from one to two months before the landslide can provide a critical clue to early warning of this kind of landslide.

## 1.

## Introduction

Loess in China occupies 0.63 million square kilometer and exists in 12 provinces, of which 0.4 million square kilometer is loess plateau. Loess plateau is the main product base of oil, gas, and coal resources. Unfortunately, it is also a geo-hazard-prone region, including collapse, landslide, earthflow, land subsidence, and ground fissures. In the past 30 years, thousands of people died by the loess geo-hazards.^{1} Currently, loess geo-hazards threaten the transportation system, local citizens, and natural resource exploration. The trigger factors of a loess landslide comprise of earthquake, strong rainfall, and seasonal irrigation.^{2}3.4.5.^{–}^{6} The characteristics of a loess landslide are special in the following aspects: group occurrence, small spatial scale, and temporal abrupt deformation for most specific landslides. Consequently, the detection of loess landslides over a large area, high-precision deformation monitoring, and early warning of the specific loess landslide are still challenging. The landslide in the southern bank of Jinghe River, Shaanxi province, China, is one of the typical cases, which is mainly triggered by farmland irrigation.^{7}8.^{–}^{9} The southern bank of Jinghe River, a small section of the subregion of Weibei loess plateau in Weihe basin, is 27.1 km in length, including three towns as Taiping, Jiangliu, and Gaozhuang and 21 villages. The heights on this terrace range from 450 to 580 m, with the slopes from 4 to 9‰. The height differences between the top and toe of the terrace range from 30 to 90 m. Since 1976, around 50 landslides over 27 sites have taken place along this region.^{7} These landslides have seasonal features, which are correlated with heavy rainfall and farmland irrigation. More than 80 percent of the landslides occurred in spring or autumn, with a certain time lag of irrigation or heavy rainfall.

Different interferometric synthetic aperture radar (InSAR) techniques have been widely applied to detect the locations of occurred and potential landslides over a large coverage region.^{10}11.^{–}^{12} Zhao et al.^{10} have reported the landslide inventory map in Northwestern US under the condition of forest, generated with L-band ALOS/PALSAR data. The length of a single landslide is $>500\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, and the deformation history is as long as several years. For the large-area loess landslide detection, previous works have been done for the postlandslides inventory mapping, investigation and mechanism analysis.^{6}^{,}^{7} Due to the particular characteristics of loess landslides, such as small spatial extent, small deformation magnitude, and abrupt occurrence in time domain,^{13} no mature InSAR techniques can be readily applied to detect and monitor them. We first point out a two-step method to detect potential loess landslide using small baseline subsets (SBAS) InSAR techniques according to the loess landslides mechanism.^{8} First, we detect land subsidence on the terrace, which can be easily realized due to its large spatial extent and long deformation history. Then we will closely search loess landslides within a given distance and direction of land subsidence center according to topography and hydrology conditions.

On the other hand, multiple InSAR techniques have been successfully employed to retrieve the prelandslide deformation history and even postlandslide surface changes.^{12}^{,}^{14}15.^{–}^{16} However, the spatial extent of a rockslide^{14} and a landslide^{15}^{,}^{16} is much larger than a loess landslide. As far as landslide monitoring is concerned, few onset times of landslides have been precisely detected and short-term early warning can be issued accordingly. Special SAR data and InSAR data processing procedures are analyzed in this paper. We take the Xingyuan landslide as an example, which took place on April 10, 2012 (Fig. 1). Meter resolution SAR data from TerraSAR-X are selected. Several specific InSAR processing measures have been taken, such as digital elevation model (DEM) error correction, phase unwrapping error detecting and repairing before the SBAS InSAR processing. Most importantly, precursory signals for Xingyuan landslide are successfully recovered, which in essence is crucial for landslide early warning and prevention.

## 2.

## Methods

Different InSAR techniques have been developed over the last 25 years, including permanent scatterers (PS) InSAR^{17}18.^{–}^{19} and SBAS techniques.^{20}^{,}^{21} For small-scale and abrupt loess landslide deformation monitoring, the PS InSAR technique can hardly be applied due to the few PS points and SAR images. On the other hand, the SBAS InSAR technique has been widely applied owing to its flexibility, high precision, and time-series deformation result recovery. However, for small-scale landslide detection and monitoring, four steps should be carefully considered, such as DEM error correction, phase unwrapping error detection and correction, high-quality interferograms selection, and downslope deformation transformation. The flowchart is shown in Fig. 2.

## 2.1.

### Digital Elevation Model Error Correction

In order to get surface deformation, external DEM is needed to subtract the topographic phase from original SAR interferogram, generated from coregistrated master and slave images. For the landslide-prone region, no suitable DEM is available due to the repeated landslide occurrences historically. That is, no real DEM can be available during the SAR acquisition periods, which will inevitably result in errors in the final deformation results. As the height is sensitive to the long perpendicular baseline component according to the height-to-phase formula, long spatial baseline and small temporal baseline InSAR pairs are used to invert the linear deformation and DEM error simultaneously.^{20} In addition, the interferograms contaminated by severe atmospheric effects should not be included. For a given interferogram $j$, after removing the flat earth and local topography, the unwrapped differential phase at pixel $x$, $r$ ($x$ and $r$ are the azimuth and range coordinates, respectively) computed from the SAR acquisitions at epoch ${t}_{A}$ (start time) and ${t}_{B}$ (end time) can be written as follows:

## (1)

$$\delta {\varphi}_{j}(x,r)=\varphi ({t}_{B},x,r)-\varphi ({t}_{A},x,r)\phantom{\rule{0ex}{0ex}}=\delta {\varphi}_{j}^{\mathrm{topo}}(x,r)+\delta {\varphi}_{j}^{\mathrm{disp}}(x,r)+\delta {\varphi}_{j}^{\text{noise}}(x,r),$$## (2)

$$\{\begin{array}{l}\delta {\varphi}_{j}^{\mathrm{topo}}(x,r)=\frac{4\pi}{\lambda}\frac{{B}_{\perp j}\mathrm{\Delta}Z}{R\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\theta}\\ \delta {\varphi}_{j}^{\mathrm{disp}}(x,r)=\frac{4\pi}{\lambda}[d({t}_{B},x,r)-d({t}_{A},x,r)]=\frac{4\pi}{\lambda}v\xb7\mathrm{\Delta}t\end{array},\phantom{\rule[-0.0ex]{1em}{0.0ex}}\forall \text{\hspace{0.17em}\hspace{0.17em}}j=1,\dots ,M$$## 2.2.

### Phase Unwrapping Detection and Correction

Phase unwrapping error often occurs in the case of decorrelation, which is always caused by surface random change, such as vegetation, plantation and irrigation, and long spatial and temporal baseline between master and slave SAR images. Generally, a minimum cost flow unwrapping algorithm with the aid of coherence-induced mask file is used to unwrap the interferogram.^{22} Different coherence thresholds can be set to get different unwrapped interferograms, and their differences mainly occur in the low-coherent areas. In case the temporally low-coherent regions are the targets of interest, two methods can be applied to correct the unwrapped errors. (1) Smaller coherence threshold is set to get denser unwrapped correct phase, which can provide more candidates for later results analysis. (2) As all interferograms are unwrapped individually, the phase contributions behave in a conservative manner. That is, ${\varphi}_{AB}+{\varphi}_{BC}-{\varphi}_{AC}=0$, where ${\varphi}_{AB}$ is the phase contribution to interferogram $AB$, which is constructed from acquisitions $A$ and $B$. If only one interferogram has unwrapping error, a phase closure technique can be used to detect the unwrapping errors.^{23} And the phase offset ($2\pi $ phase integer multiplier) can be added back to the wrongly unwrapped phase to correct it.

## 2.3.

### Small Baseline Subsets Interferometric Synthetic Aperture Radar Algorithm

Once deformation related interferograms have been successfully unwrapped, deformation time-series results can be easily retrieved under the norm of least square or singular value decomposition, which is determined by the number of subsets involved.^{20}

For a generic pixel $(x,r)$, the phase for $M$ interferograms can be given in matrix form as follows.

## (3)

$$\delta \varphi (x,r)=[A(x,r),B(x,r)]\left[\begin{array}{l}v\\ \mathrm{\Delta}Z\end{array}\right]+N(x,r),$$^{17}Then, the nonlinear deformation phase at different SAR acquisition dates will be calculated. Last, the accumulated deformation time series can be achieved by integrating the deformation in each neighboring SAR acquisition date as follows:

## (4)

$$\delta {\varphi}_{k}^{\mathrm{defo}}(x,r)=({t}_{k}-{t}_{1})v+\sum _{i=2}^{k}{\varphi}_{i}^{\mathrm{nl}\_\mathrm{defo}},$$## 2.4.

### Deformation in Slope Direction

InSAR can only measure the deformation along the line-of-sight (LOS) direction. That is, the InSAR measurement represents the component of ground deformation projected to the LOS direction. Once the landslide geometry is obtained, the LOS deformation will be projected to the downslope direction to analyze the landslide activity.

The transformation between the LOS direction and downslope direction is straightforward. The unit vector of LOS can be defined as follows:^{14}^{,}^{24}

## (5)

$$\overrightarrow{r}=\left[\begin{array}{c}{r}_{\text{east}}\\ {r}_{\text{north}}\\ {r}_{\text{up}}\end{array}\right]=\left[\begin{array}{c}-\mathrm{sin}\text{\hspace{0.17em}}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}\\ \mathrm{sin}\text{\hspace{0.17em}}\theta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}\\ \mathrm{cos}\text{\hspace{0.17em}}\theta \end{array}\right],$$## (6)

$$\overrightarrow{u}=\left[\begin{array}{c}{u}_{\text{east}}\\ {u}_{\text{north}}\\ {u}_{\text{up}}\end{array}\right]=\left[\begin{array}{c}-\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\phi \\ -\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\phi \\ \mathrm{sin}\text{\hspace{0.17em}}\phi \end{array}\right],$$## 3.

## Background and Data

## 3.1.

### Xingyuan Landslide Background

The Xingyuan landslide occurred on April 10, 2012. The perspective image and geological profile of the Xingyuan landslide are shown as Figs. 1(b) and 1(c). The landslide has a moderate slope of $\sim 40\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ and its relative height between the terrace and the bottom is $\sim 75\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. The loess layer can be divided into Malan loess, Lishi loess from the top to the bottom, among which some Paleosol imbedded. The source of landslide is 190 m in width along the edge of terrace, and 66 m along the sliding direction including slope section while the deposit is 220 m in width and 144 m in length.

## 3.2.

### Synthetic Aperture Radar Data

In order to detect and monitor small spatial loess landslide prior to the Xingyuan landslide, 11 TerraSAR-X SAR images generated in strip mode with 3 m ground pixel posting are selected, and the multilook number is set as 2 and 2 in range and azimuth directions, respectively. After setting the spatial baseline and temporal baseline thresholds, and visually checking the interferogram quality, 17 interferograms are generated for further analysis. Table 1 presents the parameters of interferometric pairs and their functions for data analysis, including DEM error calculation, phase unwrapping error correction, and time-series deformation estimation, where “a” indicates the considered interferogram. Figure 3 shows 15 interferometric SAR configurations for time-series deformation estimation, and 7 of them with longer perpendicular baselines shown in red dashed lines are also applied to estimate DEM error.

## Table 1

The parameters of interferometric pairs and their functions for later analysis.

No. | Master | Slave | Perpendicular baseline (m) | Time duration (days) | DEM error calculation (Fig. 4) | Phase unwrapping error correction | Time-series deformation estimation (Fig. 3) |
---|---|---|---|---|---|---|---|

1 | 20111107 | 20111118 | 67.0 | 11 | a | ||

2 | 20111107 | 20111129 | $-62.0$ | 22 | a | ||

3 | 20111118 | 20111129 | $-129.0$ | 11 | a | a | |

4 | 20111118 | 20111210 | $-97.6$ | 22 | a | ||

5 | 20111210 | 20111221 | 261.2 | 11 | a | a | |

6 | 20111210 | 20120101 | 84.5 | 22 | a | ||

7 | 20111210 | 20120112 | 53.2 | 33 | a | ||

8 | 20111221 | 20120101 | $-176.7$ | 11 | a | a | |

9 | 20111221 | 20120112 | $-208.0$ | 22 | a | a | |

10 | 20120101 | 20120112 | $-31.3$ | 11 | a | ||

11 | 20120101 | 20120203 | 71.8 | 33 | a | ||

12 | 20120112 | 20120203 | 103.2 | 22 | a | a | |

13 | 20120203 | 20120214 | $-156.8$ | 11 | a | a | |

14 | 20120214 | 20120225 | 174.8 | 11 | a | Fig. 6(a)a | |

15 | 20120225 | 20120307 | 14.8 | 11 | Fig. 5a | Fig. 6(b)a | |

16 | 20120307 | 20120318 | $-85.5$ | 11 | Fig. 6(c) | ||

17 | 20120318 | 20120409 | 14.4 | 22 | Fig. 6(d) |

## a

Indicates the interferogram considered

## 3.3.

### Digital Elevation Model Data

Aiming to meet the requirement of high-resolution TerraSAR data processing and avoid the interpolation error with low-resolution DEM, such as Shuttle Radar Topography Mission 3 arc sec DEM, we choose 1 to 10,000 scale topographic map, digitalize, and grid it into a DEM. As it was mapped in 1995 and several landslides had occurred in this region since then, this DEM will be corrected after calculating the DEM error with long spatial baseline InSAR pairs as shown in Table 1.

## 4.

## Results

## 4.1.

### Errors and Corrections

## 4.1.1.

#### Digital elevation model errors estimation

As the DEM applied in this research was generated 20 years ago, surface topography has been changed due to the previous landslides. In order to mitigate the deformation error resulting from low-accuracy DEM, the DEM error is first estimated and corrected. In this case, seven interferograms with long perpendicular baseline and short time duration are selected, which is shown in Fig. 4. The last frame of Fig. 4 is the estimated DEM error map, from which DEM error can be visually seen at the edge of the Jinghe River terrace including the Xingyuan landslide region. The maximum DEM error amounts to 30 m.

## 4.1.2.

#### Phase unwrapping error correction

Phase unwrapping error arises when the imaged area suffers decorrelation due to rainfall, irrigation, or large deformation gradient, which is shown as the offset of $2\pi $ phase integer multiplier. In this research, different coherence thresholds are tested to get corrected unwrapping results. Figure 5 shows the interferograms before and after phase unwrapping error correction. Figure 5(d) is the coherence map, where the Xingyuan landslide is isolated by low-coherent pixels, which results in unwrapping error by setting coherence threshold as 0.5 [see Fig. 5(b)]. After several tests, if the coherence threshold is 0.15, the corrected unwrapping phase can be obtained [see Fig. 5(a)]. As expected, the difference between correct and wrong unwrapping phase is $2\pi $ phase integer multiplier [see Fig. 5(c)]. Alternatively, robust phase unwrapping method with L0 norm method has been tested and similar results can also be obtained.^{25}

## 4.1.3.

#### Interferograms selections

As for the landslide deformation monitoring, in order to get continuous deformation results in spatial and temporal domain, especially the deformation results can be inverted as close to the landslide event as possible; four key interferograms and their coherence maps are carefully compared. From Figs. 6(a) to 6(d), we can see the decorrelated areas are altered from region I to IV, which is related to the wheat growth and irrigation in spring, 2012. So the first two interferograms are involved for SBAS calculation. Since February 25, 2012, no deformation can be obtained on the top of Xingyuan landslide (see the explanation of Fig. 10), while the last two interferograms are not considered due to their complete decorrelation.

## 4.2.

### Time Series Surface Deformation Results

Based on 15 good and corrected unwrapping interferograms, time-series deformation results are calculated with least squares method. In Fig. 7, we take November 7, 2011 as the reference date, so the cumulative deformations in LOS direction at each SAR acquisition date are calculated.

## 4.3.

### Small-Scale Loess Landslide Detection

The loess landslide is so small in spatial domain and temporally decorrelated that both PS-InSAR and large, multilooked SBAS InSAR techniques can hardly detect its location. However, as for the irrigation triggered loess landslides, surface land subsidence and underground water levels are well correlated to the deformation of landslides. Figure 8 is the cumulative surface deformation within 121 days from November 7, 2011 to March 7, 2012, where two land subsidence cones have been successfully mapped and four wells have been investigated in the field. We draw two fans extending to the cliff direction from the center of land subsidence cones, where suspicious landslides can be expected. In this case, two small landslides, namely Dongfeng and Xingyuan, are successfully detected, which both took place in 2012.

## 4.4.

### Xingyuan Landslide Deformation

We then focus on the Xingyuan landslide deformation result. According to the geometry of this landslide, that is the mean slope of 40.8 deg and the main azimuth of landslide direction of 325 deg from the north, the downslope deformation within 121 days from November 7, 2011 to March 7, 2012 is calculated and shown in Fig. 9 in three-dimensional perspective. The cumulative time-series deformation along line $AB$ in Fig. 8 from November 7, 2011 to March 7, 2012 is shown in Fig. 10.

From Fig. 9, we can clearly see the spatial deformation distribution of the Xingyuan landslide. The width of this landslide is $\sim 130\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ while its length is only 75 m along the main landslide direction, which has high consistence with the field investigation.^{13} Furthermore, from Fig. 10, we can also see the landslide length can be divided into two sections, one is on the top of the slope with 15 m in length, and the other is on the slope with around 60 m in length. Obviously, this landslide is very small in spatial domain, so that it can possibly be monitored with high-resolution SAR data and small multilook number SAR processing.

From Fig. 10, it is obvious that the deformation mainly concentrated on the slope section; the maximum deformation within 121 days is more than 4 cm located on the top of the slope as indicated with the blue column. The other deformed length is 60 m, which is highly correlated with the slope as indicated with red column. We can also see the abrupt deformation that occurred after February 14, 2012, and nonlinear deformation can be monitored from February 14 to February 25, 2012, and from February 25 to March 7, 2012. However, due to the severe decorrelation on the farmland shown in Fig. 6(a), no useful deformation results can be obtained on the start section of the profile AB since February 25, 2012. Moreover, no good interferograms can be applied to monitor the quite early deformation before April 10, 2012, when the Xingyuan landslide took place.

## 5.

## Discussion

## 5.1.

### Trigger Factor of Xingyuan Landslide

The Xingyuan landslide is one of the typical irrigation triggered loess landslides. In order to uncover the correlation between the Xingyuan landslide and farmland subsidence nearby, the landslide rate and the land subsidence rate in each 11 days are calculated before the Xingyuan landslide, which is shown in Fig. 11. Meanwhile, the values of underground water levels measured on March 16, 2012 and May 25, 2012 are also shown in Fig. 11. We can see the following features: (1) During the 154 days before the Xingyuan landslide one continuous land subsidence cone is identified at $\sim 620\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ to the south of the Xingyuan landslide, covering 570 m by 130 m. Also, the land subsidence rate increased with time. (2) The underground water level increased 5.7 m from March 16 to May 25, 2012; this shows a different mechanism from urban land subsidence. (3) As far as the Xingyuan landslide deformation is concerned, 88 days prior to the landslide is a turning point, after which the landslide rate increased dramatically. The landslide rate increased from $0.3\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}/\mathrm{d}$ 55 days prior to landslide to $1.6\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}/\mathrm{d}$ 33 days prior to the landslide, increasing five times in only 22 days. This is the first report on loess landslide onset time detection with InSAR technique. Unfortunately, due to the lack of monitoring data in 33 days prior to landslide, we cannot get the total deformation during this 33-day interval. Therefore, accelerated preslide deformation gives us an important precursor of landslides, which is important for mitigating potential hazards of landslides. In order to research the irrigation triggered landslide mechanism, physical interpretation of empirical relations should be followed.^{26}^{,}^{27}

The mechanism of irrigation triggered landslide is based on the relationship between shear strength and water content. Both rainfall and farmland irrigation make the underground water levels increase by charging surface water. Hence, the static and dynamic water pressure is increased, which makes the shear strength decrease; as a result, the landslide takes place.

## 5.2.

### Small-Scale Loess Landslide Detection

In Fig. 11, it can be clearly seen that the deformation of the Xingyuan landslide has high correlation with land subsidence, i.e., the more the land subsidence, the faster the landslide. For the Xingyuan case, we find that the land subsidence has longer deformation history but less deformation rate than that of a landslide when a landslide initiates.

From Fig. 8, we can see land subsidence is easy to be monitored with SBAS technique for its larger coverage and continuous deformation. However, it can be difficult to monitor an irrigation triggered loess landslide because of its small spatial coverage, special location, and abrupt deformation in temporal domain. Based on this research, we can first monitor the land subsidence and investigate the underground water level logs. If both continuous land subsidence and the underground water level rising can be monitored, the area that is a certain distance away from the land subsidence center, and satisfies the landslide conditions, has a high probability of being a loess landslide. Accordingly, special investigation means including field works should be taken to verify it. Moreover, the abnormal rainfall and seasonal irrigation will accelerate the landslide.

## 6.

## Conclusion

Small-scale loess landslide is very common in the loess plateau of Northwest China, which can result in severe damage to local residents. InSAR technique provides a unique way to recover the preslide time-series deformation results and the deformation related to the landslide. From this research, we can draw the following conclusions.

Small loess landslides can be possibly detected with the two-step method, i.e., first detecting large coverage land subsidence with longer deformation history and then closely searching loess landslides within given distances and direction angles according to local topography and hydrology conditions. Two loess landslides are successfully detected in this research corresponding to two land subsidence cones. The inner correlation between land subsidence and landslides must be clear when using this method.

Meter level resolution SAR data is one key factor to small-scale loess landslide monitoring. Dense SAR acquisitions and small multilook number are other factors to get detailed time-series deformation results. Reliable time-series results are achieved after taking correction of external DEM error and phase unwrapping errors, when the coherence is poor in the interested region. In addition, the onset time of a loess landslide is successfully uncovered over the Xingyuan landslide. This research method can be referred for any other similar loess landslide monitoring.

Continuous increase of preslide deformation rate is monitored for the Xingyuan landslide, and one month later, a landslide occurred, which indicates that preslide deformation rate is a crucial precursor for landslide early warning. Instability process of loess landslide must be carefully analyzed by integrating data from hydrology, geology, and irrigation.

## Acknowledgments

TerraSAR data are supplied by DLR under Grant No. LAN2471. This research is funded by the National Program on Key Basic Research Project (973 Program) (Grant No. 2014CB744703) and the Natural Science Foundation of China (Grant Nos. 41372375 and 41304016).

## References

## Biography

**Chaoying Zhao** is a professor at the Chang’an University, China. He received his master degree and his PhD in geodesy from Chang’an University in 2002 and 2009, respectively. His current research interests include different InSAR methods and their applications in geo-hazard monitoring, such as land subsidence, ground fissures, landslide, and mining-induced collapse.

**Qin Zhang** is a professor at Chang’an University. She received her PhD in geodesy from Wuhan University in 2002. Her research interests include GPS theory, data processing, and geo-hazard monitoring applications.

**Yang He** is a master student; he is mainly engaged in photogrammetry and remote sensing, and the application of InSAR.

**Jianbing Peng** is a professor at Chang’an University. He obtained his PhD in engineering geology from Chang’an University in 1999. His current research interests include the mechanism of loess landslide and ground fissure.