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.²⁰ 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:

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.²² 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.²³ 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.²⁰

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

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

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\mathrm{deg}$ and its relative height between the terrace and the bottom is $\sim 75\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.

Fig. 3

Interferometric SAR configurations for time-series deformation calculation and DEM error estimation.

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.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.

Fig. 4

Seven interferograms contaminated by heavy DEM errors, while the last one is the estimated DEM error map calculated from the other interferograms.

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.²⁵

Fig. 5

Phase unwrapping comparison with different coherence thresholds of the interferogram pair from February 25, 2012 to March 7, 2012. (a) Phase unwrapping map for the pixel coherence larger than 0.15; (b) phase unwrapping map for the pixel coherence larger than 0.5; (c) difference map between (a) and (b); (d) coherence map. On each panel, dashed rectangular indicates the landslide area.

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.

Fig. 6

Four unwrapping interferograms and their coherence maps. Left column is the coherence maps for interferograms with different start date and end date, while right column is their corresponding unwrapping maps, where Xingyuan landslide is indicated with dashed rectangular in each panel.

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.

Fig. 7

Cumulative time-series deformation maps in LOS direction at each SAR acquisition date.

Fig. 8

Cumulative surface deformation within 121 days from November 7, 2011 to March 7, 2012; the red star near Dabaozi village is the reference point, line $AB$ indicates the location of profile in Fig. 10, and three symbols within right fan show the locations of time-series results of landslide, land subsidence, and well, respectively, in Fig. 11. (a) Google image of Dongfeng landslide that occurred in 2012; (b) Google image of the Xingyuan landslide.

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.

Fig. 9

The downslope deformation of the Xingyuan landslide in three-dimensional view in 121 days.

Fig. 10

Xingyuan landslide time-series deformation profile from January 1 to March 7, 2012 along line $AB$ in Fig. 8.

From Fig. 9, we can clearly see the spatial deformation distribution of the Xingyuan landslide. The width of this landslide is $\sim 130\mathrm{m}$ while its length is only 75 m along the main landslide direction, which has high consistence with the field investigation.¹³ 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.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\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\mathrm{mm}/\mathrm{d}$ 55 days prior to landslide to $1.6\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

Fig. 11

Land subsidence and landslide velocity time series before the Xingyuan landslide and underground water-level changes before and after the Xingyuan landslide is demonstrated. The locations of land subsidence, landslide, and well are shown in Fig. 8.

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.