The land subsidence susceptibility maps of all four algorithms used in this study could be ordered according to the accuracy and time taken to build the model due to the similarities in the conditioning factors, training and testing data, and study area [28]. We used ROC curve analysis to assess the accuracy of the maps. The AUC data showed that the AdaBoost algorithm had the highest susceptibility class predictive accuracy of 81.1%, which was 1.1% higher than the multilayer perceptron algorithm, 1.7% higher than the logistic regression algorithm, and 2% higher than the LogitBoost algorithm. Since the accuracy of the algorithms used was so closely related, we analyzed the time consumption of data preparation needed to build the model. The AdaBoost algorithm needed 218.63 s to build the model, the LogitBoost algorithm needed 174.94 s to build the model, logistic regression needed 16667.96 s or 4.63 h to build the model, and multilayer perceptron needed 47,538.65 s or 13.20 h to build the model, thus identifying the LogitBoost algorithm as the fastest.

However, the LogitBoost algorithm had the lowest prediction accuracy (79.10%). The Adaboost algorithm needed 218.63 s or 3.64 min to build the model, along with having the highest prediction accuracy of 81.10% in the ROC curve analysis. Logistic regression needed 4.63 h to build the model, along with a prediction accuracy of about 79.40%, and the multilayer perceptron needed 13.20 h, along with a prediction accuracy of 80%. The maps generated using LogitBoost, logistic regression, and multilayer perceptron (Figure 7b–d, respectively) showed similar weights and indicated that the settlement area was at a high risk of land subsidence. The land-use map generated in this study indicated a relatively strong correlation between land subsidence and urban development, likely due to excessive groundwater extraction in urban areas [3,6,7]. Our map could be used by local environmental authorities and those in charge of urban development to identify areas with high subsidence risk. The method of combining SAR and GIS spatial data to generate land subsidence susceptibility maps employed in this study could be applied to other regions.

#### **6. Conclusions**

We generated a land subsidence inventory map using a time-series InSAR method based on the StaMPS algorithm from Sentinel-1 SAR datasets in ascending and descending tracks. The comparison of both tracks was conducted, finding a coefficient of correlation between the two tracks of 0.9548. The inventory map could be used as a training and testing dataset to create a land subsidence susceptibility map for Jakarta using meta-ensemble (AdaBoost and LogitBoost) and functional (logistic regression and multilayer perceptron) machine learning algorithms. We created a land subsidence inventory map through time-series InSAR analysis of Sentinel-1 SAR datasets using the StaMPS algorithm. The land subsidence susceptibility map produced by the AdaBoost machine learning algorithm had higher accuracy (AUC = 0.811) and only need 3.64 min to build the model compared with the maps created using the other algorithms (multilayer perceptron, AUC = 0.800; logistic regression, AUC = 0.794; LogitBoost, AUC = 0.791). LogitBoost was the fastest algorithm in building the model but had the lowest predictive accuracy. Logistic regression and multilayer perceptron needed 4.63 and 13.20 h, respectively, to build the model. All of the maps showed acceptable accuracy, as the AUC values were all higher than 0.5; thus, they can all be used for analyzing land subsidence susceptibility in Jakarta. Our approach based on time-series InSAR analysis, machine learning, and GIS spatial data yielded reasonable predictions of areas with high risk of land subsidence. Further research could use alternative algorithms and conditioning factors to generate land subsidence susceptibility maps in other regions.

**Author Contributions:** Conceptualization, C.-W.L., and W.L.H.; methodology, C.-W.L., A.R.A., and W.L.H.; software, A.R.A. and W.L.H.; validation, C.-W.L. and W.L.H.; formal analysis, W.L.H.; investigation, C.-W.L.; resources, C.-W.L.; data curation, C.-W.L. and W.L.H.; writing—original draft preparation, W.L.H.; writing—review and editing, C.-W.L., A.R.A., and W.L.H.; visualization, C.-W.L., A.R.A., and W.L.H.; supervision, C.-W.L.; project administration, C.-W.L.; funding acquisition, C.-W.L. All authors read and agreed to the published version of the manuscript.

**Funding:** This research was supported by a grant from the National Research Foundation of Korea provided by the government of Korea (No. 2019R1A2C1085686).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 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/).

*Article*

## **Integration of InSAR Time-Series Data and GIS to Assess Land Subsidence along Subway Lines in the Seoul Metropolitan Area, South Korea**

#### **Muhammad Fulki Fadhillah , Arief Rizqiyanto Achmad and Chang-Wook Lee \***

Division of Science Education, Kangwon National University, Gangwon-do, Chuncheon-si 24341, Korea<EMAIL_ADDRESS>(M.F.F<EMAIL_ADDRESS>(A.R.A.)

**\*** Correspondence<EMAIL_ADDRESS>
Received: 9 September 2020; Accepted: 22 October 2020; Published: 25 October 2020

**Abstract:** The aims of this research were to map and analyze the risk of land subsidence in the Seoul Metropolitan Area, South Korea using satellite interferometric synthetic aperture radar (InSAR) time-series data, and three ensemble machine-learning models, Bagging, LogitBoost, and Multiclass Classifier. Of the types of infrastructure present in the Seoul Metropolitan Area, subway lines may be vulnerable to land subsidence. In this study, we analyzed Persistent Scatterer InSAR time-series data using the Stanford Method for Persistent Scatterers (StaMPS) algorithm to generate a deformation time-series map. Subsidence occurred at four locations, with a deformation rate that ranged from 6–12 mm/year. Subsidence inventory maps were prepared using deformation time-series data from Sentinel-1. Additionally, 10 potential subsidence-related factors were selected and subjected to Geographic Information System analysis. The relationship between each factor and subsidence occurrence was analyzed by using the frequency ratio. Land subsidence susceptibility maps were generated using Bagging, Multiclass Classifier, and LogitBoost models, and map validation was carried out using the area under the curve (AUC) method. Of the three models, Bagging produced the largest AUC (0.883), with LogitBoost and Multiclass Classifier producing AUCs of 0.871 and 0.856, respectively.

**Keywords:** Seoul; synthetic aperture radar; land subsidence; GIS; machine learning; time-series

#### **1. Introduction**

Land subsidence is a threat faced by big cities with extensive development that can negatively impact the environment, social systems, and the economy [1]. Subsidence occurs due to geological causes or anthropogenic processes such as massive urban development, infrastructure development [2,3], tunneling [4–6], water extraction [7–9], and earthquakes [10]. Subsidence has been observed in several metropolitan cities, including Mexico City [11], Shanghai [12], and Jakarta [13–15]. The Seoul Metropolitan Area is the center of governance, commerce, and culture in South Korea. It has been extensively developed and is the most densely populated city in Asia [16]. Industrial development and economic growth have led to city developments such as the expansion of subway lines and the construction of many structures and buildings [17]. By examining the potential effects of land subsidence, monitoring of land deformation could be the first step of a mitigation process. Seoul, which has a high population density and extensive developments, is extremely vulnerable to land subsidence. Given the severe negative impacts of subsidence, it necessary to elucidate the factors that cause land subsidence from integrated observations, to assess damage risks and prevent damage to roads, bridge, railways, and other infrastructure.

Monitoring land deformation is essential for reducing losses due to subsidence and developing a sound mitigation plan. With advancements in knowledge and technology, monitoring techniques have improved greatly. Previously, measurements related to land subsidence could be taken using traditional geodetic [15] and leveling methods, which provide precise measurements, but these methods are inefficient and costly compared with satellite-based methods, which cover broader areas and are more cost-efficient. These days, more studies are using satellite-based synthetic aperture radar (SAR) to monitor land deformation. SAR is an advanced remote-sensing technology that has been used worldwide in several applications; it is considered active remote sensing as the sensor emits its own microwaves before recording the backscattered waves [18]. In addition, another advantage of SAR is its flexibility in acquiring data, as it can operate under any weather conditions and all day [19]. Thus, because of its operational flexibility, SAR is considered reliable for damage assessment and risk analysis.

Recently, there has been growing interest in monitoring ground deformation using differential interferometric SAR (DInSAR), which is used to estimate such small deformation. However, this method is often restrained by temporal decorrelation and atmospheric disturbance. Those problems are difficult to eliminate when estimating the representative interferograms of ground deformations [20]. To solve this problem, a long-term InSAR processing approach using a series of data, generally recognized as time-series InSAR analysis. Time-series analysis has been widely used for ground-monitoring projects, such as monitoring urban areas, groundwater, and land subsidence [21–23]. An example of InSAR time-series methods that have been developed in recent years is Persistent Scatterer InSAR (PS-InSAR). With time-series data accumulated over long periods of monitoring, deformation can be measured over broader areas, and the occurrence of deformation in an area can be studied with more precision compared to data from conventional InSAR methods. In addition, a prominent advantage of time-series analysis is a reduction in issues related to InSAR, such as temporal decorrelation and atmospheric interference.

In PS-InSAR, multiple large SAR images of the same area are processed, which can extract a number of persistent scatterers (PS) [24]. This method focuses on measuring the level of deformation associated with each of the persistent scatterers, which is a point of high density within an interferogram from a single main image. Therefore, this technique was developed for urban areas, which have more stable scatterers than mountains or forests with distributed scatterers [25]. In recent years, many PSI approaches have been developed and applied in many cases, such as PS-InSAR [26–28], Stable Point Network (SPN) [29], Interferometric Point Target Analysis (IPTA) [30–32], and Stanford Method of Persistent Scatterers (StaMPS) [33,34]. This method has been applied to measure deformation in several locations such as West Macedonia, Greece [35]; Nile Delta, Egypt [23]; Tuscany and Northern Apennines, Italy [20,36,37]; Mexico City, Mexico [11]; and Guangzhou, China [30,38].

In general, susceptibility maps are generated as part of land-subsidence mapping and modeling. These maps are used to predict and examine areas that are highly vulnerable to land subsidence and are important as part of the initial information used to prevent future land-subsidence events, assuming the same conditions will trigger land subsidence in the future [39]. Land-subsidence analysis largely focuses on elucidating the sources of subsidence, methods of evaluation, understanding subsidence events on a conceptual level, and mapping [40]. Remote sensing and Geographical Information System (GIS) analysis are commonly employed in hazard studies because of the efficiency of these techniques in data collection and analysis [41]. The utility of several GIS methods in terms of generating subsidence-susceptibility maps has been compared using statistical tools such as frequency ratio and weight of evidence [6,37]. Artificial neural networks, random forests, and fuzzy logic have also been applied to assess methods of predicting subsidence susceptibility [42–44]. Furthermore, machine-learning techniques have recently attracted attention in the environmental-modeling research community as they are highly efficient and generate improved outcomes.

Previous studies related to land subsidence in Seoul related to the risk analysis of ground subsidence around railways using ANN modeling [45]. However, studies related to land subsidence in Seoul are needed to map the potential for land subsidence and understand the causes that can lead to land subsidence. In this study, we aimed to examine ground deformation within the 2017–2020 time period in the Seoul Metropolitan Area using C-band data from the Sentinel-1 satellite. The velocity of land subsidence and the associated land deformation was investigated via Stanford Method for Persistent Scatterers (StaMPS), and a deformation time-series map was generated. The map was used as input data for a land-subsidence inventory map, which was then used to generate a land-subsidence-susceptibility map. GIS analysis was employed to complement observations of land subsidence in the urban area. Thereafter, the performance of the meta-ensemble methods in terms of land-subsidence-susceptibility modeling was assessed. The LogitBoost, Multiclass Classifier, and Bagging functional models were used as they have been shown to improve the performance of predictive models in several cases of susceptibility map [43,46,47]. Bagging more effectively reduces the bias compared to other ensembles, while LogitBoost is used to solve the overfitting problem [43,48]. Model performance in terms of predicting land-subsidence susceptibility was assessed using training and testing datasets, receiver operating characteristic (ROC), and the area under the curve (AUC). The area under the curve of the roc represents the validation of systems and its ability to predict the correct occurrence or non-occurrence of land subsidence events [49]. The ROC graph is a technique for visualizing, organizing, and selecting classifier based on their performance. Then, the area under the curve is a common method to calculate the area under the curve to compare the classifiers and convert it to the scalar values which represent the performance [50]. By linking time-series data with GIS analysis, the research findings will increase our understanding of crucial factors affecting ground surface areas. The map of ground-deformation risk in an urban area generated in this study reflects the association between subsidence risk and risk factors; this information would help identify areas of high risk and develop environmental action plans and policies.

#### **2. Materials & Methods**

#### *2.1. Study Area*

Seoul is the capital city of South Korea. It is located in the midwestern region of the Korean Peninsula at 126◦59040"E and 37◦33059"N and covers an area of 605.5 km2 [51]. The Han River, which is one of the largest rivers crossing Seoul, divides the city into north and south areas. Seoul has a population of approximately 10 million people with a density of 16,364 people/km2 , making it one of the most populous metropolitan cities in Asia [52]. The geological setting of Seoul consists of Jurassic granite, Precambrian metamorphic rocks (gneiss and schist), and Quaternary alluvium. Predominantly, coarse-grained, sandy alluvium sequence (<20 m thick) occurs along the Han River and its tributaries [53]. The alluvium is mainly distributed along the Han River and its tributaries, it is composed of coarse- to fine-grained sediments, often with high permeability. The alluvium and soil tend to be thicker close to the river, particularly its lower reaches, and thinner in mountain area [54]. In this study area, there are two types of aquifer unconsolidated alluvium aquifers and bedrock aquifers [17]. The alluvial aquifers are dominantly composed of silt and fine to coarse sands are appearance along the Han river and tributaries. The bedrock aquifers are mainly composed of fractured gneiss, schist, and granite.

As a metropolitan city that has experienced urbanization in recent years, Seoul has experienced many developments such as office, business, and residential buildings. This has an impact on increasing the density of the building in this area. The use of groundwater and other utilities in densely populated areas will have an impact on the weakening of soil conditions in these areas. This condition can indirectly lead to subsidence which can cause many losses, especially in areas with high population density such as Seoul.

Together with the increase in population, the economy has grown quickly, followed by industrial development. To meet the needs of the city inhabitants, Seoul undertook massive developments, including infrastructure, buildings, and transportation networks. At the end of 2019, a total of 23 rapid transit, light metro, commuter rail, and airport rail lines had been integrated into the Seoul Metropolitan Subway system [55]. This system operates in the Seoul Metropolitan Area, including Incheon and some

satellite cities in Gyeonggi Province. Several regional lines such as those in Chungnam and Gangwon provinces are also connected to this system. Figure 1 has shown a map of the subway line that has been operating and in the process of construction in the Seoul Metropolitan Area. New transportation routes have since been added, such as the Gimpo Gold-Line in 2019, and the Line 7 and Line 5 extensions to Hanam City is currently under construction and slated to open in 2020. To improve connectivity in this metropolitan area, several future subway lines (until 2028) are still being planned. in 2019, and the Line 7 and Line 5 extensions to Hanam City is currently under construction and slated to open in 2020. To improve connectivity in this metropolitan area, several future subway lines (until 2028) are still being planned. In densely populated areas like Seoul, ground subsidence can cause much higher casualties and property damages than in lesser populated areas. For this reason, it is of utmost importance to conduct complete monitoring on the cases of ground subsidence to prevent damages on roads, railroads, and other infrastructures.

*Remote Sens.* **2019**, *12*, 3505 4 of 27

in Chungnam and Gangwon provinces are also connected to this system. Figure 1 has shown a map

Metropolitan Area. New transportation routes have since been added, such as the Gimpo Gold-Line

**Figure 1.** Optical image of Seoul (gray) captured by Sentinel-2 on 19 June 2019 with subway lines **Figure 1.** Optical image of Seoul (gray) captured by Sentinel-2 on 19 June 2019 with subway lines superimposed.

*2.2. SAR Datasets*  In this study, SAR data from Sentinel-1B was used to generate representations of surface deformation. SAR images derived from C-band data can be used to map surface deformation over In densely populated areas like Seoul, ground subsidence can cause much higher casualties and property damages than in lesser populated areas. For this reason, it is of utmost importance to conduct complete monitoring on the cases of ground subsidence to prevent damages on roads, railroads, and other infrastructures.

#### broad areas while providing time-consistent ground-deformation data. Sentinel-1B has an acquisition *2.2. SAR Datasets*

superimposed.

cycle of 12 days. We used 93 SAR scenes from descending tracks. The descending datasets are listed in Table 1, the reference date with zero delta day and zero perpendicular baselines from the descending track is shown on October 11, 2018 as the reference date are shown in bold text. **Table 1.** Acquisition dates of data from the Sentinel-1 satellite in descending tracks, the reference date shown in bold text. **Acquisition Acquisition Acquisition**  In this study, SAR data from Sentinel-1B was used to generate representations of surface deformation. SAR images derived from C-band data can be used to map surface deformation over broad areas while providing time-consistent ground-deformation data. Sentinel-1B has an acquisition cycle of 12 days. We used 93 SAR scenes from descending tracks. The descending datasets are listed in Table 1, the reference date with zero delta day and zero perpendicular baselines from the descending track is shown on October 11, 2018 as the reference date are shown in bold text.

**Days B**⊥

ሺሻ**1 No.**

**Date (yyyymmdd)**  **Days B**⊥ ሺሻ **1**

**Date (yyyymmdd)**

#### **No. Date (yyyymmdd) Days B**⊥ *2.3. StaMPS Processing*

1 20170302 −588 78 32 20180414 −180 102 63 20190421 192 125 2 20170314 −576 128 33 20180426 −168 67 64 20190503 204 187 3 20170326 −564 128 34 20180508 −156 39 65 20190515 216 97 4 20170407 −552 91 35 20180520 −144 17 66 20190527 228 91 5 20170419 −540 61 36 20180601 −132 −2 67 20190620 252 82 6 20170501 −528 83 37 20180613 −120 117 68 20190702 264 160 7 20170513 −516 117 38 20180625 −108 87 69 20190714 276 117 One of the known methods for generating time-series data on surface deformation is StaMPS, as it can be used for analysis on the urban area like this study area. This method is commonly recognized on man-made objects such as buildings, infrastructure, and roads in urbanized areas like Seoul. Another major advantage of this method is that it does not require a prior deformation model, thus allowing analysis of different regions and several deformation causes [24].

8 20170525 −504 82 39 20180707 −96 65 70 20190807 300 41 9 20170606 −492 121 40 20180719 −84 93 71 20190819 312 100 10 20170618 −480 49 41 20180731 −72 24 72 20190831 324 87 11 20170630 −468 14 42 20180812 −60 67 73 20190912 336 99

ሺሻ **1 No.**

**Table 1.** Acquisition dates of data from the Sentinel-1 satellite in descending tracks, the reference date shown in bold text.

1 B⊥: Perpendicular Baseline.

At the start of PSI-StaMPS analysis, the interferogram process was begun to generate a couple of interferogram images from the 93 SAR scenes in descending track. Prior to interferogram generation, SAR data underwent a co-registration process, in which two SAR images were aligned to subpixel accuracy for accurate determination and noise reduction to form interferometric pairs.  When the co-registration process was complete, the co-registered images were cropped to focus on the study area before the interferogram generation process began. During the interferogram processing stage, a topographical phase was generated. Once the interferogram images were generated, the topographic phase was subtracted from the interferogram using Shuttle Radar Topography Mission digital elevation model (SRTM DEM) as the reference [56]. After the topographic phase was removed, a DInSAR phase was generated, which contained only the deformation phase.

For StaMPS processing, SAR data on 11 October 2018 was chosen as a master image and generated 92 interferograms from the descending track. After generating the interferograms, the StaMPS module was used to calculate the displacements of persistent scatterers. To begin the StaMPS process, phase stability estimation was used to select a subset of pixels based on amplitude analysis. Then, phase stability for each pixel was estimated through phase analysis [57]. Once the phase noise associated with all selected persistent scatterer (PS) pixels was estimated, the selected PS pixels were weeded out to separate the persistent points and noise, then the wrapped phase of the selected pixels was corrected for spatially-uncorrelated look-angle errors in the DEM. After correction, the corrected phase could now be unwrapped, and the PS output was generated. The parameter of the StaMPS process is shown in Table 2. Upon completion of the StaMPS process, PS results were plotted in a time-series map and a mean deformation from the line of sight (LOS) map [58,59]. The mean deformation map is converted into vertical deformation data by assuming the horizontal deformation is very small compared to the vertical deformation that causes by land subsidence [60,61]. In recent studies, the vertical deformation was used to monitor land subsidence in several places, therefore horizontal deformation in this study can be assumed as negligible and converted into vertical deformation value [13,14]. The result of vertical deformation will be assigned as a negative value from initial ground-level observation, which indicates the land subsidence measure vertically on that point. The vertical deformation can be calculated using this Equation (1) as follows [62]:

$$\mathbf{V} = \frac{\mathbf{d}\_{\text{LOS}}}{\cos \theta} \tag{1}$$

where dLOS is the deformation in line of sight and θ is the incident angle.

**Table 2.** The parameter in StaMPS processing.

#### *2.4. Generation of Susceptibility Map*

The workflow to generate land subsidence susceptibility maps, using machine learning algorithms, is illustrated in Figure 2, the summary of the methodology is as follows:

5. After the land subsidence susceptibility map was generated, all susceptibility maps were evaluated using ROC analysis.

*Remote Sens.* **2019**, *12*, 3505 7 of 27

**Figure 2.** The workflow illustration. **Figure 2.** The workflow illustration.

#### 2.4.1. Bagging 2.4.1. Bagging

Bagging is a commonly used meta-algorithm that was developed to enhance the stability and accuracy of the machine-learning algorithms used in statistical classification and regression [43]. Bagging was one of the earliest ensemble techniques that used the bootstrap sampling method [68]. The bootstrap method entails apparent random sampling with replacements to generate more than one sample that shapes a training set. Each generated subset is used to assemble a decision tree, with all trees aggregated later into the final model. This improves class accuracy by reducing the variance of class error. We used Bagging to obtain a much improved and more accurate land subsidence model because this algorithm performs well in predicting land subsidence susceptibility, as it is sensitive to small adjustments in the training data and consequently [43,46]. Bagging ensembles more effectively reduce uncertainty and bias compared to other ensembles [69]. In addition, this algorithm is capable of reflecting complex non-linear interaction between land subsidence and related factors, although it lacks a statistical significance test which can limit quantitative hypothesis testing [43]. Bagging first uses a classifier to reduce variance, then carries out classification and regression by relying on bottom-up learning. Bagging is a commonly used meta-algorithm that was developed to enhance the stability and accuracy of the machine-learning algorithms used in statistical classification and regression [43]. Bagging was one of the earliest ensemble techniques that used the bootstrap sampling method [68]. The bootstrap method entails apparent random sampling with replacements to generate more than one sample that shapes a training set. Each generated subset is used to assemble a decision tree, with all trees aggregated later into the final model. This improves class accuracy by reducing the variance of class error. We used Bagging to obtain a much improved and more accurate land subsidence model because this algorithm performs well in predicting land subsidence susceptibility, as it is sensitive to small adjustments in the training data and consequently [43,46]. Bagging ensembles more effectively reduce uncertainty and bias compared to other ensembles [69]. In addition, this algorithm is capable of reflecting complex non-linear interaction between land subsidence and related factors, although it lacks a statistical significance test which can limit quantitative hypothesis testing [43]. Bagging first uses a classifier to reduce variance, then carries out classification and regression by relying on bottom-up learning.

#### 2.4.2. LogitBoost 2.4.2. LogitBoost

follows:

LogitBoost is a boosting algorithm developed by Friedman et al., [70] to reduce bias and variance. The LogitBoost algorithm was modified from AdaBoost, which was the commonly boosting method for handling noisy data that execute additive logistic regression with least-square fits for individual class [48,71]. LogitBoost reduces training errors and enhances classification accuracy [72] by using additive logistic regression for classification with a base-learning regression scheme and an ability to perform multiclass classification. The land subsidence-inventory map was divided into two classes subsidence occurrence and subsidence non-occurrence—using the following equation [71]: LogitBoost is a boosting algorithm developed by Friedman et al., [70] to reduce bias and variance. The LogitBoost algorithm was modified from AdaBoost, which was the commonly boosting method for handling noisy data that execute additive logistic regression with least-square fits for individual class [48,71]. LogitBoost reduces training errors and enhances classification accuracy [72] by using additive logistic regression for classification with a base-learning regression scheme and an ability to perform multiclass classification. The land subsidence-inventory map was divided into two classes—subsidence occurrence and subsidence non-occurrence—using the following equation [71]:

$$\text{Lc}(\mathbf{c}) = \sum\_{i=1}^{D} \mathfrak{beta}\_i \mathbf{x}\_i + \mathfrak{k}\_0 \tag{2}$$

where D is the number of landslide-dependent factors and β*i* is the coefficient of the *i*-th component within input vector x. Probabilities were constructed using the linear logistic regression method, as follows:

$$\mathrm{p}\left(\frac{\mathsf{C}}{\mathsf{x}}\right) = \mathrm{exp}(\mathsf{Lc}(\mathsf{x})) / \sum\_{\mathsf{C}'=1}^{\mathsf{C}} \exp(\mathsf{Lc}'(\mathsf{x})) \tag{3}$$

where C is the number of classes and the least-square fit Lc(x) is resolved such P C=1 L C C (x)= 0 to set up the least number of instances per node of the logistic model trees.

#### 2.4.3. Multiclass Classifier

Multiclass Classifier is a meta-classifier that is used to process multiclass datasets with two-class classifiers. It is efficient at applying error-correcting output codes to enhance accuracy [47]. In the field of machine learning, multiclass classification can classify events into one of three or more classes. Although several classification algorithms can work with more than two classes with the aid of natural binary algorithms, the conversion to multinomial classifiers requires the use of several strategies. Multiclass classification techniques can be divided into categories such as transformation to binary, extension from binary, and hierarchical classification [73].

#### *2.5. Factors Related to Land Subsidence*

The increase in land subsidence occurrence in megacities mainly due to lowered groundwater levels and the presence of heavy buildings [17]. Excessive use of groundwater has an impact on decreasing pore pressure on the soil, coupled with the presence of heavy buildings that lead to further soil compaction [74,75]. A large amount of groundwater leaked, and dewatering accompanies a decrease in groundwater levels. Those conditions weaken the surrounding land and lead to subsidence occurrence. In addition, based on a report by the Seoul government which has conducted field investigations, one of the causes of subsidence is excessive groundwater use and damage to water utilities and sewage [3,17].

A combination of several environmental factors can influence land-subsidence susceptibility. The training data and test data point was chosen from the land subsidence inventory map as shown in Figure 3a. Here, we investigated 10 subsidence-related factors (Figure 3b–k) and evaluated the correlation between each factor and land-subsidence occurrence as shown in Table 3 below. In the previous studies, land subsidence conditioning factors such as altitude, slope, aspect, plan curvature, profile curvature, lithology, distance to the river, land use, normalized differential vegetation index (NDVI), piezometric data (groundwater drawdown) have been used with the main cause of subsidence in Iran being groundwater drawdown [49]. Another study has evaluated several factors mentioned before to identify land subsidence in the mining area in Malaysia [40]. Reclassification was employed to place subsidence related-factors into several classes using the quantile method to objectively identify and analyze the effect of each class using a specific range of values. The quantile classification method can solve unbalanced distribution by focusing on the equality of domain grids [76]. Thus, the range of each class is automatically determined based on the quantile method.

The derivative feature from the digital elevation model (DEM) contains hydrogeological and topographic conditions such as hydrological zone response, concentration, and containment of runoff volume in the landscape, which directly or indirectly affect the occurrence of land subsidence. The topography features such as elevation, slope, aspect, topographic wetness index (TWI), and profile curvature (Figure 3b–f) data were extracted using SRTM DEM 1 arc second. This feature has been widely used as conditional factors in the land subsidence susceptibility model [37,43,49]. Elevation has a role as a bridge between lithology and rain characteristics in the area. A higher area has a lower probability of additional precipitation than a lower area which has the potential to have high precipitation [49]. The elevation refers to the height of the study area which varies between 0–813 m.

occurrences.

subsidence area has correlated with the fault distance in the area between 0–1,972 m. The area with a

(**a**)

(**d**) (**e**)

**Figure 3.** *Cont*.

*Remote Sens.* **2020**, *12*, 3505 *Remote Sens.* **2019**, *12*, 3505 11 of 27

(**k**).

**Figure 3.** Train-test data point and spatial maps of 10 land-subsidence-related factor: train-test data point (**a**), elevation (**b**), slope (**c**), aspect (**d**), topographic wetness index (TWI) (**e**), profile curvature (**f**), land use (**g**), groundwater extraction (**h**), distance to river (**i**), lithology (**j**), and distance to fault **Figure 3.** Train-test data point and spatial maps of 10 land-subsidence-related factor: train-test data point (**a**), elevation (**b**), slope (**c**), aspect (**d**), topographic wetness index (TWI) (**e**), profile curvature (**f**), land use (**g**), groundwater extraction (**h**), distance to river (**i**), lithology (**j**), and distance to fault (**k**).

**Table 3.** Frequency Ratio calculation.

Slope and aspect factor as a secondary feature of DEM may relate to the land subsidence occurrence because it can affect the soil infiltration in the landscape and water utility conditions [77,78]. The groundwater or sewage infiltration through a damaged pipe may erode the soil particles [77]. That condition can indirectly affect soil conditions that lead to land subsidence occurrence.

TWI is a secondary topographic variable that specifies the degree of water accumulation in a certain location; it is commonly used to quantify topographic influence on hydrological processes. The TWI of this study area was prepared based on the DEM and categorized into five classes. Profile curvature is a geomorphic property which shows the flow intensity, the amount of sediment, and erosion [79]. The profile curvature map was categorized into three classes: convex (less than −0.01), flat (−0.01–0.01), and concave (larger than 0.01).

Land-use is related to the ecological conditions and anthropological activities of land subsidence occurrence [80,81]. Variation in land use can explain the highly dissected zones within the region and provide insight into the land subsidence activity that is likely to occur. A land-use map was created based on a digital characteristics map provided by (National Geographic Information Institute) NGII; six land-use categories were analyzed in this study and the map is shown in Figure 3g. Underground water utilities are one of the most frequently cited factors impacting land subsidence [54,82].

Groundwater extraction can correlate to land subsidence event, especially in the underground structure [49]. Groundwater extraction map of this study area was prepared from annual average groundwater outflow data measured in 231 points from Seoul Government. Prior to spatial analysis, data on groundwater extraction should be converted into raster data. The accuracy of the raster map depends on the number of data points; however, the availability of groundwater data was limited in this study. To generate raster maps from these limited data, we applied the inverse distance weighted method to make statistical inferences using observed values before interpolating to create raster maps of groundwater extraction with a 30 m × 30 m cell size as shown in Figure 3h.

In order to evaluate the relationship between groundwater conditions and the occurrence of subsidence, several factors related to groundwater conditions can be evaluated such as distance to rivers [43,83]. The distance to rivers is represented by the proximity of the rivers and drainages in the study area [40]. The distance to the river map was calculated based on the map of the river provided by the National Geographic Information Institute (NGII). Then, buffers around the river were created (measured in meters), then the raster map was divided into five classes as shown in Figure 3i.

The geological parameters of a certain area may influence the occurrence of land subsidence, which is related to the lithological and structural variation which leads to differences in strength and permeability of rocks and soil. Lithology has also been an important feature to understand the land subsidence process by describing the structure of underground materials, as most cases of subsidence occurred on landfills and alluvium layers that have natural consolidation. Besides, the groundwater withdrawal and load from the building induce the compaction rate of the alluvium [17,62]. Any fluid present in the porous medium structure is under pressure because of the weight of the structure above it. If the fluid is withdrawn from below the surface, a decrease in pore pressure can occur, resulting in the loss of the supports and possibly lead to subsidence [74]. The lithology map can be seen in Figure 3j and the description of the lithology is shown in Table 4.

**Table 4.** Description of the lithological units in the study area.

Figure 3k shows a raster map of distance to the fault which used in this research. The presence of a fault line may weaken the porous medium structure and influence the subsidence occurrence, as in the case of Las Vegas, USA [84]. In this study, we used the fault lines as one of the factors of land subsidence to consider the impact of the fault line and the ground deformation. We performed a buffering distance from the fault line with the data published by the Korean Institute of Geoscience and Mineral Resources (KIGAM) with a 1:50,000 scale then categorized into five classes.

The spatial correlation for each factor was calculated using the frequency ratio and shown in Table 3. For the classes with frequency ratio values close to one or more, it shows a high correlation between subsidence and class of those factors, and vice versa [40]. From this calculation, it is considered that the subsidence has a correlation in the area with characteristics low elevation (0–30 m) and the flat area. The three classes of slope map (0–1.8, 1.8–3.86, 3.86–7.97 degree) and four classes of aspect map (Flat, North, East, Southeast) shows a spatial correlation with the land subsidence from this frequency ratio calculation. In addition, subsidence correlated with the drying area covered by building and non-permeable surface as shown in the land-use factor. Also, the ratio of alluvium (Qa) which dominates appear around the Han river exhibit a correlation with subsidence. Besides, there are six categories from this map that correlate with this calculation of lithology factor too. The land subsidence area has correlated with the fault distance in the area between 0–1972 m. The area with a groundwater extraction rate above 350 m3 /day has a spatial correlation with the land subsidence occurrences.

#### **3. Results**

The results from the PSI-StaMPS time-series analysis on deformation and the land-subsidencesusceptibility map are presented below. The time-series results were obtained by selecting location points at which subsidence occurred.

#### *3.1. Land Subsidence Inventory Map*

The land subsidence map from the Seoul Metropolitan Area was generated via the PSI-StaMPS method, using InSAR images in descending track captured from 2017 to 2020, and is shown in Figure 4a. To enhance measurement reliability, the vertical deformation map was generated using mean line-of-sight velocity [62]. Zooming locations for each selected point are shown in Figure 4b–e and the time-series graphs of all points are shown in Figure 4f,g.

In Figure 4b, represented the subsidence map in Gimpo the western part of this study area, with the black line indicating the subway line. In this area, especially in point A, there was a subsidence of 33.5 mm recorded with a mean deformation velocity of 12.57 mm/year from 2017 to 2020. The subsidence in Gimpo mostly occurs along the subway line that was newly operational in 2019. In this location, the subsidence is associated with the compressible deposits which consist of alluvium.

As can be seen from Figure 4c, the subsidence was exhibited around the subway line where the Shincheon subway station is located. Point B was recorded the maximum subsidence of up to 29 mm from 2017–2020 with the mean deformation velocity of 7.34 mm/year. This location consists of the intersection of subway line no 5 and subway line no 2, a residential area with high-density building also appears in this area. The subsidence in point B correlates with groundwater extraction and high-density building, as those conditions influence the subsidence rate in this area.

Figure 4d shows an overview of the subsidence near the Haengsin Station, known to be a depot for the metro train. The observation in point C revealed total subsidence of 34.35 mm from 2017 to 2020 and a mean deformation velocity of subsidence of 10.25 mm/year. The location is characterized by alluvium deposit that dominantly appears around the Han river. The geological features in this area show a correlation with the subsidence.

Figure 4e shows the subsidence map in Hanam city, the eastern part of the study area, with a black line indicates the subway line. The StaMPS result in point D shows the maximum subsidence of 26.17 mm from 2017 to 2020 with a mean deformation velocity rate of 8.42 mm/year. Hanam city is an area that has many developments such as residential areas and commercial buildings; subway construction is also being carried out in this area. Land subsidence can be related to underground work and building construction which pumping a large amount the groundwater [85]. In Hanam city, the subsidence is associated with urban land use and groundwater usage of this area.

Figure 4f,g show the time series graph from four selecting points in the study area. Generally, the periodic subsidence appeared in the vertical deformation graph. A possible reason for periodic subsidence in those areas was seasonal variation in the groundwater level and surface water loading. This result occurred due to the seasonal effect of groundwater extraction, where the selected points were surrounded by high-density buildings that mostly used groundwater as a water source. During the high season of groundwater withdrawal, the groundwater level decreased. After the rainy season, the groundwater level will rise and increase the aquifer system recovery (uplift) [84]. Those conditions may influence the deformation velocity in this study area.

*Remote Sens.* **2019**, *12*, 3505 15 of 27

**Figure 4.** (**a**) Mean vertical deformation map in descending track of the Seoul Metropolitan Area generated from a Hillshade image comprising digital elevation model data from the Shuttle Radar Topography Mission. (**b**) Zooming of vertical deformation map at points A (Gimpo), (**c**) point B (Shincheon), (**d**) points C (Haengsin) and (**e**) points D (Hanam), (**f**) Vertical deformation time-series at points A (Gimpo) and B (Shincheon), and (**g**) points C (Haengsin) and D (Hanam).

#### *3.2. Land Subsidence Susceptibility Map 3.2. Land Subsidence Susceptibility Map*

Land subsidence susceptibility maps were constructed using the training dataset compiled from InSAR time-series data as land subsidence inventory map, ten land subsidence conditioning factors, and three different algorithms. Once the model training process was completed, susceptibility maps were constructed to visualize vulnerability to subsidence in the study area. In the land-subsidence-susceptibility map, each pixel in the study area was assigned a specific subsidence value using the quantile method [47]. Land subsidence susceptibility maps were constructed using the training dataset compiled from InSAR time-series data as land subsidence inventory map, ten land subsidence conditioning factors, and three different algorithms. Once the model training process was completed, susceptibility maps were constructed to visualize vulnerability to subsidence in the study area. In the land-subsidencesusceptibility map, each pixel in the study area was assigned a specific subsidence value using the quantile method [47]. Five susceptibility classes were used to reflect vulnerability to land subsidence: very low, low,

*Remote Sens.* **2019**, *12*, 3505 16 of 27

**Figure 4.** (**a**) Mean vertical deformation map in descending track of the Seoul Metropolitan Area generated from a Hillshade image comprising digital elevation model data from the Shuttle Radar Topography Mission. (**b**) Zooming of vertical deformation map at points A (Gimpo), (**c**) point B (Shincheon), (**d**) points C (Haengsin) and (**e**) points D (Hanam), (**f**) Vertical deformation time-series

Five susceptibility classes were used to reflect vulnerability to land subsidence: very low, low, moderate, high, and very high. Areas of very high susceptibility (marked red in Figure 5a–c) were most frequently found near the Han River and subway lines. The algorithms indicated that the northwestern area is very susceptible to subsidence, which may be due to several factors. For example, the geology of the northwestern area, which is near the Han River, is dominated by alluvium, which likely increases subsidence susceptibility. Most cases of observed subsidence have occurred on alluvium layers exhibiting natural consolidation; additionally, the increasing number of buildings and use of groundwater can exacerbate this condition [15,37,86]. A highly susceptible area was observed in the east, which may be associated with on-going construction in the same area. A few susceptible areas were observed along the northern Han River, mostly comprising high-density buildings and subway stations, but most of the northern area has low-to-moderate susceptibility to subsidence. Groundwater extraction in this area may increase the risk of subsidence, as some areas in which groundwater was extracted are now used for subway stations. Thus, the ground conditions in these areas might have been affected. moderate, high, and very high. Areas of very high susceptibility (marked red in Figure 5a-c) were most frequently found near the Han River and subway lines. The algorithms indicated that the northwestern area is very susceptible to subsidence, which may be due to several factors. For example, the geology of the northwestern area, which is near the Han River, is dominated by alluvium, which likely increases subsidence susceptibility. Most cases of observed subsidence have occurred on alluvium layers exhibiting natural consolidation; additionally, the increasing number of buildings and use of groundwater can exacerbate this condition [15,37,86]. A highly susceptible area was observed in the east, which may be associated with on-going construction in the same area. A few susceptible areas were observed along the northern Han River, mostly comprising high-density buildings and subway stations, but most of the northern area has low-to-moderate susceptibility to subsidence. Groundwater extraction in this area may increase the risk of subsidence, as some areas in which groundwater was extracted are now used for subway stations. Thus, the ground conditions in these areas might have been affected.

**Figure 5.** Land-subsidence-susceptibility map generated using three algorithms: the (**a**) Bagging, (**b**) LogitBoost, and (**c**) Multiclass Classifier algorithms. **Figure 5.** Land-subsidence-susceptibility map generated using three algorithms: the (**a**) Bagging, (**b**) LogitBoost, and (**c**) Multiclass Classifier algorithms.

Figure 6 shows the distribution of pixels in each susceptibility map generated by the meta-ensemble models. In the land-subsidence-susceptibility map generated using the Bagging model, 63.67% of the area exhibited very low susceptibility to subsidence, whereas 15.53%, 7.55%, 6.42%, and 7.04% of the area exhibited low, moderate, high, and very high susceptibility to subsidence, respectively. In the map constructed using the Multiclass Classifier model, 63.66% and 18.31% of the area exhibited very low and low susceptibility to subsidence, respectively, whereas 10.98%, 1.62%, and 5.42% of the area exhibited moderate, high, and very high susceptibility to subsidence, respectively. Lastly, based on the LogitBoost model, most of the area was not very susceptible to subsidence, with 64.44%, 18.07%, 6.62%, 4.41%, and 6.66% of the map classified as areas of very low, low, moderate, high, and very high susceptibility, respectively. The distribution of pixels in very low class and very high class in Figure 6 has a similar pattern between each algorithm. A very high class can be considered as the subsidence area. Meanwhile, medium and high classes are considered as areas of future land subsidence and very low and low classes are areas with the lowest probability of land subsidence in the future. With this description, it is possible to know the area and the extent of the potential for subsidence that will occur in the future. Generally, the consistency of this model can be evaluated based on the presence of past land subsidence in land subsidence susceptibility classes. The existence of a higher percentage of land subsidence pixels in a higher degree of susceptibility classes indicates higher consistency and vice-versa. 7.04% of the area exhibited low, moderate, high, and very high susceptibility to subsidence, respectively. In the map constructed using the Multiclass Classifier model, 63.66% and 18.31% of the area exhibited very low and low susceptibility to subsidence, respectively, whereas 10.98%, 1.62%, and 5.42% of the area exhibited moderate, high, and very high susceptibility to subsidence, respectively. Lastly, based on the LogitBoost model, most of the area was not very susceptible to subsidence, with 64.44%, 18.07%, 6.62%, 4.41%, and 6.66% of the map classified as areas of very low, low, moderate, high, and very high susceptibility, respectively. The distribution of pixels in very low class and very high class in Figure 6 has a similar pattern between each algorithm. A very high class can be considered as the subsidence area. Meanwhile, medium and high classes are considered as areas of future land subsidence and very low and low classes are areas with the lowest probability of land subsidence in the future. With this description, it is possible to know the area and the extent of the potential for subsidence that will occur in the future. Generally, the consistency of this model can be evaluated based on the presence of past land subsidence in land subsidence susceptibility classes. The existence of a higher percentage of land subsidence pixels in a higher degree of susceptibility classes indicates higher consistency and vice-versa.

*Remote Sens.* **2019**, *12*, 3505 17 of 27

63.67% of the area exhibited very low susceptibility to subsidence, whereas 15.53%, 7.55%, 6.42%, and

Figure 6 shows the distribution of pixels in each susceptibility map generated by the meta-

**Figure 6.** Distribution of pixels classified as areas of very low, low, moderate, high, and very high susceptibility in the land-subsidence-susceptibility maps generated by three machine-learning **Figure 6.** Distribution of pixels classified as areas of very low, low, moderate, high, and very high susceptibility in the land-subsidence-susceptibility maps generated by three machine-learning algorithms.

#### algorithms. *3.3. Model Validation*

subsidence susceptibility in the study area.

*3.3. Model Validation*  A good land-subsidence-susceptibility map should be able to predict future subsidence in the target area and provide initial information for preventative actions. To validate our susceptibility maps, the accuracy of all used algorithms in this study was evaluated by ROC curve analysis. ROC curve analysis has been used as a standard way of validating the probability models used to generate land subsidence susceptibility maps, according to the area under the curve (AUC) [6]. The AUC, which ranges from 0.5 to 1, was used to assess model accuracy. An AUC value near 0.5 indicates that a model is inaccurate, whereas a value near 1.0 indicates an ideal model with a good fit [50]. AUCs were calculated to compare model performance, with the model with the highest AUC value was taken to be the best model. The Bagging model produced the largest AUC (0.883), followed by the LogitBoost model (0.871) and the Multiclass Classifier model (0.856) as shown in Figure 7. Thus, the Bagging model generated the best subsidence-susceptibility map in this study. However, all models produced good AUC values, indicating that they all performed well in terms of predicting land-A good land-subsidence-susceptibility map should be able to predict future subsidence in the target area and provide initial information for preventative actions. To validate our susceptibility maps, the accuracy of all used algorithms in this study was evaluated by ROC curve analysis. ROC curve analysis has been used as a standard way of validating the probability models used to generate land subsidence susceptibility maps, according to the area under the curve (AUC) [6]. The AUC, which ranges from 0.5 to 1, was used to assess model accuracy. An AUC value near 0.5 indicates that a model is inaccurate, whereas a value near 1.0 indicates an ideal model with a good fit [50]. AUCs were calculated to compare model performance, with the model with the highest AUC value was taken to be the best model. The Bagging model produced the largest AUC (0.883), followed by the LogitBoost model (0.871) and the Multiclass Classifier model (0.856) as shown in Figure 7. Thus, the Bagging model generated the best subsidence-susceptibility map in this study. However, all models produced good AUC values, indicating that they all performed well in terms of predicting land-subsidence susceptibility in the study area.

**Figure 7.** Receiver operating characteristic (ROC) curves associated with land-subsidencesusceptibility maps generated by three meta-ensemble algorithms (Bagging, LogitBoost, and **Figure 7.** Receiver operating characteristic (ROC) curves associated with land-subsidence-susceptibility maps generated by three meta-ensemble algorithms (Bagging, LogitBoost, and Multiclass Classifier).

#### Multiclass Classifier). **4. Discussion**

#### **4. Discussion**  *4.1. Land Subsidence Inventory Map*

*4.1. Land Subsidence Inventory Map*  StaMPS was employed to analyze land subsidence in the Seoul Metropolitan Area, with a deformation time-series map generated for all-terrain in the area based on descending-track data StaMPS was employed to analyze land subsidence in the Seoul Metropolitan Area, with a deformation time-series map generated for all-terrain in the area based on descending-track data acquired from March 2017 to May 2020. Subsequently, a vertical deformation map was generated from the time-series analysis.

acquired from March 2017 to May 2020. Subsequently, a vertical deformation map was generated from the time-series analysis. The results indicate that occurrences of subsidence were distributed over several locations in the study area, such as at Gimpo City (Figure 4a, point A). At this location, subsidence occurred near a new subway line that opened in 2019. Subsidence also occurred near a station with two subway lines near the southern Han River; high-density buildings are also found in this area (Figure 4a, point B). Other areas where subsidence occurred include a Seoul metro depot for several subway lines (Figure 4a, point C) and a newly developed area with several recently constructed residential and commercial buildings (Figure 4a, point D). Cases of land subsidence in the Seoul Metropolitan Area almost mostly occurred in the vicinity of subway lines or where the ground was weak. In particular, all areas of subsidence were located near subway lines and stations, implying that subway operation may be associated with subsidence in the study area [38,87]. Additionally, subway-tunnel excavations might have impacted the surrounding soil and the environment. During construction, how underground The results indicate that occurrences of subsidence were distributed over several locations in the study area, such as at Gimpo City (Figure 4a, point A). At this location, subsidence occurred near a new subway line that opened in 2019. Subsidence also occurred near a station with two subway lines near the southern Han River; high-density buildings are also found in this area (Figure 4a, point B). Other areas where subsidence occurred include a Seoul metro depot for several subway lines (Figure 4a, point C) and a newly developed area with several recently constructed residential and commercial buildings (Figure 4a, point D). Cases of land subsidence in the Seoul Metropolitan Area almost mostly occurred in the vicinity of subway lines or where the ground was weak. In particular, all areas of subsidence were located near subway lines and stations, implying that subway operation may be associated with subsidence in the study area [38,87]. Additionally, subway-tunnel excavations might have impacted the surrounding soil and the environment. During construction, how underground water is discharged, and the excavation method should be taken into consideration. Further analysis is needed to examine the impacts of construction on land subsidence in this area.

water is discharged, and the excavation method should be taken into consideration. Further analysis is needed to examine the impacts of construction on land subsidence in this area. In terms of urbanization, the construction of buildings near subway lines and stations may add to the load on the soil and increase the risk of subsidence. High demand for transportation and urbanization increases the intensity of building construction and the amount of groundwater extracted. Besides, the groundwater extraction to fulfill social demand can influence the subsidence rate in the study area [3,88]. From the time-series deformation graph in Figure 4f, we can see a seasonal variation of subsidence rate in the study area between 2017–2020. It can also be noted that the high subsidence may appear in summer seasons. On the other hand, the subsidence rates lower after the rainy season which appears in July–August [3], the aquifer conditions after the rainy season are expected to have some influence on the subsidence rate. More study based on the water-level data analysis is required to better assess the possibility of the deformation, and further details of the In terms of urbanization, the construction of buildings near subway lines and stations may add to the load on the soil and increase the risk of subsidence. High demand for transportation and urbanization increases the intensity of building construction and the amount of groundwater extracted. Besides, the groundwater extraction to fulfill social demand can influence the subsidence rate in the study area [3,88]. From the time-series deformation graph in Figure 4f, we can see a seasonal variation of subsidence rate in the study area between 2017–2020. It can also be noted that the high subsidence may appear in summer seasons. On the other hand, the subsidence rates lower after the rainy season which appears in July–August [3], the aquifer conditions after the rainy season are expected to have some influence on the subsidence rate. More study based on the water-level data analysis is required to better assess the possibility of the deformation, and further details of the structure and hydrologic parameters of groundwater should be resolved [89]. A combination of InSAR time-series analysis and

analysis of hydrology of the area subsidence and the geomechanical parameters of the underlying aquifer structure area is a potential research topic to find out the cause of subsidence.

However, information extraction from the StaMPS technique is sometimes difficult due to a large number of PSs, thereby long interpretation times. The large amounts of PSs may cause several deficiencies in the analysis process, such as a reduction in the extraction of useful information from the dataset. Also, to obtain better time-series analysis results, several optimization steps can be taken for persistent scatterer points [20]. Several optimization methods for selecting PS points have been carried out [20,90]. This, in turn, could serve as a reference for future work on land subsidence studies that could increase efficiency and potentially lead to better deformation analyzes. Further studies should investigate other factors related to land subsidence. The results of GIS analysis are discussed below.

#### *4.2. Land Subsidence Susceptibility Maps*

Land subsidence should be mapped accurately to prepare subsidence inventories for target areas as an essential part of susceptibility analysis. To generate a subsidence-inventory map, we used InSAR remote-sensing data, which covered a broad area and were collected efficiently. Moreover, the InSAR time-series (StaMPS) method allows measuring of land subsidence and its rate to be measured to millimeter-level [37]. ArcGIS software was used for database construction, coordinate conversion, overlay analysis, and susceptibility modeling. Subsidence-related factors were identified based on information derived from the literature before susceptibility maps were generated [40,49]. Meta-ensemble machine learning was applied to estimate land-subsidence susceptibility, using three algorithms—Bagging, LogitBoost, and Multiclass Classifier.

Land subsidence susceptibility maps revealed that the northwestern and eastern areas, as well as a small area in the center, were most susceptible to land subsidence. We analyzed subsidence-related factors by comparing general patterns of subsidence with factor maps. The results revealed that most cases of subsidence occurred in areas where the ground consisted of alluvial layers, especially for subsidence that occurred near the Han River. However, there is a potential for subsidence in central areas that have different geological conditions from these two regions which are dominated by the alluvium layer. In this case, there are other factor influences besides geological factors in this subsidence modeling. The central area was assessed as moderately to highly susceptible to subsidence. In this area, there is a high density of buildings, and groundwater had been extracted at several spots near a subway station. Accordingly, groundwater outflow during subway operation could be another cause of land subsidence in this area [91]. If large amounts of groundwater are extracted, the surrounding soil structure may be affected. Thus, the weakened soil may be less able to withstand the pressure from aboveground buildings. Based on spatial distribution analysis, land use and groundwater extraction most strongly influence subsidence. In this study, the groundwater-extraction map was obtained via the interpolation of data on groundwater extraction near the subway infrastructure in Seoul, which might have generated errors in spatial distribution. Access to groundwater data for areas outside Seoul would allow for more accurate analyses of land subsidence.

In addition, several other factors that have not been identified in this study can be evaluated. The selection of these factors is based on the previous literature which may have some differences such as the condition of the area study and the subsidence mechanism. For this reason, additional analysis of factors related to the subsidence mechanism is needed to adjust to several concepts or assumptions of the subsidence mechanism that could potentially occur. However, adding some details such as aquifer conditions and groundwater levels can help evaluate the correlation of these factors and have the potential to improve the land subsidence susceptibility maps.

Susceptibility maps were validated based on ROC curves and AUC values and by comparing map predictions with testing data, which comprised 50% of the total dataset. The results indicate that the meta-ensemble approach performed better than the other approaches. A traditional model based on frequency ratio produced an AUC of 0.844, whereas the AUCs produced by the meta-ensemble models Multiclass Classifier, LogitBoost, and Bagging were 1.2%, 2.7%, and 2.95% larger, respectively. Therefore, these techniques can reduce bias and account for factor weights to improve the accuracy of predictions.

Although all models produced good AUC values and thus performed well in terms of predicting land subsidence in the study area, the susceptibility map constructed using the Bagging model was the most accurate. These results agree well with previous findings that model performance in terms of predicting subsidence improves with the use of machine learning [92]. Therefore, the Bagging model should be used for the susceptibility map. In fact, the Bagging model uses more recently well-organized techniques in soft computing modeling that not only enable improvement of a single classifier but can also deal with complex and high-dimensional modeling problems. Given the complexity of land subsidence and the interaction of several related factors, novel combinations of model-method can considerably improve the accuracy of land subsidence prediction.

#### **5. Conclusions**

This study aimed to assess and map land-subsidence susceptibility in the Seoul Metropolitan Area using InSAR data from Sentinel-1 acquired between 2017 and 2020. A deformation time-series map was generated using StaMPS, which revealed that land subsidence occurred in four areas (Figure 4a), with subsidence rates of 6–12 mm/year. Subsidence mostly occurred near subway lines and where a new subway line was being constructed. Besides, the subsidence occurrence in areas with high-density building and heavy groundwater extraction may lead to weakening of the ground.

To identify the factors influencing land subsidence, 10 potential subsidence-related factors were analyzed. Factor maps were overlaid with subsidence maps and each pixel within a layer was evaluated in a GIS environment. The training and testing datasets were prepared from time-series InSAR from the Sentinel-1 SAR dataset using the StaMPS method. Then, the spatial correlation for each factor was calculated using the frequency ratio. Meta-ensemble algorithms (Bagging, LogitBoost, and Multiclass Classifier) were employed to generate land-subsidence-susceptibility maps, and model performance in terms of reliability and prediction accuracy was compared using ROC analysis.

The land-subsidence-susceptibility maps revealed that the northwestern and eastern areas, as well as a small central area, were most vulnerable to land subsidence. The susceptibility of the northwestern and eastern areas most appear in the geological condition which is dominated by alluvium. By contrast, in the central area, which is moderate to highly susceptible, land use and groundwater extraction are the main factors influencing subsidence risk. From the ROC analysis, the AUC produced by each model was computed. All models performed well (AUC > 0.8). Bagging produced the largest AUC of 0.883, followed by LogitBoost (0.871) and Multiclass Classifier (0.856). Compared with the frequency-ratio method, machine-learning models produced more accurate predictions and are thus more appropriate for subsidence analysis in this study area.

Accurate predictions are essential for environmental planning to control and mitigate the impacts of land subsidence. Land-subsidence-susceptibility mapping is a valuable method for identifying areas with a high risk of land subsidence. Despite limitations associated with the datasets used in this study, we demonstrated that the analysis of remote-sensing and GIS spatial data via the machine-learning approach generates reliable and accurate predictions of land subsidence. Further research is needed to determine the effect of aquifer conditions, subway construction and operation on land subsidence. A large dataset of PS points may influence a deficiency in extracting useful information. The optimization approaches for selecting PS points must be proposed to overcome those limitations in future work such as optimization hotspot analysis and other statistic methods [20,90]. Furthermore, with the high complexity of the relationship between land subsidence and other factors, a novel combination of a machine learning and meta-heuristic algorithm as a hybrid method can improve the results of the land subsidence susceptibility map.

**Author Contributions:** Conceptualization, C.-W.L.; methodology, C.-W.L., A.R.A., and M.F.F.; software, A.R.A. and M.F.F.; validation, C.-W.L. and M.F.F.; formal analysis, M.F.F.; investigations, C.-W.L.; resources, C.-W.L.; data curation, C.-W.L. and M.F.F.; writing—original draft preparation, M.F.F.; writing—review and editing, C.-W.L.,

A.R.A., and M.F.F.; visualization, C.-W.L. and M.F.F.; supervision, C.-W.L.; project administration, C.-W.L.; funding acquisition, C.-W.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by a grant from the National Research Foundation of Korea (No. 2019R1A2C1085686), which is funded by the Government of Korea.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 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/).

### *Article* **Mapping Urban Green Spaces at the Metropolitan Level Using Very High Resolution Satellite Imagery and Deep Learning Techniques for Semantic Segmentation**

**Roberto E. Huerta 1 , Fabiola D. Yépez 1,\* , Diego F. Lozano-García 2 , Víctor H. Guerra Cobián 1 , Adrián L. Ferriño Fierro 1 , Héctor de León Gómez 1 , Ricardo A. Cavazos González 1 and Adriana Vargas-Martínez 2**

**Citation:** Huerta, R.E.; Yépez, F.D.; Lozano-García, D.F.; Guerra Cobián, V.H.; Ferriño Fierro, A.L.; de León Gómez, H.; Cavazos González, R.A.; Vargas-Martínez, A. Mapping Urban Green Spaces at the Metropolitan Level Using Very High Resolution Satellite Imagery and Deep Learning Techniques for Semantic Segmentation. *Remote Sens.* **2021**, *13*, 2031. https://doi.org/10.3390/ rs13112031

Academic Editors: Chang-Wook Lee, Hyangsun Han, Hoonyol Lee and Yu-Chul Park

Received: 10 March 2021 Accepted: 27 April 2021 Published: 21 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** Urban green spaces (UGSs) provide essential environmental services for the well-being of ecosystems and society. Due to the constant environmental, social, and economic transformations of cities, UGSs pose new challenges for management, particularly in fast-growing metropolitan areas. With technological advancement and the evolution of deep learning, it is possible to optimize the acquisition of UGS inventories through the detection of geometric patterns present in satellite imagery. This research evaluates two deep learning model techniques for semantic segmentation of UGS polygons with the use of different convolutional neural network encoders on the U-Net architecture and very high resolution (VHR) imagery to obtain updated information on UGS polygons at the metropolitan area level. The best model yielded a Dice coefficient of 0.57, IoU of 0.75, recall of 0.80, and kappa coefficient of 0.94 with an overall accuracy of 0.97, which reflects a reliable performance of the network in detecting patterns that make up the varied geometry of UGSs. A complete database of UGS polygons was quantified and categorized by types with location and delimited by municipality, allowing for the standardization of the information at the metropolitan level, which will be useful for comparative analysis with a homogenized and updated database. This is of particular interest to urban planners and UGS decision-makers.

**Keywords:** neural networks; urban vegetation; urban open spaces; Monterrey Metropolitan Area; sustainable development

#### **1. Introduction**

Urban green spaces (UGSs) face significant challenges due to rapid urbanization and climate change [1]. UGSs are crucial in order to safeguard the quality of urban life [2]. City managers are urged to integrate UGSs in urban development plans [3]. The conservation of ecosystem services of UGSs can mitigate the impacts of urban development; can reduce ecological debts; and is the simplest, fastest, and most effective way to ameliorate important challenges in cities, such as urban heat islands and air pollution [4–6]. UGSs enhance resilience, health, and quality of life of citizens, especially benefiting those with high accessibility. UGS accessibility is a crucial aspect of sustainable urban planning [7] and social justice [8].

UGS survey data are not commonly updated or freely accessible to local users. There is a need for uniform and spatially explicit inventories of existing UGSs [9] and the quantification of their proximity services [10]. One of the most critical functions of the UGSs within cities is the provision of essential environmental services, as it is related to human wellbeing [11,12]. UGSs contribute to the reduction of harmful effects that cause cardiovascular,

respiratory, and metabolic diseases [13], and they mitigate the stress caused by increases in temperature and noise levels [14]. Additionally, UGSs promote physical activity and social interaction, improving the physical and mental health of residents who use these facilities [15].

The constant environmental, social, and economic transformations of cities make the management of UGSs a major challenge for government administrations in metropolitan areas with extensive and rapid urban development [16,17]. Depending on the management of the UGSs, negative or positive effects can be promoted, from the destruction of these spaces [18] to promoting adequate conditions for management and maintenance [19]. Inventories of UGSs allow the monitoring of their status and help provide guidelines for the development of adequate management strategies [20].

UGSs exist in diverse shapes, sizes, vegetation covers, and types (i.e., park, residential garden, median strip, square, and roundabout) [21,22]. Traditional methods to obtain polygons of UGSs have relied on the visual interpretation of aerial imagery, remotely sensed data interpretation, and manual digitalization [23]. Similar to tree inventories and other natural elements in UGSs, data collection methods are intensive and involve manual measurements of dendometric parameters in the field. These methods are timeconsuming and costly [24]. The integration of remote sensing and geographic information systems for mapping and monitoring UGSs has been advantageous, as this reduces the resources required by traditional methods [25–27]. Moreover, the use of these techniques and frameworks allows the computation of vegetation indices that highlight vegetation properties such as vegetation cover and vigor [28].

Methods based on different machine learning algorithms, including decision tree [29], maximum likelihood classification [30], random forest, and support vector regression [31], have been used to map UGSs. Other authors propose the use of object-based image analysis, which takes advantage of both the spectral and contextual information of the classifying objects [32]. With technological advancement and the evolution of deep learning, optimization of the acquisition of UGS inventories is possible through the detection of spectral and geometric patterns available in satellite imagery [33]. Convolutional neural networks (CNNs) have performed well at high-level vision tasks, such as image classification, object detection, and semantic segmentation [34]. A combination of multitemporal MODIS and Landsat-7 imagery was used to classify UGSs in Mumbai metropolitan area in India [35]. The results indicate that for over 15 years, the overall UGSs were reduced to 50%. Other authors analyzed four different methods of classifying UGSs: support vector machine, random forest, artificial neural networks, and naïve Bayes classifier [36]. They found that support vector machines produce higher accuracy classifications in a short amount of time. Multitemporal high-resolution imagery was employed to map open spaces in Kampala, Uganda, with the use of a cloud computation method and machine learning that combined nine base classifiers [37]. The results produced a map of open spaces with an 88% classification accuracy. A deep learning classification based upon a high-resolution network (HRNet) method of high-resolution GaoFen-2 imagery was used for the city of Beijing, China, indicating that the HRNet combined with phenological analysis significantly improved the classification of UGSs [38].

CNNs are convenient models for semantic segmentation because they produce hierarchies that help determine low-, medium-, and high-level characteristics [39,40]. These models are automatically trained using previously labeled input information, and they produce class identification results [41]. With the use of labeled samples, a network can update its weights until it obtains a proper mapping of the inputs and a minimal loss [42].

Due to the absence of a dense layer, the use of fully convolutional networks (FCNs) allows the generation of outputs in which each pixel has a classification according to the input information [43]. Based on the FCN model, the U-Net architecture uses the same principle and considers a symmetric encoder–decoder composition. This process first reduces the size and increases the number of bands of the training images and their activation maps generated in each layer of the network to subsequently carry out the opposite

process considering information from the encoder in the segmentation of fine details [44]. These types of networks have achieved wide success with state-of-the-art results for a wide variety of problems from medical applications [45,46] to their employment in remote sensing for road [47] and building extractions [48], as well as land cover classification [49], but they have not been used to make many advances in the UGS area.

Detailed geometric information on UGSs is typically presented as a shapefile that is not updated frequently; therefore, it does not reflect changes occurring due to rapid urban development processes. Additionally, spatial data or information about the availability of UGSs is not generally accessible to urban residents [50], restricting their use. The need to improve and make available geospatial data of green and public spaces is recognized by the United Nations sustainable development agenda as it helps to create more inclusive, safe, resilient, sustainable cities [51]. The generation of UGS inventories in conjunction with other public space inventories aid in the calculation of Sustainable Development Goal (SDG) 11.7.1, i.e., quantifying the average share of green and public spaces in cities. This allows for the obtainment of SGD 11.7, which is to "provide universal access to safe, inclusive and accessible, green and public spaces". It is therefore necessary to later relate the information available in digitization, vectorization, and computation to demographic data to generate accessibility maps [52,53]. This study evaluates two deep learning model techniques for semantic segmentation of UGS polygons. The process involves different convolutional neural network encoders on the U-Net architecture with the use of threeband compositions of very high resolution (VHR) satellite imagery channels and vegetation indices as input data. This precise and updated data collection and new UGS cartography at the metropolitan level would improve the understanding of connectivity and accessibility of UGSs as a basis for management and decision-making for land use in urban areas.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area chosen to test this method was the Monterrey Metropolitan Area (MMA) (Figure 1), located at the coordinates 25◦4000000 N 100◦1800000 W. It has a total area of 6687.10 km2 of which 27.57% is built-up area. The MMA is comprised of Monterrey, the capital of the State of Nuevo Leon, and 11 surrounding municipalities [54]. Its population, as of 2015, was 4.7 million inhabitants [55]. Within its orography, the Sierra Madre Oriental, Sierra San Miguel, the hills of Topo Chico, La Silla, and Las Mitras are prominent.

**Figure 1.** Study area location. (**A**) State of Nuevo Leon within the Mexican Republic; (**B**) Monterrey Metropolitan Area (MMA) within the state of Nuevo Leon; (**C**) orthomosaic of WorldView-2 coverage of the MMA.

The methodological workflow is shown in Figure 2, where the methodology is divided into three sections: data preprocessing, CNN model implementation, and evaluation of semantic segmentation of UGSs. For this study, we used a UGS definition based on the "Regulation of Environmental Protection and Urban Aesthetics of Monterrey" [56]. This document describes UGSs as land surfaces containing vegetation, gardens, groves, and complementary minor buildings for public use within the urban area or its periphery. Input label polygons for the CNN models were obtained from three sources, the UGS database of the National Geostatistical Framework (2010) of the National Institute of Statistics and Geography (INEGI), the collaborative Open Street Maps (OSM) project [57], and the database of median strips of the MMA arranged by the Department of Geomatics of the Institute of Civil Engineering of the UANL [58].

**Figure 2.** Summary of the current methodological process for semantic segmentation of UGSs using deep learning. Data preprocessing in blue, CNN model implementation in yellow, and evaluation of semantic segmentation of UGSs in green. Input data in gray, and all the processes appear in light blue, yellow, and green. Abbrevation: Urban Green Spaces (UGS), Convolutional Neural Networks (CNN).

#### *2.2. Input Data*

An orthomosaic of the MMA with a 0.5 m pixel resolution was used for the classification of UGSs (Figure 1). It was generated from nine WorldView-2 (WV2) satellite images

obtained between June and October 2017. The spectral information it contained includes the red, green, blue (RGB) and near-infrared (NIR) bands in the ranges of 630–690 nm, 510–580 nm, 450–510 nm, and 770–895 nm, respectively.

#### *2.3. Data Preprocessing*

Information of the three original databases was reclassified based on the function and geometric structure of the UGSs to generate a common UGS database. The database includes polygons representing (1) median strips along streets and avenues, which are characterized by their elongated and narrow shapes; (2) residential gardens, which have pixels that correspond to vegetation managed by the municipality; (3) roundabouts, which have a round shape; (4) squares, which are spaces mostly used for recreation that maintain a symmetry and lack elements related to sports; (5) parks, which are embedded in residential areas, used for recreation and sports, and tend to be asymmetrical. Original classifications (some of them in Spanish) and their equivalent names after the reclassification process are shown in Table 1.

**Table 1.** Label reassignment on the reclassification process of the three databases.

1 Original data in Spanish.

Polygons that presented overlap were discarded with the employment of ArcMap "select by location" tool. The three reclassified databases were merged to produce the input shapefile for a rasterization process. The resulting product had a resolution of 0.5 m and was carried out for the generation of the final label raster. The pixel values determined the presence or absence of UGSs corresponding to median strips, roundabouts, parks, squares, and residential gardens. An additional sixth class named non-UGS was added to cover background pixels.

With the use of the green, red, and NIR bands, normalized difference vegetation index (NDVI) [59,60], two-band enhanced vegetation index (EVI2) [61,62], and normalized difference water index (NDWI) [63] were calculated using Equations (1)–(3), respectively.

$$\text{NDVI} = \frac{(NIR - Red)}{(NIR + Red)} \tag{1}$$

where *NIR* represents the near-infrared channel and *Red* is the red channel.

$$\text{EVI2} = 2.5 \frac{(NIR - Red)}{(NIR + 2.4 \ast Red + 1.0)} \tag{2}$$

where *NIR* represents the near-infrared channel and *Red* is the red channel.

$$\text{NDWI} = \frac{(Green - NIR)}{(Green + NIR)} \tag{3}$$

where *Green* represents the green channel and *NIR* is the near-infrared channel.

This study used 12 three-band compositions to determine their potential for the segmentation of UGSs. The bands used for the combinations were the produced indices (NDVI, NDWI, EVI2) and the spectral single bands NIR, red, green, and blue obtained from the original WV2 data (Figure 1).

As part of the process, 24,667 orthomosaics with dimensions of 256 × 256 pixels were produced from the label raster and each of the 12 three-band compositions was generated for the MMA. To obtain those orthomosaics, the original compositions and their respective label rasters were clipped using first a 2 × 2 mosaic fishnet (Clip1) with resulting orthomosaics that cover an area of 1336.64 km2 . Using the results of Clip1, a second fishnet of 8 × 8 mosaics (Clip2) was applied to each section to obtain 135 segments of 167.08 km2 , with 50 produced for the quadrant of the cardinal NE position, 43 for the NW, 9 for the SE, and 33 for the SW. Both fishnets developed for the generation of training samples from the MMA orthomosaic and UGS labels are shown in Figure 3A. Subsequent split raster process was performed in ArcMap, ArcGIS v10.8.1 software, for each of the quadrants previously generated, and over 24,000 orthomosaics were obtained for each of the three-band combinations as well as their equivalent ground truths. All the data had a spatial resolution of 0.5 m. The data were divided in a proportion of 85% for training and 14% for validation [64]. An additional 1% of the information was used for the evaluation of the model. The results were hosted on Google Drive's cloud storage service for later use through the Google Colab platform, which provided a Tesla P100 PCIe 16 GB GPU.

**Figure 3.** Production of training samples. (**A**) A 1:350,000 scale map with the result of data homogenization and rasterization of UGSs (light green polygons) with two raster extraction processes shown in red mosaics fishnet (Clip 1) and yellow fishnet (Clip 2). (**B**) Irregular geometry 1:7000 scale map of UGSs (light green polygons) extracted from Clip 2 and (**C**) Irregular geometry 1:13,000 scale map of UGSs (light green polygons) extracted from Clip 2.

#### *2.4. CNN Model Implementation*

Twenty-four semantic segmentation models were implemented via CNN, two for each band composition generated. ResNet-34 and ResNet50 encoders pre-trained by the

ImageNet dataset [65] were used in a dynamic U-Net architecture implemented in the fastai deep learning library [66]. This library works with the Python language and the PyTorch library [67] as a backend. Figure 4A shows the architecture for the model based on the pre-trained ResNet-34 encoder. Each model received 256 × 256 pixel images as input. ResNet network architectures integrate connection jumps (Figure 4B) that avoid the leak gradient problem present in other types of networks [68]. This helps to maintain performance and precision despite increases in the number of training layers [69]. At the point of the greatest compression in the FCN, a decoder was attached that follows the principle of the U-Net architecture to finally obtain an output equal in size to the input images.

**Figure 4.** Dynamic U-Net model used for semantic segmentation. (**A**) The encoder consists of a ResNet-34 into which the orthomosaic and UGS labeling information is integrated. The numbers on the left of the graphs are the sizes of the input and output images, and the sizes of the activation maps. The numbers on the right are the channels/filters. (**B**) Building blocks used in the encoder section of the model.

Preliminary tests consisted of trial and error based on the limited literature related to UGS segmentation using deep learning models [38,70]. According to the capabilities of the system, a batch size of 16 samples was assigned. Data augmentation included (1) transformations with image turns at different angles, finding a 50% probability of being horizontal or vertical; (2) random symmetrical deformations with values of 0.1 magnitude; (3) random rotations with angles of 20◦ ; (4) changes in the focus of images, up to 200%; and (5) changes in light and contrast by a factor of 0.3. These techniques generated transformations for each epoch within the models and increased the size of the training samples by 60 times. Examples of image transformations can be observed in Figure 5.

**Figure 5.** Example of images produced by data augmentation. (**A**) Horizontal flip; (**B**) random rotations; (**C**) vertical flip; (**D**) change in focus.

The accuracy score is an evaluation metric that quantifies the percentage of correctly classified pixels made by the predictions of the model [71]. It is calculated by Equation (4).

$$\text{Accuracy Score} = \frac{TP + TN}{\text{TP} + \text{TN} + \text{FP} + \text{FN}} \tag{4}$$

where *TP* represents the true positives, *TN* is the true negatives, FP is the false positives, and FN is the false negatives.

In semantic segmentation, the loss function metric is an algorithm used to evaluate the difference between training results and labeled data. To determine the most appropriate loss function metric for our data, we considered their spatial characteristics. As UGS represents a small portion of pixels at the metropolitan level, most cities are covered by built-up areas occupied by streets, buildings, and other impervious surfaces. This configuration causes an imbalance of classes that could produce errors and bias towards the background class that covers most of the area of interest. Semantic segmentation studies that used deep learning [72–74] have proved that the Dice coefficient or F1 score [75] is a loss function adequate for these kinds of problems. The Dice coefficient is calculated according to Equation (5).

$$\text{Dice Coefficient} = \frac{2|I\_{GT} \cap O\_{SEG}|}{|I\_{GT}| + |O\_{SEG}|} \tag{5}$$

where *IGT* is the input ground truth and *OSEG* is the output segmentation.

A Dice coefficient of 0 indicates that there is no overlap between the data, whereas a value of 1 means that the data has total overlap [76]. Because the input data consisted of three-band composed images, the Dice coefficient was computed for each class and then averaged via arithmetic mean through the fastai implementation [77]. The optimal learning rate for each model was defined using the *learn.lr\_find()* method present in the fastai library [66]. This hyperparameter increases the learning rate from an exceptionally low value to the point where the loss gradient decreases [78]. Each ResNet34 model had 100 epochs to test the functionality of the implementations, this value was established according to the literature [79,80]. ResNet50 was implemented using 10 epochs due to reported Google Colab limitations on GPU, RAM, and session time availability. Mean, standard deviation, and confidence intervals of 95% were produced for the models to determine the statistically significant differences between the Dice coefficient calculated for each three-band combination.

#### *2.5. Evaluation of Semantic Segmentation of UGSs*

The model with the highest Dice coefficient was evaluated using the additional 1% testing subset taken from the original information. This subset, which was produced using the semantic segmentation, results in footprints used to obtain vector information corresponding to UGS polygons. Vectorization was implemented using Solaris library (https://solaris.readthedocs.io/, 26 November 2020) on Google Colab. The generated polygons were downloaded and then analyzed in ArcGIS to evaluate the model. The evaluation data contained 68 mosaic samples (Figure 6) with a total coverage of 32.08 hectares (ha). These samples permitted the evaluation of the effectiveness of the CNN with images different from those of the training and validation sets.

**Figure 6.** The ground truth coverage and the testing subset with 68 mosaics used to evaluate the CNN model.

As part of the evaluation, the intersection over union (IoU) was calculated. This metric computes the amount of overlap between the predicted polygons and the ground truth data [81] (Equation (6)).

$$\text{Interest over Union} = \frac{|I\_{GT} \cap O\_{SEG}|}{|I\_{GT} \cup O\_{SEG}|} \tag{6}$$

where *IGT* is the input ground truth and *OSEG* is the output segmentation.

The recall analysis was obtained for the evaluation. This metric calculates the proportion of positives identified correctly as shown in Equation (7).

$$\text{Recall} = \frac{TP}{TP + FN} \tag{7}$$

where *TP* represents the true positives and *FN* is the false negatives.

The overall user and producer accuracies and the kappa coefficient were processed for the accuracy assessment of the evaluation data. This index of agreement is obtained through the computation of a confusion matrix with errors of omission and commission between classified maps and ground truth data; a kappa coefficient of 1 represents a perfect agreement, and a value of 0 indicates that the agreement is not as it was expected by chance [82,83]. For the calculation of the index, 2000 random points were generated by stratified random sampling and were labeled and verified against the reference data. The kappa coefficient (Equation (8)) is computed as follows:

$$\text{Kappa Coefficient} = \frac{N \sum\_{i=1}^{r} \boldsymbol{\chi\_{i}} - \sum\_{i=1}^{r} (\boldsymbol{\chi\_{i+}} \ast \boldsymbol{\chi\_{+i}})}{N^2 - \sum\_{i=1}^{r} (\boldsymbol{\chi\_{i+}} \ast \boldsymbol{\chi\_{+i}})} \tag{8}$$

where *r* represents the number of rows and columns in the error matrix, *N* is the total number of pixels, *xii* is the observation in row *i* and column *i*, *xi+* is the marginal total of row *i*, and *x+i* is the marginal total of column *i*.

#### **3. Results**

#### *3.1. Data Preprocessing*

The final dataset used as model input is shown in Table 2. It presents numbers of polygons, the area covered by the UGS classes within the MMA, and the proportions of area. The UGS pixel ratio was 3.003%. The data augmentation technique helped to improve the amount of supporting information for the training process by increasing the number of mosaics from 24,667 to 1,480,020 mosaics. With this increase, the CNN improved the learning process.

**Table 2.** Results of the homogenization of UGS databases.

Parks represented the most prominent type of UGS with 1.8% coverage of the MMA (Table 2). Median strips represented 0.84% cover. Residential gardens represented 0.34% cover. The classes with the lowest coverage were squares and roundabouts with 0.01% and 0.001%, respectively.

#### *3.2. Semantic Segmentation of UGSs*

The highest Dice coefficient and accuracy results of the semantic segmentation for each of the 12 three-band compositions are presented in Table 3. NDVI–red–NIR composition achieved the best results using ResNet34 encoder with a Dice coefficient value of 0.5748 and an accuracy of 0.9503. Red–green–blue composition achieved the best results using ResNet50 with a Dice coefficient of 0.4378 and an accuracy of 0.9839. In contrast, EVI2– NDWI–NIR composition had the lowest values for both encoders. For the ResNet34 encoder model, the mean Dice coefficient was 0.49, the standard deviation was 0.09, and the statistical significance using 95% confidence intervals ranged from 0.42 to 0.55. For the ResNet50 encoder model, the mean Dice coefficient was 0.42, the standard deviation was 0.58, and the statistical significance using 95% confidence intervals ranged from 0.28 to 0.36.

Figure 7 illustrates the behavior of the training and validation process for the highest Dice coefficient for both encoders. As observed, the ResNet34 learning process extended to the 100 epochs (fluctuating between 0.45 and 0.63) and presented its peak at the 83rd epoch, reaching a Dice coefficient of 0.5748. In contrast, the ResNet50 learning process activity occurred during the first 10 epochs, reaching the best Dice coefficient of 0.4378 on the 4th epoch.

**Table 3.** Semantic segmentation model validation results for UGSs in VHR satellite images.

**Figure 7.** Plot of training loss, validation loss, and Dice coefficient for both encoders.

The best segmentation process is represented by the lowest loss value. An example of this is shown in Figure 8A–D, where NDVI–red–NIR composition using ResNet34 reflects how the learning process increases from a loss of 1.14 to 0.77. This behavior is also observed for the RGB combination using ResNet50 where the loss was from 1.47 to 0.84 (Figure 8E–H).

IoU metric was 0.75 for the evaluation of the NDVI–red–NIR composition. This was calculated using the polygons presented in Figure 9A. The recall analysis revealed that the ground truth data had an overlap of 96.07% with the predicted data and the proportion of the overlapping polygons corresponding to the predicted data was 80.04% (Figure 9B).

Results of the confusion matrix and kappa coefficient produced for the evaluation dataset are shown in Table 4. Both the ground truth and the predicted data contained polygons corresponding to parks and median strips classes. The kappa coefficient was 0.94, and the overall accuracy calculated was 0.97. The user accuracy was 1 for the parks and 0.96 for median strips, and the producer accuracy was 0.92 for the parks and 1 for median strips.

**Figure 8.** Progression of the segmentation process of UGSs for both encoders. (**A**–**D**) Samples obtained using an NDVI–red–NIR composition from ResNet34 encoder. (**E**–**H**) Samples obtained using a red–green–blue composition from ResNet50 encoder. The information is presented at the same scale; each square surface is 1.63 hectares (ha).

**Figure 9.** Model assessment. (**A**) The overlap coverage between the polygons produced from the CNN model and the ground truth data. (**B**) Recall analysis. T is the total coverage. Ground truth vs predicted columns show the area in ha and proportion of positives identified correctly.

**Table 4.** Confusion matrix and kappa coefficient.

#### **4. Discussion**

The two automated methods tested for the identification of individual UGS polygons at the metropolitan level generated new databases that provide useful information, including geometry, condition, and spatial attributes, for the decision-making process regarding these important open public spaces. Updated databases improve national inventories with detailed geospatial and geometric information of UGSs, which is important for assessing their distribution and management. Updated and accurate UGS information, however, is difficult to acquire or access, especially in developing countries with no access to VHR imagery. Latin American countries lack this kind of information, where the UGS inventories are typically based on photo-interpretation techniques and depend on the user experience.

The information produced with this methodology can be used in conjunction with demographic information to analyze the accessibility and connectivity of UGSs at the metropolitan level. While this is a portion of the information needed for the quantification of the achievement of the SDG 11.7, the generation of a similar database considering nongreen public spaces should also be contemplated to cover the analysis of both elements.

The methodology of segmentation of UGS polygons at the city level proposed in this work will allow effectively updating this information for urban spaces in Mexico. The method includes the typical UGS classes, such as median strips, roundabouts, parks, squares, and residential gardens present in every metropolitan area. OSM and INEGI open access databases used in this research are available for all of Mexico. This information was complemented by information produced by the local university through a metropolitan project funded by the state government, proving the importance of integrating multilevel governance (or institutions) to enhance and update geospatial data such as UGS inventories.

Methods to increase data representation are a necessity when VHR imagery is limited. In this work, data augmentation techniques using simple strategies showed their effectiveness by providing over 1 million additional orthomosaics. The variations that occurred in each of these image transformations helped to reach a more complete training set representing the complexity occurring in the study area due to temporal or environmental disparities.

The prospected combinations produced by four-band VHR imagery and their implementation using two encoders allowed the assessment of 24 segmentation models. This kind of modeling is only possible with a high computation capability. When that is not available, other options such as deep learning processing cloud services (e.g., Google Colab, Amazon AWS, Microsoft Azure) can be implemented.

According to the semantic segmentation model validation results and its statistical analysis, there is a significant difference (95%) between the dice coefficient of the different band combinations in both models. The best three-band combination for the semantic segmentation of UGSs is NDVI–red–NIR when using ResNet34 encoder and red–green– blue when using ResNet50 encoder. A future analysis regarding the learning process could help to identify the learning patterns and the influence of each band within the models by using interpretability, representation learning, and visualization methods [84].

A large difference between the validation accuracy and the dice coefficient was observed in Table 3. This variance is associated with the data imbalance caused by background non-UGS pixels. As the non-UGS class has the highest number of pixels in the analyzed orthomosaics, the accuracy of the model is high as it is quantifying a high percentage of correctly classified pixels for the entire area. The dice coefficient is a more reliable

parameter in semantic segmentation processes with data imbalances because it reflects a metric based only on the segmented classes.

Semantic segmentation studies using similar approaches to map urban tree coverage, buildings, and roads [85–87] reported dice coefficients of 0.94, 0.84, and 0.87, respectively. The result in this study is lower (0.57), which may be related to the complexity of mapping UGS polygons. The referenced studies focus on the segmentation of classes that represent the city coverage; however, this research seeks to segment the pixel class and also several types of geometry. Additionally, UGS polygons are composed of a mix of pixels representing not only vegetation but also other kinds of infrastructure, such as sidewalks and playgrounds, which decrease the certainty for the segmentation process.

The kappa coefficient produced in the accuracy assessment of the evaluation data indicated a strong agreement between the predicted polygons and the ground truth data. With a value of 0.94, the kappa coefficient was similar to high accuracy results obtained in recent studies related to UGS mapping methods [88,89]. This indicates that the methodology used in this study is accurate for extracting and updating geometrical UGS databases at the metropolitan level.

#### **5. Conclusions**

This study evaluated two deep learning model techniques for semantic segmentation of UGS polygons with the use of different CNN encoders on the U-Net architecture to improve the methodology of UGS cartography. The models have the capability to detect patterns for all types of UGSs reported in Mexico, even with a high variation in shape or size, and to segment hundreds of thousands of polygons that represented 3% of the total MMA.

Results demonstrate that this methodology is an accurate digital tool for extracting and updating geometrical UGS databases at the metropolitan level (Dice coefficient of 0.57, recall of 0.8, IoU of 0.75, and kappa coefficient of 0.94). The implementation of these models could update UGS inventories necessary to assess urban management as cities grow or change. This methodology produces UGS geospatial data that are essential for quantifying the accomplishment of the SDG 11.7 regarding green spaces. This information in combination with demographic data could be used to elaborate UGS accessibility maps necessary to assess UGS accessibility. This new cartography may improve urban management for the conservation of natural resources and the environmental services they provide, as well as making their maps more accessible to urban residents and decision-makers.

**Author Contributions:** Conceptualization, F.D.Y.; data curation, R.E.H.; formal analysis, R.E.H. and F.D.Y.; funding acquisition, F.D.Y.; investigation, R.E.H. and F.D.Y.; methodology, R.H.G. and F.D.Y.; project administration, F.D.Y.; software, R.E.H.; supervision, F.D.Y., R.A.C.G., D.F.L.-G., V.H.G.C., A.L.F.F., and H.d.L.G.; validation, F.D.Y. and D.F.L.-G.; visualization, R.H.G. and F.D.Y.; writing original draft, R.H.G. and F.D.Y.; writing—review and editing, D.F.L.-G., V.H.G.C., A.L.F.F., H.d.L.G., R.A.C.G. and A.V.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Council of Science and Technology, grant number 559018 and grant for international research (Foreign Mobility 2019-1).

**Acknowledgments:** The authors extend their sincere thanks to Gustau Camps-Valls, researcher of the Image Processing Lab of the University of Valencia, who guided Roberto Huerta in neural network subjects through a stay in the Image Processing Laboratory. Map data are copyrighted by OpenStreetMap contributors and available from https://www.openstreetmap.org, (accessed on 26 November 2020). The authors thank William D. Eldridge for his support and comments.

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

#### **References**

## *Article* **Susceptibility Analysis of the Mt. Umyeon Landslide Area Using a Physical Slope Model and Probabilistic Method**

#### **Sunmin Lee 1,2 , Jungyoon Jang 3 , Yunjee Kim 2 , Namwook Cho 2 and Moung-Jin Lee 4,\***

Received: 7 July 2020; Accepted: 13 August 2020; Published: 18 August 2020

**Abstract:** Every year, many countries carry out landslide susceptibility analyses to establish and manage countermeasures and reduce the damage caused by landslides. Because increases in the areas of landslides lead to new landslides, there is a growing need for landslide prediction to reduce such damage. Among the various methods for landslide susceptibility analysis, statistical methods require information about the landslide occurrence point. Meanwhile, analysis based on physical slope models can estimate stability by considering the slope characteristics, which can be applied based on information about the locations of landslides. Therefore, in this study, a probabilistic method based on a physical slope model was developed to analyze landslide susceptibility. To this end, an infinite slope model was used as the physical slope model, and Monte Carlo simulation was applied based on landslide inventory including landslide locations, elevation, slope gradient, specific catchment area (SCA), soil thickness, unit weight, cohesion, friction angle, hydraulic conductivity, and rainfall intensity; deterministic analysis was also performed for the comparison. The Mt. Umyeon area, a representative case for urban landslides in South Korea where large scale human damage occurred in 2011, was selected for a case study. The landslide prediction rate and receiver operating characteristic (ROC) curve were used to estimate the prediction accuracy so that we could compare our approach to the deterministic analysis. The landslide prediction rate of the deterministic analysis was 81.55%; in the case of the Monte Carlo simulation, when the failure probabilities were set to 1%, 5%, and 10%, the landslide prediction rates were 95.15%, 91.26%, and 90.29%, respectively, which were higher than the rate of the deterministic analysis. Finally, according to the area under the curve of the ROC curve, the prediction accuracy of the probabilistic model was 73.32%, likely due to the variability and uncertainty in the input variables.

**Keywords:** probabilistic method; Monte Carlo simulation; physical slope model; Mt. Umyeon landslides

#### **1. Introduction**

Landslides cause substantial economic and social losses, especially in urban areas where many people live. Landslides are destructive and represent the most frequent risk factors in mountainous areas; especially in urban areas where damage to forests and infrastructure, such as buildings and

roads, can lead to soil erosion and thus land deformation. Therefore, analysis of mitigation risks and countermeasures for mitigation are essential steps in reducing natural disasters.

Geographic information system (GIS)-based landslide susceptibility analysis has been conducted to predict areas with high probabilities of landslides. Based on spatial data constructed from factors that influence landslide occurrence, such as topography, hydrology, forests, and geology, studies have been conducted to estimate landslide susceptibility [1–5]. Many traditional statistical methodologies, such as evident belief function, frequency ratio, analytical hierarchy process, and logistic regression have been applied [5–9]. Recently, machine learning methodologies have been applied to estimate the spatial uncertainty of landslides [1,4,10,11]. Typically, decision tree-based models [12], such as random forest and boosted tree models [13,14], have been applied. Support vector machines [15] and artificial neural networks [16,17] are other commonly applied machine learning methods for landslide susceptibility analysis. Various studies analyzed the relationship between landslide location data and the factors that cause landslides and calculate their effects on landslides. There is a disadvantage to this method because it is not possible to conduct landslide susceptibility analysis on target areas before landslides occur when there is no information on the locations of landslides. The landslide damage area in South Korea is expanding every year; thus, in addition to areas where landslides have already occurred, the frequency of new landslides is increasing [18]. Therefore, it is necessary to predict landslide susceptibilities and prepare countermeasures for areas without prior information.

Among the methods for evaluating landslide susceptibilities, physical-model based analysis estimates the stability by treating a slope as a specific physical model and inputting slope information [19,20]. This method enables susceptibility analysis regardless of the information about landslide occurrence location information, so it is possible to analyze the susceptibility of an area before landslides occur [21]. Therefore, this analysis method has the advantage of considering the occurrence mechanism and process of landslides and has been reported as one of the most effective techniques of landslide vulnerability and risk analysis [22]. In particular, among the physical slope models, the infinite slope model, which has a similar form of landslide fracture surface and is easier to analyze than other models, is most commonly used [23,24]. In addition, probabilistic techniques are used to effectively deal with the spatial variability in the geotechnical properties used as an input into physical slope models and inaccurate results due to complex geological conditions [25–33].

Landslide susceptibility analysis using physical models uses input data related to the topographical and geological characteristics of the slope [34]. In the process of obtaining the geotechnical characteristics, uncertainty occurs due to the spatial variability of the ground and complex geological conditions, which increases the possibility of obtaining incorrect analysis results. Therefore, probabilistic techniques such as Monte Carlo simulation have been used to quantify uncertainty [26,30,33]. In addition, studies are being conducted on various hydrogeological models capable of estimating pore water pressure due to rainfall infiltrating underground for the calculation of pore water pressure affecting stability by reducing shear strength of slope materials [35,36]

Thus, this study applied a physical model-based method to analyze landslide susceptibility before landslide occurrence. In this study, an infinite slope model was used as the physical slope model, and we used Monte Carlo (MC) simulation for the probabilistic analysis. Information on the actual occurrence of landslides in Mt. Umyeon was used to validate the accuracy of our techniques [37]. Figure 1 shows the detailed workflow used in this study. We used Mt. Umyeon as the study area because an urban landslide occurred here in July 2011. The Mt. Umyeon landslide is a representative example of serious human injury because it was located in the center of the metropolitan area of Seoul.

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 3 of 19

**Figure 1.** The workflow of this study. **Figure 1.** The workflow of this study.

#### **2. Study Area 2. Study Area**

The study area, Mt. Umyeon, is located in Seoul, the capital of the Republic of Korea, within 126°59′–127°01′ E, 37°27′–37°28′ N (Figure 2). Mt. Umyeon is located at the center of a densely populated area, with highways to the east, rivers and parks to the south, and major cultural facilities around the Mt. Umyeon area. The study area, which was selected based on the type of watershed affected by rainfall, is approximately 22.18 km2 (width: 5.075 km, length: 4.37 km). The maximum height of Mt. Umyeon is 293 m; the southern slope has a large slope and a valley, while the northern slope is gentle. Mt. Umyeon is a relatively low mountainous region of gneiss formed by retardation and weathering. The terrain is vulnerable to landslides because the gneiss in the bedrock is distributed with severe weathering and many faults. In addition, the dark veins are partially infiltrated, and overall, weathering is severe and the outcrop is poor [38]. The study area, Mt. Umyeon, is located in Seoul, the capital of the Republic of Korea, within 126◦590–127◦010E, 37◦270–37◦280N (Figure 2). Mt. Umyeon is located at the center of a densely populated area, with highways to the east, rivers and parks to the south, and major cultural facilities around the Mt. Umyeon area. The study area, which was selected based on the type of watershed affected by rainfall, is approximately 22.18 km2 (width: 5.075 km, length: 4.37 km). The maximum height of Mt. Umyeon is 293 m; the southern slope has a large slope and a valley, while the northern slope is gentle. Mt. Umyeon is a relatively low mountainous region of gneiss formed by retardation and weathering. The terrain is vulnerable to landslides because the gneiss in the bedrock is distributed with severe weathering and many faults. In addition, the dark veins are partially infiltrated, and overall, weathering is severe and the outcrop is poor [38].

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 4 of 19

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 4 of 19

**Figure 2.** Study area: (**a**) Seoul; (**b**) Mt. Umyeon area. **Figure 2.** Study area: (**a**) Seoul; (**b**) Mt. Umyeon area.

In this region, a number of landslides occurred due to heavy rains, with a cumulative rainfall of 587.5 mm for 3 days from 26 to 28 July, 2011 [37]. About 150 large and small landslides occurred [38], and the area of debris flow was very wide compared to about 11 square kilometers selected as the radius of the study area (Figure 3). Landslides were presumed to result from heavy rains over a period of about 1 h following the weakening of ground due to previous heavy rains of 230.0–266.5 mm from about 15 h before the landslide [38]; the estimated time of the landslide was 09:00 on 27 July, 2011. A number of landslides in the form of debris flow have been reported and landslide occurrence locations were collected based on field investigation and visual analysis of aerial photographs before and after the landslides [38] by points [3,5,13]. In this region, a number of landslides occurred due to heavy rains, with a cumulative rainfall of 587.5 mm for 3 days from 26 to 28 July, 2011 [37]. About 150 large and small landslides occurred [38], and the area of debris flow was very wide compared to about 11 square kilometers selected as the radius of the study area (Figure 3). Landslides were presumed to result from heavy rains over a period of about 1 h following the weakening of ground due to previous heavy rains of 230.0–266.5 mm from about 15 h before the landslide [38]; the estimated time of the landslide was 09:00 on 27 July, 2011. A number of landslides in the form of debris flow have been reported and landslide occurrence locations were collected based on field investigation and visual analysis of aerial photographs before and after the landslides [38] by points [3,5,13]. **Figure 2.** Study area: (**a**) Seoul; (**b**) Mt. Umyeon area. In this region, a number of landslides occurred due to heavy rains, with a cumulative rainfall of 587.5 mm for 3 days from 26 to 28 July, 2011 [37]. About 150 large and small landslides occurred [38], and the area of debris flow was very wide compared to about 11 square kilometers selected as the radius of the study area (Figure 3). Landslides were presumed to result from heavy rains over a period of about 1 h following the weakening of ground due to previous heavy rains of 230.0–266.5 mm from about 15 h before the landslide [38]; the estimated time of the landslide was 09:00 on 27 July, 2011. A number of landslides in the form of debris flow have been reported and landslide occurrence locations were collected based on field investigation and visual analysis of aerial photographs before and after the landslides [38] by points [3,5,13].

Table 1 shows the spatial datasets used in this study. To build the input dataset, all input data **Figure 3.** Landslide status on the northern side of Mt. Umyeon in 2011 [39,40]. **Figure 3.** Landslide status on the northern side of Mt. Umyeon in 2011 [39,40].

#### were constructed in a raster format with a 5 × 5 m grid structure using the ArcGIS 10.3 software. For the analysis of landslide susceptibilities, relevant factors were selected through literature review of **3. Spatial Datasets 3. Spatial Datasets**

previous studies [41–46]. First, a 1: 5000 topographic map was obtained, from which we collated data such as elevation, slope, and specific catchment area. A digital elevation model was constructed by extracting a contour vector layer, including elevation attributes, from a digital Table 1 shows the spatial datasets used in this study. To build the input dataset, all input data were constructed in a raster format with a 5 × 5 m grid structure using the ArcGIS 10.3 software. For the analysis of landslide susceptibilities, relevant factors were selected through literature review of previous studies [41–46]. First, a 1: 5000 topographic map was obtained, from which we collated data such as elevation, slope, and specific catchment area. A digital elevation model was constructed by extracting a contour vector layer, including elevation attributes, from a digital Table 1 shows the spatial datasets used in this study. To build the input dataset, all input data were constructed in a raster format with a 5 × 5 m grid structure using the ArcGIS 10.3 software. For the analysis of landslide susceptibilities, relevant factors were selected through literature review of previous studies [41–46]. First, a 1: 5000 topographic map was obtained, from which we collated data such as elevation, slope, and specific catchment area. A digital elevation model was constructed by extracting a contour vector layer, including elevation attributes, from a digital topographic map

(Figure 4a). Finally, the slope and specific catchment area (SCA) were constructed based on the digital elevation model (Figure 4b,c).

**Table 1.** Landslide-related factors used to construct spatial database.

a Aerial photograph before and after Mt.Umyeon landslides from Kakaomap [47]. b Topographical factors were extracted from digital topographic map by National Geographic Information Institute. c The detailed soil map produced by Rural Development Administration. d Field investigation data produced by Korean Geotechnical Society [38]. e The 16-h accumulated precipitation of from seven Automatic Weather System (AWS) observatories by Korea Meteorological Administration (KMA).

Second, the z-model was used to build the depth data of the fracture surface. The z-model is a submarine model that reflects the topographic characteristics of the slope by calculating its depth with respect to the elevation [48]. According to the results of ground drilling and seismic surveys by the Korean Society of Civil Engineers [37], the distribution of the soil thickness in the Mt. Umyeon area is 1.18–10.13 m. Therefore, the minimum and maximum depths of the soil in the z-model were set to 1 and 10 m, respectively, and soil thickness was calculated (Figure 4d) [37].

Third, unit weight, cohesion, friction angle, and hydraulic conductivity were obtained through field surveys with indoor experiments. Direct shear tests were conducted on four sampling points (SPs) on the obtained samples of cohesion and friction angle, and borehole shear tests and in-situ permeability tests were conducted at eight drilling points (DPs) in the study area (Figure 5) [38]. Table 2 summarizes the results obtained from this process. In order to interpolate and analyze geological characteristics of the entire study area based on the obtained data [32,41], kriging spatial interpolation analysis from ArcGIS, a spatial processing interpolation technique that reflects the correlation between the distance from the surrounding value and the value located around it, was performed to construct data in raster format (Figure 6).

**Table 2.** Geotechnical properties of unit weight, cohesion, friction angle, and hydraulic conductivity.

topographic map (Figure 4a). Finally, the slope and specific catchment area (SCA) were constructed

**Table 1.** Landslide-related factors used to construct spatial database.

Elevation [m] Slope gradient [˚] Specific catchment area (SCA)

Unit weight [kN/m3] Cohesion [kPa] Friction angle [degree] Hydraulic conductivity [cm/s]

Precipitation e Rainfall intensity Point a Aerial photograph before and after Mt.Umyeon landslides from Kakaomap [47]. b Topographical factors were extracted from digital topographic map by National Geographic Information Institute. c

Automatic Weather System (AWS) observatories by Korea Meteorological Administration (KMA).

**Data Source. Factors Data Type Scale**  Aerial photograph a Landslide location Point -

Soil map c Soil thickness [m] Polygon 1:25,000

GRID 1:5000

Point -

based on the digital elevation model (Figure 4b,c).

Topographical map b

Field Investigation d

**Figure 4.** Map showing the input parameters constructed using topographic maps: (**a**) elevation; (**b**) slope; (**c**) specific catchment area (SCA); (**d**) soil thickness. *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 7 of 19

**Figure 5.** Location map of sampling point for acquisition of geotechnical properties. **Figure 5.** Location map of sampling point for acquisition of geotechnical properties.

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 8 of 19