Landslide Movement Monitoring with InSAR Technologies

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.

1. Introduction

Landslides represent a major geological hazard causing injury or death and causing economic loss [1, 2, 3, 4, 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, 8, 9, 10, 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, 13, 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, 18, 19]. With advances in algorithm development, InSAR has been widely used to monitor slope movement [20, 21, 22], update landslide inventory maps [11, 23], and generate landslide risk maps [24, 25]. Multi-source space-air-ground 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 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.

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

Daguan County Town is built on an ancient landslide (Figure 2(a)). The altitude ranges from $500m$ to $1500m$ and the angle of slope from ${10}^{\circ }$ to ${60}^{\circ }$. 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 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.

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

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, 37, 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]:

$$\mathit{y}=\mathit{A\gamma }$$E1

where $\mathit{y}={\left[y_1,\dots ,y_N\right]}^{\mathrm{T}}$ (${\left(\cdot \right)}^{\mathrm{T}}$ is the transpose operation) represents the complex values of interferograms after removing the APS, $\mathit{\gamma }$ is the reflectivity vector and $\mathit{A}$ is the sensing matrix containing the steering vector $\mathit{a}\left(\mathrm{\Delta }h,\mathrm{\Delta }v\right)$ as its columns:

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

where ${\xi }_i=2b_1/\left(\lambda r_{0}sin\beta \right)$ ($b_i$ is the perpendicular baseline, $\lambda$ is the wavelength, $r_0$ is the slant range and $\beta$ is the incidence angle) and ${\eta }_i=2t_i/\lambda$ ($t_i$ is the temporal baseline). $\mathrm{\Delta }h$ and $\mathrm{\Delta }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 \left(\mathrm{\Delta }h,\mathrm{\Delta }v\right)=\frac{\left|a{\left(\mathrm{\Delta }h,\mathrm{\Delta }v\right)}^{\mathrm{H}}y\right|}{{‖a\left(\mathrm{\Delta }h,\mathrm{\Delta }v\right)‖}_2{‖y‖}_2}$$E3

where $\left|\cdot \right|$ is modulus operation for each element, ${‖\cdot ‖}_2$ is 2-norm, and ${\left(\cdot \right)}^{\mathrm{H}}$ stands for the conjugate-transpose operation. The maximum value of $\gamma$ is calculated by sampling $\mathrm{\Delta }h$ and $\mathrm{\Delta }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:

$$\mathrm{\Delta \Phi }=\mathit{DJ}=\left[\begin{array}{cc}2{\mathit{\pi \xi }}_1& 2{\mathit{\pi \eta }}_1\\ ⋮& ⋮\\ 2{\mathit{\pi \xi }}_N& 2{\mathit{\pi \eta }}_N\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }h\\ \mathrm{\Delta }v\end{array}\right]$$E4

where $\mathrm{\Delta }\mathtt{\Phi }$ 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]:

$$\mathit{J}={\left({\mathit{D}}^{\mathrm{T}}{\mathit{W}}_M\mathit{D}\right)}^{-1}{\mathit{D}}^{\mathrm{T}}{\mathit{W}}_M\mathrm{\Delta }\mathtt{\Phi }$$E5

where ${\mathit{W}}_M$ 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 $\sigma \mathit{I}$ ($\mathit{I}$ is the identity matrix) in the inversion:

$$\mathit{X}={\left({\mathit{G}}^{\mathrm{T}}{\mathit{W}}_R+\sigma \mathit{I}\right)}^{-1}{\mathit{G}}^{\mathrm{T}}{\mathit{W}}_R\mathit{H}$$E6

where $\mathit{X}$ contains the absolute estimates of PS, $\mathit{G}$ is the conditioning matrix consisting of −1,0 and 1 (−1 and 1 represent the start and end of the arc, respectively), ${\mathit{W}}_R$ is a diagonally weighted matrix containing the temporal coherence, and $\mathit{H}$ contains the arc on relative estimates. The adjustment parameter $\sigma$ 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) $\mathbf{C}$ could be calculated. Then the inversion of $\mathbf{C}$ could be used to change the optical phase by the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. It should be noted that $\mathbf{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]:

$$\mathit{\theta }=arg\underset{\theta }{max}\left({\mathtt{\Lambda }}^{\mathrm{H}}\left(\mathit{C}\circ \mathtt{\Psi }\right)\mathtt{\Lambda }\right)=arg\underset{\theta }{max}\left({\mathtt{\Lambda }}^{\mathrm{H}}\mathit{C}\mathtt{\Lambda }\right)$$E7

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 Y-axis is north, and the Z-axis is up. The deformation of a target on the earth’s surface is:

$$U=U_x{s}_x+U_y{s}_y+U_z{s}_z$$E8

where $U_x$, $U_y$ and $U_z$ are the eastern, northern, and vertical components of $U$ (the real movement of a landslide), and $s_x$, $s_y$, and $s_z$ 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 $U_y$ is negligible. $U_A$ and $U_D$ were used to represent the LOS deformation velocity for the ascending and descending, respectively. It should be noted that we assumed that the 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:

$$\left\{\begin{array}{c}{U}_a\approx U_x{a}_x+U_z{a}_z\\ U_d\approx U_x{a}_x+U_z{a}_z\end{array}\right.$$E9

where $a_x$ and $a_z$ 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\left(p_1,p_2,\dots ,p_m\right)=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{U}_i^{\ast }-U_i}$$E10

where $\left(p_1,p_2,\dots ,p_m\right)$ is the material properties set to be defined, $U_i^{\ast }$ and $U_i$ are simulated and InSAR measurement results, respectively. $\left|\cdot \right|$ 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, 54, 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.

Figure 3.

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

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 precipitation-induced 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 $20m\times 20m$ 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 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}^{\circ }$ 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{\mathit{km}}^2$) except S21. There are some landslides without sufficient points or without significant movements (larger than $5\mathit{mm}/\mathit{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].

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

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.

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

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.

Inclinometer tubes were installed to measure the deformation at different depths. The oscillation range of these two tubes was $\pm {15}^{\circ }$. We can obtain the deformation of I1 at depth of $1m$, $15m$ and $30m$, the deformation of I2 at depth of $1m$, $13m$ and $25m$. The measurement result provided two direction movement, x-direction 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 $1m$ to validate InSAR measurement results. The InSAR results agreed well with I1, but not with I2. The difference is $4.9\mathit{mm}/\mathit{yr}$ and $31.2\mathit{mm}/\mathit{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 $15m$. 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\mathit{mm}/\mathit{yr}$ and the corresponding material parameters were assumed to be optimal.

Material property	E-Modulus Max depth (m)	K-Modulus	N-exponent	k(0)	Poisson’s ratio	Pore water pressure (KPa)	Cohesion (MPa)	Friction angle (degree)
Range	[11, 20]	[160, 250]	[0.1, 0.5]	[1.1, 2]	[0.1, 0.4]	[−600, 0]	0.2	20
Precision	1	10	0.1	0.1	0.05	100	0	0
Result	15	200	0.2	1.2	0.2	−200	0.2	20

Table 1.

Main material properties estimated by the GA.

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\mathit{mm}/\mathit{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 Technology building, Gengyun Building, and Zizhi Building) along the cross-section 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\mathit{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.

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.

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, 61, 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 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\mathit{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.

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

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\mathit{mm}/\mathit{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\mathit{mm}$. The movements became small at the landslide toe, because retaining walls were used to maintain stability of Zhaoma Highway therein.

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.

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

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

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:

1. The scatterers detected from SAR images are mainly located in the buit-up areas of Daguan County. Because of the geometric distortion and decorrelation, the scatterers were relatively sparse at lower and higher elevations. 24 of 31 landslides were identified successfully based on the InSAR measurement results. Besides, there was good agreement between InSAR and GNSS and inclinometer data, and the few large errors may have been selected for comparison due to inconsistencies between them. In general, the InSAR measurement results could reflect the landslide movement of this area accurately.

2. By coupling InSAR measurement, we used GA to calculate the optimal material properties. The objective function was minimized to $1.3\mathit{mm}/\mathit{yr}$. The numerical modeling described the full-scale landslide motion from surface to bottom, which could help us to understand the physical process of landslide. The simulated results implicated that the the head of S12 suffered a maximum cumulative deformation because of the external loading caused by buildings of Daguan No.1 Middle School. The maximum deformation of S15 occurred at the toe. This may be due to the reduction of the effective modulus of elasticity of the stratum after soil washing.

3. In the numerical modeling, three boundary conditions are set as external loading. It becomes feasible to evaluate the independent effects of each boundary condition. The results showed that the main load dominates the 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.

Acknowledgments

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