*Article* **Use of Mamdani Fuzzy Algorithm for Multi-Hazard Susceptibility Assessment in a Developing Urban Settlement (Mamak, Ankara, Turkey)**

#### **Tugce Yanar 1, Sultan Kocaman 1,\* and Candan Gokceoglu <sup>2</sup>**


Received: 7 December 2019; Accepted: 17 February 2020; Published: 19 February 2020

**Abstract:** Urban areas may be affected by multiple hazards, and integrated hazard susceptibility maps are needed for suitable site selection and planning. Furthermore, geological–geotechnical parameters, construction costs, and the spatial distribution of existing infrastructure should be taken into account for this purpose. Up-to-date land-use and land-cover (LULC) maps, as well as natural hazard susceptibility maps, can be frequently obtained from high-resolution satellite sensors. In this study, an integrated hazard susceptibility assessment was performed for a developing urban settlement (Mamak District of Ankara City, Turkey) considering landslide and flood potential. The flood susceptibility map of Ankara City was produced in a previous study using modified analytical hierarchical process (M-AHP) approach. The landslide susceptibility map was produced using the logistic regression technique in this study. Sentinel-2 images were employed for generating LULC data with the random forest classification method. Topographical derivatives obtained from a high-resolution digital elevation model and lithological parameters were employed for the production of landslide susceptibility maps. For the integrated hazard susceptibility assessment, the Mamdani fuzzy algorithm was considered, and the results are discussed in the present study. The results demonstrate that multi-hazard susceptibility assessment maps for urban planning can be obtained by combining a set of expert-based and ensemble learning methods.

**Keywords:** multi-hazard; susceptibility mapping; developing urban settlements; landslide; flood; logistic regression; Mamdani fuzzy algorithm; M-AHP

#### **1. Introduction**

Improved disaster management is an important focus locally and globally to reduce the losses caused by natural disasters [1]. The harmful effects of natural hazards on human lives and economies increase with the inadequate land-use planning in developing countries. Actual infrastructure and land use should be taken into account in urban planning, which are often neglected [2]. Urban planning is a complex procedure that needs to consider existing infrastructure, human use, and natural hazards. The recent developments in geoinformation technologies for data collection and analysis, such as photogrammetry, remote sensing, three-dimensional (3D) geographical information systems (GIS), Web-GIS, volunteered geographical information (VGI), and advanced spatial analysis methods, can support this procedure and provide the essential tools to develop a combined approach.

Among the geology- and climate-related natural hazards (i.e., landslides, floods, earthquakes, droughts, wildfires, tornados, volcanic eruptions, and avalanches) [3], urban areas are mostly affected by landslides and floods. Landslide is one of the most common natural hazards with global spread, and it may damage buildings, infrastructure, and other facilities in urban areas [4]. Between 1950 and 2018, 23,041 landslides were observed in Turkey [5]. Although there is a vast amount of research on landslide susceptibility assessment in the literature (e.g., References [6–10]), most studies were on open lands and forests. Landslide susceptibility assessment in urban areas is difficult due to dense construction and buildings that modify and largely cover the topography. In addition, further research is needed due to the complexity of the problem, as well as incomplete and temporally inaccurate landslide inventories [11]. Production of regional landslide susceptibility maps can be difficult due to the requirement of actual data. Such maps are essential for urban planning and disaster mitigation efforts carried out by governments.

Although landslide conditioning parameters can be manifold, the main limitation for the number of factors employed to produce landslide susceptibility maps is the data availability. In general, geomorphological (e.g., topographical, hydrological, etc.), geological (e.g., lithology), and land-use and land-cover (LULC) parameters must be considered for this purpose [4]. A dense digital terrain model (DTM) can be used to characterize the geomorphology in detail in an urban area. In comparison to LULC, the geological and geomorphological characteristics change slowly. Since the LULC data are very important for landslide susceptibility assessment [12], actual data to demonstrate the LULC are required for obtaining high accuracy. In addition to the conditioning parameters, existing landslide inventories are employed in the susceptibility assessment models. On the other hand, landslide inventory extraction in settlement areas and heavy construction sites is also difficult due to covered or modified topography that considerably obstructs the visibility of landslides.

Floods also constitute one of the most commonly occurring and destructive natural hazards in the world, which also cause the highest number of fatalities [13,14]. According to the statistics [15], increasing numbers of flood events are being observed in Turkey. Depending on climate change and rapid urbanization, the number will continue to rise in the next decade in Turkey. Flood susceptibility maps also contain essential information for mitigation efforts, and they must be taken into account by decision-makers [16–18]. Flood susceptibility assessments for Ankara city were produced by Sozer et al. [19,20] using an expert-based decision support system called the modified analytical hierarchy process (M-AHP) [21]. A small part of the flood susceptibility map that covers the study area was employed here for multi-hazard susceptibility assessment. It should be noted that the flood susceptibility must be evaluated at a regional level, which includes the extent of the basin, since using the data of only a small area with limited altitude extent would lead to incomplete data and misinterpretation.

On the other hand, landslides and floods are often effective over the same regions since they occur in areas with similar geomorphological and climate conditions. Heavy precipitation triggers both floods and landslides, which occur one after the other. Therefore, areas that are prone to floods and landslides need to be assessed together to understand the combined effects of both.

The main aim of this study was to develop a spatial analysis methodology to predict the combined flood and landslide susceptibility level in a developing urban settlement, which can then be used as a basis for land-use planning. A part of Mamak District in Ankara, Turkey was selected as the study area for this purpose, because Mamak District is an unplanned settlement area of Ankara.

In the initial phase of the study, the landslide susceptibility map was produced by extracting and using the up-to-date LULC data from Sentinel-2 satellite images, high-resolution DTM, and lithology data obtained from existing geodatabases [4]. Since the flood susceptibility maps must be produced at regional level, i.e., basins, a part of the map produced by Sozer at al. [19,20] was reclassified and used for the purposes of this study. To assess the multi-hazard susceptibility level (MHSL), a Mamdani fuzzy algorithm was developed, and the results are presented and discussed in the later sections. This is the first application of the Mamdani fuzzy algorithm to combine two susceptibility maps in the international literature.

#### **2. Background on Multi-Hazard Assessment**

Many natural hazards types affect urban settlements. To mitigate the effects of future natural hazards, all possible risks in an area should be assessed. Multi-hazard assessment models can be generated by integrating multiple susceptibility assessment maps belonging to different types of natural hazards for a specific area. Moreover, using multi-layer information is a more reliable and effective approach to disaster prevention [22]. Different types of susceptibility maps created by using various parameters and factors can be combined with different methods, such as AHP, which is based on expert opinion. Using AHP, a suitability map can be determined by the weight coefficients and uncertainties for each hazard [23–25]. Furlan et al. [26] provided a gradual analysis of all components contributing to the risk at a particular site. This method involved the assessments of hazard, exposure, vulnerability, and risk. In this analysis, since the input dataset was large and heterogeneous, the multi-criteria decision analysis (MCDA) method was used to evaluate the parameters [27]. In all these studies, scores and weights provided by the experts had great importance, and it was stated by the researchers that they affected the accuracy of the results. Chen et al. [28] considered debris flows and river and flash flooding to be common in one area. These hazard types were examined in four scenarios (major, moderate, minor, and frequent events). In this approach, the losses caused by each hazard were firstly calculated individually. Afterward, the spatial probability of the element at risk, the physical vulnerability, and the quantification of the exposed elements at risk were multiplied. The effects of hazard types were compared based on the results.

It is also important to include social and economic dimensions in multi-hazard assessment methods [29]. China's disaster risk index was calculated for 31 provinces by using four types of factors: exposure (population exposed to earthquakes, floods, droughts, low temperatures/snow, and gale/hail), susceptibility (based on public infrastructure, income health, and economic status), coping capacity (based on governance, medical care, and material security), and adaptive capacity (related to future natural events) [30]. In addition to these factors, the exposure parameter was analyzed as hazard exposure and hazard loss. Hazard exposure referred to the presence of assets and values that may be adversely affected in hazardous areas. The hazard loss was defined by the extent of physical damage, monetary loss, human loss, and economic deterioration. At the same time, the scope of the other parameters was extended [31]. Using a multiple linear regression method, influencing factors of community resilience were calculated. It was also stated that this model can be developed for employing a different weighting scheme by using expert knowledge and the entropy. Moreover, for holistic multi-hazard assessment methodologies, it was proposed that anthropogenic processes should be used in relation to natural disasters, and it was mentioned that the environment in which natural hazards are experienced is shaped by human activities. The main idea was to include this relationship in the multi-hazard assessment process [32,33]. Additionally, Gallina et al. [33] and Basheer Ahammed and Pandey [22] pointed out that the climate change perspective is a forgotten piece that should be evaluated in a multi-risk assessment of natural and anthropogenic systems.

Barrantes [34] proposed a natural multi-hazard assessment model that can be used when working with limited data. This model proposed an algorithm that takes the spatial overlap of the values of each risk and the potential interaction between different natural risks and the temporal frequencies into account. In this algorithm, all risk combinations were evaluated and potential interactions between natural hazards matrix were formed. Potential multi-hazard risks arose with the set of intersections of all combinations. Liu et al. [35] used the MmhRisk-HI (Model for multi-hazard Risk assessment with a consideration of Hazard Interaction) method for multi-hazard risk assessment. The method had two main components. The first component analyzed the relationship between the hazardous environment and the hazards, showing the probability of multiple hazard occurrence. In this component, the probabilities were calculated with functions used according to the relationship levels (independent, mutex, parallel, series) between the natural hazards. The second component calculated the possible damages and loss rate by employing a Bayesian network. In Bernal et al. [36], a fully probabilistic multi-hazard risk model was assessed for hazard, exposure, vulnerability, and

loss using specific stochastic event sets. The average and maximum yearly losses were calculated with this model. Quantitative approaches were also used in some regional risk assessments [37]. The study of quantitative multi-hazard risk evaluation using elements at risk, their exposure, and their vulnerability can be examples for regional multi-hazard assessment. In another quantification approach, the interrelation between hazards was examined in two ways: by cascading hazards and by compound hazards. These relationships were evaluated with three hazard interrelation modeling approaches (stochastic, empirical, mechanistic) [38]. When the results were evaluated using different model parameters in each approach, it was seen that the extreme copulas method in the stochastic approach, the linear regression method in the empirical approach, and hydrodynamic models in the mechanistic approach were most prevalent. As mentioned, the use of mechanistic and stochastic methods in multi-parameter (more than two) hazards imposes certain restrictions, such as uncertainties caused by statistical assumptions (e.g., distribution, dependence model selection) or the effect of the data quality used for validating the mechanistic models [38].

Assessing multiple hazards in urban areas and predicting future risks can help decision-makers to prioritize actions and manage the risks [32,39–41]. Although the effects of natural risks on the urban area were examined individually and the results were compared visually, a quantitative approach of the co-evaluation process of hazards was lacking in Chang et al. [39] and Jacobs et al. [42]. A combined and quantitative assessment of hazards provides more accurate results than individual assessments and visual comparison. However, the choice of parameters and their weights, as well as the quality of data, are also important. The choice of using a qualitative, semi-quantitative, or quantitative approach may vary depending on the target [43]. In Omidvar and Karimi [44], a method developed using the theory of probability and Boolean logic was used. This study was conducted for multi-hazard reliability measurements according to available urban data. Multi-hazard reliability emerged with operations on different combinations of hazard risks.

Machine learning methods are also available for multi-hazard risk assessment. In Reference [45], risk analysis for each natural hazard type was carried out by modifying the weights of each parameter according to the hazard type using the random forest (RF) method. Afterward, a multi-hazard map was produced by combining the results. The model accuracy of 96.70% supported the usefulness of machine learning in multi-hazard risk assessments. Mirzaei et al. [46] used Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) as a multi-criteria decision-making model. It worked with similarity indexes of hazard maps and showed multi-hazard assessment. In Sheikh et al. [47], TOPSIS–Mahalanobis distance, TOPSIS, and simple additive weight (SAW) methods were combined. Although the TOPSIS method was criticized for using only a geometric distance, it was mentioned that it gives more clear results for natural hazards.

There are also other studies in the literature in which different multi-criteria decision-making methods were used in risk assessment. In Pourghasemi et al. [48], risk maps were created by using the stepwise weight assessment ratio analysis (SWARA) method which is an expert-oriented method, as well as the adaptive neuro-fuzzy inference system (ANFIS) method which involved an FIS (fuzzy inference system) and gray wolf optimization (GWO) to find the optimal solution. In multi-hazard analysis, the occurrence rate of hazard risk combinations was determined using weighted overlay analysis of risk maps by Mukhopadhyay et al. [27]. Moreover, with the use of fuzzy modeling, direct standardization of multiple indicators, aggregation, and deriving the impact are possible. By using the gamma fuzzy overlay model, the relationship between multiple input criteria was explored [49]. Kappes et al. [50] emphasized the importance of multi-hazard assessment and compiled difficulties when analyzing multi-hazards. Consequently, the present study introduces the production of a multi-hazard map for a settlement area by employing the Mamdani fuzzy inference algorithm.

#### **3. Materials and Methods**

Mamak District is a rapidly developing area located in the eastern part of Ankara, Turkey. Similar to many other large cities in Turkey, Ankara is affected by urban sprawl, and Mamak is one of the development centers for this sprawl. Approximately 640 thousand people live in the area. Since Mamak is on the route to Eastern Turkey, there is substantial transportation infrastructure. Furthermore, Bayindir Dam, which is currently used as a recreational area, is also located here. A part of Mamak District which is prone to both landslides and flooding was selected as the study area since it features continuous urban expansion and infrastructure (Figure 1). The area is ca. 30 km2 and the minimum and maximum altitudes are 924 m and 1284 m, respectively.

**Figure 1.** The location of the study area and an overview of the Sentinel-2 red–green–blue (RGB) image used in the study (upper left coordinates: 32◦56 51.372" E, 39◦56 27.108" N; lower right coordinates: 33◦0 57.578" E, 39◦53 41.689" N).

The landslide susceptibility map of the study area was produced by applying logistic regression (LR) to LULC data (sourced from Sentinel-2 imagery), the geomorphological features (sourced from DTM), and the lithology data [4]. Sentinel-2 images are distributed freely by the European Space Agency (ESA). It was found that actual land-use data can be produced from Sentinel-2 images, which are geometrically corrected (i.e., orthorectified, L2A) and obtained regularly over a large geographical extent by ESA [51]; they also provide multi-band (13 bands in total) data at spatial resolutions of 10 m, 20 m, and 60 m. The satellite constellation also has high transmission frequency [52] and, thus, is widely used in natural hazard assessments [53]. Sentinel-2 imagery can be easily employed by non-experts via the Sentinel Application Platform (SNAP) Tool from ESA. The LR is a rather simple method and can take non-numerical parameters such as different LULC types into account. The output probabilities were classified as low, moderate, and high. The flood susceptibility map of Ankara was produced in a previous study by Sozer et al. [19] for Ankara city and reclassified here into the same three classes (i.e., low, moderate, and high) to be used in the multi-hazard susceptibility assessment model. The study workflow is depicted in Figure 2. More details on the input data and methodology are provided in the sub-sections below.

**Figure 2.** Overall workflow of the study.

#### *3.1. Input Datasets for Landslide Susceptibility Map Production*

A DTM with 5-m resolution was obtained from the General Directorate of Mapping (GDM), Turkey. The geomorphological parameters such as altitude, slope, general curvature, plan and profile curvatures, topographic wetness index (TWI), stream power index (SPI), distance to channel networks, and ridgelines were derived from the DTM. The Sentinel-2 satellite imagery from 23 March 2019 were employed for classifying the LULC using the RF method. The lithology data were digitized into vector form using the WebGIS portal data of the General Directorate of Mineral Research and Exploration (GDMRE/MTA), Turkey [54], and converted into raster data with 5-m grid spacing. The LULC map was also resampled to 5-m grid data to perform the LR technique. The data sources and the spatial resolutions are summarized in Table 1.

**Table 1.** Properties of input parameters. LULC—land use and land cover; DTM—digital terrain model; GDMRE/MTA—General Directorate of Mineral Research and Exploration.


For training the LR method, boundary polygons of eight landslides covering ca. 2000 grid points were manually delineated by the expert using the DTM and the satellite images. The altitude range map and the landslide inventory represented by the red polygons are depicted in Figure 3. In order to utilize in the LR model estimation process, the vector landslide inventory map was rasterized with 5-m grid spacing. The landslides detected in the study area occurred in schists and volcanic units. These units are highly susceptible to weathering and landslides. The main characteristics of the landslides were circular, and the depth of failure surfaces was controlled by the thickness of weathered zones. In addition, the damage on the buildings in the study area was observed, and it is well known that the main cause of damage is the landslides in the area. However, it is impossible to draw the borders of the landslides due to urbanization on the slopes.

**Figure 3.** The elevation map and the manually delineated landslides (red polygons).

#### *3.2. Geomorphological Characteristics of the Study Area*

In order to understand the topography, some scalars such as primary and secondary derivatives of a DTM can be used. Primary features (e.g., slope, aspect, curvature) are computed from elevations, whereas secondary features (e.g. SPI, TWI) are obtained from the second derivatives of elevations [55]. These derivatives can be computed by using spatial analysis tools and software. Here, SAGA GIS from the SAGA User Group Association, Germany [56] and ArcGIS from ESRI Inc., Redlands, CA, USA [57] were used for this purpose. A statistical summary of the elevation data and the derivatives are given in Table 2. Similar statistics were derived for the landslide positive samples (~2000 grid points) and are provided in Table 3.


**Table 2.** Statistics of topographic attributes. SPI—stream power index; TWI—topographic wetness index.


**Table 3.** Statistics of topographic attributes in the landslides.

Altitude shows the elevation measures of an area [58]. The slope gradient represents the variation in elevations [58] (Figure 4). Here, the slopes were employed to relate the topographical changes to landslide formation. The aspect was computed as the angle from the north to depict the direction of the slope. The aspect parameter was used to understand which slopes (north, south, etc.) would affect the landslides more [58]. The aspect values are also depicted in Figure 4.

**Figure 4.** The slope gradient (**left**) and aspect (**right**) maps of the study area.

The curvature calculated from the DTM shows the changes in slope and aspect (Figure 5). The curvature parameter can be classified as plan and profile, and it needs to be analyzed separately. The planimetric component is based on the rate of slope variation along the contour lines, and the profile component is computed along the slope to determine the rate of slope gradient change [58]. Negative, positive, and zero curvatures reflect concave, convex, and flat surfaces, respectively [59]. The plan and profile curvatures of the study area are given in Figure 6.

**Figure 5.** General curvature map of the study area.

**Figure 6.** Plan (**left**) and profile (**right**) curvature maps of the study area.

The SPI is an indicator for erosive power of flowing water [60]. SPI is effective on landslides and denotes potential erosion energy [61]. The TWI shows the positions and size of the water-saturated regions [55] (Figure 7). The distances to channels were used for understanding the effect of the drainage network on landslides. The vertical distances to the channels were computed from the elevation data and, thus, the network was formed [62]. The landslide probability and the distances to the ridges were negatively correlated [59]. The ridges were computed using the DTM, and the distance values to both the ridges and the channel networks are provided in Figure 8.

**Figure 7.** SPI (stream power index, on the left) and TWI (topographic wetness index, on the right) maps of the study area.

**Figure 8.** Distances to channels (**left**) and to ridgelines (**right**).

#### *3.3. Land-Use and Land-Cover Extraction from Sentinel-2 Imagery*

The use of up-to-date LULC data is essential for natural hazard assessments and disaster mitigation efforts. In this study, the LULC data were extracted from Sentinel-2 optical satellite imagery with 10-m spatial resolution. Seven LULC classes as shown in Figure 9 were extracted from red–green–blue (RGB) bands using SNAP software. The RF method was applied for the classification by collecting training data on the images. Using the RF method, classification is made with an ensemble of decision trees created by using training samples and variables [63]. In decision trees, the most rated pixels from all trees in the forest are classified. Because of the higher accuracy compared to other machine learning methods, it is widely used for image classification [64–67]. As seen in the study of Lim et al. [68] on Sentinel-2 images, the RF method can work with high accuracy in image classification. The RF classifier has two main parameters: the number of trees (T) and the number of variables (M) [69].

**Figure 9.** Land-use and land-cover map of the study area.

In this study, the classification process was completed by creating 10 trees using four bands, i.e., red, green, blue, and the gray-level co-occurrence matrix angular second moment (GLCM-ASM) parameter using 2077 training samples. Since it was difficult to separate the industrial units from

the roads only with RGB information, the GLCM of the images was also computed and added to the classification as proposed by Stumpf and Kerle [70], and the GLCM-ASM parameter was found to be particularly useful in the present study. The distribution of the training samples to the classes is shown in Table 4. When the accuracy of classification was calculated by the cross-validation method, the correct prediction percentage was 93.73%, the classification precision was 92.01%, and the kappa value was 97.13% (Table 4). As can be seen in Table 4, the classification accuracies of all bands were improved by including the GLCM-ASM parameter, except for the industrial units, which remained high for both versions. In addition, it was visually verified that this parameter was especially useful to separate the road and industrial unit classes.


**Table 4.** Classification accuracy of the random forest method. ASM—angular second moment.

#### *3.4. Lithological Characteristics of the Study Area*

In addition to the LULC and elevation data, the lithology type is extremely important for natural hazard assessments [71]. Lithology type and the structural differences generally affect the robustness and permeability of rocks and soils [72]. The lithology map was obtained from the geosciences portal (Yer Bilimleri Portali) of GDMRE/MTA. The lithology map is shown in Figure 10, and detailed descriptions are provided in Table 5. This vector map was preprocessed for conversion to raster data with 5-m grid spacing.

**Table 5.** Age and general descriptions of the lithologies in the study area [54].


**Figure 10.** Lithology map of the study area [54].

#### *3.5. Landslide Susceptibility Map Production with Logistic Regression Method*

For the generation of landslide susceptibility maps, various mathematical and machine learning methodologies can be applied. Many studies on the assessment of landslide susceptibility using logistic regression were published in the literature (e.g., References [73–78]). In this study, multivariate LR was employed to derive the landslide susceptibility distribution of the area. LR is a statistical model, and it was used here to predict the potential landslide areas since it is fast and accurate for landslide susceptibility assessment purposes [9,59,79]. The LR is a supervised method and uses dependent (i.e., landslide conditioning factors) and independent (i.e., actual landslide inventory) variables. The dependent variable is a binary value which depicts the occurrence/non-occurrence of the event [62]. Independent variables were the 11 conditioning factors used as input data layers here (e.g., elevation data, slope, aspect, LULC, etc.). In the model estimation stage, the relationship between the variables was analyzed using the landslide positive samples (i.e., inventory data) and a number of randomly chosen non-landslide samples. Equations (1) and (2) were used for computing the logistic regression method.

$$\mathbf{Y}\_{\mathrm{i}} = \beta\_0 + \beta\_1 \mathbf{X}\_{\mathrm{i}} \tag{1}$$

$$P\_i = (\mathbf{Y} = 1 | \mathbf{X}\_i) = 1/(1 + \mathbf{e}^\circ(-\mathbf{Y}\_i)) \tag{2}$$

where Yi represents the dependent variables, xi represents the independent variables, β<sup>0</sup> is a constant, β<sup>i</sup> represents the *i*-th regression coefficient, and P is the probability of the existence of landslides [59]. In Vorpahl et al.'s [80] and Park et al.'s [81] studies, one of the most accurate results among the landslide susceptibility maps created with similar parameters used in this study was achieved using the LR method compared to other methods. After calculating the LR model parameters, the landslide susceptibility map was produced for the whole area. The ratio of the landslide positive and negative (non-landslide) samples was 1:2.

#### *3.6. Flood Susceptibility Map of the Study Area*

The flood susceptibility map of Ankara City was obtained from a previous study [19,20] (Figure 11) with the M-AHP method [21] using flow accumulation, slope, topographic altitude, distance to permanent river, distance to dry drainage, land cover, topographic wetness index, and lithology parameters. Each parameter was weighted with the M-AHP method according to expert opinion, and the highest score of each parameter was defined a priori. The output flood susceptibility map was reclassified by dividing the resulting probability values into three classes with equal intervals, and they were clipped for the study area (Figure 12). The original flood susceptibility map was divided into five equal classes by Sozer et al. [19]. However, in accordance with the purpose of the study, three

classes of the flood susceptibility map were used in the present study. The histograms prepared for the five classes and three classes are given in Figure 13. The output susceptibility classes were categorized as low, moderate, and high, and they were used for the multi-hazard susceptibility assessment with the Mamdani fuzzy method.

**Figure 11.** Flood susceptibility map produced by Sozer et al. [19] (the rectangular area to the east is the selected study area).

**Figure 12.** Flood susceptibility map of the study area (modified after Reference [19]).

*ISPRS Int. J. Geo-Inf.* **2020**, *9*, 114

**Figure 13.** Histograms of flood susceptibility classes for five classes (**left**) and three classes (**right**).

#### *3.7. Multi-Hazard Susceptibility Assesment with Mamdani Fuzzy Method*

The multi-hazard susceptibility assessment map was derived with a combined assessment of flood and landslide susceptibility results using the Mamdani Fuzzy Method, which was first developed by Mamdani and Assilian [82]. This method is able to reduce uncertainties while solving complex problems using "if–then" rules. The stages of a Mamdani FIS are fuzzification, rule evaluation, aggregation, and defuzzification [82]. Fuzzy inference systems (FIS) produce a crisp output for supplied crisp inputs by using fuzzy set theory [83]. The general structure of a Mamdani FIS can be found in several books and publications (i.e., References [84–87]). Osna et al. [87] developed an integrated tool for construction of a Mamdani FIS for Netcad Software, Netcad, Ankara, Turkey. With the Mamdani fuzzy logic [88] operator in Netcad, landslide and flood susceptibilities could be evaluated together in the study area. The Mamdani fuzzy algorithm was previously used for landslide susceptibility mapping [86,87,89], but the present study is the first attempt at the combination of two susceptibility maps to obtain a multi-hazard susceptibility map.

The inputs of the Mamdani fuzzy model constructed in the study were landslide susceptibility and flood susceptibility maps, while the output was the multi-hazard susceptibility level (MHSL). Traditionally, a fuzzy model is built using expert knowledge in the form of linguistic rules. Three membership functions, i.e., low, moderate, and high, were defined for each input and output in the Mamdani FIS implemented here. The membership functions are shown in Figure 14, which also constitute the fuzzification stage of the system. In Figure 14, the vertical axes of the graphs denote the membership degree, and the horizontal axes represent the susceptibility levels, which range from 0–1 for landslide and 8–66 for flood. In the literature, many methods, such as intuition, rank ordering, angular fuzzy sets, genetic algorithms, inductive reasoning, soft partitioning, etc., exist for membership value assignment (e.g., References [90–92]). Although it is possible to select membership functions in a site-specific or target-oriented manner, or non-linearly, a generic approach was preferred here to prove the usability of the approach. In this study, the constructed fuzzy model employed two inputs and one output using three membership functions, and the fuzzification of crisp numbers and degree of membership of each crisp input were calculated at this stage. One of the fundamental features of a Mamdani FIS constitutes the linguistic if–then rules, namely, rule evaluation. In the present study, the if–then rules were generated by the expert (last author), which prevented an exhaustive data analysis process. The total number of linguistic rules generated by the expert was nine (Table 6). The final fuzzy output of the model was produced by aggregation of all local results from fuzzy rules triggered in the rule evaluation phase [87]. In the Mamdani FIS constructed here (Figure 15), the maximum operator was considered for aggregation, as suggested by Reference [87]. Finally, the center of gravity was used for defuzzification. Employing the Mamdani FIS constructed, the landslide and flood susceptibility maps were used as inputs, and the MHSL map was produced as presented in the section below.

**Figure 14.** The membership functions of each input. The vertical axes in both graphs represent the degree of membership, while the horizontal axes reflect the susceptibility level range for landslide (**left**) and flood (**right**).

**Table 6.** *If–then* fuzzy rules used for the multi-hazard susceptibility level (MHSL) assessment in the study area.


**Figure 15.** The general structure of the Mamdani fuzzy inference system (FIS) constructed.

#### **4. Results and Discussion**

#### *4.1. The Landslide Susceptibility and MHSL Maps*

The output landslide susceptibility map is shown in Figure 16. Although the landslide occurrence probability values ranged from 1%–99% (i.e., from 0–1, as shown in the landslide susceptibility membership graph in Figure 14), these values were reclassified into three categories (low, moderate, and high) by using equal interval classification for easy interpretation. These values were also used in the fuzzy assessment model in the next step. The map demonstrates the existing landslide hazard potential, especially in the western parts of the area. The field observations of the expert support the findings of the results obtained from the study.

**Figure 16.** Landslide susceptibility map of the study area.

The accuracy of the output map was evaluated to understand the quality of the results. The factors that affect the accuracy are the data quality, the applied methods, the number of input parameters used in the process, and the approach for map production [72]. The ROC (receiver operating characteristic) curve, which is a measure of the capability of the current model in classification [93], was used to evaluate the accuracy (Figure 17). The figure shows that the areas with and without landslide susceptibility were classified with 96% accuracy.

In the present study, a plausible and practical methodology to combine different susceptibility maps was proposed. The Mamdani fuzzy algorithm was used and the MHSL map was obtained (Figure 18). The results show that some of the slopes and valleys have high multi-hazard potential. To minimize the losses caused by landslide and flood, the high multi-hazard susceptibility zones shown in Figure 18 must be investigated carefully, and necessary engineering measures must be provided at the construction stage.

**Figure 17.** Receiver operating characteristic (ROC) curves of landslide susceptibility map.

**Figure 18.** Multi-hazard susceptibility level map of the study area.

#### *4.2. Discussion*

The study area is subject to urban transformation projects due to unplanned settlements. The Mamak Urban Transformation Project was implemented in a part of the study area, which covers ca. 7.4 km2. The project location is shown in Figure 19. It was divided into 11 stages and carried out by TOKI (Toplu Konut Idaresi Baskanligi), which is a state organization carrying out large-scale construction works for new houses, in the Ankara Metropolitan Municipality and Mamak Municipality. The main aim of the project was to transform the slums, i.e., unplanned settlement areas with insufficient facilities and infrastructure, and modernize these areas [94]. While the initial number of slums was 13,662, the number of slums destroyed as of 2019 was 8389 throughout the project. In total, 30,000 dwellings are planned to be constructed in the next phase of the project. Considering the natural hazard potential and related risks within the study area, the MHSL map could be used to analyze the vulnerability of future urban development and transformation plans.

**Figure 19.** A part of urban transformation project within the study area [91].

The DTM textured with the Sentinel-2 RGB image (Figure 20), the landslide susceptibility map (Figure 21), the flood susceptibility map (Figure 22), and the MHSL map (Figure 23) of the study area were visualized in 3D with the QT Modeler software from Applied Imagery, Chevy Chase, MD, USA [95] for interpretation of the results.

**Figure 20.** The DTM (digital terrain model) of the study area textured with the Sentinel-2 image.

**Figure 21.** The DTM of the study area textured with the landslide susceptibility map (output of logistic regression).

**Figure 22.** The DTM of the study area textured with the flood susceptibility map (modified into three classes after Sozer et al. [19]).

**Figure 23.** The DTM of the study area textured with the MHSL (multi-hazard susceptibility level) map. The circles denote important focal areas for city planning in northwest and south Mamak as mentioned above.

Figure 23 demonstrates the locations of the existing and the planned urban transformation sites. The existing project area (purple circle in Figure 23) includes areas with multi-hazard risk. When the risks are evaluated individually, there is low landslide risk in this area, but the excess of flood risk is also remarkable. In this respect, it can be considered to review the resilience of the project against natural hazard risks.

There are slum areas (blue circles in Figure 23) which represent potential urban transformation sites. Land-use decisions in these areas should be prepared elaborately. It is more appropriate to evaluate these areas at high susceptibility and utilize them as urban green areas by establishing agreements with the property owners. There are also newly constructed buildings (dark-green circle in Figure 23) in the west of the Mamak region. As can be seen, there are multiple areas with multi-hazard risk. The status of new buildings in these areas should be examined. The construction quality and the landslide resistance of these structures can be used as criteria to measure the accuracy of the previous partial urban transformation in this area.

With regard to the methodological approach employed here, the accuracy of the LR method was found sufficient for the purposes of the study. Although the number of the training samples in manually delineated landslide areas was low, using the 1:2 ratio for landslide/non-landslide samples worked efficiently. The produced MHSL map was representative for the purposes of the study, and it can be used as base map for urban planning and transformation purposes. The ensemble (i.e., LR) and expert-based (i.e., M-AHP) methods employed for the production of the individual hazard susceptibility maps were useful, and the Mamdani Fuzzy algorithm was able to handle the complexity of the problem.

#### **5. Conclusions**

In this study, a fuzzy model for integrated multi-hazard susceptibility assessment was developed. In addition, the usability of Sentinel-2 images in obtaining up-to-date LULC data in urban development areas for the production of landslide susceptibility maps was evaluated. The RF classification method was employed for producing LULC classes from Sentinel-2 RGB images, and the GLCM-ASM parameter was added to the classification to distinguish the industrial areas from roads. A high-resolution DTM and the lithology data were integrated into the landslide susceptibility map, and the LR method was applied for this purpose. An MHSL map was produced using the Mamdani fuzzy algorithm. The produced MHSL map can be used as essential data for urban development and transformation plans. Further analysis and planning can be carried out for this purpose.

The main difficulty encountered during the study was the preparation of a fully completed landslide inventory map due to urbanization. However, a combined methodology to obtain the MHSL map was described and applied successfully. The high multi-hazard susceptibility zones must be investigated carefully before construction; alternatively, if possible, these zones should be avoided with regard to construction purposes to minimize losses sourced from natural hazards. These maps are highly useful for planning stages, and, if these maps are considered during the planning stage, serious benefits can be obtained.

The methodology introduced in the present study for producing multi-hazard susceptibility represents a first in international literature. The use of an expert-based fuzzy inference system for the combination of two susceptibility maps portraying different natural hazards yielded very promising results. The study showed that the Mamdani type fuzzy inference system is a suitable approach for producing multi-hazard susceptibility mapping. the use of this approach for the combination of several multi-hazard susceptibility maps may provide new possibilities for suitable site selection efforts. However, there is no methodology available for the accurate assessment of a multi-hazard susceptibility map. For this reason, the study area selected was relatively small and from a well-known area. Consequently, the final output map was assessed with field observations.

As a final concluding remark, the fuzzy algorithm proposed for combining different natural hazards is a flexible and transparent modeling approach; hence, the model can be tuned or re-constructed easily when new information is obtained. As a future recommendation, some attempts at performance assessments of the multi-hazard susceptibility maps should be carried out, and some numerical indices for accuracy and performance assessments should be developed.

**Author Contributions:** Conceptualization and methodology, Candan Gokceoglu, Sultan Kocaman, and Tugce Yanar; investigation, validation, formal analysis, writing—original draft preparation, and visualization, Sultan Kocaman and Tugce Yanar; supervision, project administration, funding acquisition, and resources, Sultan Kocaman; data curation, Tugce Yanar; writing—review and editing, Candan Gokceoglu All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is part of the M.Sc. thesis of the first author, and it was funded by the Hacettepe University Scientific Research Projects Coordination Unit (Grant FYL-2019-18017).

**Acknowledgments:** The authors sincerely thank the General Directorate of Mapping and the General Directorate of Mineral Research and Explorations for the provision of data, as well as Recep Can, Burhan Sozer, Orhan Firat, and Hakan A. Nefeslioglu from Hacettepe University for their continuous support. The authors would also like to thank the anonymous reviewers for their constructive comments and suggestions, which greatly helped to improve the manuscript.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
