Landslide susceptibility mapping using MT-InSAR and AHP enabled GIS-based multi-criteria decision analysis

Abstract Landslide susceptibility maps (LSMs) are generally prepared by integrating multiple prominent thematic layers, including DEM derived products (elevation, slope, and aspect), and other parameters such as lithology, geomorphology, LULC, etc. These parameters can be assigned optimum weights using the analytic hierarchy process (AHP) method, followed by a GIS-based weighted overlay analysis. In recent years, multi-temporal interferometric synthetic aperture radar (MT-InSAR) techniques have been rigorously explored, for land deformation detection and monitoring, by extracting highly stable measurement pixels using tens of SAR acquisitions simultaneously. In this research work, a GIS-based multi-criteria decision analysis to prepare LSMs is proposed, with MT-InSAR derived displacement estimates used as a critical input parameter. An LSM is generated by processing 20 ERS-1/2 and Envisat ASAR images, acquired over ∼120 sq. km wide river basin, located in Uttarakhand, India. The generated LSM is found to be congruent with the susceptible maps made available by the Geological Survey of India (GSI) under the National Landslide Susceptibility Mapping (NLSM) program. Preliminary results indicate that the majority of the unstable zones along the Alaknanda River are correctly identified. The approach is further implemented to generate an updated susceptibility map using 60 scenes of freely available Sentinel-1A dataset, followed by validation through actual field survey. This resulted in the generation of an updated susceptibility map, which helped in the identification of 44.5% new landslide susceptible zones (LSZs). Furthermore, the status of previously identified zones is also quantified. The performance of the proposed approach suggests its usability in generating and updating near-real-time LSMs.


Introduction
Landslides are one of the natural hazards which have caused significant socio-economic damages over the years. There is a need for extensive research to reduce the loss and effectively manage timely rescue efforts. Every year, the road infrastructure development takes place in hilly regions of India (road widening, new route constructions, bypass, road cuttings, etc.), which leads to the development of many steep slopes. This further leads to the destabilization of many landslides, which are susceptible to be triggered by rainfall or earthquake shocks. The above factors motivate for detailed investigation and development of effective methods to map deforming regions and update the unstable zones. The heavily trafficked Char-Dham route, Uttarakhand, India is selected for detailed analysis. Landslide susceptibility maps (LSMs) are generally prepared for landslide-prone areas, categorizing susceptible regions based on input parameters. This further provides a vital source of information tools for geologists, disaster management agencies, and disaster response forces (Nadim et al. 2006;Hong et al. 2007). Factors and information considered in generating LSMs play a critical role in making it more precise and informative. However, careful investigation of the techniques used to provide reliable results with enough accuracy so that accurate LSM could be generated. SAR interferometry has been intensively explored for the mapping of slow-moving landslides. Slow-moving landslides are more prone to become hazardous if some of the triggering factors become active (Gonz alez et al. 2015;Polcari et al. 2016;Zhang et al. 2018;Tiwari et al. 2020). MT-InSAR techniques involve the processing of multi-temporal SAR image acquisitions. The processing reduces spatially uncorrelated phase components and generates a time-series of 1D LOS displacements. The generated displacement maps can be utilized for early warnings (Hooper et al. 2012;Crosetto et al. 2016;Intrieri et al. 2018). The first known MT-InSAR technique is based on the selection of permanent scatterers (Ferretti et al. 2000), PSInSAR TM exploits the amplitude information from SAR acquisitions but requires deformation models. The technique succeeded in estimating deformation trends for urban regions based on its PS selection criteria. However, the method failed in the case of nonurban regions (Ferretti et al. 2000(Ferretti et al. , 2001. Some other variants to this method were also proposed by researchers by utilizing amplitude dispersion and a variety of coherence standards (Berardino et al. 2002;Mora et al. 2003;Werner et al. 2003). For a review of Persistent Scatterer Interferometry (PSI) methods, please refer to Crosetto et al. (2016). Hooper et al. (2004) proposed a unique method, named Stanford Method for Persistent Scatterer interferometry (StaMPS), which utilizes phase information for selection of persistent scatterers (PS) pixels based on phase noise. The StaMPS-based MT-InSAR has been widely used for deformation studies (Hooper et al. 2004;Hung et al. 2011;Cuenca et al. 2013;Motagh et al. 2014). This approach generates displacement estimates for the PS pixels selected based on stable phase history over the observation period (Hooper 2006, Hooper et al. 2007, 2012. Researchers have proposed various methods for the preparation of LSMs. Further, they have examined various causative factors for landslide susceptibility, such as elevation, aspect, slope, lithology, curvature, distance to roads, distance to rivers, distance to faults, stream power index, LULC, lithology, geomorphology, precipitation, soil type, and texture, etc. (Hong et al. 2007;Ahmed 2015;Xiao et al. 2020). However, it has been observed that these methods get biased while assigning the weights to the input parameters. Therefore, to overcome this issue, the analytical hierarchy process (AHP) was introduced (Kavzoglu et al. 2014;El Jazouli et al. 2019;Semlali et al. 2019). Several researchers have also explored AHP fuzzy logic to assign weight priorities (Tazik et al. 2014;Noorollahi et al. 2018).
It has also been noticed that the surface displacement estimates are not considered while preparing LSMs. Researchers have independently used MT-InSAR for landslide inventory preparation. However, no efforts have been made towards exploring the contribution of the MT-InSAR estimates in the LSM generation, along with other factors. Solari et al. (2019) utilized hotspot analysis based on a 1D LOS velocity map generated from MT-InSAR processing, to detect unstable zones. The authors recommended the requirement of more factors, such as in situ measurements, optical images, etc. to categorize the moving regions as LSZs. Zhao et al. (2019) applied an SBAS-InSAR approach for LSM optimization results using the random forest method of machine learning. Liu et al. (2018) proposed a method for landslide inventory mapping by utilizing MT-InSAR results obtained from multi-frequency SAR datasets, SAR intensity images, and DEM derived maps. Solari et al. (2020) proposed a method for the detection of active deformation areas (ADA) based on landslide terminology given by Cruden and Varnes (1996). Since MT-InSAR is considered as a powerful tool to investigate slow-moving landslides, its estimates can be a critical input for detecting susceptible zones when used along with landslide influencing and triggering factors, such as geomorphology, lithology, rainfall, LULC, slope, aspect, etc.
The objective of this research work is to generate an updated LSM for the Srinagar-Rudraprayag region, Uttarakhand using MT-InSAR and AHP enabled multicriteria decision analysis. To accomplish this, the proposed methodology is initially validated using landslide points with identified susceptible landslide zones, obtained from GSI (Geological Survey of India 2020). Similar approach is implemented on a stack of 60 Sentinel-1A images, to generate an updated LSM. To highlight the efficacy of the proposed method, a comparison with LSM prepared with and without the inclusion of the MT-InSAR results is presented. Augmenting the previous research works on LSM generation, this study uses multi-constellation SAR satellites for deformation analysis and landslide inventory preparation based on multi-criteria decision analysis. The proposed method uses MT-InSAR derived displacement estimates along with commonly used different thematic layers, such as geomorphology, lithology, elevation, slope, aspect, and LULC.

Study area
Uttarakhand has witnessed a large number of landslides each year, mostly caused by anthropogenic activities creating unstable slopes along dozer roads. The landslides are generally triggered by intense rainfall and underlying earthquakes. GSI has declared most of the regions of the study area shown in Figure 1 as landslide-prone regions, on account of frequent landslides and the presence of faults such as Main Central Trust (MCT) and North Almora Thrust (NAT). This study focuses on landslides present between the Srinagar and Rudraprayag area along the National Highway 58 (NH-58). To prepare an LSM, multi-SAR constellations are utilized, as given in Table  1. The thematic layers used in this work are listed in Table 2.

MT-InSAR processing
The MT-InSAR processing chain utilizes temporal SAR acquisitions for the generation of 1D-LOS displacement estimates. The chain can be implemented using an integrated open-source platform, comprising the Sentinels Application Platform (SNAP) and StaMPS, apart from other proprietary software. A stack of 20 Envisat and ERS1/2 scenes and 60 Sentinel-1A scenes, covering periods of 1996-2010 and 2017-2019, respectively, are initially pre-processed in SNAP. Later, these are exported to and processed in the StaMPS environment to generate displacement maps. The step-by-step processing flowchart is given in Figure 2. Interested readers can follow Hooper et al. (2004Hooper et al. ( , 2007 for detailed understanding. SAR acquisitions of ERS1/2 and Envisat consist of dissimilar geometries. Therefore, it is not possible to co-register scenes in the same way as for the scenes  obtained by the same sensor or a sensor of similar characteristics. The tool available in SNAP called 'Cross Interferometry Resampling Operator' has been used to coregister multi-geometry acquisitions and to generate single master cross-  interferograms to a master scene. The master is selected in such a way that it minimizes spatial, temporal, Doppler, and thermal decorrelation. Sentinel-1A SAR acquisitions collected in TOPS mode (Terrain Observation with Progressive Scan) with VV (Vertical-Vertical) polarization have been utilized. SAR acquisitions in this mode are distributed in three sub-swaths (IW1, IW2, and IW3) and multiple bursts. Therefore, pre-processing of Sentinel-1A imagery includes two additional steps, splitting (selection of sub-swath which contains study region) and de-bursting (to merge multiple bursts). The steps included in the pre-processing of data are splitting, orbital correction, image stacking, co-registration, and single master interferometric stack generation.
For a non-zero baseline, the topographic phase includes the phase introduced due to deformation which makes the extraction of deformation signal complex (Hanssen 2001). Therefore, differential interferograms have been generated by subtraction of 1 arc-second SRTM digital elevation model (DEM) from single master interferograms to remove the topographic signal. Further, the topographic phase-corrected single master interferograms are exported in gamma format for processing in StaMPS environment (Dwivedi et al. 2015;Tiwari et al. 2016).
The StaMPS technique can be summarized as a four-segment process: (i) Interferogram formation to reduce decorrelation, (ii) PS identification for the initial selection of PS candidates based on amplitude dispersion index, (iii) PS selection based on PS phase stability, and (iv) Displacement estimation (Hooper 2006). In MT-InSAR processing, differential interferograms generated by SNAP-StaMPS integration, are exported to StaMPS followed by PS Identification, PS selection, and Displacement estimation.
Initially, a set of PS candidates (PSC) is selected on the basis of amplitude dispersion index ðD A Þ, which is defined as where r A and l A represent standard deviation and mean of amplitude values correspondingly. Further, the phase history of the PSC is analysed and a probability of being a PS is calculated for every PSC. Subsequently, PS pixels are selected based on PS probabilities, and pixels which are not persistent throughout all the interferograms, or pixels showing dominant behaviour by the scattering of adjacent pixels are rejected. For selected PS pixels, various phase noises correlated and uncorrelated (in time and space both), are removed to make accurate unwrapping. Further, the deformation signal is isolated by the 3D unwrapping of the interferometric phase for all the PS pixels. Lastly, 1D LOS displacement estimates are generated (Hooper et al. 2007). For better interpretation of displacements, it will be advantageous to project 1D LOS displacements (V l ) to downward slope direction (V s ).
V s ¼ V l sin b cos h sin k cos g À sin b sin h cos k cos g þ cos b sin k where g denotes slope, k is aspect, b represents incidence angle, and h is satellite heading with respect to true north, respectively (Zhao et al. 2019).

Analytic hierarchy process (AHP)
Decision-making problems require assistance from many potential factors. Optimal assignment of priority to the influencing factors can lead to better decisions. Multicriteria analysis, which includes a list of potential factors, may lead to optimal results if optimal weights are assigned. AHP is a structure-based decision making procedure proposed by Saaty (1980), which solves complex decision-making problems based on a simple mathematical model. Researchers have appreciated the capability of AHP in assigning optimal weights to input factors and obtaining suitable results such as landslide susceptibility mapping (Moradi and Rezaei 2014;Tazik et al. 2014;Kumar and Anbalagan 2016;Noorollahi et al. 2018;El Jazouli et al. 2019;Nguyen and Liu 2019;Semlali et al. 2019). The concept of AHP relies on the relative importance of factors defined on a scale of importance. Because of its unique consistency test, this method assigns reasonable weights to factors under evaluation. The following section presents a brief description of the AHP method.
Successive steps of AHP implementation include development of hierarchical structure, relative comparison of factors on scale of importance contrary to objectives, computation of weights for each factor, calculation of consistency vector (CV), computation of average value of CV and Consistency index (CI), selection of Random Consistency Index (RI) with respect to number of factors, calculation of Consistency ratio (CR), comparison of CR with CR threshold to conclude the suitability of relative importance decided between factors and finally, decision making.
Initially, a pairwise comparison matrix (PCM) is prepared on the basis of importance of each input factor with respect to other potential factors (Table 3).
For the PCM of n factors represented by Equation (2), the weights of each factors could be calculated as shown in Equation (3). For the evaluation of proposed PCM, a consistency check is performed by CV (Equation (4)).  where RI is calculated from the randomly generated PCM and CI represents deviation in consistency which is computed using k (the average of the elements of CV).
CI ¼ k À n n À 1 If CR (a measure of acceptability of any proposed PCM) exceeds a specified threshold, the proposed PCM is rejected, and the relative importance of the included factors is reconsidered. If CR is below the threshold, as mentioned in Table 4, the assigned weights are considered reliable and acceptable.
For the implementation of AHP, a good literature survey is required for the selection of factors to be included in susceptibility mapping and for defining priority in between factors i.e. relative importance. There are greater chances of getting inconsistent results. Hence, more than one trial is recommended for reaching the CR (threshold), thereby achieving reasonable results. Hence, an expert suggestion could make this process easier. The possibility of achieving different results by different Table 4. RI and CR values based on numbers of included n factors (Saaty 1994 groups using the same factors for a study site could not be ignored. This is because the threshold value allows CR to be less than 10%, which is possible by altered PCMs. In this research, AHP has played a crucial role in weight assignment of thematic layers, provided that a new factor (MT-InSAR result) is included.

Spatial data integration within GIS
Spatial data integration deals with operations on different types of geospatial datasets (thematic layers) for visualization or querying purposes. Spatial analysis can be performed on the detection of landslide-prone regions, based on influencing factors. The 1D LOS displacement results can be integrated with slope, aspect, lithology, geomorphology, rainfall, etc. to highlight LSZs. For the development of the updated LSZs using overlay analysis, weightage for all thematic layers is adapted from the AHP method, as shown in Table 5. The methodology of spatial data integration is shown in Figure 3.

MT-InSAR results
To obtain LSZs, the two Single Look Complex (SLC) image stacks (Stack-1: combined 20 scenes of Envisat ASAR and ERS1/2, Stack-2: 60 scenes of Sentinel-1A TOPS) have been processed using the MT-InSAR technique. Master selection of both the stacks has been done by computing the highest sum-correlation of all possible single master interferometric pairs. As a result, 19 interferograms from Stack-1 and 59 interferograms from Stack-2 are generated. Further the topographic signal is removed by using 1 arc-second SRTM DEM (30 m). Through interferometric processing, $11,500 and $3600 PS pixels are detected from the two stacks, respectively. 1D LOS mean displacement map from the ERS1/2 and Envisat dataset along with displacement rate projected along slope, coherence map and master atmospheric error are shown in Figure 4(a-d) with a range of deformation rate from À4.35 to 8.61 mm/yr. PS pixels having a deformation rate between À1.5 mm/yr and 1.5 mm/yr are considered as stable. PS pixels of values less than À1.5 mm/yr, with a mean deformation rate of À2.53 mm/yr, can be considered as unstable. An initial investigation of the MT-InSAR results with pre-existing LSM prepared by NLSM (GSI), confirms the presence of unstable regions in the proximity to LSZs. However, only based on this analysis, the identified moving regions cannot be declared as LSZs. To do so, other essential factors such as optical imagery datasets, geomorphology, lithology, LULC maps and DEM derivatives, presence of faults, rainfall, distance from roads, etc. have to be considered. Figure 4(e-h) shows the 1D LOS mean displacement map, displacement rate projected along slope, coherence map and master atmospheric error obtained from MT-InSAR processing of the Sentinel-1A dataset. The deformation rate varies from À11.34 to 11.88 mm/yr. High negative deformation rates from À11.34 mm/yr to À1.5 mm/yr with mean value of À4.42 mm/yr emphasize availability of unstable regions. Analysis of the 1D-LOS mean velocities of both the results indicates the existence of subsidence and uplift in the study area. Further, PS pixels detected in the nearby regions of Srinagar city show movement towards the satellite LOS. PS pixels of the South-East region in the study area show movement away from the satellite LOS. These detected movements indicate instability in the region.

AHP weight analysis
The pairwise comparison matrix has been prepared according to the significance of thematic layers. The methodology for calculating weight for all factors is also presented. LSM is prepared for two cases: with and without including MT-InSAR displacement maps. The idea is to see the impact of MT-InSAR derived displacements on the preparation of LSM. Tables 5 and 6 present priorities given to factors considered for the generation of LSMs. Consistency ratio (CR) of 0.046 and 0.057 (less than 0.1 for n ¼ 7 and 6, respectively) confirm the reliability of weights assigned to each factor. Assigned weights for both the cases are shown in Equations (7) and (8).   each factor is weighted, considering their effect in the decision making procedure. For example, high slopes and high elevation regions are considered more prone to landslide than other respective regions present in both the layers. SE and NW directions are preferred in aspect, as suggested in a previous research work by Khanduri et al. (2018), who detected NW-SE deformation trend in the study region. From the lithology point of view, limestone, dolomite, quartzite, and phyllites are abundant in the study region and possess more number of landslides. Therefore, the above-mentioned sub-factors are given priority in re-classification. Similarly, piedmont alluvial plain and highly dissected structural hills and valleys in geomorphology, and mixed forest, forest, and savannas are given priority in the LULC layer.
For spatial integration analysis, a 1D LOS displacement map of ERS1/2 and Envisat and all thematic layers, mentioned above, are used as input layers to detect susceptible landslide zones. Figure 6 shows the output of the weighted overlay  operation, i.e. LSZs in the study area. The generated results are congruent to the high LSZs provided by the GSI, wherein some of the prominent landslide zones are highlighted. For example, the Sirobagarh landslide, and few slides along the Alaknanda River are identified. Further, based on the analysis of the results, a total of 11.76 and 13.42 square km of the study area is highlighted as high and medium susceptible zones, respectively, by the GSI. On the other hand, the proposed approach emphasizes to 2.92 and 20.45 square km area as high and medium susceptible zones, respectively, and identifies 138 high and 623 medium LSZs.

Updated LSM
It is observed that the proposed approach efficiently detects the affected landslide zones. Further, the susceptible landslide zones generated using Sentinel-1A SAR acquisitions of time-interval 2017-19 improve the understanding of the current scenario of landslide activities in the study site. Figure 7 shows high and medium LSZs and highlights 4.54 and 25.62 square km area as updated high and medium susceptible zones. Comparing the present LSM with previously obtained LSM generated using Envisat and ERS1/2, 131 high and 502 medium LSZs have been identified in which 27% (11.39 square km) are active. Figure 8 shows the current status of LSZs  (refer to Table 7). Based on the density of the newly identified LSZs, it is also observed that the upper part of the NAT fault line appears more active than in the past. Some of the newly identified susceptible zones are observed along the Alaknanda River. This may be due to regular highway construction activities along National Highway (NH-58), passing alongside the river. Figure 9 shows the proportion of LSZs in all unstable zones. In the field visit conducted for validation, we noticed that new road constructions and their widening have also resulted in the development of many unstable slopes along the way, causing fear of failures that might be triggered in monsoon. Figure 10 shows some of the landslide sites observed during field visit in Novemebr-2019 which have been used for validation of the generated LSM. Further, to present the potential effect of MT-InSAR results in LSM preparation, a comparison is made between LSM generated with and without including MT-InSAR results ( Figure 11). LSM generated without including MT-InSAR  contains a large number of LSZs than that available as ground truth. Hence, it is evident that the inclusion of MT-InSAR displacement maps, as input in LSM preparation, substantially improved our understanding in terms of terrain susceptibility for a landslide.

Conclusions
MT-InSAR displacement results are ambiguous to detect unstable zones. The proposed approach of utilizing MT-InSAR derived displacements along with other influencing parameters has substantially improved the mapping of susceptible landslide zones. The proposed approach performed efficiently well in identifying the activity of landslide susceptibility zones for 10 years. The approach can be further used for updating existing maps and for generating new LSMs. Among the total zones which were identified from the updated LSM, 44.5% (16.35 square km) are newly identified, 28.5% (28.70 square km) are stabilized and 27% (11.39 square km) zones are still active. A comparison between LSMs with and without MT-InSAR results is also presented to evaluate its contribution to the overall LSM preparation. The overall performance of this approach suggests its usability in generating and updating LSMs developed by state-of-the-art techniques. The methodology, added with quick data processing mechanisms, can also be utilized for generating near real-time LSM in the future.