An Application of Sentinel-1, Sentinel-2, and GNSS Data for Landslide Susceptibility Mapping

Abstract

In this study, we used Sentinel-1 and Sentinel-2 data to delineate post-earthquake landslides within an object-based image analysis (OBIA). We used our resulting landslide inventory map for training the data-driven model of the frequency ratio (FR) for landslide susceptibility modelling and mapping considering eleven conditioning factors of soil type, slope angle, distance to roads, distance to rivers, rainfall, normalised difference vegetation index (NDVI), aspect, altitude, distance to faults, land cover, and lithology. A fuzzy analytic hierarchy process (FAHP) also was used for the susceptibility mapping using expert knowledge. Then, we integrated the data-driven model of the FR with the knowledge-based model of the FAHP to reduce the associated uncertainty in each approach. We validated our resulting landslide inventory map based on 30% of the global positioning system (GPS) points of an extensive field survey in the study area. The remaining 70% of the GPS points were used to validate the performance of the applied models and the resulting landslide susceptibility maps using the receiver operating characteristic (ROC) curves. Our resulting landslide inventory map got a precision of 94% and the AUCs (area under the curve) of the susceptibility maps showed 83%, 89%, and 96% for the F-AHP, FR, and the integrated model, respectively. The introduced methodology in this study can be used in the application of remote sensing data for landslide inventory and susceptibility mapping in other areas where earthquakes are considered as the main landslide-triggered factor.

1. Introduction

Natural mass-movement hazards, such as landslides, cause extensive damage to the environment and social infrastructure, severely affecting urban development and land use. Landslides, threatening the lives and property of local people, are defined as one of the most harmful geological hazards that frequently occur in mountainous areas. Landslides are known as the downslope mass movements of earth, debris, and rock, due to gravity impact [1]. They always happen in variety of sizes and shapes and are triggered by different criteria, such as rainfall, earthquakes, slope erosion, and volcanic eruptions [2]. In this study area, a severe earthquake was the main triggering factor for both large and small landslides.

Susceptibility mapping is considered as the primary approach for any landslide analysis. Thus, researchers have applied landslide susceptibility mapping (LSM) as a practical analysis approach to obtain better insights into this phenomenon [3]. LSM is also the best way to identify areas with potential slope instability. In this regard, remote sensing and geographic information system (GIS) data were extensively used for LSM studies [4].

Differential synthetic aperture radar interferometry (DInSAR) data are suitable for monitoring mass movements, especially for substantial surface deformations caused by earthquakes. Several previous studies used DInSAR data to analyse surface deformations [5]. For instance, Casagli et al. (2016) [5,6] used interferometric synthetic aperture radar (InSAR) and object-based image analysis (OBIA) and discussed instance cases in Italy, Austria, Slovakia and Taiwan. Dabiri et al. (2020) [7] applied Sentinel-1 and Sentinel-2 data to assess the landslide-induced geomorphological changes in Hítardalur Valley, Iceland. In one of the most recent works studying the Kermanshah earthquake, Goorabi (2019) [8] evaluated Sentinel-1A, B, C-band coherence maps to detect landslides caused by the earthquake. There is also a wide range of previous studies that used OBIA for landslide detection, such as Hölbling et al. (2008) [9] who applied an object-based approach to detect mass movements. The results of their study were validated using persistent scatterer interferometry measures based on ERS1/2, RADARSAT-1 data, ENVISAT ASAR, and a landslide inventory map derived from two different sources. Shirvani et al. (2019) [10] applied OBIA along with a machine learning model of random forest for detection of old and new landslides using Sentinel-1 and -2 data. Tavakkoli Piralilou et al. (2019) [11] also integrated OBIA with several machine learning models for landslide detection within a multi-scale image analysis. They studied the affected area of the Gorkha earthquake (magnitude of 7.8 M, April 2015) in the higher Himalayas. Subsequently, they integrated their resulting landslide detection maps of different scales using the Dempster–Shafer theory to reduce detection errors. Recent studies have also applied deep learning approaches and convolutional neural networks, in particular, to landslide detection [12,13]. However, deep learning approaches need massive training samples of the labelled remote sensing images to show their priority to the other approaches. Moreover, the performance of the integration of DInSAR data along with optical data in deep learning approaches is still unclear. Based on our literature review, we can conclude that both DInSAR and OBIA can be considered as sophisticated approaches of spatial geodesy to detect land deformation hazards and landslides in particular. The wide spatial coverage of DInSAR data and its sensitivity to surface deformations in both vertical and horizontal directions, make it a practical method for detecting huge mass movements [14].

For the susceptibility modelling of any land deformation hazard, and LSM in particular, well-prepared input data sets, such as DInSAR data, are necessary for spatial modelling in GIS. The integration of different GIS spatial analysis models for identifying areas of possible future landslides through GIS-based LSM is an acknowledged approach to hazard mitigation [15]. Such approaches lead to increased accuracy of the final results, which is essential for any environmental planning manager to mitigate adverse consequences [16]. LSM is a popular approach to identify landslide-prone regions. Integrated models combining a selection of related landslide factors have produced acceptable LSM results [17]. The simultaneous use of DInSAR and GIS data, along with an appropriate expert knowledge-based or data-driven model, or an integration of the two, has been shown to produce very successful results for monitoring land deformations [18]. Several recent LSM studies have been conducted using knowledge-based models. Technically speaking, these studies implemented a range of models, such as the analytic hierarchy process (AHP). The AHP has produced relatively good generic LSMs, especially on small and medium scales, where we do not have access to an appropriate inventory data set [19]. The analytic network process (ANP) model is a general form of the AHP that was developed in the context of interdependent networks to be used for more complex spatial problems and to minimise the inaccuracies of traditional models [20]. Although the ANP is useful for solving complex spatial problems, it requires expert knowledge to organise the pairwise comparison matrices, similar to the AHP method. Thus, the LSMs resulting from both the ANP and the AHP models are strongly dependent on expert knowledge and preferences [21]. Therefore, several studies have attempted to improve the accuracy of spatial decision support systems by integrating them with other models, including fuzzy logic [20,22], Boolean intersection methods, such as order weighted average (OWA) and weighted linear combination (WLC) [23,24], uncertainty and sensitivity analyses [25,26], interval calculus [27,28], and also some data-driven models, such as the frequency ratio (FR) [29].

A literature review shows that some studies focused on landslide detection using different data sources and approaches. Moreover, some other studies applied different models for LSM based on expert knowledge and/or landslide inventory datasets. Current LSM studies represent that several different mathematical concepts have been effectively combined with knowledge-based models to optimise the model and enhance the overall accuracy of natural hazard susceptibility mapping. Combining various models is a common approach to obtain more reliable LSM results. However, landslide detection approaches and LSM models were not applied in one comprehensive study within a synergistic combination. To this end, we aim to use data from different sources of Sentinel-1A and -2B, a field survey of the global positioning system (GPS) data, and Google Earth™ images, along with DInSAR and OBIA for detecting the occurrence of landslides in our case study area. In addition, we want to use the fuzzy AHP and FR models separately and in integration to create more accurate LSMs. Thus, in this study, first an earthquake-induced landslide inventory at a regional scale was generated using Sentinel 1 and 2 images along with OBIA. The integration of optical and radar images for landslide detection application is very rarely done for regional-scale landslide mapping. Using the generated inventory, we trained a novel landslide susceptibility model by combining data-driven FR, and knowledge-based fuzzy AHP approaches. These two approaches are already well established as independent methods. The use of ensemble to improve results are also widely used across disciplines. We showed that our integrated model outperformed the base models, and the results look convincing. All of our outputs were validated with points generated from a GPS survey.

2. Overview of the Study Area

The Kermanshah earthquake was an oblique overthrust earthquake with a magnitude of Mw 7.3, which happened on 12 November 2017. The epicentre of the Kermanshah earthquake was located at 34.905° N latitude, and 45.956° E longitude [30] in the northern border region of Iran and Iraq, about 350~km north of Baghdad, at a depth of 33.9 km. The earthquake was responsible for over 600 fatalities, and 8000 people were injured in Iran alone [8]. The earthquake was widely felt in areas as far away as the most western part of Iran and in other countries such as Turkey and Kuwait [31]. Geologically, the earthquake was located within the High Zagros Fault (HZF) and the Mountain Front Fault (MFF) northwest (NW) of the Zagros Mountains [32], close to the boundary between the Arabian and Eurasian plates. The thrust faulting in the foothills of the north-western Zagros Mountains [33] was the main cause of the Kermanshah earthquake. The Arabian Plate moves about 2.5 cm/yr⁻¹ in a north/north-easterly direction towards the Eurasian Plate [31]. The tectonic setting of this area was fully described by McQuarrie (2004) [34].

Mass movements, especially landslides, caused immediately by the quake and its significant aftershocks, affected all the mountainous areas of this province [8]. This area is bound to the west by the border to Iraq (Figure 1), however, as we do not have access to the required data from the area located in Iraq, we cut the image frame along the Iranian border. The study area contains six cities and 470 villages and has a population of around 269,588 (Iranian Census and Centre 2016). Elevations in this area range from 237 m in plains to 2566 m in mountainous regions. The slopes in the study area range from less than 2% gradient in the western parts to more than 70% gradient in the mountainous regions. The study area is highly susceptible to landslide hazards [8] which result in considerable economic losses when they affect infrastructure. Figure 1 depicts several samples of the landslides observed during the field observation.

3. Workflow

3.1. Overall Methodology

The workflow of the present study is as follows:

▪

Landslide inventory generation: we applied Sentinel-1A and Sentinel-2B images to detect the earthquake-induced landslides using DInSAR interference and OBIA.

▪

Landslide susceptibility mapping: a) we collected and prepared all conditioning factors for further susceptibility analyses; b) we applied the FR and fuzzy AHP (FAHP) models for LSM using the resulting landslide inventory and expert knowledge, respectively; c) we integrated the resulting weights of the FAHP model for each conditioning factor with the FR model for each class of each conditioning factor to produce a new LSM.

▪

Accuracy assessment: we collected GPS points in an extensive field survey and randomly divided the points into 30% and 70% for validating the resulting landslide inventory dataset and the performance of each model, respectively.

Figure 2 schematically depicts this workflow, and the following sections contain detailed explanations and further descriptions. Finally, the discussion section outlines the experimental results and discussions regarding the impact of integrating of FR and FAHP models.

3.2. Landslide Inventory Dataset Generation

Landslide inventory data deliver key geoinformation about the location and distribution of landslides and affected areas in landslide-prone regions. This inventory data is critical for validating GIS-based landslide susceptibility maps [25]. Reliable landslide inventory maps provide geoinformation about slope failures and displacements affecting mass movements. Landslide inventory data sets usually include different type of ground movements showing how the landslide material (rock or soil or both) is displaced such as slides (including rockslides and landslides), failures, falls spreads, topples, and flows. These different landslide types can be characterised based on their shape. For the type of fall, for instance, the accumulation part of the landslide is usually cone-shaped, and in the scarp part, the bedrock is evident, while back scarps are the common characteristic of slides. The landslide type of flows, and especially debris flows, mostly occur across streams [35]. However, landslides usually form within more than one type movement, such as rockslide-debris flow or fall-slide. In our case study area, eight massive surface displacements were recorded based on our field survey. The biggest one resulted in a 90 cm vertical and 50 cm westward movement with a slip of around 3 m, reported by [36,37]. Although there were different ground movement types, mostly combinations of landslide types were detected. The observed combination forms of the Kermanshah earthquake were mostly flow-slide along with some small fall-slide in steep areas, based on the landslide classification of the International Consortium on Landslides (ICL) guide (http://icl.iplhq.org). As differentiating the reactivation of different post-earthquake types of landslides was a problematic task, our landslide inventory data set was prepared based on mostly huge and obvious ground movements. Nowadays, it is possible to monitor mass-movement hazards, such as landslides, using remote sensing data. Section 3.2.1 and Section 3.2.2 summarise the process of generating the landslide inventory map.

3.2.1. The Differential Synthetic Aperture Radar Interferometry (DInSAR)

Synthetic aperture radar (SAR) data was applied to detect landslides that occurred after the earthquake. These images are sensitive to changes in vegetation cover, which is why a pair of SAR images acquired at times of varying vegetation cover can be incoherent, even when acquired within one month [18]. However, this does not pose a severe problem for this study as the vegetation cover is scarce in almost the whole study area. The landslide inventory database of the areas affected by the earthquake was generated using DInSAR interference. In this study, the interferometric wide swath mode of Sentinel-1A images was used for further processing. The European environmental monitoring program Copernicus launched the satellite in 2014. For our study, we selected the frames 108 and 109 of Sentinel-1A that cover a broad region affected by the mass movements. In this study, we used three sets of repeat-pass SAR data with a minimal perpendicular baseline, acquired by Sentinel-1A, as the main input for the registration and unwrapping operations to produce land surface displacement maps. Two of the images were acquired before the earthquake (30 October and 11 November 2017), and the third was acquired less than one week after the event (17 November 2017). All pixels of the slave image were moved to align with all pixels of the master image [38]. In this step, we co-registered the precise orbit state vectors and digital elevation model (DEM), and converted them to the same coordinate system as that of the image [39]. Small differences within overlap bursts help extract the horizontal motion of the land surface parallel to the satellite path [40]. The de-bursting process was used to remove the demarcation area between two bursts. We then used the resulting de-burst stacks from the previous process to subtract the topographic phase.

The signal-to-noise ratio of the interferogram image can be significantly improved by preserving the signal features and reducing the number of residual points by convolving the interferogram image with several filter coefficients (e.g., Boxcar, GAMMA, Gaussian and Goldstein). In this study, we convolved the results from the previous stage with the Goldstein filter [41] to improve the signal-to-noise ratio. This filter is a non-linear frequency weighting method that has been applied in several studies [37]. The Goldstein filter is a low-pass filtering approach that can quickly smooth the intensity of Fourier-transformed patches of the overlapping interferogram image [42]. Thus, this filter is practical for reducing phase noise before the unwrapping process. The main error source affecting the interferometric phase in such a large deformation in our case study is decorrelation error. The interferometric phase in the landslide-affected areas significantly reduced due to decorrelation error. The application of the interferometric measurements derived from Sentinel-1A acquisitions and relative problems are fully described by Vajedian et al. (2018) [32] for the same earthquake, northern of our case study area, in the region between Iran–Iraq borders.

Besides, coherence information resulting from SAR coherence Sentinel-1A acquisitions is used as one of the input data sources for landslide inventory mapping in this study. Low coherence between two images in the landslide-affected areas along with optical and topographical data has a positive effect on the annotation of these features. The loss of coherence is considered as auxiliary information mostly for the landslide-affected areas that were not detected using other input data, e.g., optical data. The problem of existing low coherence with vegetated areas is also solved using appropriate slope thresholds and normalised difference vegetation index (NDVI) data in the OBIA procedure (see Section 3.2.2).

3.2.2. Object-Based Image Analysis (OBIA)

Earth observation is the most critical factor in landslide detection and monitoring to assist hazard analyses [9]. In this study, OBIA was used to support the DInSAR approach in landslide detection and inventory preparation (see Figure 3). In this study, multiresolution segmentation (MRS) was applied for the segmentation of our input data. We developed a semi-automated approach based on the spectral and spatial characteristics of the landslides. In this context, we applied segmentation to develop initial segments for object-based landslide mapping.

We developed the segmentation with a scale of 15, a shape index of 0.6, and a compactness of 0.4, using the estimation of scale parameters (ESP2) tool and cross validation. Accordingly, efficient object-based features were obtained based on previous research (e.g., [10,43]). These earlier works determined that the NDVI, brightness, shape index, mean flow direction, and curvature—together with texture analysis based on grey level co-occurrence matrix (GLCM) values (e.g., GLCM_correlation, GLCM_mean, GLCM_{standard deviation}, GLCM_contrast, and GLCM_entropy)—are the most efficient object-based features for detecting and delineating landslides. In this study, the only pre-processing of Sentinel-2B was atmospheric correction. Then, we applied the object-based features to detect and delineate landslides using bands 3, 4, and 2 of Sentinel-2B and bands 4 and 8 for the NDVI, SAR coherence from Sentinel-1A, and a DEM, as well as its derivations, including, aspect, slope, flow direction, and curvature. Figure 4 shows some of the object-based features that we used for landslide detection and delineation. For this goal, these object-based features and the obtained threshold values were applied for each feature by considering five known landslide points as sample data (see

Table 1. The resulting threshold values for five samples of known landslides.

	Brightness	NDVI	Mean Curvature	Mean Flow Direction	GLCM Contracts	GLCM Correlation	GLCM Entropy	GLCM Mean
Sample 1	222	0.014	−0.1	0.1	372	5.10	5.64	0.8
Sample 2	230	0.02	−0.18	0.2	278	8.3	6.7	1.3
Sample 3	214	0.018	−0.14	0.4	379	9.7	5.8	2.8
Sample 4	189	0.021	−0.19	0.3	289	8.6	6.1	1.03
Sample 5	205	0.016	−0.12	0.8	357	7.89	7.68	2.04

). Therefore, we used spectral and topographical properties, as well as DInSAR productions, for landslide detection in the object-based environment. Thus, the first step was to classify the low coherence areas as landslide areas at a primary segmentation level (see Figure 4). Accordingly, we employed the obtained threshold values for each feature (Table 1) as the basis for the fuzzy nearest neighbour classification. Subsequently, the outcomes of both the DInSAR and OBIA approaches were manually enhanced using Google Earth™ images to generate a reliable landslide inventory dataset for further analyses. The resolution of the generated inventory landslide created in this study was based on the resolution of Sentinel-2 NIR and Red bands of 10 m.

Since the literature recognises accuracy assessment as a critical step for evaluating any classified map based on remote sensing satellite imagery [44], we aimed to validate the reliability of our results using known landslide points collected by a GPS in a field operation. Based on the known landslide points, the accuracy assessment was derived using the kappa coefficient and error matrix. Our computed kappa coefficient was 0.97, and the overall obtained accuracy was 94%, which indicated a very valid outcome for the detected landslides.

3.3. Landslide Susceptibility Mapping

3.3.1. Conditioning Factors

We selected eleven landslide conditioning factors in our study area, including topographical factors, such as elevation, slope angle, and aspect. Thereby, we used a digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) (https://www2.jpl.nasa.gov/srtm/) with 30 m spatial resolution as the basis of the aspect and slope information layers. In previous research, there is a consensus that slope is the most crucial factor for landslide susceptibility assessment in any study of LSM [15]. Degradation is much greater in mountainous areas characterised by steep topography, particularly at the top of hills and peaks. This phenomenon is very obvious in the inventory dataset depicted in Figure 5. Based on previous studies in a similar area, we categorised slope layers into five classes. Furthermore, previous studies suggest that elevation is another essential factor for LSM. It is one of the effective conditioning factors in earthquake-induced landslides. Thereby, landslide susceptibility strongly correlates with the valleys and peaks of mountainous regions. Geomorphological and geological processes are also affected by elevation [45], which we used to define the local relief of the case study area. Relief classes represent the variation in elevation values of our case study area, ranging from the lowest elevation in riverbeds to the highest elevation in mountainous areas [46]. To investigate how many landslide occurrences in different relief classes, we prepared five elevation classes. The slope aspect is considered in the literature as one of the crucial conditioning factors for any landslide event. Some meteorological factors, including the average amount of sunshine and the direction of rainfall, also determine landslide susceptibility [47]. The conditioning factor of aspect indirectly impacts the slope steadiness through hydrologic processes, such as precipitation, and impacts the weathering procedure [48]. In this study, we categorised the aspect layer into nine classes, including flat (0°), north (0°–22.5°; 337.5°–360°), northeast (22.5°–67.5°), east (67.5°–112.5°), southeast (112.5°–157.5°), south (157.5°–202.5°), southwest (202.5°–247.5°), west (247.5°–292.5°), and northwest (292.5°–337.5°) (see Figure 5).

Geological factors, including lithology, distance to faults, and soil type, are represented in Figure 6. We derived the applied geological conditioning factors from the available 1:100,000 geological maps. Different lithological units and soil types play a different role in controlling landslides. Badlands constitute more than 60% of the soil in the study area. The more lithological units and formations are composed of harder minerals, weaker shock waves at the surface, and less severe damages caused by an earthquake in a given area. Additionally, the gaps and cracks present in the geological formations of an area are where shock waves are transferred from the epicentre to the surface [49]. Distance to faults is important, especially for our study where regions closer to faults are more prone to the effects of earthquakes [50]. Fault distance plays a significant role in an area’s resilience to the effects of the destruction by earthquake wave transmission power since proximity to faults causes high seismic risk and destruction. A greater distance from faults can diminish the risk. Although some previous studies, such as [51], found the factor of distance to faults to be ineffective in susceptibility modelling, we still considered it in our LSM procedures. In the present study, the lithology elements were grouped into four main groups, as described in

Table 2. Brief description of the lithology of the study area.

Groups	Geo-Units	Descriptions	Age
Group 1	$K{E}_{gu}^{pd}$	Grey and brown, medium - bedded to massive fossiliferous limestone.	Late Cretaceous
Group 2	$K^{bgp}$	Undivided Bangestan Group, mainly limestone and shale, Albian to Companian.	Cretaceous
Group 3	$Q_1^{ft}$	High level piedmont fan and vally terrace deposits	Quaternary
Group 4	$Q_2^{ft}$	Low level piedment fan and vally terrace deposits	Quaternary

. Most GPS points of the landslide inventory dataset fall within groups A and B. Group A contains grey and brown medium-bedded to massive fossiliferous limestone. Group B includes the undivided Bangestan group, mainly limestone and shale, Albian to Companian.

Hydrological factors include distance to streams and precipitation. Streams and rivers are also considered as landslide conditioning factors by erosion, as they saturate the materials in the lower part of valleys and erode the slopes. Streams generate post-glacial alluvial sand, which is one of the leading causes of surface vibration and lubrication [52]. The rainfall layer is known as the main triggering factor for slope failures and represents an important factor in LSM in some mountainous areas, such as Austria [48]. The amount and spatial variance of rainfall can vary by seasons and topographical features. In this study, we categorised the factor of annual precipitation information into five classes, namely, 0–424 mm, 424–451 mm, 451–476 mm, 476–502 mm, and >502 mm.

Anthropological factors include land use classes and distance to roads. Road construction is the leading cause of diminution in the load on the topography [53,54]. Landsat satellite data (https://www.usgs.gov/) were used as the basis of our land use layer, which has a spatial resolution of 30 m. For generating the land use map, the supervised classification and maximum-likelihood algorithm were applied [55] and five land use classes were prepared, such as Irrigated agriculture Orchard/Forests, and Grassland. Finally, we considered the normalised difference vegetation index (NDVI) as our vegetation factor. We calculated the NDVI value using the formula NDVI = (NIR − Red)/(NIR + Red), where NIR and Red are the near-infrared and red bands of the electromagnetic spectrum of Landsat 8 (launched on February 11 2013). Since our study area has sparse vegetation cover, we classified this factor into only four classes (see Figure 7). Considering the influence of individual conditioning layers, we prepared a map with a pixel size of 30 m for further analysis. The number of classes in each factor and the intervals of distances are based on the studies mentioned in the literature review in the introduction section [1,4,17,22,24,25].

3.3.2. Frequency Ratio (FR)

The FR model is a useful bivariate statistical approach that can be effectively applied to geospatial problems. The results of the FR model represent the likelihood of the presence or absence of landslides, and their distribution, for any class of landslide conditioning factors [56]. The considered classes for each factor can be easily weighted and ranked based on the ratio of landslides in each class, divided by the total number of landslides in the whole study area [52]. The FR model is an effective geospatial assessment tool to find spatial correlations among an inventory map and each class of considered factors [24]. Indeed, the FR model is defined as the ratio of the area where landslides occur (generated landslide inventory dataset) to the total study area. In this study, we computed the FR weights by calculating the ratios of each class of considered factors relative to the study area using our landslide inventory dataset, which we generated using DInSAR data. First, we rasterised the inventory dataset, including landslide polygons, and then overlaid the output with the conditioning factors layer. We then produced the final susceptibility map using a simple linear combination formula (see Equation (1)) of the resulting FR weights for all considered LSM factors.

LSM_FR = FR₁ + FR₂ + … + FR_n

(1)

where n represents the number of factors, and FR is the rating of each factor. The LSM generated according to the FR weights is presented in the results section.

3.3.3. Fuzzy AHP

Fuzzy logic mimics human reasoning in the way it handles uncertainty and uses approximate information to make decisions [57,58]. Compared to Boolean logic, where each element is considered as either true or false, the elements of fuzzy logic are considered as both the true and false but to varying degrees. This makes it possible for any element of a set to belong to its complement as well, but only to a certain degree [43,59]. Zadeh (1965) [60] introduced fuzzy sets to handle vague concepts in an obvious way. Since then, fuzzy sets have been effectively used to solve complex spatial problems that are vaguely modelled and have a degree of uncertainty. As almost all geographic phenomena are represented by vague definitions, fuzzy sets are extensively combined with different knowledge-based (e.g., fuzzy AHP) and machine learning models (e.g., adaptive neuro-fuzzy inference system) to model and map spatial problems [61]. The combination of fuzzy logic and AHP has been widely used in landslide susceptibility modelling and mapping [62,63].

In this study, in order to calculate the criteria weights using the fuzzy AHP, questionnaire forms were prepared to gather expert knowledge for further processes. Therefore, five geographers and two geologists from the University of Tabriz (https://tabrizu.ac.ir/), four local experts (during a field survey), a geologist and three landslide experts from the University of Salzburg (https://www.uni-salzburg.at), two GIS experts from the University of Vienna (https://www.univie.ac.at/) and (https://en.um.ac.ir/), a member of the Faculty of Environment and Labour Safety from the Ton Duc Thang University (https://www.tdtu.edu.vn/), and an earthquake engineer and synthetic aperture radar expert from the Tokyo Institute of Technology (https://www.titech.ac.jp/) were asked to fill in the questionnaire forms, and a total of 19 forms were obtained for weighting the condition factors based on expert knowledge. Accordingly, we used the fuzzy set theory for the decision-making processes of the AHP, based on the judgment of our target decision-makers. Thereby, we used the values resulting from the decision-makers’ factor comparisons to rank the factors according to their importance regarding the considered susceptibility. These comparisons are weighted qualitatively using linguistic variables. Therefore, fuzzy sets are needed to quantify the qualitative comparisons applying appropriate membership functions [64]. We used a fuzzy triangular set (see Figure 8) to convert the qualitative comparisons to the quantitative values for weighting our factors and two different triangular fuzzy numbers (see Figure 9) for the possible degree. The relationship between the qualitative comparisons using linguistic variables and the quantitative values is presented in

Table 3. The applied triangular fuzzy number of linguistic variables.

Linguistic Variables	Triangular Fuzzy Numbers	Reciprocal Triangular Fuzzy Numbers
Extremely strong	(9,9,9)	(1/9, 1/9, 1/9)
Very strong	(6,7,8)	(1/8, 1/7, 1/6)
Strong	(4,5,6)	(1/6, 1/5, 1/4)
Moderately strong	(2,3,4)	(1/4, 1/3, 1/2)
Equally strong	(1,1,1)	(1,1,1)
Intermediates	(7,8,9), (5,6,7), (3,4,5), (1,2,3)	(1/9, 1/8, 1/7), (1/7, 1/6, 1/5), (1/5, 1/4, 1/3), (1/3, 1/2, 1)

.

Then, the pairwise comparison matrices of the AHP were integrated with the fuzzy sets for further analyses and calculations and to reduce the subjectiveness in the criteria weightings. Incidentally, fuzzy numbers were defined to make pairwise comparisons and to structure the corresponding matrices.

3.3.4. Integrated Model

In the integrated model, the resulting FR weights were integrated with the weights derived from the FAHP to reduce uncertainties in our multi-criteria decision-making model. Therefore, we generated separate LSMs based on only the FR and FAHP models and another LSM based on the integrated model. We used the FAHP model to weight the applied conditioning factors and the FR model to weight the classes of each conditioning factor.

4. Results

4.1. Resulting Landslide Inventory Dataset

The resulting interferogram and coherence for before and after the earthquake are presented in Figure 10 and Figure 11. We integrated the results of the DInSAR method to the OBIA process as auxiliary information for landslide detection. Figure 12 presents an enlarged landslide area based on the coherence and interferogram results before and after the earthquake.

4.2. Resulting Landslide Susceptibility Maps

In the present study, we considered three models to generate LSMs and classify regions with a high susceptibility to landslide occurrence. The results of all criteria are shown in

Table 4. The final weight of each criterion from the fuzzy analytic hierarchy process (FAHP).

Criteria	Rank of Weight	Weight
NDVI	6	0.095
Soil type	8	0.069
Altitude	1	0.218
Land cover	7	0.073
Distance to streams	10	0.006
Distance to roads	11	0.003
Slope	2	0.170
Lithology	4	0.119
Distance to fault	9	0.012
Rainfall	5	0.111
Aspect	3	0.119

:

The resulting weight of each class (see

Table 5. The normalised weights of each conditioning factor based on each model.

Conditioning Factors	Classes	Percentage of Domain (%)	Percentage of Landslide (%)	FR	Normalised FR	FAHP	Normalised FAHP
Slope angle (°)	0–5	32.04	0	0	0	0.1	0.24
	5–10	21.31	0	0	0	0.07	0.17
	10–15	14.39	25	1.73	0.51	0.25	0.6
	15–20	11.26	37.5	3.33	1	0.41	1
	20<	20.97	37.5	1.78	0.53	0.15	0.36
Altitude (m)	625–700	34.87	0	0	0	0.9	1
	700–1200	19.77	37.5	1.89	0.85	0.29	0.32
	1200–1500	17.01	37.5	2.2	1	0.3	0.33
	1500–2000	18.71	12.5	0.66	0.3	0.12	0.13
	2000–2566	9.6	12.5	1.3	0.59	0.19	0.21
Distance to faults (m)	<1000	7.96	25	3.14	1	0.52	1
	1000–3000	15.06	25	1.66	0.52	0.32	0.61
	3000–9000	40.93	37.5	0.91	0.28	0.09	0.17
	9000	36.02	12.5	0.34	0.1	0.06	0.11
Distance to roads (m)	<700	52.6	62.5	1.88	1	0.4	1
	1400–700	28.27	12.5	0.44	0.23	0.35	0.87
	2800–1400	15.28	25	1.63	0.86	0.18	0.45
	2800<	3.86	0	0	0	0.07	0.17
Distance to rivers (m)	<700	40.12	37.5	0.93	0.28	0.64	1
	700–1400	31.9	12.5	0.39	0.12	0.21	0.32
	1400–2800	20.21	25	1.23	0.38	0.09	0.14
	>2800	7.75	25	3.22	1	0.05	0.078
Rainfall (mm/yr)	424 >	45.98	37.5	0.81	0.25	0.33	1
	451–424	11.92	37.5	3.14	1	0.29	0.52
	476–451	12.38	25	2.01	0.64	0.23	0.21
	502–476	10.92	0	0	0	0.08	0.11
	>502	18.8	0	0	0	0.06	0.078
NDVI	0.034 >	0.23	0	0	0	0.21	0.58
	0.17–0.034	77.99	62.5	0.8	0.45	0.36	1
	0.37–0.17	21.43	37.5	1.74	1	0.22	0.61
	>0.37	0.35	0	0	0	0.2	0.55
Lithology	Group (A)	48.21	12.5	0.26	0.1	0.52	1
	Group (B)	36.39	87.5	2.4	1	0.32	0.61
	Group (C)	0.55	0	0	0	0.09	0.17
	Group (D)	14.38	0	0	0	0.06	0.11
Aspect	Flat	18.46	0	0	0	0.02	0.08
	North	6.37	0	0	0	0.11	0.47
	Northeast	10.61	0	0	0	0.07	0.3
	East	5.42	25	4.61	1	0.23	1
	Southeast	5.25	0	0	0	0.1	0.43
	South	14.78	37.5	2.53	0.54	0.11	0.47
	Southwest	10.68	25	2.34	0.5	0.15	0.65
	West	15.46	12.5	0.8	0.17	0.09	0.39
	Northwest	7.27	0	0	0	0.03	0.13
	North	5.64	0	0	0	0.02	0.086
Land cover	Settlement	0.17	0	0	0	0.05	0.15
	Irrigated agriculture	24.45	12.5	0.49	0.04	0.33	1
	Orchard and Forests	2.15	25	11.62	1	0.32	0.96
	Grassland	21.17	37.5	1.77	0.15	0.23	0.69
	Bare soil and rock bodies	51.5	25	0.48	0.04	0.06	0.18
Soil type	Inceptisols/Vertisols	20.71	50	2.41	0.6	0.08	0.16
	Bad Lands	60.64	0	0	0	0.11	0.22
	Rock Outcrops/Entisols	9.23	12.5	1.35	0.33	0.31	0.63
	Rock	9.4	37.5	3.98	1	0.49	1

) was used to rank our considered landslide factors and classes using a GIS environment. If the resulting FR value is higher than 1, this indicates a robust relationship between the polygon of the landslides in the inventory map used for training the models and the total considered area. This, in turn, means there is high landslide susceptibility, whereas FR values lower than 1 refer to low landslide susceptibility. The same holds true for the delivered factor weights based on the FAHP. Thus, higher FAHP weights indicate high landslide susceptibility. To produce the weightings of the FR model, we utilised the prepared inventory map, whereas to derive the weights of the FAHP model, we considered our target expert comparisons in the pairwise matrices. The process of weighting and ranking the sub-factor classes based on both the FR and FAHP models is presented in Figure 13. Although, interestingly, most sub-factor classes displayed the same pattern, there were some significant differences in the first class of altitude, rainfall, and lithology factors. Based on the expert knowledge of our decision-makers, group (A), containing the grey and brown fossiliferous limestone, has the highest importance for the landslide phenomenon, but this class received lower FR quantities. However, the remaining classes were almost the same.

For our integrated model, we used factor weights from the AHP with a combination of the FR values, which we derived for the defined classes of our factors. Thus, the resulting LSM based on the integrated model considered the results of both the FAHP and the FR models. We applied the natural breaks classification technique to classify regions that obtained similar values, distinguished by breakpoints. This technique is widely used for classifying the results of LSMs. Therefore, the classification was based on the resulting values of the same categories, such as values between ‘‘High’’ and ‘‘Moderate’’ susceptibility. Thus, we classified our LSMs into five classes of landslide susceptibility using the natural breaks classification technique (see Figure 14).

5. Accuracy Assessment

5.1. Validation of the Resulting Landslide Inventory Dataset

We conducted a field survey in April 2017 to check the accuracy of the integrated model of DInSAR and OBIA for landslide detection. Accuracy was important for our study because we used these results to train our two LSM models. Therefore, we also used field survey GPS locations to generate a reliable landslide inventory dataset to validate both the performance of the inventory dataset and the resulting LSMs. To validate the resulting displacement map from the interferometry and OBIA approaches, we compared it to the randomly selected 30% (27 GPS points) of landslides in the validation inventory dataset. We reserved the remaining 70% (62 GPS points) of landslides to validate the LSMs resulting from the FAHP, FR, and the integrated models.

To measure the accuracy of the extracted landslide polygons, we applied the validation metric of precision [18,65]. Therefore, we considered the number of landslide GPS locations that match the extracted landslide polygons as true-positives (TP), and the rest as false-positives (FP). Then, we calculated the precision from the following Equation (2):

$$\mathrm{Precision} = \frac{\mathrm{TPs}}{\left(TPs + FPs\right)}$$

(2)

Based on this calculation, our map achieved a precision score of over 94%. The resulting accuracy shows that the displacement map is reliable and can be used as inventory data.

5.2. Validation of Resulting LSMs

Validation is a crucial step in modelling and generating any natural hazard susceptibility map [53,66]. As previously mentioned, we reserved 62 GPS points from the validation inventory dataset specifically for this task. To validate different resulting LSMs within the validation data set, we used receiver operating characteristics (ROC) curves [67]. ROC curves show a full scene of a trade-off between the true-positive rate (TPR) and the false-positive rate (FPR) in the LSMs. The vertical axis of the ROC curve refers to the TPR, whereas the horizontal axis denotes the FPR, represented in Figure 15 [68]. TPRs are the pixels that were correctly classified as landslide susceptible areas, whereas FPRs refer to the pixels that wrongly classified as landslide susceptible areas. To generate ROC curves, we must then obtain the correctly and incorrectly labelled pixels in the inventory data set and plot these as TPRs vs. FPRs. The area under the curve (AUC) is a common accuracy measure to evaluate LSMs. The resulting AUC values designate the ratio between the number of correctly labelled pixels and incorrectly labelled pixels [18]. As such, higher AUC values show that the resulting susceptibility map has a higher accuracy. AUC values close to 1 show a perfect susceptibility assessment, whereas a value of 0.5 indicates that the resulting susceptibility was generated by chance [69]. The AUC values of our study for the FAHP, FR, and the integrated models were 0.83, 0.89, and 0.96, respectively.

6. Discussion

Based on the resulting accuracies, our integrated model produces more accurate susceptibility categories than the FAHP and FR models. Thus, for our specific case, we can conclude that a combination of knowledge-based and data-based models generates a more accurate result than the models used separately. Moreover, the resulting LSM created by the FR model obtained a higher accuracy than that of the FAHP. The higher accuracy of the FR model indicates that our data-based model achieved a better performance than our knowledge-based model. Nevertheless, the results of our knowledge-based models may improve if we used more experts to evaluate the pairwise comparison matrices or used other models, such as the conventional AHP. The accuracy of the maps resulting from the data-based models is highly dependent on the accuracy and amount of the training data, and the accuracy will change when using another section of the inventory dataset. We obtained almost similar results with some related works. Mondal and Maiti (2013) [2], for instance, obtained an overall classification accuracy of 0.92 by integrating the AHP and FR models for landslide susceptibility modelling and mapping of the Shivkhola Watershed, Darjeeling Himalaya.

Table 6. Area of each category and the number of observed landslide events.

Landslide Intensity Zone/Class	Percent of Study Area			No. of Settlements			Landslide GPS Points
	FR	FAHP	FR-FAHP	FR	FAHP	FR-FAHP	FR	FAHP	FR-FAHP
Very high	10.24	10.18	10.7	40	42	50	9	11	24
High	19.91	22.08	19.79	77	70	78	22	18	9
Medium	29.37	26.33	24.94	138	122	100	24	25	19
Low	2584	19.47	24.17	127	92	107	8	10	9
Very low	14.61	22.1	20.38	88	150	141	1	0	3

shows the observed landslides and the number of settlements in each class of the final LSMs. The resulting LSM based on the integrated model contains 24 landslide areas in the resulting very high susceptibility class. However, the LSMs produced by the FAHP and FR models resulted in 11 and 9 very high susceptibility areas, respectively.

From a lithological point of view, groups A and B account for 12.5 and 87.5 percent of the landslide area in our case study area, respectively. Both of these groups contain geo-units of limestone, which is much more susceptible compared to the other geo-units. The last class (20° <) of slope factor showed the highest susceptibility weight and had a significant role on the LSM for our study area; in contrast, the last class of the distance to roads (2800m <) factor showed the lowest impact in on the resulting LSMs. In some more detail, the fourth class of the precipitation criterion (502 mm/yr <), weighting approximately 0.08, had almost no importance for the LSM. The first class of distance to fault criterion (0–1000 m) was found to be the most significant compared to the other classes in this factor. Our results of conditioning weighting and ranking using the AHP also are comparable with similar studies. Pourghasemi and Pradhan (2012), for instance, find the slope, altitude, and aspect factors as the most three important landslide susceptibility conditioning factors, which is completely the same as our results in this study. However, slope degree was identified as the most important factor for their study area of the Haraz watershed, Iran. But in our case, altitude is ranked as the first, based on the overall opinions of our questioned experts.

In this study, we applied an object-based method to extract the landslide areas confirmed by GPS data from an extensive field survey. The extracted landslide polygons were used as the inventory dataset for training and validating the applied approaches. Several studies attempted to use different available methods and techniques to extract landslide areas, such as object-based methods [9], various machine learning techniques [11,70], and deep learning models [12]. Based on the data we used to extract landslide areas in this study, we found a knowledge- and object-based method to be the most practical way to deal with both optical and SAR imageries from Sentinel-2B and Sentinel-1A, respectively. To our knowledge, our study is one of the rare studies that attempted to generate a landslide inventory dataset from satellite imagery and use it for susceptibility modelling and mapping. However, several studies evaluated different approaches for landslide susceptibility modelling only. The FAHP and FR methods have been used separately in several LSM studies. The fuzzy theory approach has been integrated with some multi-criteria decision-making methods, such as the AHP and ANP, to deal with the associated uncertainty and enhance the precision of the resulting susceptibility maps [22,62]. Although the AHP method itself was integrated with the FR method for different susceptibility models by [29], the present study is unique in that it attempts to integrate the capabilities of fuzzy theory, multi-criteria decision making, and FR for LSM. The advantages of our proposed workflow can be summarized as its integration of patterns based on data-driven methods with expert knowledge for improving the accuracy of the susceptibility mapping. This integration could deal with the associated uncertainty within the pairwise comparisons, which is a common issue for spatial decision support systems, such as AHP. At the same time, we showed that the results of a data-driven method can be only based on randomness. Thus, the integration of a data-driven method and expert knowledge could significantly enhance the accuracy of landslide susceptibility mapping. However, this approach cannot be considered as an automated or even a semi-automated method for LSM.

We had two major limitations in our current study. First, we did not have access to some auxiliary data, such as hanging-wall information, for our study area. This is considered a limitation as some previous studies have found that the hanging-wall effect plays a significant role in the damage pattern of landslides and has more impact than the footwall of the thrust [71]. Our second limitation was that we did not have any information about the correlation between the dip of the geo-units and that of the slope, since the agreement between these parameters within an area shows a higher landslide susceptibility compared to other areas with less or no agreement.

7. Conclusions

In this study, we used an integrated model to map the earthquake-induced landslide susceptibility of Kermanshah, Iran. To generate an earthquake-induced landslide inventory dataset, we integrated the DInSAR and OBIA approaches and manually enhanced the dataset using Google Earth™ (Version 7.3.3.7673) images. We then used this resulting inventory dataset to train various models that served as the basis for generating LSMs of the area affected by the Kermanshah earthquake. We also carried out fieldwork and prepared a GPS-point dataset of landslide occurrences, which we randomly divided into two sections to validate the training dataset produced by the integrated model and the final resulting LSMs. We used the following methodology to produce the LSM in this study: first, we used an FR model to calculate the weights of classes that are necessary to aggregate LSM conditioning factors. Next, we applied a FAHP to weight both the factors and their classes. We then followed the standard approach for using the FAHP for GIS- and multi-criteria decision making. This approach is based on a multi-criteria evaluation and a fuzzy set theory, which finds the areas susceptible to landslide phenomena by considering expert knowledge. It should be noted that we first used the FR model to generate an LSM based on the spatial correlation between the results of the interferometry and the class of related factors. This study proposes to evaluate the FR model when integrated with a GIS-based multi-criteria decision-making method for LSM. By applying a validation methodology, we determined the improved certainty of multi-criteria decision making when using an FR model. To validate the resulting LSMs of each model and the integrated model, we applied the ROC method using our landslide inventory data set of GPS points.

The accuracy of the resulting LSM based on the FR model was higher than that of the FAHP model, which shows that the accuracy of the results strongly depends on the quality of the training dataset. The training dataset generated by the integrated model provided somewhat appropriate data for LSM using the FR model. On the other hand, the higher accuracy of the integrated model illustrates that the integration of an optimised decision support system, such as FAHP, with a high-quality training dataset led to an even better LSM result. In summary, the AHP model is a technical model which is sensitive to the inherent uncertainty to expert decision making. Even with the integration of fuzzy rule sets, the FAHP model still relies on pairwise comparisons based on experts’ decisions. Moreover, the accuracy of the FR model directly depends on the quality of the training data. Therefore, the integrated model using a reliable training dataset is less dependent on user decisions and expert knowledge. The drawback of the integrated model is that it requires far more calculations compared to using only one model. The methodologies used in this study can be easily transferred to other areas that are prone to earthquakes and landslides. The results of this study can be useful for prioritising efforts to moderate the adverse consequences of mass movements after earthquakes. Our future work will focus on reducing the damages caused by earthquake-induced landslides in this study area, thereby creating a vulnerability map considering different factors that can show which vulnerable areas coincide with areas highly susceptible to landslides.