Application of New Technologies

#### **Chapter 10**

## Landslide Movement Monitoring with InSAR Technologies

*Peifeng Ma, Yifei Cui, Weixi Wang, Hui Lin, Yuanzhi Zhang and Yi Zheng*

#### **Abstract**

Synthetic aperture radar interferometry (InSAR) is a technology that has been widely used in many areas, such as topographic mapping, land and resource survey, geological exploration, disaster prevention and mitigation, volcanic and seismic monitor and so on. Landslide, as a representative geohazard, include a wide range of phenomena involving downhill ground movement. InSAR, a technology which can measure surface deformation at the millimeter level over serveral days or years, is suitable to detect landslides with chronical and widespread movements. In this chapter, we introduce main process methods of InSAR data, including Persistent Scatter Interferometry (PSInSAR) and Distributed Scatter Interferometry (DSInSAR). A study area, Daguan County Town, one of the most landslide-prone areas in China is induced to demonstrate the practicability of InSAR in detecting landslides. Combined InSAR results with geological, geotechnical and meterological data, the distribution of landslide in Daguan County in spatial and temporal dimensions would be displayed. We also coupling numerical modeling and InSAR for characterizing landslide movements under multiple loads. The numerical results revealed that body loads dominated the cumulative downhill movements by squeezing water and air from voids, and precipitation caused seasonal movements with the direction perpendicular to the slope surface.

**Keywords:** InSAR, landslide movement, numerical modeling, spatial–temporal analysis, human activities

#### **1. Introduction**

Landslides represent a major geological hazard causing injury or death and causing economic loss [1–5]. More than 20,000 hazards associated with landslides occurred in China from 2013 to 2014, totally causing 10 billion CNY loss approximately [6]. Landslides tend to occurred in areas where natural or human activities are more frequent [7–11]. Exposing the causes of landslides is essential for characterizing slope mobility and hazard mitigation. The traditional method of landslide investigation is to conduct field surveys [10], but it is laborious and difficult to find the landslide boundaries. Synthetic aperture radar interferometry (InSAR) monitoring, based on satellite, is now an effective way for landslide deformation monitoring from regional to local scales [12–14]. In particular, the launch of Sentinel-1 satellites opens a new era of global coverage monitoring with revisit times of 12 days using one satellite and 6 days using two, making regular landslide monitoring feasible [15, 16]. Advanced

multitemporal InSAR methods have been proposed to achieve centimeter-to-millimeter deformation monitoring using multi-baseline SAR images by detecting temporally coherent targets [17–19]. With advances in algorithm development, InSAR has been widely used to monitor slope movement [20–22], update landslide inventory maps [11, 23], and generate landslide risk maps [24, 25]. Multi-source space-airground monitoring technologies are also combined with InSAR to study landslide movements to overcome the limitations of single technology [11, 26].

Landslides in built-up areas should first be considered for monitoring, comparing with landslides in non-built-up areas. Human activities, e.g., land development, are frequent in these areas, so monitoring is necessary to prevent greater economic losses and casualities. In earthquake area zone of southeast China, with scarce land resources, the construction of new cities is often carried out on steep slopes. In these areas, infrastructures and buildings are at high risk since the land urbanization may destroy the original ecological environment [7]. InSAR can measurement surface deformation and the landslides boundaries can be confirmed through vertical deformation velocity. However, this method can only monitor the movement on surface. To reveal the physical mechanisms of landslide, the numerical modeling should be introduced to simulate physical process underground. Numerical modeling is purely mathematical and thus is very different from physical modeling in the laboratory or field modeling [27]. The models are an abstract of the real geology object, so the geometries are usually simplified. Besides, the material properties of the simulated strata are retrieved through back analysis method based on InSAR measurements. The phenomena of landslides could be simulated after the definition of model, using Finite Element Analysis (FEA) method [28, 29].

The coupling analysis of numerical modeling and InSAR measurements have been proved an effective method to deduce the deformation of dams [30]. In this chapter, we review an experimental study [31] to show the performance of coupling numerical modeling and InSAR for monitoring slope. Daguan County Town, a built-up area located on a hillside, is considered as our study area. An advanced multi-temporal InSAR method would be used to measure surface movements with robust detection of persistent scatterers (PSs) and distributed scatterers (DSs). Then, the detection results would be validated by ground data from Global Navigation Satellite System (GNSS) stations and inclinometer tubes. The coupling analysis of numerical modeling and InSAR measurement results would reveal the mechanisms of landslides, and influence of three main factors (e.g., precipitation, body loads and construction-induced loads) would be discussed.

#### **2. Study area and datasets**

Daguan County, located in the northeast of Yunnan and southwest of the Sichuan earthquake fault zone (**Figure 1**), is one of the most landslide-prone areas in China [6]. From 1844 to 1974, nine earthquakes with magnitude greater than 5 occurred in Daguan County. Earthquake and fault activity have led to a wide distribution of ancient landslides in this area [32, 33]. The region has a subtropical monsoon climate, with the rainy season normally lasting from June to September [6]. The historical geological evolution of the county consists of three main periods of movements (i.e., the Jinning Movement, the Yanshan Movement and Jialing River Movement), which formed folds and fractures. The geological strata range from the oldest to the youngest include the Silurian (consisting of the Lower Huanggexi, Middle Shichuankan and Upper Caidianwan formations), the Devonian (consisting of the Lower Cuifengshan and Middle Qingmen formations) and the Quaternary (consisting of artificial, remnant slopes, fallen rocks and alluvial

#### *Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

#### **Figure 1.**

*(a) Study area and SAR images coverage. Rectangle a and D indicate the coverage of the ascending and descending Sentinel-1 images, respectively. The black rectangle indicates the location of Daguan. (b) the geological environment of Daguan County town. (c) Landslide inventory map. S1-S22 are the 22 landslides, D1-D5 are the 5 debris flows and R1-R4 indicate 4 rockfalls. The polygon colors indicate: (1) large and high risk; (2) large and medium risk; (3) medium and high risk; (4) medium and medium risk; (5) small and high risk; (6) small and medium risk landslides; (7) high risk and (8) medium risk debris flows; and (9) high and (10) medium risk rockfalls [31].*

deposits) [33]. Accompanied by frequent crustal movements, the rocks in the area are extremely fragmented and therefore easily infiltrated by rain and groundwater. Prolonged rainfall can saturate, soften, and erode the soil, presumably leading to slope failure. The main soft and hard interaction surfaces that are prone to landslides include artificial and rockfall layers, joint fractures and relief joints, overlying soft soils and underlying bedrock, soft and hard rocks, and the interaction between different weathering surfaces. Totally speaking, the steep topography, tectonic action and stratigraphic lithology together shape the landslide-prone environment.

Daguan County Town is built on an ancient landslide (**Figure 2(a)**). The altitude ranges from 500*m* to 1500*m* and the angle of slope from 10<sup>∘</sup> to 60<sup>∘</sup> . Its main aspect of the slope is to the west. Human activity has significantly altered the landscape and topography, with buildings and infrastructure densely located in this area. In order to promote sustainable development, the Daguan County government began rezoning and revising the land use in 2010, and the land has been rapidly developed over the past few years. Construction works have triggered instability and response of ancient landslides by altering above-ground loads and subsurface geological and hydrological conditions [34, 35]. The Yunnan Geotechnical Engineering Survey and

**Figure 2.**

*(a) Photo of the Daguan County town. (b) Stair buckling in the region of S12. (c) Wall cracking in the region of S15. Stratigraphic cross-sections of S12 (d) and S15 (e) [31].*

Design Institute conducted a field survey of the geological environment and distribution of geological hazards in the area in 2011 and generated an inventory map based on its identification of 31 landslides, including 22 landslides, 5 debris flows and 4 rockfalls (**Figure 1(c)**). **Figure 2(d)** and **(e)** are stratigraphic cross-sections of the two most severe landslides (e.g., S12 and S15). The geological setting consists of Silurian marl remants (*S*2*d*) and artificial fill, and the sliding soils are loose Quaternary gravelly soils. The slides of S12 and S15 caused buckling of the stairs and cracking of the walls. Daguan No.1 Middle School is located close to the top of the landslide in S12 (**Figure 2(d)**), and a building (e.g., Genyun Building) here was rebuilt in May 2018. In S15, near the toe of the slope, is the retaining wall (**Figure 2(e)**), whose construction began in June 2018, to dampen movement by applying lateral foces. These construction works may affect the mobility of the landslide.

The SAR images is from Sentinel-1 in the wide swath imaging mode. The monitoring coverage is as shown in **Figure 1(a)**. A total of 60 ascending images from February 12, 2017 to February 26, 2019 and 56 descending images from February 19, 2017 to March 5, 2019.

#### **3. Methods**

#### **3.1 Deformation estimation of PS and DS points in a two-tier network**

#### *3.1.1 Detection of the most reliable PS in the first-tier network*

The problem of decorrelation usually occur in landslide areas. To solve this problem, we used a two-tier network to detect PSs and DSs with robust deformation estimation [36–38], which can help improve the measurement points density. The first-tier network was constructed to find the most reliable PSs over the study area. Depend on amplitude dispersion, the primary PSs candidates (PPSC) could be selected firstly [39]. A Delaunay Triangulation Network (DTN) could then be constructed using PPSC. Some additional redundant arcs should be added to improve the density of the DTN, because of the sparse PPSC [37]. After removing the atmospheric phase screen (APS), the signal model of *N* observations was written like [40, 41]:

$$
\gamma = A\gamma \tag{1}
$$

where *y* ¼ *y*1, … , *yN* � �<sup>T</sup> (ð Þ� <sup>T</sup> is the transpose operation) represents the complex values of interferograms after removing the APS, *γ* is the reflectivity vector and *A* is the sensing matrix containing the steering vector *a*ð Þ Δ*h*, Δ*v* as its columns:

$$\mathfrak{a}(\Delta h, \Delta v) = \left[ \exp \left( j2\pi (\xi\_1 \Delta h + \eta\_1 \Delta v) \right), \dots, \exp \left( j2\pi (\xi\_N \Delta h + \eta\_N \Delta v) \right) \right]^T \tag{2}$$

where *ξ<sup>i</sup>* ¼ 2*b*1*=*ð Þ *λr*<sup>0</sup> sin *β* (*bi* is the perpendicular baseline, *λ* is the wavelength, *r*<sup>0</sup> is the slant range and *β* is the incidence angle) and *η<sup>i</sup>* ¼ 2*ti=λ* (*ti* is the temporal baseline). Δ*h* and Δ*v* are relative height and mean deformation velocity to be determined, respectively. We used beamforming and an M-estimator to calculate the relative estimates [36]. The beamforming-based inversion is as follows:

$$\gamma(\Delta h, \Delta v) = \frac{\left| a(\Delta h, \Delta v)^{\mathsf{H}} \mathcal{Y} \right|}{||a(\Delta h, \Delta v)||\_2 ||\mathcal{Y}||\_2} \tag{3}$$

where j j� is modulus operation for each element, k k� <sup>2</sup> is 2-norm, and ð Þ� <sup>H</sup> stands for the conjugate-transpose operation. The maximum value of *γ* is calculated by sampling Δ*h* and Δ*v* to describe the temporal coherence of points. The arc between points would retain if the corresponding temporal coherence is larger than a given threshold; otherwise, the arc would be removed. The threshold we set here is 0.72, typical for millimeter-level deformation estimation [42]. For the preserved arcs, the preliminary estimates were used to unwrap the temporal phase. Then the inversion problem is transformed to:

$$
\Delta \Phi = \mathbf{D} \mathbf{J} = \begin{bmatrix} 2\pi \xi\_1 & 2\pi \eta\_1 \\ \vdots & \vdots \\ 2\pi \xi\_N & 2\pi \eta\_N \end{bmatrix} \begin{bmatrix} \Delta h \\ \Delta \nu \end{bmatrix} \tag{4}$$

where ΔΦ is the unwrapped phase. In the presence of low-quality images, the preliminary relative estimates may be biased. To address this issue, we introduced an M-estimator to re-calculate the estimates using the unwrapped phase [43]:

$$\mathbf{J} = \left(\mathbf{D}^{\mathrm{T}} \mathbf{W}\_{M} \mathbf{D}\right)^{-1} \mathbf{D}^{\mathrm{T}} \mathbf{W}\_{M} \Delta \boldsymbol{\Phi} \tag{5}$$

where *W<sup>M</sup>* is an iteratively computed weight matrix that represents the phase quality. Compared to the unweighted least square estimate, the M-estimate assigns smaller weights to the phase outliers, thus improving the robustness of the estimate. After solving the relative estimates, we integrate them with the network adjustment method. When the adjustment matrix is poorly conditioned, its inversion may be unstable. To solve this problem, we apply a ridge estimator to the network

adjustment [44]. It introduces a conditioning matrix *σI* (*I* is the identity matrix) in the inversion:

$$\mathbf{X} = \left(\mathbf{G}^{\mathrm{T}}\mathbf{W}\_{R} + \sigma\mathbf{I}\right)^{-1}\mathbf{G}^{\mathrm{T}}\mathbf{W}\_{R}\mathbf{H} \tag{6}$$

where *X* contains the absolute estimates of PS, *G* is the conditioning matrix consisting of �1,0 and 1 (�1 and 1 represent the start and end of the arc, respectively), *W<sup>R</sup>* is a diagonally weighted matrix containing the temporal coherence, and *H* contains the arc on relative estimates. The adjustment parameter *σ* is determined according to the L-curve method [44]. The ridge estimator outperforms the traditional least-square estimator in regulating possible ill-conditioned problems. The effectiveness of the M-estimator and the ridge estimator has been evaluated in [36] and is omitted here for simplicity.

#### *3.1.2 Detection of the remaining PS and all the PS in the second-tier network*

The PSs detected in the first-tier network were regarded as reference points to build the second-tier network. Identifying statistically homogeneous pixels (SHPs) by Kolmogorov–Smirnov (KS), the complex covariance matrix (CCM) **C** could be calculated. Then the inversion of **C** could be used to change the optical phase by the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. It should be noted that **C** is rank deficient, and the inversion is not stable when the amount of SHPs is less than *N*. To solve this problem, we have revised it as following [19, 45]:

$$\boldsymbol{\theta} = \arg\max\_{\boldsymbol{\theta}} \left( \boldsymbol{\Lambda}^{\mathrm{H}} (\mathbf{C} \circ \boldsymbol{\Psi}) \boldsymbol{\Lambda} \right) = \arg\max\_{\boldsymbol{\theta}} \left( \boldsymbol{\Lambda}^{\mathrm{H}} \mathbf{C} \boldsymbol{\Lambda} \right) \tag{7}$$

This expression improves the robustness of estimation by assigning a larger weight to the higher-coherence phase and avoiding matrix inversion. Then, we employed a more efficient phase-linking method to obtain the optimal phase [46]. This method is called a coherence-weighted phase-linking (CWPL) method [45]. Finally, we used the reconstructed optimal phase to identify whether the DSC is a true DS using temporal coherence (T\_DS) thresholding (0.65 in this case). The workflow we can see on the top left of **Figure 3**.

#### **3.2 2D deformation velocity decomposition**

Landslides are always described as a movement with a predominantly vertical orientation. However, a single track only provides deformation in the LOS direction. To describe the movement of landslides properly, we appreciated that the up-down and east–west components of the deformation could be calculated using the observations obtained from both ascending and descending orbits [47]. Suppose that in the Cartesian coordinate system, the direction of the X-axis is east, the direction of the Yaxis is north, and the Z-axis is up. The deformation of a target on the earth's surface is:

$$U = U\_x \mathfrak{s}\_x + U\_y \mathfrak{s}\_y + U\_x \mathfrak{s}\_x \tag{8}$$

where *Ux*, *Uy* and *Uz* are the eastern, northern, and vertical components of *U* (the real movement of a landslide), and *sx*, *sy*, and *sz* are unite vectors in the respective directions [48]. Because of the polar-orbit of Sentinel-1, the LOS deformation is insensitive to movement in the north–south direction, and the *Uy* is negligible. *UA* and *UD* were used to represent the LOS deformation velocity for the ascending and descending, respectively. It should be noted that we assumed that the *Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

**Figure 3.** *Workflow of PS and DS detection, 2D deformation decomposition, and numerical modeling [31].*

mean velocity of the landslide motion is the same for both types of orbit images, since their time spans are similar. Therefore, the projections of the east–west and up-down motions in the LOS direction could be written as:

$$\begin{cases} U\_d \approx U\_x a\_x + U\_x a\_x \\ U\_d \approx U\_x a\_x + U\_x a\_x \end{cases} \tag{9}$$

where *ax* and *az* represent the unit LOS vectors obtained from the orbit parameters (**Figure 3**).

#### **3.3 Numerical modeling of landslide movements using finite element analysis**

In order to describe the underground process of landslide movement, we performed numerical model of landslide movement in GeoStudio software [49, 50]. SIGMA/W is a tool in GeoStudio that can perform stress and deformation analyses of earth structures, we can use it to simulate the physical process of ground volume change in response to self or external loading. Before solving, three components should be identified, which are geometry, material properties and boundary conditions. The geometry is the cross-section of a slope, it can be defined by simplifying the stratigraphic data.

The settings of material properties is crucial to make the model we defined close to the reality. We decided the material properties by combining the empirical knowledge and back analysis method [30, 51]. The purpose of Back analysis is to make the deformation of numerical modeling and InSAR measurement results consistent by modifying material properties. We first defined an objective function using the 2-norm of the difference between the modeled and measured deformation:

$$f(p\_1, p\_2, \dots, p\_m) = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} U\_i^\* - U\_i} \tag{10}$$

where *p*1, *p*2, … , *pm* � � is the material properties set to be defined, *U* <sup>∗</sup> *<sup>i</sup>* and *Ui* are simulated and InSAR measurement results, respectively. j j� is the absolute value operation and *n* means the number of InSAR scatterers located at the landslide. We then used the genetic algorithm (GA) to iteratively search for the optimal material properties [52]. The range and precision of soil properties were initialized in the GA based on a empirical knowledge [53–55]. The optimal set of property is the combination of material parameters by minimizing the objective function. If the objective function was not minimized, the GA is subjected to natural selection and mutation with crossover and mutation probability set to 0.9 and 0.1, respectively. Considering the case of non-convergence, we set an upper limit on the number of iterations to 50, which is also considered optimal when it is reached.

The boundary conditions were used to replicate real-world loading in deformation simulation. Considering the actual forces of the study area, three kinds of boundary condition are introduced in modeling. The first one is body load, which represents the body gravity related to the volume of the element and material, and it can be calculated by the unit weight and volume of soils. The second one is precipitation-induced hydraulic condition. The precipitationinduced hydraulic condition was regarded as a function varied with time. The input precipitation data was collected during the acquisition time of Sentinel-1 images. Amplitude of InSAR seasonal deformation varies with altitude under the influence of precipitation accumulation and stratified atmospheric delay [56, 57]. To calibrate this effect, we set different rainfall intensity at different elevations. The third boundary condition is construction-induced force (e.g., building loads). The building pressure was estimated by the concrete-steel density, thickness of walls and number of floors in the buildings obtained from field survey. After setting boundary conditions, we used the linear elastic model to simulate the deformation [58]. Linear elastic model assumes that stress is proportional to strain, and the load–displacement response is likely linear elastic along the lower initial portion of the stress–strain curve (**Figure 3**).

#### **4. Experimental results**

#### **4.1 InSAR results and validation**

#### *4.1.1 2-D deformation mapping*

Before decomposition of deformation, we identified the common points detected from ascending and descending tracks. The InSAR results are converted to raster data at 20*m* � 20*m* resolution, the value of each grid cell was taken as the mean value of all points within it, and the position was in the center of each grid cell. **Figure 4(a)** and **Figure 4(b)** present the east–west and up-down direction deformation velocity map, respectively. The positive and negative velocities in **Figure 4(a)** indicate the eastward and westward motion, respectively, whereas the positive and negative velocities in **Figure 4(b)** indicate vertical uplift and subsidence, respectively. Sparse points were observed in the area with low elevation close to the Daguan River, caused mainly by slope-induced shadow and layover problems. High-elevation area also displayed a sparse distribution of points, caused

#### *Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

**Figure 4.**

*Spatial pattern of InSAR measurements. (a) and (b) deformation velocitity in east–west and vertical direction, respectively. The black triangles indicate the locations of the reference point and the white dashed line marks the boundary of the moving and stable areas. The polygons in the map indicate the landslide boundaries from the inventory map, and two landslides with blue boundaries are S12 and S15 [31].*

mainly by vegetation-induced decorrelation effects. The results demonstrated that the moving direction was generally westward and downward, following the downhill direction. The westward movement was more significant than the downward movement because the slope angles of most positions were less than 45<sup>∘</sup> and the horizontal projection of downhill movements was larger than the vertical projection. The overlap between the landslide polygons and 2D deformation suggests that the movements generally agreed well with the landslide areas. However, the spatial extents of landslides and monitored movements differed, indicating a possible change in landslide mobility from the time of the field survey (2011) to the acquisition time of the Sentinel-1 images (2017–2019). There was a sharp boundary between the moving and stable areas in the north of the study area (**Figure 4**). We, therefore, infer a fault here. Considering that Fenping Fault crossed the study area in **Figure 1(b)**, we conclude that the sharp boundary indicates the location of Fenping Fault. The Fault acts as barriers to groundwater flow, resulting in large differential deformation between the two sides of the fault. Considering the types of landslide may influence the interpretation of measurement results and subsequent numerical modeling, we confirmed that all landslides in this area are translational slides from field surveys. The polygon with more than 10 InSAR points within the boundary would be considered to be detected successfully. So, 24 of 31 landslides in the inventory map were successfully detected. The undetected landslides are almost small-sized (smaller than 0*:*01*km*<sup>2</sup> ) except S21. There are some landslides without sufficient points or without significant movements (larger than 5*mm=yr*), but we cannot guarantee that they are stable because rapidly moving targets may be shadowed or exceed the detection limitation of InSAR methods [59].

#### *4.1.2 Validation of InSAR measurements using ground data*

As the yellow labels shown in **Figure 5**, there are three global navigation satellite system (GNSS) stations (G1, G2 and G3) and two inclinometer tubes (I1 and I2) in the study area. We selected neighboring points from InSAR measurement results and ground monitor stations for comparison.

The measurements from GNSS are in three directions. We selected GNSS measurements along east–west and up-down directions to compare with InSAR 2-D deformation results. The linear deformation velocity was derived by linear fitting function. It should be noticed that there are missing data at G2 and G3, so we extend the time span of fitted data to July 2019 for them. The results implicated that the movements in horizontal were more significant than that in vertical movements. We can see from the mean deformation velocity that most InSAR results agreed well with the GNSS results except the horizontal movement of G2. That may be due to the fact that the InSAR point selected for comparison was located on a relatively stable structure rather than moving ground.

**Figure 5.**

*Comparison between InSAR and ground data. (left top) location of GNSS stations and inclinometer tubes. (right top) comparison between InSAR and inclinometer data. (bottom) comparison between InSAR and GNSS data, horizontal movement in the upper row and vertical movement in the lower row [31].*

*Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

Inclinometer tubes were installed to measure the deformation at different depths. The oscillation range of these two tubes was <sup>15</sup><sup>∘</sup> . We can obtain the deformation of I1 at depth of 1*m*, 15*m* and 30*m*, the deformation of I2 at depth of 1*m*, 13*m* and 25*m*. The measurement result provided two direction movement, xdirection is the tangential direction of the ground, y-direction is the vertical direction. We assumed that x-direction movement is along western because it is the same as the slope orientation, and the x-direction movement of I1 and I2 at different depth from August 2018 to February 2019 are shown on **Figure 5** (Right top). We use the deformation velocity at the depth of 1*m* to validate InSAR measurement results. The InSAR results agreed well with I1, but not with I2. The difference is 4*:*9*mm=yr* and 31*:*2*mm=yr*, respectively. The reason is that the InSAR scatterer selected for comparison may be located on a relatively stable structure. After validation, we can reasonably integrate them into numerical modeling for geological parameter retrieval.

#### **4.2 Numerical modeling results**

Two cross-sections of slopes (S12 and S15) are selected for numerical modeling analysis since they are typical slopes in this area and had been mapped by field survey. The size of each element of models was set to 15*m*. The InSAR measurement results of S12 and S15 were coupled to estimate the optimal material properties using GA. **Table 1** shows the range, precision and the optimal result of eight main parameters we used. The external stress loading was set as a function with four variables, and they are Max depth, K-Modulus, N-exponent and k(0). Cohesion and friction angle were set according to empirical knowledge. The minimized deformation is 1*:*3*mm=yr* and the corresponding material parameters were assumed to be optimal.

#### *4.2.1 Numerical modeling of landslide S12*

We combined the 2-D deformation to calculate the downhill movements of S12 (**Figure 6(a)**). The maximum combined movement was 23*:*2*mm=yr* and the direction was generally consistent with the downhill direction. Numerical modeling was conducted to derive the movements of S12. By iteratively searching for the optimal soil properties using GA, surface deformation by numerical modeling became consistent with measured movements by InSAR, and then we assumed that the material property set was valid. The cumulative deformation by numerical modeling is shown in **Figure 6(c)**. InSAR can measure only surface deformation, whereas numerical modeling depicted full-scale movements of S12 from the surface to the bottom. Deformation at the surface was generally more significant than that at the bottom, consistent with inclinometer data. Daguan No. 1 Middle School is close to the head of the slope, and the weights of three main buildings (Science and


#### **Table 1.**

*Main material properties estimated by the GA.*

#### **Figure 6.**

*(a) Combined deformation velocity of the landslide S12 by InSAR. The background is optical image from Google earth. P1 and P2 are selected for time series analysis in figure 7. (b) Relationship between the thickness of gravel soils and deformation by InSAR and numerical modeling. Two photos show the Gengyun building before and after reconstruction, respectively. (c) Simulated deformation of the sliding layer in numerical modeling. The black and purple arrows indicate the surface downhill movements by InSAR and numerical modeling, respectively. The rectangles b1, b2 and b3 indicate the locations of science and technology building, Gengyun building, and Zizhi building, respectively. The black solid rectangles indicate the retaining walls. P3 is selected for time series analysis in figure 7.*

Technology building, Gengyun Building, and Zizhi Building) along the crosssection were estimated to impose building-induced loads, aggravating landslide movements. In particular, Genyun Building was rebuilt during this time, and led to the maximum cumulative deformation of 34*:*9*mm* in numerical modeling. In contrast, there are no buildings or infrastructures at the toe of S12, and the retaining walls mitigate downhill movements by imposing lateral loading. That makes the movement here less than that of the head, despite the relatively equal thickness of the soil layers (i.e., similar body load). In the central position of S12, the gravel soil is thin due to the location of Cuiping Road, which reduces the body load, and the movement here is relatively small.

#### *Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

Two moving scatterers (P1 and P2 in **Figure 6(a)**) are selected to study the temporal evolution (**Figure 7**). The independent deformation induced by three boundary conditions is also demonstrated. In the ascending track of Sentinel-1, seasonal movement is highly correlated with precipitation with a 1 to 3 month delay. The boundary conditions we defined for Daguan County are body loads, precipitation-induced hydraulic changes and construction-induced forces. Since the model has shown to be close to the measurements, we control for a single boundary condition to discuss the effect of each boundary condition on the landslide and the results are shown in **Figure 7**. The simulated time series deformation under precipitation-induced hydraulic condition showed similar trend with InSAR results (**Figure 7(a)** and **(c)**). We conclude that the seasonal trend is caused by precipitation, which drives movements by changing pore pressures that decrease when surface water cannot infiltrate the landslide body, and increase when surface water infiltrates the landslide body [60–62]. The delay is related to the pore pressure diffusion time since the onset of intense precipitation [63]. Interestingly, seasonal movement was not distinct from the descending track. This is because the LOS direction of the descending track is generally parallel to the slope surface. The simulated results showed that the direction of seasonal deformation was generally perpendicular to the surface (**Figure 7(f)**). Consequently, the LOS deformation from the descending track is insensitive to the seasonal rebound and subsidence. In the descending track, P1 and P2 showed continuous movements away from the sensor, which indicated a continuous downhill movement. Gravel soils were consolidated in response to body loads due to squeezing of water and air from the voids. The resulted movements were continuous and showed a decelerating trend with the increased consolidation (**Figure 7(d)**). Body loads caused larger cumulative deformation than the other two boundary conditions, indicating that it dominates

#### **Figure 7.**

*Time series deformation of P1 and P2 in the (a) ascending (P1-a and P2-a) and (b) descending (P1-D and P2-D) images. Linear fitting is conducted before and after may 2018 for P1-D and P2-D. monthly rainfall data is collected from the National Meteorological Information Center. (c–e) Are time series deformation of P3 numerically modeled by precipitation-induced hydraulic change, body loads, and construction-induced force, respectively. (f–h) Are cumulative deformation numerically modeled by precipitation-induced hydraulic change, body loads, and construction-induced force, respectively [31].*

#### **Figure 8.**

*(a) Combined deformation velocity of the landslide S15 by InSAR. The background is optical image from Google earth. P1 and P2 are selected for time series analysis in figure 9. The blue lines indicate the gullies. (b) Relationship between the thickness of gravel soils and deformation by InSAR and numerical modeling. The photo shows the retaining walls. (c) Simulated deformation of the sliding layer in numerical modeling. The thickness of cross-section has been increased by 4 times for better visualization. The black and purple arrows indicate the surface downhill movements by InSAR and numerical modeling, respectively. The black solid rectangles indicate the retaining walls. P3 is selected for time series analysis in figure 9.*

cumulative downhill movements. The continuous movements of P1-D and P2-D implicated a constant velocity before May 2018 and a subsequent acceleration, suggesting different causes. As described above, the Gengyun Building started reconstruction in May 2018. Construction works caused additional loading associated with softened soils and hydrological change [64]. To model it, we defined that Gengyun Building-induced force started to be effective from May 2018 and the pressure increase from 0 to 48*kPa* gradually based on the calculated building weight. In this sense, Gengyun Building-induced force became a transient boundary condition during the monitoring period. The simulated results showed that P3 was relatively stable at the beginning and moved significantly after May 2018 due to reconstruction works of Gengyun Building (**Figure 7(e)**). This component was added to the movements caused by body loads and yielded an acceleration in time series deformation. Compared with seasonal fluctuation caused by precipitation, construction works caused permanent change of the time series trend.

#### *4.2.2 Numerical modeling of landslide S15*

Landslide S15 is located lower than S12 and has caused wall cracking and fracturing at Daguan Vocational School (**Figure 2(c)**). InSAR points were only present in the upper part of S15 because the lower part is covered by vegetation (**Figure 8(a)**). The maximum combined movement velocity was 21*:*4*mm=yr*. Two buildings were located close to the landslide head. However, because they were constructed many years ago, gravel soils have been consolidated and measured movements were not significant at the head. The relationship between the thickness of gravel soils and deformation implicated that the movements were more significant in the lower part. The effective modulus of elasticity in the lower part of S15 calculated by GA was smaller than that in the upper part. Two gullies formed the landslide boundary (**Figure 8(c)**). In rainy seasons, the rainfall converged to the gullies, ingressing and washing away the soils and decreasing slope stiffness [65, 66]. That decreased the effective modulus of elasticity of the sliding layer. The movements were therefore more significant in the lower part of S15. The largest cumulative deformation is 54*:*7*mm*. The movements became small at the landslide toe, because retaining walls were used to maintain stability of Zhaoma Highway therein.

Similar to S12, the seasonal movements were significant in the ascending Sentinel-1 images and were less distinct in the descending images (**Figure 9**). The seasonal trend of S15 showed a rebound in winter and subsidence in summer, which is opposite to the seasonal variance of S12 and precipitation. That may be caused by the different positions of the reference points when processing the data. The phenomenon that amplitude of seasonal trend is different at different elevations has been studied in [56, 57]. Hu et al. [57] attributes it to different accumulated precipitation at different elevations. The total precipitation accumulated during the wet seasons in the mountains at higher elevation is larger than that in the valleys/basins at lower

#### **Figure 9.**

*Time series deformation of P1 and P2 in the (a) ascending (P1-a and P2-a) and (b) descending (P1-D and P2-D) images. Linear fitting is conducted before and after July 2018 for P1-D and P2-D. monthly rainfall data is collected from the National Meteorological Information Center. (c–e) Are time series deformation of P3 numerically modeled by precipitation-induced hydraulic change, body loads, and construction-induced force, respectively. (f–h) Are cumulative deformation numerically modeled by precipitation-induced hydraulic change, body loads, and construction-induced force, respectively. The thickness of cross-section has been increased by 4 times for better visualization [31].*

elevation. Dong et al. [56] suggests that different amplitude of seasonal trend is caused by stratified atmospheric delay at different elevations. These two factors may both exist and influence the seasonal amplitude. To facilitate numerical modeling, we set different rainfall intensity at different elevations to calibrate the stratification effect of seasonal deformation. The simulated seasonal movement of S15 is present in **Figure 9(c)**. In this study, the elevation of S15 was lower than the reference point (**Figure 4**) and thus, the seasonal movements showed an opposite variance. Body loads dominated cumulative downhill deformation among the three conditions (**Figure 9(g)**), and the time series deformation showed a continuous trend (**Figure 9(d)**). Contrary to the accelerated deformation of S12, P1 and P2 in S15 showed a decelerated trend from July 2018, suggesting that the slope tends to be stable. This is the result of constructing retaining walls (**Figure 8(b)**). In numerical modeling, we configured the retaining walls from July 2018. The retaining walls imposed force opposite to body loads and prevented landslide movements [67]. Theoretically, retaining walls cannot cause deformation. To facilitate modeling, we assumed that it caused deformation opposite to downhill deformation. In this sense, retaining wall-induced force became a transient boundary condition. For S12, the new building construction aggravated downhill movements in time series. For S15, the deformation induced by new retaining walls counteracted downhill deformation induced by body loads, and thus, the total time series deformation tended to be stable after July 2018.

## **5. Conclusions**

On slopes in urban areas, landslides usually have complex loading characteristics, including those influenced by natural and human factors. In this chapter, we reveal the landslide motion in Daguan County through the coupling of InSAR and numerical modeling methods. The coupled approach makes it possible for us to derive full-scale motions. We summarize the main findings as follows:


*Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

> cumulative deformation of the landslide by squeezing the water and air in the void. Precipitation induced seasonal motion in a direction perpendicular to the slope surface, and this seasonal feature only appeared in the ascending SAR images. The human activities led to the permanent changes in the time series trend, both positive and negative effects. In particular, the reconstruction projection of Genyun Building accelerated the downhill movement in S12 after May 2018. In S15, the construction of the retaining walls applying a force opposite to the body load, which reduced the downhill movement after July 2018.

#### **Acknowledgements**

The authors would like to thank European Space Agency (ESA) for providing free and open Sentinel-1 data. This work was supported by the National Natural Science Foundation of China (41971278, 41601356, 42077238, 41941019), the Research Grants Council (RGC) of Hong Kong (CUHK14504219) and the Second Tibetan Plateau Scientific Expedition and Research Program (2019QZKK0903).

#### **Author details**

Peifeng Ma1 \*, Yifei Cui<sup>2</sup> , Weixi Wang3 , Hui Lin<sup>4</sup> , Yuanzhi Zhang<sup>5</sup> and Yi Zheng<sup>1</sup>

1 Institute of Space and Earth Information Science, The Chinese University of Hong Kong, Hong Kong, China

2 State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing, China

3 Research Institute for Smart Cities, School of Architecture and Urban Planning, Shenzhen University, Shenzhen, China

4 School of Geology and Environment, Jiangxi Normal University, Nanchang, China

5 University of Chinese Academy of Sciences, Beijing, China

\*Address all correspondence to: peifengma@cuhk.edu.hk

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **References**

[1] Haque U, Blum P, Da Silva PF, Andersen P, Pilz J, Chalov SR, et al. Fatal landslides in Europe. Landslides. 2016;**13**(6):1545-1554

[2] Intarawichian N, Dasananda S. Frequency ratio model based landslide susceptibility mapping in lower Mae Chaem watershed, Northern Thailand. Environmental Earth Sciences. 2011; **64**(8):2271-2285

[3] Liu D, Cui Y, Guo J, Yu Z, Chan D, Lei M. Investigating the effects of clay/ sand content on depositional mechanisms of submarine debris flows through physical and numerical modeling. Landslides. 2020;**17**(8): 1863-1880

[4] Pazzi V, Morelli S, Fanti R. A review of the advantages and limitations of geophysical investigations in landslide studies. International Journal of Geophysics. 2019;**2019**

[5] Strozzi T, Klimeš J, Frey H, Caduff R, Huggel C, Wegmüller U, et al. Satellite SAR interferometry for the improved assessment of the state of activity of landslides: A case study from the cordilleras of Peru. Remote Sensing of Environment. 2018;**217**:111-125

[6] Wang Q, Li W, Yan S, Wu Y, Pei Y. GIS based frequency ratio and index of entropy models to landslide susceptibility mapping (Daguan, China). Environmental Earth Sciences. 2016;**75**(9):780

[7] Cui Y, Cheng D, Choi CE, Jin W, Lei Y, Kargel JS. The cost of rapid and haphazard urbanization: Lessons learned from the Freetown landslide disaster. Landslides. 2019;**16**(6): 1167-1176

[8] Fan X, Scaringi G, Korup O, West AJ, van Westen CJ, Tanyas H, et al. Earthquake-induced chains of geologic

hazards: Patterns, mechanisms, and impacts. Reviews of Geophysics. 2019; **57**(2):421-503

[9] Guo C, Cui Y. Pore structure characteristics of debris flow source material in the Wenchuan earthquake area. Engineering Geology. 2020;**267**: 105499

[10] Guzzetti F, Mondini AC, Cardinali M, Fiorucci F, Santangelo M, Chang KT. Landslide inventory maps: New tools for an old problem. Earth-Science Reviews. 2012;**112**(1–2):42-66

[11] Rosi A, Tofani V, Tanteri L, Stefanelli CT, Agostini A, Catani F, et al. The new landslide inventory of Tuscany (Italy) updated with PS-InSAR: Geomorphological features and landslide distribution. Landslides. 2018;**15**(1):5-19

[12] Cigna F, Bianchini S, Casagli N. How to assess landslide activity and intensity with persistent Scatterer interferometry (PSI): The PSI-based matrix approach. Landslides. 2013; **10**(3):267-283

[13] Crosetto M, Monserrat O, Cuevas-González M, Devanthéry N, Crippa B. Persistent scatterer interferometry: A review. ISPRS Journal of Photogrammetry and Remote Sensing. 2016;**115**:78-89

[14] Massonnet D, Rossi M, Carmona C, Adragna F, Peltzer G, Feigl K, et al. The displacement field of the landers earthquake mapped by radar interferometry. Nature. 1993; **364**(6433):138-142

[15] Bonì R, Bordoni M, Vivaldi V, Troisi C, Tararbra M, Lanteri L, et al. Assessment of the Sentinel-1 based ground motion data feasibility for large scale landslide monitoring. Landslides. 2020;**17**:2287-2299

*Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

[16] Rucci A, Ferretti A, Guarnieri AM, Rocca F. Sentinel 1 SAR interferometry applications: The outlook for sub millimeter measurements. Remote Sensing of Environment. 2012;**120**:156-163

[17] Ferretti A, Fumagalli A, Novali F, Prati C, Rocca F, Rucci A. A new algorithm for processing interferometric data-stacks: SqueeSAR. IEEE Transactions on Geoscience and Remote Sensing. 2011;**49**(9):3460-3470

[18] Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing. 2001; **39**(1):8-20

[19] Ma P, Wang W, Zhang B, Wang J, Shi G, Huang G, et al. Remotely sensing large-and small-scale ground subsidence: A case study of the Guangdong–Hong Kong–Macao Greater Bay Area of China. Remote Sensing of Environment. 2019;**232**:111282

[20] Bianchini S, Raspini F, Solari L, Del Soldato M, Ciampalini A, Rosi A, et al. From picture to movie: Twenty years of ground deformation recording over Tuscany region (Italy) with satellite InSAR. Frontiers in Earth Science. 2018; **6**:177

[21] Hilley GE, Bürgmann R, Ferretti A, Novali F, Rocca F. Dynamics of slowmoving landslides from permanent scatterer analysis. Science. 2004; **304**(5679):1952-1955

[22] Intrieri E, Raspini F, Fumagalli A, Lu P, Del Conte S, Farina P, et al. The Maoxian landslide as seen from space: Detecting precursors of failure with Sentinel-1 data. Landslides. 2018;**15**(1): 123-133

[23] Righini G, Pancioli V, Casagli N. Updating landslide inventory maps using persistent Scatterer interferometry (PSI). International Journal of Remote Sensing. 2012;**33**(7):2068-2096

[24] Singh LP, Van Westen C, Ray PC, Pasquali P. Accuracy assessment of InSAR derived input maps for landslide susceptibility analysis: A case study from the Swiss Alps. Landslides. 2005;**2**(3): 221-228

[25] Solari L, Bianchini S, Franceschini R, Barra A, Monserrat O, Thuegaz P, et al. Satellite interferometric data for landslide intensity evaluation in mountainous regions. International Journal of Applied Earth Observation and Geoinformation. 2020;**87**:102028

[26] Hu X, Bürgmann R, Schulz WH, Fielding EJ. Four-dimensional surface motions of the Slumgullion landslide and quantification of hydrometeorological forcing. Nature Communications. 2020;**11**(1):1-9

[27] Kavanagh JL, Engwell SL, Martin SA. A review of laboratory and numerical modelling in volcanology. Solid Earth. 2018;**9**(2):531-571

[28] Smith IM, Griffiths DV, Margetts L. Programming the Finite Element Method. Chennai, India: John Wiley & Sons; 2013

[29] Zienkiewicz OC, Taylor RL, Taylor RL, Taylor RL. The Finite Element Method: Solid Mechanics. Vol. 2. Barcelona, Spain: Butterworth-Heinemann; 2000

[30] Zhou W, Li S, Zhou Z, Chang X. Insar observation and numerical modeling of the earth-dam displacement of shuibuya dam (China). Remote Sensing. 2016;**8**(10):877

[31] Ma P, Cui Y, Wang W, Lin H, Zhang Y. Coupling InSAR and numerical modeling for characterizing landslide movements under complex loads in urbanized hillslopes. Landslides. 2021;**18**(5):1611-1623

[32] Chen X, Zhou Q, Ran H, Dong R. Earthquake-triggered landslides in

Southwest China. Natural Hazards and Earth System Sciences. 2012;**12**(2):351-363

[33] Liu X, Wang S, Zhang X. Influence of geologic factors on landslides in Zhaotong, Yunnan province, China. Environmental Geology and Water Sciences. 1992;**19**(1):17-20

[34] Siles G, Trudel M, Peters DL, Leconte R. Hydrological monitoring of high-latitude shallow water bodies from high-resolution space-borne D-InSAR. Remote Sensing of Environment. 2020; **236**:111444

[35] Yan G, Yin Y, Huang B, Zhang Z, Zhu S. Formation mechanism and characteristics of the Jinjiling landslide in Wushan in the three gorges reservoir region, China. Landslides. 2019;**16**(11): 2087-2101

[36] Ma P, Lin H. Robust detection of single and double persistent scatterers in urban built environments. IEEE Transactions on Geoscience and Remote Sensing. 2015;**54**(4):2124-2139

[37] Ma P, Liu Y, Wang W, Lin H. Optimization of PSInSAR networks with application to TomoSAR for full detection of single and double persistent scatterers. Remote Sensing Letters. 2019;**10**(8):717-725

[38] Shi G, Lin H, Bürgmann R, Ma P, Wang J, Liu Y. Early soil consolidation from magnetic extensometers and full resolution SAR interferometry over highly decorrelated reclaimed lands. Remote Sensing of Environment. 2019; **231**:111231

[39] Ferretti A, Prati C, Rocca F. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing. 2000;**38**(5):2202-2212

[40] Stoica P, Moses RL. Introduction to spectral analysis. Pearson. Education. 1997 [41] Zhu XX, Bamler R. Tomographic SAR inversion by *L*\_{1}-norm regularization—The compressive sensing approach. IEEE Transactions on Geoscience and Remote Sensing. 2010; **48**(10):3839-3846

[42] Colesanti C, Ferretti A, Novali F, Prati C, Rocca F. SAR monitoring of progressive and seasonal ground deformation using the permanent scatterers technique. IEEE Transactions on Geoscience and Remote Sensing. 2003;**41**(7):1685-1701

[43] Huber PJ. Robust estimation of a location parameter. In: Breakthroughs in Statistics. New York, NY: Springer; 1992. pp. 492-518

[44] Hansen PC, O'Leary DP. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM Journal on Scientific Computing. 1993; **14**(6):1487-1503

[45] Zhang B, Wang R, Deng Y, Ma P, Lin H, Wang J. Mapping the Yellow River Delta land subsidence with multitemporal SAR interferometry by exploiting both persistent and distributed scatterers. ISPRS Journal of Photogrammetry and Remote Sensing. 2019;**148**:157-173

[46] Guarnieri AM, Tebaldini S. On the exploitation of target statistics for SAR interferometry applications. IEEE Transactions on Geoscience and Remote Sensing. 2008;**46**(11):3436-3443

[47] Rucci A, Vasco D, Novali F. Monitoring the geologic storage of carbon dioxide using multicomponent SAR interferometry. Geophysical Journal International. 2013;**193**(1):197-208

[48] Hanssen RF. Radar Interferometry: Data Interpretation and Error Analysis. Vol. 2. Dordrecht: Springer Science & Business Media; 2001

[49] Krahn J. Stress and Deformation Modeling with SIGMA/W. Alberta,

*Landslide Movement Monitoring with InSAR Technologies DOI: http://dx.doi.org/10.5772/intechopen.105058*

Canada: GEO–SLOPE International, Ltd; 2004

[50] Segerlind LJ. Applied Finite Element Analysis. Vol. 316. New York: Wiley; 1976

[51] Zhang J, Tang WH, Zhang L. Efficient probabilistic back-analysis of slope stability model parameters. Journal of Geotechnical and Geoenvironmental Engineering. 2010; **136**(1):99-109

[52] Houck CR, Joines J, Kay MG. A genetic algorithm for function optimization: A Matlab implementation. Ncsu-ie tr. 1995;**95**(09):1-10

[53] Kulhawy FH, Mayne PW. Manual on Estimating Soil Properties for Foundation Design. Cornell Univ., Ithaca … : Electric Power Research Inst., Palo Alto, CA (USA); 1990

[54] Rawls WJ, Brakensiek D. Estimating soil water retention from soil properties. Journal of the Irrigation and Drainage Division. 1982;**108**(2):166-171

[55] Ritter A, Hupet F, Muñoz-Carpena R, Lambot S, Vanclooster M. Using inverse methods for estimating soil hydraulic properties from field data as an alternative to direct methods. Agricultural Water Management. 2003; **59**(2):77-96

[56] Dong J, Zhang L, Liao M, Gong J. Improved correction of seasonal tropospheric delay in InSAR observations for landslide deformation monitoring. Remote Sensing of Environment. 2019;**233**:111370

[57] Hu X, Wang T, Pierson TC, Lu Z, Kim J, Cecere TH. Detecting seasonal landslide movement within the Cascade landslide complex (Washington) using time-series SAR imagery. Remote Sensing of Environment. 2016;**187**:49-61

[58] Hashash Y, Jung S, Ghaboussi J. Numerical implementation of a neural network based material model in finite element analysis. International Journal for Numerical Methods in Engineering. 2004;**59**(7):989-1005

[59] Bonì R, Bordoni M, Colombo A, Lanteri L, Meisina C. Landslide state of activity maps by combining multitemporal A-DInSAR (LAMBDA). Remote Sensing of Environment. 2018; **217**:172-190

[60] Coe JA, Ellis WL, Godt JW, Savage WZ, Savage JE, Michael J, et al. Seasonal movement of the Slumgullion landslide determined from global positioning system surveys and field instrumentation, July 1998–march 2002. Engineering Geology. 2003;**68**(1–2): 67-101

[61] Cui Y, Chan D, Nouri A. Coupling of solid deformation and pore pressure for undrained deformation—A discrete element method approach. International Journal for Numerical and Analytical Methods in Geomechanics. 2017;**41**(18): 1943-1961

[62] Handwerger AL, Roering JJ, Schmidt DA. Controls on the seasonal deformation of slow-moving landslides. Earth and Planetary Science Letters. 2013;**377**:239-247

[63] Zhao C, Lu Z, Zhang Q, de La Fuente J. Large-area landslide detection and monitoring with ALOS/PALSAR imagery data over northern California and southern Oregon, USA. Remote Sensing of Environment. 2012;**124**: 348-359

[64] Attanayake PM, Waterman MK. Identifying environmental impacts of underground construction. Hydrogeology Journal. 2006;**14**(7): 1160-1170

[65] Thomson S, Tiedemann C. A review of factors affecting landslides in urban areas. Bulletin of the Association of

Engineering Geologists. 1982;**19**(1): 55-65

[66] Wong H, Ho K. The 23 July 1994 landslide at Kwun lung Lau, Hong Kong. Canadian Geotechnical Journal. 1997; **34**(6):825-840

[67] Yingren Z, Shangyi Z. Calculation of inner force of support structure for landslide/slope by using strength reduction FEM [J]. Chinese Journal of Rock Mechanics and Engineering. 2004;**20**

Section 6
