Measuring landslide vulnerability status of Chukha, Bhutan using deep learning algorithms

Landslides are major natural hazards that have a wide impact on human life, property, and natural environment. This study is intended to provide an improved framework for the assessment of landslide vulnerability mapping (LVM) in Chukha Dzongkhags (district) of Bhutan. Both physical (22 nos.) and social (9 nos.) conditioning factors were considered to model vulnerability using deep learning neural network (DLNN), artificial neural network (ANN) and convolution neural network (CNN) approaches. Selection of the factors was conceded by the collinearity test and information gain ratio. Using Google Earth images, official data, and field inquiry a total of 350 (present and historical) landslides were recorded and training and validation sets were prepared following the 70:30 ratio. Nine LVMs were produced i.e. a landslide susceptibility (LS), one social vulnerability (SV) and a relative vulnerability (RLV) map for each model. The performance of the models was evaluated by area under curve (AUC) of receiver operating characteristics (ROC), relative landslide density index (R-index) and different statistical measures. The combined vulnerability map of social and physical factors using CNN (CNN-RLV) had the highest goodness-of-fit and excellent performance (AUC = 0.921, 0.928) followed by DLNN and ANN models. This approach of combined physical and social factors create an appropriate and more accurate LVM that may—support landslide prediction and management.


Materials and methods
Study area. The Chukha Dzongkhag (district) lies in the southern foothills of the Bhutan Himalayas (Fig. 1).
It shares border with the neighbouring Indian state of West Bengal and thus serves as a vital route for mutual trading between the two countries. The major border town of Phuentsholing is the gateway city that connects India to Western Bhutan. Some of the important hydropower plants such as Chukha and Tala are located in the Chukha district. Because of proximity to India and ease of accessibility to large Indian market, many of the country's industrial infrastructures are located along the foothill flat areas of Chukha district. From the statistic report of Bhutan Govt., an estimated 4.22% of the entire district (1879.77 km 2 ) comprises of fields cultivated with oranges and potatoes that constitutes the main source of rural cash income.
In the study area, the monsoon season extends from June-September of each year and about 78% of the annual rainfall occurs during this period, making the area highly vulnerable to landslide hazards and flooding. In addition, the rugged mountainous terrain with steep slopes further exacerbates the condition that triggers the landslides, exposing the area to the vagaries of weather events. The landslides obstruct the national highway (Phuentsholing-Thimphu Highway) which is the main trade route of the country, causing heavy losses to lives and infrastructures ( Table 1). One of the major landslides in the region was along the Phuentsholing-Thimphu highway (26.85 N, 89.33 E to 27.15 N, 89.55 E) just after the 2016 monsoon. The intensity and frequency of such events are expected to increase with climate change, resulting in an increasing risk to the residents and also affecting the economy of the country.
Methodology. The steps followed in this work are presented in Fig. 2.
(1) Prepared landslide inventory map based on landslide locations identified through field investigation. (2) Collected data for several geo-environmental and socio-economic factors.
(3) Selected factors using collinearity test and information gain ratio (IGR) methods.

Preparation of landslide inventory map (LIM).
Mapping of the past and present records of landslide events is called an inventory map 24 . In this study, these records were collected from the Project Dantak, Border Road Organisation (BRO), Govt. of India, and from two Royal Government of Bhutan organizations-Department of Geology and Mines (DGM) and National Centre for Hydrology and Meteorology (NCHM). For location verification, field survey was carried out using Global Positioning System devices in January 2020 (Fig. 3). A total of 350 landslides were mapped for modelling the relative landslide vulnerability, of which 50.4% are rock falls, 40.2% are debris slides and 9.4% are rotational slides. Correspondingly, the same amount of non-landslide points was randomly selected for training and validating the models. Out of the total landslide location, 70% was used as training dataset and the remaining 30% as testing dataset 24 .

Preparation of the landslide vulnerability conditioning factors (LVCFs). In this study, 22
geo-environment and 9 socio-economic factors were used for landslide vulnerability modelling which were obtained from the high-resolution (12.5 × 12.5 m) PALSAR DEM of Alaska DEM facility, Landsat 8OLI/TIRS of 30 m × 30 m resolution from USGS Earth Explorer, the demographic data from the National Statistics Bureau, Bhutan, rainfall data of last five years from National Centre for Hydrology and Meteorology, and geological map of 1:500,000 scale from Department of Geology and Mines, Royal Government of Bhutan. The soil resistivity and textural classes (clay, silt and sand) data were extracted from laboratory tests of collected samples from the study area. Data processing and modelling were performed using software such as SPSS, MS-Excel, QGIS 3.16, and R studio.
Geo-environmental factors. The elevations of the study area vary from a maximum of 4411 m to as low as 98 m (Fig. 4a) and the slope ranges from 0° to 80° (Fig. 4b). Aspect of slope is classified into nine categories namely, Flat, North, South, West, East, Northeast, Southeast, Northwest, Southwest (Fig. 4c). The spatial distribution of plan curvature and profile curvature ranges from − 32.89 to 32.12 (Fig. 4d) and from − 45.17 to 43.52 (Fig. 4e). The spatial value of convergence index (CI) ranges from − 92.25 to 89.53 (Fig. 4f) and topographical position index varies between − 17.59 and 26.35 (Fig. 4g). The value of terrain ruggedness ranges from 0 to 40 (Fig. 4h). The spatial value of topographical wetness index varies between 0.68 and 25.41 ( Fig. 4i) and valley depth varies from 0 to 817 (Fig. 4j). The length of slope and value of relative slope position of this study area ranges between 0 and 119 ( Fig. 4k) and from 0 to 1 (Fig. 4l), respectively. Rainfall map was prepared based on the last 5 years' average annual precipitation data of the different stations using the Inverse Distance Weighted interpolation method. The average maximum and minimum precipitations for this study area are 4130 mm and 1357 mm, respectively (Fig. 4n). The land use/land cover (LU/LC) map of the district has been prepared from Landsat 8 OLI/TIRS satellite imagery following the supervised maximum likelihood classification method. The LU/LC categories are broad leave forest, build-up areas, mixed forest, grassland, miscellaneous, conifer forest, and agriculture respectively (Fig. 4p). The normalized differential vegetation index (NDVI) was also generated using RED and NIR band of Landsat 8OLI/TIRS imagery and values differs from − 0.33 to 0.81 (Fig. 4q). The positive values of NDVI indicate the dense vegetation cover areas, while the negative values indicate low vegetation cover areas. The geology of the study can be categorized into Buxa Group, Daling-Shumar Group, the Lesser Himalayan zone, Paro formation, and the structurally lower Greater Himalayan zone (Fig. 4r). The soil resistivity values were determined by surface electrical resistivity method using soil samples and the map depicting resistivity www.nature.com/scientificreports/ distribution was prepared applying Inverse Distance Weighted (IDW) interpolation method. The soil resistivity spatially varies from 276 to 799 Ω-m (Fig. 4v). Following the same method sand, clay, and silt maps were also prepared. In different parts of the district the percentage of sand, silt, and clay was found to be 34-50%, 12-30%, and 28-45%, respectively (Fig. 4s, 4t, 4u). Drainage of the study area has been extracted from the open series topographical maps and DEM. Buffering tool of QGIS 3.16 was used to prepare the distance from the river map ( Fig. 4m) and distance from lineament map (Fig. 4o).
Socio-economic factors. To identify the social vulnerability of landslides, socio-economic factors were chosen. The distance from road was prepared using the Euclidian distance buffering tool (Fig. 5f). The demographic data was derived from the district headquarter of Chukha district, Bhutan. The thematic maps of the population density, old population density, house frequency, literacy rate, medical facility, disability prevalence rate, number of households that require at least 30 min to reach the road head ( Fig. 5a-i) were prepared in GIS platform based on the block-wise data.   26 . Based on the diagnosis, a total of 22 geoenvironment and 9 socio-economic factors have been selected to proceed with the modelling.
Information gain ratio (IGR) method. The information gain ratio (IGR) is one of the machine learning techniques 27 which evaluates the correlations of landslide occurrence with landslide conditioning factors and the role of these factors in the frequency of their correlations 28 . The IGR is employed to reduce a bias to multi-value property, taking into account the size and number of sections, when selecting a function. The IGR for element TWI, for example, is estimated as:   www.nature.com/scientificreports/ www.nature.com/scientificreports/

Applied deep learning and benchmark machine learning approaches
Artificial neural network (ANN). One of the most popular ANN used for landslide susceptibility is the Multi-layer perceptron Neural Network (MLP-NN) 29 . Basically, three layers which include an input layer, one hidden layer, and an output layer were involved in the topology of MLP-NN models 30 . Although their output is regulated by their structure, the activation functions and the updating of the link weights between the processing components 31 . The MLP-NN was introduced in the current study based on the "RSNNS" R package 32 and by using a grid search technique, the number of the hidden processing elements was tuned.

Convolution neural network (CNN).
DLAs are modelled on the structure of the human brain and based on an ANN. CNN, introduced by LeCun et al. 33 is a well-known deep learning algorithms (DLA). Recently, a number of disciplines, including earth science, have been increasingly using CNN for classification and prediction 34 . The utilisation of several layers, pooling, local connections, and mutual weighting distinguishes CNN from traditional neural networks. The central concept behind CNN is that images are used as input parameters. The set of indicators can be greatly decreased and processing can be accelerated. The convolution layers (CLs), pooling layers (PLs), and linear rectified unit layers are involved in the typical structure of CNN model (Fig. 6). CLs deliver the best data classification results by learning the convolutions. By decreasing the quantity of convolution architectures, PLs control overfitting that allows consistent conversion and enhances computational efficiency 35 . The ReLU improves the network's nonlinear capabilities by ReLU activation. Researchers have created and applied numerous analysis structures based on data form, image, and purpose, including GoogleNet 36 , ZFNet 37 , and VGGNet 38 among the more common structures. Several articles have clarified these layer forms, their basic learning criteria, and how the CNN processes the data for training 33 . The CNN-2D structure is used here because it is relevant for earth science studies 39 . The input data in CNN must be 1-D images: a 1-D input grid cell (vector) with different features which be translated into a 2-D input grid cell to ensure optimum initialization efficiency and this technique is used for mapping landslide vulnerability. In current research, we compare the number of LCVFs to the number of attribute values for each factor, and the larger of the two numbers is used to determine the size of the related two-dimensional matrix. The research region, for example, contains 25 geological categories, which is more than the number of LCVFs. As a result, for each grid cell, we created a 25 × 25 matrix. Since there is no sorting of data and analysis is continuous, the images are large.  www.nature.com/scientificreports/ Deep learning neural networks (DLNN). The major advantage of DLNN model is that it uses raw dataset to create a high-level function 18 . DLNN comprises of three layers-an input layer, followed by hidden layers which lead to output layers 40 . Figure 7 demonstrates the conceptual setup of the DLNN model included in this analysis. The overall pattern of DLNN model is to work in such a way that the input layer delivers the signals that are diverse landslide factors, and after processing and interpretation of this information in multiple hidden layers, the impacts are shown in the model's last year, the output layer. The output layer has two probable labels, i.e. a negative label (non-landslides) and a positive label (landslides). From the last hidden layer, these classification results are collected and displayed in the output layer 41 . DLNN has unique advantages over the conventional ML algorithm, and therefore much more attention has been put on the use of the DLNN model in the area of prediction analysis. DLNN outperforms all other ML models in several ways, by making optimum use of unstructured data by specific observations to recognize the training dataset, being versatile enough as to identify new data, and being able to create new learning models by introducing more layers to the neural network architecture. Mathematical equation mentioned below was applied in DLNN as per the Kim 40 : where x is an input signal, and h is an activation function. The following can be represented on the basis of the ReLU activation function as: The cost function is the distinction between class outcomes that are experiential and expected. The loss function (L) of a cross-entropy is used to detect patterns and is given by: where, the number of the training data sets is expressed by N D , T indicates the class outputs detected and Y displays the class outputs expected.

Methods used for validating the models
Receiver operating characteristics (ROC) curve. For any two distinct vectors where the first vector describes the binary presence-absence state of a particular action and the second vector gives the related probability predictions, the ROC curve can be prepared which is a common cut-off dependent diagnosis 42,43 . The overall model success can be recognized based on AUC values, as per the classification provided by Hosmer and Lemeshow 44 . There are four major elements of ROC curve: e true positive (P), false positive (Q), false negative (R) and true negative. The following measures including sensitivity or true positive rate (TPR), false positive rate (FPR), true negative rate or specificity (TNR), false negative rate or miss rate (FNR), efficiency, precision, negative predictive value (NPV), Matthews correlation coefficient (MCC) and Cohen's Kappa have been computed from these elements for validating the models: where N is sample size, P and P is predicted and observed values of dependent variable, respectively.

MAE.
Mean absolute error (MAE) is measured as the sum of the differences between the expected value and the actual value of the model, without taking their direction into account: where N is sample size, P and P are predicted and observed values of dependent variable, respectively.
Relative landslide density (R-Index). R-index was used to evaluate the vulnerability maps of landslides.
R-index is calculated as following 45 .
where xi is the percentage of the area that is vulnerable to landslides in each vulnerability class, and Xi is the percentage of landslides in each vulnerability class. The maximum values of these vulnerability classes indicate the highest goodness-of-fit and excellent accuracy 46 .

Results
Multicollinearity assessment. The findings of collinearity results indicate that no linearity exists among the LVCFs because the Tolerance and VIF values of these factors do not surpass their limits (Table 2) which suggests the aptness for inclusion in the modelling.

Results of IGR.
The result of IGR method is shown in Fig. 8. The IGR values of the geo-environmental LVCFs for landslide vulnerability prediction were higher for geology, soil resistivity, rainfall, and elevation (Fig. 8a) and for socio-economic LVCFs, it is higher for distance to road and population density (Fig. 8b).
Landslide vulnerability analysis. Using the three models i.e. ANN, CNN, and DLNN and considering three aspects of vulnerability i.e. geo-environmental, socio-economic, and relative or overall vulnerability, of area as most vulnerable to landslides. The very high vulnerable zones have mostly occurred in the southern, south-west and south-east portion of the district and very low vulnerability is found in northern, north-east, and north-west regions. These results appear to be associated with the presence of weak geology, heavy rainfall during monsoon, and Phuentsholing-Thimphu national highway that passes through the first three portions of the district. Apart from this, rapid population growth in the above areas due to the fast development of Phuentsholing, the business city of Bhutan, resulted in new anthropogenic activities and thus weakening the soil (Fig. 9). The spatial extension of other vulnerability classes is given in Fig. 10. Table 3 shows the comprehensive results of the ROC curve and other validation measures. The AUC computed using train data indicates success rate and AUC computed using test or validation data denotes prediction rate of the models 47 Table 3). The result of reliability measures such as MAE and RMSE value is lowest for CNN-RLV followed by DLNN-RLV and ANN-RLV model. The R-index values for each model using geo-environmental aspect, socio-economic aspect or overall factors, is highest for very high vulnerability class followed by high vulnerability class (Table 4; Fig. 11). www.nature.com/scientificreports/ Factors importance analysis. The importance of the LVCFs in predicting landslide vulnerability was evaluated by random forest (RF) and chi-square attribute evolution technique (CSAE). The outcome of RF showed that the geology, rainfall, soil resistivity, and elevation were the most predictive factors for landslide vulnerability modelling in this research, followed by the sand and silt distribution and other geo-environment LVCFs (Fig. 12a). Among the socio-economic LVCFs, distance to road, population density factors have the most importance (Fig. 12c). The CSAE method yielded similar results to the RF model (Fig. 12b,d).

Discussions
Identifying or mapping the areas with future landslide potential is one of the most useful document for appropriate land use planning and mitigation decision making 13,47,48 . Evaluation of landslide susceptibility and vulnerability of a hilly area (as mountainous areas are subjected to the landslides) is indeed, necessary as it will serve as an essential dataset which will be used for identifying the points and places of relatively high landslide susceptibility and vulnerability which can be utilized for efficient and safe planning, design as well as the construction activities in an identified landslide zone. In the studied region, as well as other similar locations, there has been a significant growth in human socio-economic activities as well as other geo-environmental variables in recent decades. As a result, the frequency with which landslides occur in certain locations has risen. Such occurrences result in a significant loss of human life and property, as well as it results in a negative impact on different sectors such as tourism and infrastructure development 49 . A better working LVM has always been a subject of significant relevance within the community of numerous scholars across the world who are researching and investigating landslides. This is due to the fact that the methodology and conditioning factors used can have a substantial impact on the models' predictive performance 48 . It is very relevant to employ the mechanisms of factor selection to maximize the efficiency of landslide models by eliminating unwanted or trivial variables before training those 50 . CA and IGR method was adopted to accomplish this. In this study, CA revealed that all of the conditioning variables were effective and independent. In order to pick the most appropriate conditioning factor, there must be a non-collinear link between the LVCFs, which may be done via multicollinearity analysis 51 . Applying CA and IGR analysis, this study took into account a total of twenty-two geo-environmental and nine socio-economic factors, and the incorporation of these LVCFs has been justified in several studies 4,5,21 . It is very important to select LVCFs associated with past and present landslide events in vulnerability modelling. Although www.nature.com/scientificreports/ landslides are quasi-natural hazard but when it results in destruction and harm to the economy and society, it indicates the association of socio-economic causes in the event happened. Therefore, to enhance the performance of modelling, this research was conducted taking into consideration both physical and social aspects; these factors are found to be very significant, especially in hillslope regions 52 .
The study was approached towards mapping the landslide vulnerability considering geo-environmental and socio-economic aspects to improve the classification accuracy using deep learning and benchmark machine learning methods viz. CNN, DLNN and ANN. In all cases of vulnerability mapping (Physical, socio-economic In this study, ANN-RLV model has assured 88.7% success rate by the application of the AUC measure which is similar to the aforementioned study. In landslide susceptibility assessments, Bui et al. 9 used a DLNN model and compared its predictive efficiency with state-of-the-art machine learning models in Kon Tum province, Vietnam. Using ROC curve, the performance of the models revealed that the DLNN model had the highest goodness-of-fit and outperforming ANN, SVM model. Relatively better performance of DLNN than of ANN was also achieved in landslide vulnerability mapping for the present research. The efficiency of deep learning models compared to ML models was found to be better as per the study of Yao et al. 54 where authors have developed the deep neural network model based on semi-supervised analysis (SSL-DNN) for the landslide susceptibility estimation. For comparison, supervised models were introduced, including deep neural network (DNN), SVM, and logistic regression (LR). The result revealed that all comparable models were surpassed by the proposed SSL-DNN (AUC = 0.898) which is greatly supportive of outcome of the deliberating research. Application and enhanced competence of DLNN model can also be found in other hazard vulnerability and susceptibility mapping as Band et al. 55 proposed a DLNN model and an ensemble particle swarm optimization (PSO) algorithm with DLNN (PSO-DLNN), for gully erosion susceptibility mapping. These models were compared with ANN and SVM model. The PSO-DLNN model has the highest efficiency followed by the ANN and SVM. Therefore, the derived outcome of this research has a similar covenant to the very previous studies delegating relevance of adopted methodological scheme. Above all, r, the convolution neural network (CNN) model has provided superior results outperforming DLNN and ANN models in all of the social vulnerability, landslide susceptibility and relative or combined vulnerability assessment as exposed by all the adopted validation measures both in training set and testing set based analysis (Table 4). Therefore, accuracy of the DLNN model was better than the www.nature.com/scientificreports/ conventional ANN machine learning technique. This is because greater number of samples and huge data could be handled by this model and the outcomes can be estimated with greater precision. The key benefit of deep learning is its formal system of self-governing DLNN layer organization learning. The conventional machine learning is incapable of processing such a vast number of inputs, and the result is less optimal in comparison to the deep learning model 55 . The accuracy of machine learning for different purposes was substantially enhanced in deep learning systems 55 . Yi et al. 50 developed a convolutional neural network (CNN) model for the spatial prediction of landslides. The result of CNN was compared with three conventional ML algorithms, i.e., logistic regression, multilayer perceptron (MLP) neural network and radial basis function (RBF) neural network which found CNN as the best fitted and excellent predictive model, followed by the MLP, logistic regression RBF. Wang et al. 39 applied the deep learning and ML models such as logistic regression, SVM, RF models for LS assessment. The result of that research also proved that CNN had the highest performance of predictive modelling followed by the ML models. With the agreement of the results of these studies, the present work also confirms the relatively highest adaptability of CNN model in deriving LVMs as reflected in the produced results of validation and accuracy measures. ReLU activation function was applied in the present study considering the aforesaid literature. The advantages of CNN model are that it considers all the neighbourhood information and can determine manifold stages of representations from input data 56. It maintains the association of pixels using several factors and identifying internal elements 59 . Following Jenk's algorithm of natural breaks classification, LVMs were divided into five groups of susceptibility classes. This strategy of clustering data helps to reduce the mean-variance of each class from the mean within class range and to increase the discrepancy between each class from the means of the other classes 57 . From the analysis of the LVMs of this study, the very high vulnerability zone of landslide is found in the southern, southwestern and south-eastern portion where the soil resistivity and geology is very weak. The amount of annual average precipitation is maximum in this part of the district. With very high elevation and steep slope, certain geological configurations influenced by socio-economic aspects tend to become unstable causing landslide. Phuentsholing is a highly urbanized centre located in the south upland slope with dense road network built through modification of the general slope that accelerates landslide processes.
Landslides are caused by a number of factors in a given area, but not all of them are equally responsible. During the field investigation, it was found that the landslides in the study region are caused by both natural causes (geological structure, heavy rainfall and very steep slope) as well as by human interferences (such as slope cutting for the construction of road, deforestation for expansion of agricultural land) which makes the area vulnerable to landslide. Relative landslide density index (R), in this work, helped to analyse the association Table 3. Results of validation measures based on training and validation datasets.  www.nature.com/scientificreports/ between produced vulnerability classes of landslide models and percentage of inventory landslides. The highest R-index value can be found in very high vulnerability class of each model followed by high vulnerability class which is positive for validating models. The factor importance using RF (mean decrease Gini) and CSAE (average merit) represented the most important predictive LVCFs which are fittingly prominent in the southern part of the district (Fig. 12). Among the physical conditioning factors, both the RF and CSAE based evaluation has identified elevation, rainfall, geology, soil resistivity, soil clay and silt percentage etc. as the foremost persuasive factors for land sliding. The alike importance of the factors has reflected in several pieces of research such as elevation has been identified as a critical LVCF in various literature since most landslides occur in mountainous regions with a specific gradient 4,7,10 . Landslides are influenced by the structure, ordination, age, and exposure of the underlying surface 52 . The soil Table 4. Relative landslide density index of susceptibility classes of LSM models. www.nature.com/scientificreports/ resistivity parameter has a positive connection with landslide susceptibility, showing that reducing soil resistance upturns the chance of a landslide, especially at higher elevations. Increased soil clay and silt content at medium to high altitude and upslope areas remain more unstable due to lower integrity than rocky sections, making it more vulnerable to seepage erosion, liquefaction, and fluidization. However, the type of soil and the amount of vegetation cover plays an important influence in this 46 . Among the socio-economic factors, RF and CSAE method have recognised distance to road, population density, agricultural density, house frequency factor as crucial www.nature.com/scientificreports/ for making the area landslide vulnerable. A negative correlation exists for the distance to road LVF, indicating that the severity of possible landslide events increases in places as roads get closer, and vice versa. A substantial amount of work by Chan et al. 58 , Weigand et al. 59 has established the role of these factors in triggering landslides.

Conclusion
Landslides have been seen in recent decades as the most critical natural risk that poses serious threat to both life and property all over the world. Thus short-term and long-term solutions are considered necessary to confront these daunting challenges. The landslide vulnerability map has recently become an important means of delineating landslide-prone regions and management. With the aid of sophisticated methods, proper data, and integration of remote sensing and a geographical information system, this can be accomplished. DLNN, ANN, and convolution neural network (CNN) models were used and r all aspects of a landslide event were considered which are novel approaches to perform the landslide vulnerability mapping of the district. Therefore, along with geo-environmental data, potential socio-economic factors were also considered as LVCFs using advanced factor selection techniques. CNN model achieved highest accuracy in modelling the vulnerability in the study area. As per the finding of the models, lower part of the district is highly susceptible and it needs immediate measures for managing. For the landslide risk supervision in the present and also for the future, the LVM can suggest implementing different management strategies like afforestation, barrier construction, and proper land use planning.