**4. Ultra-rapid freezing**

Ultra-rapid freezing can be considered a midway technique between slow freezing and Vitrification. It is quicker than the slow-freezing technique, does not involve the use of programmable machines and requires lower concentrations of cryoprotectant agents (CPA) than those used in vitrification.

Experimental results demonstrate that this technique has lower performances than slow freezing's and vitrification's ones (AbdelHafez et al., 2010).


**3. Data and Methods**

The methodology applied for landslide susceptibility, hazard, and risk assessment is shown as a flowchart in Figure 2. A detailed explanation is provided in the following subsections.

**Figure 2.** The flowchart of landslide risk assessment in the present study.

#### *3.1. Landslide Inventory*

A landslide inventory map is the first step to assess susceptibility. It shows information on all historical and active landslides. Combined with field surveys, relevant literature records and news reports of landslide records were used in this study to verify the spatial distribution of landslides using Google Earth high-resolution images and an unmanned aerial vehicle's digital orthophoto map (DOM) (0.1 m) provided by the Xi'an Center of China Geological Survey for visual interpretation.

Generally, slope units are defined by topographic characteristic lines (such as ridgelines and gully lines) and waterway paths, which are closely related to the DEM of the mountain area [36]. Therefore, these topographic characteristic lines and waterway paths are also the basic means to determine slope units in this study. Additionally, to better reflect the terrain of landslides or unstable slopes, we divided slope units according to topographic and geomorphic characteristics in the detailed field surveys. The basic requirement in a field survey is that every gully and slope must be investigated and the results are presented in the form of slope units. In this study, a total of 1841 slope units covering the whole study area were surveyed and the location and boundary of discernible slope units were obtained by the geographic information system (GIS) (Figure 3a). According to the current morphological characteristics and active state of the slopes in the field investigation, the slope units were divided into loess landslides, unstable slopes, and slopes to be evaluated. Loess landslides are the main historical landslides in the study area; unstable slopes show some deformation signs, such as creep slip, collapse, toppling, etc., and are developing toward becoming potential landslides. Finally, the landslide inventory map of the study area was aggregated and shown in Figure 3b, including 344 loess landslides and 411 unstable slopes, detailed in Section 4.1.

**Figure 3.** The inventory map of the present study. (**a**) Slope unit mapping and road elements based on topography and field surveys. (**b**) The distribution of landslides and unstable slopes in the study area.

### *3.2. Factors Influencing Landslides*

Selecting appropriate environmental factors and inducing factors is the basis of risk assessment, which depends on data availability, scale, and study area, and affects future predictions [21,48]. Based on the field survey and previous work on landslides in the study area, we considered eight factors: slope, profile curvature, relief, the normalized difference vegetation index (NDVI), landslide density, building density, the thickness of loess, and the thickness of exposed bedrock. They are described below.

• Slope and profile curvature: A slope gradient is the measurement of the steepness of a surface. If the slope is too low, the gravitational potential energy is insufficient, and if the slope is too high, the material accumulation cannot provide the material basis for landslides. A profile curvature is used to describe the complexity of the terrain, which is divided into convex, straight, and concave profiles, and reflects convergent and divergent drainages in addition to variations in erosion rate [49]. In the study, the slope and profile curvature were calculated using ArcGIS and a DEM with a spatial resolution of 2 m (Figure 4a,b).

**Figure 4.** The conditioning factors used in this study: (**a**) slope; (**b**) profile curvature; (**c**) relief; (**d**) NDVI; (**e**) landslide density; (**f**) building density; (**g**) loess thickness; (**h**) thickness of exposed bedrock; and (**i**) deformation. The above parameters have been normalized.

opment of landslides [50]. The NDVI value was calculated using the expression NDVI = (NIR − R)/(NIR + R) from Landsat-8 images, where NIR is the reflectivity of the near-infrared portion of the electromagnetic spectrum and R is the reflectivity of the red portion of the electromagnetic spectrum (Figure 4d).

• Landslide density and building density: The landslide density directly reflects the development quantity of disasters in an area. In urban areas, construction activity is one of the most dominant human activities that cause slope instability. Settlement along slopes in urban areas is an important factor in slope failure. Therefore, we used building density to reflect the effect of human activities on slope stability. Landslide density and building density were calculated for the slope units by vectorizing landslides and building contours and then interpolating them into grid data (Figures 3 and 4e,f).

The thicknesses of the loess and the exposed bedrock were measured during field investigations (Figure 4g,h). The loess thickness on a hillslope, which coincides with the failure depth, is a critical parameter in performing the slope stability analysis. The overlying loess thickness plays an important role in hydrological effects, such as the ratio of the saturated depth to the losses [51]. The thickness of the exposed bedrock of the slope has a great impact on the landslide scale, landslide type, and slope deformation [52]. Due to the undeveloped tectonic activity in the study area, the effect of earthquakes and faults has not been considered in this study. In addition, because precipitation within the relatively small study area is mostly unvarying, precipitation data were excluded from the analysis processes.

A statistical description of the influencing factors is shown in Table 1. To eliminate the dimensional influence of factors, the minimax normalization method was applied [53]. The continuous factor values of each factor were normalized, so all the values fall in the (0, 1) interval, where the normalized data were calculated following the equation below:

$$X\_i^{\prime} = \frac{\mathbf{x}\_i - \mathbf{x}\_{\min}}{\mathbf{x}\_{\max} - \mathbf{x}\_{\min}} \tag{1}$$

where *Xi* is the normalized input and *xi*, *xmin*, and *xmax* are the actual, minimum, and maximum input data, respectively. The results of normalized factors are shown in Figure 4.

**Table 1.** Statistical description of the influencing factors.

### *3.3. DInSAR*

In general, the displacement of a pixel is calculated using the interference phase difference between two SAR images by using the pixel product of a reference image and slave image—this is the basic principle of InSAR [54]. DInSAR is applied to the removal of the topographic phase contribution from the interferogram deformation phase using a two-pass, three-pass, or four-pass technique; however, it is worth noting that the twopass technique, which imports an external DEM, yields a more reliable and operational outcome [54,55]. Furthermore, several limitations of InSAR technology must be considered at the beginning of use. One limitation is geometric distortion caused by topography, especially in mountainous areas with high elevations, which is affected by the look side of

radar observation modes [56]. Another limitation is the poor coherence, even incoherence of interferograms caused by diffuse vegetation, which is very obvious in the C-band Sentinel-1 images of the study area [57].

In this research, two ascending SAR images acquired from ALOS-2 on 5 November 2018 and 20 May 2019 were selected for interference calculations (Figure 5a and Table 2). Due to the relatively flat terrain of the study area, the SAR images from a single orbit can be used to detect and monitor the deformation of most of the slopes. The sensor of ALOS-2 can transmit and receive the L-band with strong penetrating ability and can capture the ground deformation under the dense vegetation. The external DEM for removing the topographic phase and geocoding is the 1-arc-second (~30 m) Shuttle Radar Topography Mission (SRTM) data from NASA. Ground deformation along the LOS (light-of-sight) of the Yan'an City area was obtained after registration and resampling, differential interference, coherence calculation, filtering and phase unwrapping, orbit refining, and reflattening in addition to geocoding. All of these were processed with the DInSAR tool of the SARScape software (Figure 5b). The normalized deformation factor image is shown in Figure 4i.

**Figure 5.** Area covered by two ascending SAR images from ALOS-2 (blue rectangle) in (**a**), and image of the LOS displacement of the study area (red polygon) using DInSAR technology in (**b**).

**Table 2.** Details of two SAR images and displacement from ALOS-2.

#### *3.4. Random Forest*

Random forest has been widely used in data classification and management and has excellent performance in landslide susceptibility mapping. Random forest is an ensemble machine learning algorithm based on a decision tree. The classifier is a recursive process from root nodes to child nodes, which is similar to the combination of a decision tree and the flowchart of a tree structure [58]. The bootstrap method is used to extract multiple samples from the original samples. Starting from the nodes of a tree, the optimal features among different internal nodes are selected, and the corresponding branches are determined based on the test output. Finally, the results are obtained from the leaf nodes of the decision tree.

Random forest has strong generalization ability and can deal with multi-dimensional and large learning sets. Compared with other statistical learning models, random forest does not easily generate overfitting. It improves prediction accuracy without significantly increasing the amount of calculation. It has a higher tolerance for outliers and noise, resulting in data loss and imbalance. In this study, a random forest module was built based on the R language. Before running the random forest module to perform the landslide susceptibility assessment, the training and validation datasets must be selected. In the study, landslide inventories including stable slopes, unstable slopes, stable landslides, and unstable landslides were selected as training sites. Using ArcGIS, these slope units were converted into 109,981 vector points, where 3000 points were selected randomly at each of the landslide sites and non-landslide sites to train and test the classifier. Finally, its performance was evaluated with the ROC curve and confusion matrix.

#### **4. Results**

#### *4.1. Characteristics of Landslides*

Based on the field survey, the geohazards of Yan'an City were counted as 334 landslides, 411 unstable slopes, and land subsidence locations. Their depths of the sliding surface are mainly shallow (less than 30 m) [47]. The landslides can be classified as loess landslides and loess-bedrock interface landslides because most landslides occur in the loess layer or on the top of the bedrock (Figure 6a,b). The geometric morphology and characteristics of the loess landslides and the unstable slopes including the types, lengths, widths, height, slope angles, and others, such as the longitudinal shape and depth of slide surfaces, were mapped using the GIS and field investigations with high-resolution DEM (~2 m) (Table 3). The length and width of landslides are mainly in the range of 50 to 200 m. The height and slope angle are also condition factors of loess landslides. A higher or steeper slope has a higher degree of stress concentration and tensile stress range so it is more prone to failure and sliding. The study area is located in the loess hilly gully region with dense gullies, and the relative height differences of 60 to 150 m leads to the height of landslides usually being less than 120 m. Since the late Cenozoic era, the Loess Plateau has been in a

state of intermittent uplift, with rivers cutting sharply and ravines crisscrossing, creating topographic conditions for loess landslides. The longitudinal shape can control the values and positions of the stress inside the slope body and plays a key role in the stability. For example, flat and convex slopes tend to be more easily destroyed under stress, suggesting an unstable evolutionary trend, whereas sunken and stepped slopes tend to be more stable with less stress concentration. Therefore, the longitudinal shapes of landslides in this study are mainly flat and convex slopes.

**Figure 6.** Examples of geohazard types mapped in the study area: (**a**) loess–bedrock interface landslide; (**b**) loess landslide; (**c**) soil–bedrock unstable slope; (**d**) soil unstable slope; and (**e**,**f**) cracks and damages in the ground and buildings due to land subsidence. Arrows indicate the direction of the slide and the location of the cracks.

Unstable slopes refer to a slope with creep slip, collapse, toppling, lateral tensile fracture, and other deformation characteristics or trends, and that is regarded as a potential geohazard. The 411 unstable slopes from the field survey were divided into three types according to their material composition: (i) soil unstable slope, (ii) rock unstable slope, and (iii) soil–bedrock unstable slope. There were about 285 soil unstable slopes in the study area, accounting for 69% of the total number of unstable slopes. Soil–bedrock unstable slopes and rock unstable slopes are fewer, numbering 125 and five, respectively, accounting for 31% of the total number of unstable slopes and mainly occurring in the Quaternary loess and the Jurassic sandstone strata. The unstable slopes have similar characteristics to landslides in their ranges of length, width, height, and area but the slope angles of unstable slopes are relatively larger. The characteristics of landslides and unstable slopes in the present study are summarized in Table 3.

In addition, under the pressure of population growth and development as well as the preservation of historical and cultural sites, in 2012 the government built a new district called Yan'an New District by cutting mountains and filling ditches. However, because of the special microstructure and complex engineering-geological conditions of the loess, land subsidence in Yan'an New District has become one of the geohazards that requires much attention. The surface deformation along the radar LOS calculated by the DInSAR technique was very similar to that of the small baseline subsets InSAR (SBAS-InSAR) approach from Sentinel-1 images (see [59,60] for more details). DInSAR and field surveys show that the land subsidence area is the ribbon (Figure 7), mainly concentrated in the filling area manifested as wall cracking or collapse, ground subsidence, and cracks (Figure 6e,f).

**Figure 7.** Surface deformation image of Yan'an New District from 2018 to 2019 calculated using DInSAR. Positive values indicate that the ground object deformation is close to the radar along the radar LOS, and negative values indicate that the ground object deformation is far away from the radar along the radar LOS.

#### *4.2. Landslide Susceptibility Mapping*

The normalized factors were used as the input, and the landslide susceptibility index was the output data. The mean decrease accuracy and mean decrease Gini coefficient of the random forest are used to order the eight variables (Figure 8). The vertical axis represents the eight variables, with the mean decrease accuracy and mean decrease Gini coefficient decreasing from top to bottom. It shows that the importance of the hazard density is the highest and that the build density, thickness of exposed bedrock, loess thickness, relief, and the NDVI are the next most-important. The ROC curve is widely used to evaluate the classification results of the random forest classifier through the area under the ROC curve (AUC) [35]. The vertical axis and horizontal axis represent the true positive rate (TPR) and false-positive rate (FPR) using the random forest classifier, respectively. TPR and FPR, also called the sensitivity and specificity, are the ratio of the landslide sample points correctly detected by the classifier and the ratio of the non-landslide sample points incorrectly classified as landslide sample points, respectively [61]. The larger the AUC, that is, the closer the vertex of the curve is to the upper left corner, the better the classifier's test capability. In this research, the AUC is 0.96, which indicates excellent classification results of the random forest classifier (Figure 9). In addition, the confusion matrix shows that the overall accuracy of the random forest classifier is 0.903 and that the predicted precision of non-landslides and landslides is 0.927 and 0.881, respectively, which is a good method to analyze the prediction accuracy (Table 4). Four levels of susceptibility, i.e., very high (>0.711), high (0.711–0.458), moderate (0.458–0.231), and low (<0.231), were categorized based on the natural breaks classification conducted using the ArcGIS software (Figure 10). The natural breaks classification was determined based on natural groupings inherent in the data. Then, the classification interval was identified to provide an optimum grouping of similar values and maximize the differences between classes [62,63]. Additionally, the distribution of the landslide susceptibility index using the natural breaks is shown in Figure 10b.

**Figure 8.** Mean decrease accuracy and mean decrease Gini of variables assigned by the random forest classifier. The vertical axis is the inducing factor variable. hazard—hazard density; build—build density; lithe—thickness of exposed bedrock; loess—loess thickness; height—relief; ndvi—NDVI; curv—profile curvature.

The results show that the distribution of landslides and unstable slopes in the study area is closely related to the susceptibility partitioning (Table 5). Over one-third of the slope units in the study area are in the high- and very high-susceptibility areas, accounting for 21% and 16% of the total, respectively, with a total area of 10.1 km2; the remaining slope units are in the moderate- and low-susceptibility areas. Moreover, 35% of the landslides and the unstable slopes are located in the very high-susceptibility area, accounting for 16% of the total number of the slope units; 33% of the landslides and the unstable slopes are located in the high-susceptibility areas, accounting for 21% of the total number of the slopes; 21% of the landslides and the unstable slopes are located in the moderate-susceptibility area, accounting for 27% of the total number of the slopes; and 11% of the landslides and the unstable slopes are located in the low-susceptibility areas, which account for 36% of the total number of slopes.

**Figure 9.** ROC curves and AUC value for evaluating the classification results of landslide susceptibility using the random forest classifier.

**Table 4.** Confusion matrix of the random forest classifier.

The results of the landslide susceptibility assessments show that two regions are highly prone to landslides, one being the urban center of Baota District, and the other Nanniwan Airport (Nnwa) and its surrounding areas. In the urban center of Baota District, including Yangjialing Village (Yjl), Nanshi Street (Ns), Baiping Village (Bp), Hutoumao Village (Htm), Zezigou Village (Zzg), Nanzhaibian Village (Nzb), Majiawan Village (Mjw), Huanghaowa Village (Hhw), Mata Village (Mt), Erzhuangke Village (Ezk), and Shanlangcha Village (Slc), where landslides occur frequently, the landslide susceptibility is high and very high because of the very high-density population and frequent human activities (Figure 10). Due to the effects of road construction, domestic water discharge, crop planting, slope toe excavation, and other activities, landslides, including rock falls, slope failures, unstable slopes, and creep, occur frequently, which poses a great threat to the lives and properties of the local residents. The other highly landslide-prone areas are Nanniwan Airport (Nnwa) and its surrounding areas, including Yangjiawan Village (Yjw), Maozegou Village (Mzg), Sanshipu Village (Ssp), and Yejiagou Village (Yjg). Nnwa is the area of mountain excavation and valley infilling on the Loess Plateau and its construction destroys the stability of the surrounding slopes, resulting in the development of landslides and unstable slopes in the surrounding areas.

**Figure 10.** Landslide susceptibility map of the study area in (**a**), and the distribution of the landslide susceptibility index using the natural breaks in (**b**).

**Table 5.** Slope unit statistics based on landslide susceptibility, hazard, and risk zone.

#### *4.3. Landslide Hazard Mapping*

Reliable landslide hazard mapping is crucial for hazard mitigation and risk management. In this study, InSAR technology was used to obtain the landslide hazard assessment, aiming for an ongoing and quantitative practice [63]. The fieldwork showed that the landslide type in the study area was relatively single, mainly loess landslides, and that the geological environment and inducing conditions, such as rainfall, are similar in the small study area, leading to the relatively simple mechanism of loess landslide activity. Therefore, DInSAR was used as a comprehensive indicator to reflect slope displacement, whether caused by rainfall or human activities, in landslide hazard assessment. The spatial probability of landslides (landslide susceptibility) and the intensity of ground surface deformation were used in the weighted overlay model parameters to calculate landslide hazard [64]. The weighted overlay technique is defined to develop a map using the overlays of several raster layers by giving weight to each raster layer according to expert opinions [65]. The weighted overlay analysis was applied to obtain the landslide hazard assessment using the following equation:

$$\, \, WX \, \_i = \sum\_{j=1}^{m} R(j \, ) \, \times \, \, X(i, j) \tag{2}$$

where *m* is the total number of factors to assess, *WXi* is the hazard index of the assessment units, *R*(*j*) is the weight value of each factor, and *X*(*i*,*j*) is the value of the assessment factors. In this study, *X*(*i*,*j*) is the landslide susceptibility index obtained from the random forest and the ground–surface deformation intensity that was defined using the normalized ground deformation data obtained from DInSAR during the monitored time; the weight values of both were set at 0.5 after analyzing the geological environment and inducing conditions of landslides in the study area, respectively. Finally, the hazard indexes of slope units were calculated by summing the product of assessment factors and corresponding weight values. Four levels of hazard, i.e., very high (>0.594), high (0.594–0.416), moderate (0.416–0.269), and low (<0.269), were categorized based on the natural breaks classification, and the LOS displacement in different hazard levels were counted, which are illustrated in Figures 11 and 12, respectively.

The number and LOS displacement values of slope units in different hazard levels are illustrated in Figure 11. The distribution histogram shows the maximum and mean displacement values, as well as the number of slope units in different hazard levels. It shows that the displacement values of slope units are distributed in a normal curve and that the higher hazard of slope units presents a higher displacement value than the lower hazard on the whole. The results show that 40% (732) of the slope units in the study area are in the high- and very high-hazard areas for landslides, accounting for 24% and 16% of the total, respectively, with a total area of 11.1 km2 (Figure 12). There was a small increase in the number and distribution of the hazard zones in the urban areas compared with the susceptibility map. About 34% (254) of the landslides and the unstable slopes are located in the very high-hazard areas, accounting for 14% of the total number of slopes; 32% (234) of the landslides and the unstable slopes are located in the high-hazard areas, accounting for 13% of the total number of slopes; 23% (170) of the landslides and the unstable slopes are located in the moderate-hazard areas, accounting for 9% of the total number of slopes; and 11% (83) of the landslides and the unstable slopes are located in the low-hazard areas, which are 5% of the total number of slopes. The results show that the spatial distribution of landslide hazard areas was consistent with the field investigations.

**Figure 11.** Distribution histogram of the number and LOS displacement values of slope units in different hazard levels: (**a**) low hazard; (**b**) moderate hazard; (**c**) high hazard; and (**d**) very high hazard. The red lines represent the average displacement values in different hazard levels.

#### *4.4. Landslide Risk Mapping*

The JTC-1 Joint Technical Committee on Landslides and Engineered Slopes noted that landslide risk is a measure of the probability and severity of the adverse effects of landslides on health or property, which must consider the hazard mapping and vulnerability of landslides [5]. Vulnerability assessment is a fundamental element in the evaluation of landslide risks [66]. Vulnerability to landslides is expressed in economic (monetary, quantitative) and heuristic (qualitative) scales. When using economic measurements, vulnerability is commonly expressed in the element values, such as monetary, intrinsic, and utilitarian values [67]. Due to a lack of information about properties and population distribution data, the Kriging interpolation of building distribution and building density was used for the vulnerability assessment in this study. The location and spatial distribution of buildings reflect the distance between buildings and slope units, which indirectly indicates the extent to which buildings and populations are threatened by landslides. Additionally, the building density can also indicate the properties and populations. Of course, this assumes that the sizes and values of the buildings are similar and that the differences in the populations attached to the different buildings are slight. The equation for landslide risk calculation is expressed as follows:

$$R = H\_L \times V\_L \tag{3}$$

where *HL* and *VL* represent the landslide hazard and vulnerability, respectively. The landslide risk index obtained from Equation (3) is divided into four levels according to the natural breaks method after normalization, namely, very high-risk (>0.406), high-risk (0.406–0.223), middle-risk (0.223–0.101), and low-risk (<0.101). The results of the risk assessment zones and statistics are shown in Figure 13 and Table 5.

**Figure 12.** Landslide hazard map of the study area in (**a**), and the distribution of landslide hazard index using the natural breaks in (**b**).

The landslide risk assessment map shows that the risk in the urban center is higher than that in the suburban areas, where the risk decreases with increasing distance from the urban center (Figure 13). Additionally, a total of 20 and 167 extra slopes are in very high and high-hazard zones besides landslides and unstable slopes (Table 5). About 6% (116) of the slope units are located in very high-risk zones with a total area of about 2 km2; 20% (377) of the slope units are located in high-risk zones with a total area of about 4.4 km2, which are mainly distributed in concentrated areas (i.e., Yjl-Sy-Slc-Hhw) (Figure 13 and Table 5). The building and population densities in these areas are high, which may lead to significant economic losses and casualties, so it is necessary to pay more attention and conduct landslide risk management to mitigate the landslide risks. Compared with the landslide susceptibility and hazard maps, Nnwa and its surrounding areas are classified into moderate- and low-risk areas because of the low population density in the areas. In addition, many engineering solutions, including slope geometry modification, underground drainage systems, gravity retaining walls, and anti-slide piles, have been applied to stabilize the slopes. Therefore, the slopes, which are originally very highly prone to landslides, are classified as low-risk zones for landslides.

**Figure 13.** Landslide risk map of the study area in (**a**), and the distribution of the landslide risk index using the natural breaks in (**b**).

#### **5. Discussion**

Landslide risk assessment has been attracting the attention of researchers and governments in order to effectively deal with landslides in the study area. For this purpose, machine learning and DInSAR technology were used to evaluate the landslide susceptibility in the main urban parts of Yan'an City. Currently, the susceptibility, hazard, and risk of landslides in the whole of Yan'an City have been determined in existing studies. In terms of the methodology, they can be divided into the heuristic model and generalized objective functions based on experts' knowledge scoring [44,68,69], a quantitative model of evidence weight [42], and a physically deterministic model [43]. In terms of map units and scopes, they can be divided into grid cells (25 m or 30 m), catchment basin units [42–44,68,69], administrative boundaries [42–44], and watershed boundaries [68,69]. However, those results are not only subjective but also can only meet the needs of a wide range of risk management options and cannot truly reflect the geomorphic characteristics of the slope in the study area, which can be useful for risk management in a large administrative area. The research reviewed that an inventory including detailed landslide information and a reasonable mapping unit as well as model type is a prerequisite for obtaining highly accurate assessment results [20]. Firstly, a landslide inventory with more detailed information can provide more input data to the model to analyze the relationship between landslides and geological environment factors to obtain a comprehensive landslide susceptibility map [70]. Some previous studies in the study area may not have shown this. Therefore, a complete landslide inventory, mapping all landslides in the study area, was determined through the interpretation of UAV images and site-by-site investigations. Secondly, field surveys and risk assessment of slope units on a large scale in small areas can provide planners with an adequate and applicable landslide risk map, especially in areas of critical concern such as urban centers. Research shows that grid cells or pixels are still the most commonly used map units in current papers on landslide assessment, and only a few papers have used slope units [20,21]. To reflect the geomorphic characteristics of slopes, slope units were used as map units in this study. Thirdly, the selection of a model is also an important factor affecting the accuracy of landslide susceptibility assessment. There are more and more models and methods developed for landslide susceptibility assessment, among which machine learning with good performance can be used to solve the nonlinear relationship between landslides and geological environment factors [16,34]. For this reason, random forest was selected to predict landslide susceptibility in the study area, with good proven performance [34,71]. Therefore, the accuracy of the landslide susceptibility assessment in this study was improved by using a machine learning model and slope units.

In addition, DInSAR technology was introduced in the process of hazard assessment to calculate slope deformation, and the hazard was calculated by giving the same weight coefficient of susceptibility and slope deformation. InSAR technology was used to perform a time-effective analysis, and the results can present the active state of slopes directly to predict the failure time and assess the hazard class of landslides [72,73]. At present, combining the ground deformation products from the InSAR technology with a landslide risk assessment map has become a concern in the relevant research [63,74]. However, many existing studies focus on the early identification and long-term monitoring of temporal and spatial evolution using InSAR technology [75,76] but there is still insufficient attention paid to risk assessment. The application of InSAR data in landslide risk assessment can improve the reliability of landslide predictions and make a reliable landslide risk map [45].

Finally, landslide risk was calculated by multiplying the hazard with the vulnerability composed of the spatial distribution and density of the buildings. The susceptibility, hazard, and risk assessment in this paper have a similar trend to the previous paper on Yan'an City in the area of different levels and spatial distributions [42–44], such as the area and percentage of low-susceptibility or hazard zones being greater than that of the highersusceptibility or hazard zones, and the high-risk zones in the spatial distribution patterns are similar, which can also imply the accuracy of this work to a certain extent. Moreover, the risk assessment in this study can provide more specific guidance for risk management and prevention in practice.

There are some limitations that need to be considered in future research. Firstly, due to the lack of detailed population and property data, only the spatial location information of buildings was used for the vulnerability. The precondition for this question is to assume that the values of the buildings in the study area are the same, which would lead to certain information loss for vulnerability. Secondly, influenced by the observation mode of the radar satellites, the deformations obtained by the SAR images are ultimately along the LOS. However, the deformation rate along the slope (Vslope) can more intuitively reflect the real motion of the slopes, which can be transformed through the spatial geometric relationship between the radar LOS and slope. Due to the limitations of the image numbers in the SAR dataset from the study area, we have to use the ascending SAR images for InSAR processing, with which it is difficult to form an effective complement for the descending data. Therefore, we will also try to transform the LOS displacement into the slope direction displacement in future research. PSI or SBAS algorithms can be selected to obtain long-term ground deformation products if the SAR datasets have sufficient and long-term images, which can reflect the long-term movement status and trend of slopes [63,74]. It is expected that more accurate risk assessment maps for the study or elsewhere could be produced by improving the above limitations.

#### **6. Conclusions**

Quantitative risk assessment is very effective for landslide risk management and urban areas need more detailed investigations and assessments. In this study, the quantitative landslide risk assessment was based on susceptibility and hazard assessments. The random forest classifier and eight environmental factors influencing landslides, including slope, profile curvature, relief, NDVI, landslides density, building density, the thickness of loess and the thickness of exposed bedrock were used to examine landslide susceptibility in Baota District, Yan'an City. Combined with DInSAR technology, landslide hazard mapping was developed to reflect the hazards quantitatively. Surface deformation, which can be caused by many factors (e.g., precipitation, slope groundwater, and engineering), can be detected by DInSAR technology with centimeter precision. Finally, the landslide risk map was obtained by being combined with the landslide susceptibility and hazard assessment and divided into very high-risk, high-risk, middle-risk, and low-risk areas according to the natural breaks method.

In this study, a total of 1841 slope units were mapped in the study area, including 334 landslides and 411 unstable slopes determined by field investigations, in which the main material of landslides and unstable slopes is loess and only a few of them contain bedrock. The length and width of landslides and unstable slopes are mainly between 50 m and 150 m, the slope angles are mainly between 20◦ and 50◦, and the heights are predominantly between 30 and 90 m, where the slope angles and heights of most of the unstable slopes are larger than those of the landslides. The areas are usually less than <sup>20</sup> × <sup>10</sup><sup>3</sup> <sup>m</sup>2. Reliable risk assessment was achieved using 1841 slope units, which were divided based on the terrain, optical images, and DEM. Remote sensing InSAR technology was applied to determine the quantitative landslide hazard zones. The classification results of the random forest classifier were evaluated with the receiver operating characteristics (ROC) curve and confusion matrix. The confusion matrix shows that the overall accuracy of the random forest classifier is 0.903 and that the AUC value is 0.96, with good prediction accuracy and classified ability of landslide susceptibility. The results of the landslide risk assessment indicate the risk level and the corresponding quantity of the slope units and total areas. Approximately 6% of the slope units located in the very high-risk zones and 20% of the slope units located in the high-risk zones must receive more attention to monitor the dynamics.

The present research has significant implications for landslide risk mitigation in Baota District, Yan'an City. Our scientific landslide risk map is expected to promote landslide prevention based on a zoning strategy and provide a valuable decision to support the local and regional government for disaster prevention, mitigation, and management, which eventually can effectively reduce the impacts of geohazards.

**Author Contributions:** Conceptualization, Y.Z.; methodology, W.L., Y.Z. and Y.L. (Yiwen Liang); software, W.L. and Y.L. (Yiwen Liang); validation, P.S., W.L., Y.Z. and Y.L. (Yiwen Liang); formal analysis, W.L. and Y.L. (Yiwen Liang); investigation, Y.L. (Yuanxi Li), X.S. and A.W.; resources, Y.Z.; data curation, Y.Z.; writing—original draft preparation, W.L. and Y.L. (Yiwen Liang); writing—review and editing, Y.Z. and X.M.; visualization, W.L. and Y.L. (Yiwen Liang); supervision, Y.Z. and X.M.; project administration, Y.Z.; funding acquisition, Y.Z. and X.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the National Natural Science Foundation of China (Grant No. 42007232, 42130709), the Natural Science Foundation for Young Scientists of Gansu Province (Grant No. 20JR5RA223), the Science and Technology Project of Gansu Province (Grant No. 18JR2JA006), and the Fundamental Research Funds for the Central Universities (Grant No. lzujbky-2021-ey05), Geological Survey Project of China (Grant No. DD20189270).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The ALOS-2 images were provided by the Japan Aerospace Exploration Agency (JAXA). We are grateful to the editor and four anonymous reviewers for their constructive suggestions and comments to improve the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.