# **Artificial Intelligence Methods Applied to Urban Remote Sensing and GIS**

Edited by Chang-Wook Lee, Hyangsun Han, Hoonyol Lee and Yu-Chul Park Printed Edition of the Special Issue Published in *Remote Sensing*

www.mdpi.com/journal/remotesensing

## **Artificial Intelligence Methods Applied to Urban Remote Sensing and GIS**

## **Artificial Intelligence Methods Applied to Urban Remote Sensing and GIS**

Editors

**Chang-Wook Lee Hyangsun Han Hoonyol Lee Yu-Chul Park**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin


Yu-Chul Park Geophysics Kangwon National University Chuncheon Korea, South

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Remote Sensing* (ISSN 2072-4292) (available at: www.mdpi.com/journal/remotesensing/special issues/AI urban remote sensing).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1604-2 (Hbk) ISBN 978-3-0365-1603-5 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editors**

#### **Chang-Wook Lee**

Changwook Lee received his B.S. degree from Kangwon National University, Chuncheon, Korea, and M.S. and Ph.D. degrees from Yonsei University, Seoul, Korea in 1999 and 2009, respectively. He held a postdoctoral position in InSAR for the ARTS contract with the US Geological Survey EROS data center, and all work was performed as determined by the USGS project manager working in coordination with SSSC management. The work was performed at the USGS Cascades Volcano Observatory in Vancouver, Washington, by National Aeronautics and Space Administration (NASA) project supports from 2009 to 2011. He is an associate professor in the department of science of education, Kangwon National University, in Chuncheon, Korea. His research interests include SAR, InSAR and time-series processing technique development for natural disaster monitoring and resource characterization. He has authored more than 100 articles in these research fields.

#### **Hyangsun Han**

Hyansun Han received his B.S. and M.S. degrees in geophysics, and Ph.D. degree in remote sensing from the Department of Geophysics, Kangwon National University, South Korea, in 2006, 2008 and 2013, respectively. From March 2013 to September 2014, he was a postdoctoral researcher with the Research Institute for Earth Resource at Kangwon National University. From October 2014 to August 2015, he was a research professor with the Ulsan National Institute of Science and Technology, South Korea. He was a senior research scientist with the Korea Polar Research Institute (KOPRI), South Korea, from September 2015 to February 2020. He is an assistant professor in the Department of Geophysics, Kangwon National University. His research interests include satellite remote sensing of cryosphere, land and ocean, and various applications of remotely sensed data. He received Student Paper Awards at the Korea Remote Sensing Symposia in 2006 and 2007, Excellence Research Awards at Kangwon National University in 2013, and KOPRI's Scientist Awards in 2016, 2017 and 2018.

#### **Hoonyol Lee**

Hoonyol Lee received his B.S. degree in geology and M.S. degree in geophysics from the Department of Geological Sciences, Seoul National University, Seoul, South Korea, in 1995 and 1997, respectively, and his Ph.D. degree in radar remote sensing from the Department of Earth Sciences and Engineering, Imperial College London, University of London, London, U.K., in 2001. From 2001 to 2003, he was a postdoctoral research associate with Imperial College London. From 2003 to 2004, he was a senior researcher with the Korea Institute of Geoscience and Mineral Resources (KIGAM), Daejeon, South Korea. He has been with the Department of Geophysics, Kangwon National University since 2004. He was a Visiting Scholar with the Department of Geological Sciences, University of Oregon, Eugene, OR, USA, from August 2008 to July 2009, and with the Department of Civil, Environmental, and Construction Engineering, University of Central Florida, Orlando, FL, USA, from August 2014 to July 2015. He is currently the president of the Korean Society of Remote Sensing since 2021.

#### **Yu-Chul Park**

Yu-Chul Park received his Bachelor of Science degree in geology in 1992 from Seoul National University, his master of science degree in geophysics in 1995 from Seoul National University, and his doctor of science degree in mathematical geology from Purdue University in 2000, USA.

## **Preface to "Artificial Intelligence Methods Applied to Urban Remote Sensing and GIS"**

Recently, remote sensing and GIS techniques have gained increasing importance for rapid urbanization, the expansion of urban growth, and the enlargement of populations, due to the application of artificial intelligence, machine learning, and deep learning algorithms. This Special Issue aims to present the state-of-the-art research on optics, SAR, hyperspectral images, and GIS techniques for monitoring urban area environments corresponding to changes in times using publicly available and commercial datasets such as satellite and UAV data.

Given the information above, the aim of this Special Issue is to present the observed urban area and monitor the surrounding urban area in "Artificial Intelligence Methods Applied to Urban Remote Sensing and GIS". This research paper of *Remote Sensing* will cover a wide range of fields, including GISs, remote sensing, earth science, computer science, and environmental science, to analyze the urbanization phenomenon along with theoretical research and practical developments. Some of the prospective/encouraged topics for this Special Issue include:


#### **Chang-Wook Lee, Hyangsun Han, Hoonyol Lee, Yu-Chul Park**

*Editors*

### *Article* **Improvement of Earthquake Risk Awareness and Seismic Literacy of Korean Citizens through Earthquake Vulnerability Map from the 2017 Pohang Earthquake, South Korea**

**Ju Han <sup>1</sup> , Arip Syaripudin Nur <sup>2</sup> , Mutiara Syifa <sup>3</sup> , Minsu Ha <sup>3</sup> , Chang-Wook Lee 2,3 and Ki-Young Lee 3,\***


**Abstract:** Earthquake activities in and around the Korean Peninsula are relatively low in number and intensity compared with neighboring countries such as Japan and China. However, recent seismic activity caused great alarm and concern among citizens and government authorities, and uncovered the level of preparedness toward earthquake disasters. A survey has been conducted on 1256 participants to investigate the seismic literacy of Korean citizens, including seismic knowledge, awareness and management using a questionnaire of citizen earthquake literacy (CEL). The results declared that the citizens had low awareness and literacy, which means that they are not properly prepared for earthquake hazards. To develop an earthquake risk reduction plan and program efficiently and effectively, not only must it appropriately characterize the target audience, but also indicate high potential earthquake zones and potential earthquake damage. Therefore, this study mapped and analyzed the seismic vulnerability in southeast Korea using LogitBoost, logistic model tree (LMT), and logistic regression (LR) machine learning algorithms based on a building damage inventory map. The damaged buildings' locations were generated after the 2017 Pohang earthquake using the damage proxy map (DPM) method from the Sentinel-1 synthetic aperture radar (SAR) data. DPMs detected coherence loss, which indicates damaged buildings in urban areas in the Pohang earthquake and shows a good correlation with the Korea Meteorological Administration (KMA) report with modified Mercalli intensity (MMI) scale values of more than VII (seven). The damage locations were randomly divided into two datasets: 50% for training the vulnerability models and 50% for validating the models in terms of accuracy and reliability. Fifteen seismic-related factors were used to construct a model of each algorithm. Model validation based on the area under the receiver operating curve (AUC) was used to determine model accuracy. The AUC values of seismic vulnerability maps using the LogitBoost, LMT, and LR algorithms were 0.769, 0.851, and 0.749, respectively. We suggest that earthquake preparedness efforts should focus on reconstruction, retrofitting, renovation, and seismic education in areas with high seismic vulnerability in South Korea. The results of this study are expected to be beneficial for engineers and policymakers aiming at developing disaster risk reduction plans, policies, and programs due to future seismic activity in South Korea.

**Keywords:** seismic vulnerability map; DPM method; Sentinel-1; machine learning; seismic literacy

#### **1. Introduction**

Earthquake activities in and around the Korean Peninsula are relatively low in number and intensity compared with neighboring countries such as Japan and China, because it is located within the Eurasian intracontinental region [1]. However, seismographs often

**Citation:** Han, J.; Nur, A.S.; Syifa, M.; Ha, M.; Lee, C.-W.; Lee, K.-Y. Improvement of Earthquake Risk Awareness and Seismic Literacy of Korean Citizens through Earthquake Vulnerability Map from the 2017 Pohang Earthquake, South Korea. *Remote Sens.* **2021**, *13*, 1365. https:// doi.org/10.3390/rs13071365

Academic Editor: Peter Hofmann

Received: 4 March 2021 Accepted: 27 March 2021 Published: 2 April 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/).

record sudden occurrences of moderate earthquakes; historical documents show that several damaging earthquakes happened in the country [2], indicating that the Korean Peninsula is not completely safe from earthquake disasters.

On 15 November 2017, an M<sup>L</sup> 5.4 earthquake occurred in Pohang, South Korea at 05:29:31 UTC [3], causing widespread damage in and around the city [4]. The earthquake was the second largest to occur in the Korean Peninsula since earthquake monitoring was initiated by the Korea Meteorological Administration (KMA) in 1978 [4,5]. In terms of the magnitude, the Pohang earthquake was not larger than the Gyeongju earthquake. However, the damage of the Pohang earthquake was much more than that of the Gyeongju earthquake. The Pohang earthquake caused more than USD 75 M of indirect damage to over 57,000 structures and over USD 300 M of total economic impact, as estimated by the Bank of Korea, injured 135 residents, and displaced more than 1700 people into emergency housing [6]. Meanwhile, the Gyeongju earthquake resulted in approximately USD 9.5 M in damage to 5368 properties, and 23 injured people [7]. More than 100 heritage buildings and monuments sustained damage from the earthquakes [8]. Twenty-one kilometers from Pohang, there is the Gyeongju Historic Area that was registered as a UNESCO World Cultural Heritage Site in November 2000, an area that embodies the time-honored history and culture of Gyeongju, the ancient capital of the Silla Kingdom (57 BC–935 AD). Some damage was found in this area, such as to Dabotap Pagoda (dislocated banister), Cheomseongdae Observatory (shifted and tilted), and Gyeongju Gyochon Traditional Village (cracked walls) [9]. Therefore, we need to conduct research about it.

Moreover, satellite images could detect surface deformation after the earthquake in Pohang, and therefore, for the first time, surface deformation measurements for an earthquake in Korea historically. A radar interferometry image was taken by the Sentinel-1 satellite, highlighting a deformation of −5 to 5 cm (blue to red) that occurred near Pohang city center. This image was obtained from a two-pass interferogram created using GAMMA software [10]; in this process, two synthetic aperture radar (SAR) images (4 and 16 November 2017) were co-registered to form an interferometric pair, which were then cropped to the area of interest. Topographic interferometric synthetic aperture radar (InSAR) images were then produced through interferogram generation. These images were derived from global 1 arcsecond Shuttle Radar Topography Mission (SRTM) data (30 m resolution), followed by topographic phase removal and differential synthetic aperture radar interferometry (DInSAR) phase generation, leaving only the deformation phase image [11]. A phase unwrapping procedure was then applied to generate an unwrapped DInSAR image, for conversion into displacement values (cm). Considering the uniqueness of this earthquake, we choose Pohang earthquake as our study case.

Additionally, we also need to know citizens' risk awareness about earthquakes in Korea as they had never experienced or felt the degree of natural hazard caused by the Pohang earthquake directly before. However, many people were not aware of earthquakes risk from the survey data. A survey was conducted with 1256 Korean citizens during spring 2020. Figure 1 shows the survey results that only 6% of the participants were aware that where their live is "absolutely not at all" safe from earthquakes, 29% of the participants thought "hardly not", 42% thought "normal", 20% thought "to some extent", and 3% thought they were "absolutely" safe from earthquakes. These results are an important warning sign for regulators and authorities, given the recent earthquakes that caused great human and material losses.

*Remote Sens.* **2021**, *13*, x FOR PEER REVIEW 3 of 26

Yes, absolutely 3%

**Figure 1.** Chart of earthquake risk awareness of Korean citizens. **Figure 1.** Chart of earthquake risk awareness of Korean citizens.

**Figure 1.** Chart of earthquake risk awareness of Korean citizens. However, a small number of people in the earthquake area (close distance) were aware of the extreme fear and danger of earthquakes compared to the majority of people according to the survey results. Figure 2 shows that people who live near the epicenter of the Pohang earthquake have higher awareness of earthquakes because they were directly affected by the damage from the Pohang earthquake. Meanwhile, people who live far from the epicenter, such as people living in Seoul, have a lower awareness level of earthquakes because they only felt a slight shock from the Pohang earthquake, which was not fatal. Therefore, we need to inform many people of the dangers caused by such an earth-However, a small number of people in the earthquake area (close distance) were aware of the extreme fear and danger of earthquakes compared to the majority of people according to the survey results. Figure 2 shows that people who live near the epicenter of the Pohang earthquake have higher awareness of earthquakes because they were directly affected by the damage from the Pohang earthquake. Meanwhile, people who live far from the epicenter, such as people living in Seoul, have a lower awareness level of earthquakes because they only felt a slight shock from the Pohang earthquake, which was not fatal. Therefore, we need to inform many people of the dangers caused by such an earthquake, and when earthquakes occur in another area in the future, we need to be able to recognize the earthquake and respond to it. However, a small number of people in the earthquake area (close distance) were aware of the extreme fear and danger of earthquakes compared to the majority of people according to the survey results. Figure 2 shows that people who live near the epicenter of the Pohang earthquake have higher awareness of earthquakes because they were directly affected by the damage from the Pohang earthquake. Meanwhile, people who live far from the epicenter, such as people living in Seoul, have a lower awareness level of earthquakes because they only felt a slight shock from the Pohang earthquake, which was not fatal. Therefore, we need to inform many people of the dangers caused by such an earthquake, and when earthquakes occur in another area in the future, we need to be able to recognize the earthquake and respond to it.

No, not at all 6%

> I think the place where I live is safe from earthquakes

quake, and when earthquakes occur in another area in the future, we need to be able to

**Figure 2.** Map of citizens' awareness risk awareness score throughout Korea. **Figure 2.** Map of citizens' awareness risk awareness score throughout Korea.

**Figure 2.** Map of citizens' awareness risk awareness score throughout Korea. To reduce the effects of earthquake disasters and develop an earthquake risk reduction plan and program efficiently and effectively, not only is performing sustainable preparation, such as seismic literacy of citizens, necessary, but so is indicating high potential To reduce the effects of earthquake disasters and develop an earthquake risk reduction plan and program efficiently and effectively, not only is performing sustainable preparation, such as seismic literacy of citizens, necessary, but so is indicating high potential To reduce the effects of earthquake disasters and develop an earthquake risk reduction plan and program efficiently and effectively, not only is performing sustainable preparation, such as seismic literacy of citizens, necessary, but so is indicating high potential earthquake zones and potential earthquake damage by producing seismic vulnerability maps. Seismic vulnerability assessment involves the comprehensive evaluation of factors that affect risks

associated with earthquakes within predefined areas. Urban areas are at higher risk of seismic disasters than outlying areas due to their higher building and infrastructure density and larger population. Therefore, in assessing seismic vulnerability, it is essential to select suitable influential factors and methods for the area of interest. Several methodologies have been applied for seismic vulnerability assessment and mapping during the past few decades [7,12–15].

Seismic vulnerability assessment studies commonly analyze case studies using a combination of multicriteria decision making (MCDM) and geographic information system (GIS) approaches [12,13,16,17]. Among these, the analytical hierarchy process (AHP) is one of the most widely known MCDM methodologies; it stratifies and quantifies the importance of each applied influential factor to determine its relative importance, and assesses vulnerability by applying weights to all factors [12,16,18]. However, this method can be subjective because the opinion of the researcher can affect the weight assignment process; therefore, it is somewhat unsuitable for objective assessment. To address this problem, recent studies have applied hybrid models that combine various methodologies [14,19].

Many recent studies related to seismic vulnerability assessment and mapping have been conducted using machine learning techniques [12,16,18–21]. For example, Han et al. (2019) [20] used a logistic regression (LR) model and applied the support vector machine (SVM) methodology to four kernel models (linear, polynomial, radial basis function, and sigmoid) to derive a suitable model for seismic vulnerability assessment; this study was notable in that the results of several seismic vulnerability models were compared analytically; such analyses are rarely conducted in this field, despite the broad application of machine learning techniques in recent years.

Providing training data plays an important role in the accuracy of the vulnerability map. Here, we used the damage proxy map (DPM) method to extract a building damage map for training and testing datasets. The DPM method is part of an ongoing collaborative effort between the Jet Propulsion Laboratory (JPL) and the California Institute of Technology, called the Advanced Rapid Imaging and Analysis (ARIA) project. The DPM method, using synthetic aperture radar (SAR) satellite data, has been shown to be useful for damage mapping following an earthquake and other natural disaster events, including the 2015 MW 7.8 Gorkha, Nepal earthquake using COSMO-SkyMed and ALOS-2 satellites [20], the 2019 typhoon Hagibis, Japan using Sentinel-1 satellites [21], the 2019 MW 7.1 Ridgecrest earthquake in California using Sentinel-1 satellites [22], and the 2014 eruption of Kelud volcano (Indonesia) using COSMO SkyMED satellites [23].

This study aims to improve the risk awareness and seismic literacy of Korean citizens through an earthquake vulnerability map of all buildings in southeast Korea. To produce the earthquake vulnerability map, we generated a damage proxy map (DPM) after the 2017 Pohang earthquake from the Sentinel-1 SAR dataset as a dependent variable, then applied machine learning to construct models using 15 seismic-related factors as independent variables. Model performances were verified using a receiver operating characteristic (ROC) curve. Finally, dangerous and safe areas were identified in southeast Korea by creating maps based on the model with the highest accuracy for each methodology, and the result were assessed. The results of this study should improve citizen earthquake risk awareness and seismic literacy, especially in high seismic vulnerability areas, and facilitate the construction of seismic vulnerability models that will be useful to reduce future losses due to earthquakes in South Korea.

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

#### *2.1. Study Area*

The study area covers three metropolitan cities; Busan, Daegu, and Ulsan, and two provinces; Geyongsangbuk and Gyeongsangnam, surrounding the epicenter of the Pohang earthquake. For simplicity, we refer to it as southeast Korea. The blue line in Figure 3 shows the border of the study area. In total, southeast Korea is home to 12,961,687 people and has an area of 32,285 km<sup>2</sup> [24]. Within the total area, urban areas account for 5.49%, followed

by agriculture at 21.88%, forestry (66.71%), and other areas (5.92%). The proportions of males and females in these areas are 50.21% and 49.79%, respectively [25]. occurred in southeast Korea [30]. Among the 88 earthquakes of magnitude 2.0 or higher, 23 (26.17%) occurred in the same area.

The study area covers three metropolitan cities; Busan, Daegu, and Ulsan, and two provinces; Geyongsangbuk and Gyeongsangnam, surrounding the epicenter of the Pohang earthquake. For simplicity, we refer to it as southeast Korea. The blue line in Figure 3 shows the border of the study area. In total, southeast Korea is home to 12,961,687 people and has an area of 32,285 km2 [24]. Within the total area, urban areas account for 5.49%, followed by agriculture at 21.88%, forestry (66.71%), and other areas (5.92%). The proportions of males and females in these areas are 50.21% and 49.79%, respectively [25].

The study area was chosen considering the high seismic activity in these areas. A series of aftershocks of the Pohang earthquake was observed to have occurred with a magnitude of 3.0 or more until 31 May 2018. A magnitude of 4.3 occurred near the epicenter, 2 h after the mainshock, and a magnitude 4.6 earthquake occurred 4 km to the southwest at around 20:03:04 UTC on 11 February 2018 [26]. One year prior to the Pohang earthquake, an ML 5.8 earthquake shocked Gyeongju at 11:32:55 UTC on 12 September 2016. The epicenter was 8.7 km from the south of the city and 15 km beneath the surface [27]. The Gyeongju earthquake was just 40 km from the site of the Pohang earthquake. The Gyeongju earthquake was accompanied by 600 aftershocks, including an ML 5.1 foreshock that occurred near the mainshock at 10:44:32 UTC [28]. In 2018, 115 earthquakes with magnitudes of more than 2 occurred in the Korean Peninsula; among these, 36 earthquakes (31.13%) occurred in southeast Korea [29]. In 2019, 957 earthquakes with a magnitude of less than 2.0 occurred in the Korean Peninsula; among these, 294 earthquakes (30.72%)

*Remote Sens.* **2021**, *13*, x FOR PEER REVIEW 5 of 26

**2. Materials and Methods** 

*2.1. Study Area* 

**Figure 3.** Study area of this study, indicated by a blue line. Location of the 2017 Pohang earthquake and surface deformation generated from Sentinel-1 synthetic aperture radar (SAR) data acquired on 4 and 16 November 2017. **Figure 3.** Study area of this study, indicated by a blue line. Location of the 2017 Pohang earthquake and surface deformation generated from Sentinel-1 synthetic aperture radar (SAR) data acquired on 4 and 16 November 2017.

The Korean Peninsula lies at the eastern margin of the Eurasian Plate. About 30–15 million years ago, north-northeast (NNE)-striking strike–slip faults and NNE- to NE-striking normal faults settled predominantly in southeastern Korea and adjacent offshore areas The study area was chosen considering the high seismic activity in these areas. A series of aftershocks of the Pohang earthquake was observed to have occurred with a magnitude of 3.0 or more until 31 May 2018. A magnitude of 4.3 occurred near the epicenter, 2 h after the mainshock, and a magnitude 4.6 earthquake occurred 4 km to the southwest at around 20:03:04 UTC on 11 February 2018 [26]. One year prior to the Pohang earthquake, an M<sup>L</sup> 5.8 earthquake shocked Gyeongju at 11:32:55 UTC on 12 September 2016. The epicenter was 8.7 km from the south of the city and 15 km beneath the surface [27]. The Gyeongju earthquake was just 40 km from the site of the Pohang earthquake. The Gyeongju earthquake was accompanied by 600 aftershocks, including an M<sup>L</sup> 5.1 foreshock that occurred near the mainshock at 10:44:32 UTC [28]. In 2018, 115 earthquakes with magnitudes of more than 2 occurred in the Korean Peninsula; among these, 36 earthquakes (31.13%) occurred in southeast Korea [29]. In 2019, 957 earthquakes with a magnitude of less than 2.0 occurred in the Korean Peninsula; among these, 294 earthquakes (30.72%) occurred in southeast Korea [30]. Among the 88 earthquakes of magnitude 2.0 or higher, 23 (26.17%) occurred in the same area.

The Korean Peninsula lies at the eastern margin of the Eurasian Plate. About 30–15 million years ago, north-northeast (NNE)-striking strike–slip faults and NNE- to NE-striking normal faults settled predominantly in southeastern Korea and adjacent offshore areas when the East Sea opened in the early to middle Tertiary as a back-arc basin; smaller-scale coetaneous basins also formed, including the Pohang Basin. Although this region is about 400–500 km in length, its seismicity is affected by complex interactions of the Indo-Australian and Eurasian plates, as well as by subduction of Philippine Sea plates beneath the Japan and Ryuku trenches [31]. Several faults are distributed within the study area, including Dongrae, Moryang, Miryang, Ulsan, Wangsan, and Yangsan [32]. Due to seismic history and geographic characteristics, the probability of earthquake occurrence in southeast Korea is considered relatively high, and secondary damage in the event of an earthquake with a medium or higher magnitude constitutes

an unusually high risk. Therefore, sustainable preparation and management planning for such events is required.

#### *2.2. SAR Datasets*

A building damage inventory map for producing the seismic vulnerability map in southeast Korea was generated using Sentinel-1 synthetic aperture radar (SAR) C-band data (5.5 cm wavelength) provided by the European Space Agency (ESA). Pohang is located in between two frames (470 and 475), path 61 of the Sentinel-1 imagery; therefore, we obtained six Sentinel-1 single look complex (SLC) images with vertical transmission and vertical return (VV) polarization for the Pohang earthquake, with four scenes prior to the event on 23 October 2017, and 4 November 2017, and two scenes after the event on 16 November 2017. Furthermore, we needed to merge the scenes before processing. The images were co-registered, with the 4 November 2017 scenes as a reference.

#### *2.3. Damage Proxy Map (DPM)*

DPMs generated from the comparison of pre- and co-event SAR images can help identify damage caused by earthquakes using remote sensing imagery [20]. The method relies on the reduction in the coherence of the radar echoes between satellite-based SAR images taken before and after the earthquake to identify anomalous changes in ground surface properties. Coherence measures the change in radar backscatter from the ground, a proxy for the ground-surface property changes. A low coherence implies a large change to the ground surface that reflected the SAR radiation [33]. Changes can be caused by damage to the ground itself or damage to structures.

The process started with image co-registration, with the 4 November 2017 scenes as a reference. This co-registration process was done with sub-pixel accuracy to match scenes to one another. We used the complex pixel value, c, of the pre-processed SLCs for the change detection analysis, where damage is inferred from loss of coherence or decorrelation between SAR images [20]. We computed the pre—and co-event interferometric coherences, γ (Equation (1)), from a pair of SLCs before the event and another pair spanning across the event, respectively [21].

$$\gamma = \frac{|\langle \mathbf{c\_1c\_2^\*} \rangle|}{\sqrt{\langle \mathbf{c\_1c\_1^\*} \rangle \langle \mathbf{c\_1c\_2^\*} \rangle}}, 0 \le \gamma \le 1 \tag{1}$$

where c<sup>1</sup> and c<sup>2</sup> are complex pixel values of two co-registered SAR images and \* denotes the complex conjugate. The resulting coherence ranges from 0 (incoherent) to 1 (coherent). The coherence is equal to 1 if the observation is identical in the two images because of the stable object-like buildings in the scene. The pre-event coherence represented change unrelated to the event and was assumed to be the background value. Then, we obtained a coherence difference (COD) by subtracting γpre-event from γco-event. Therefore, the process could generate a COD ranging from −1 to 1. A negative COD (or coherence gain) usually indicates surface changes occurring between the pre-event scenes and is associated with changes not related to the event. Coherence gain could happen in an agricultural area when the fields are full of crops. Then, when harvested during the period time of the pre-event interferometric pair, the area has low coherence. After harvesting, leaving an empty field, the area has a greater coherence in the co-event interferometric pair (γco-event > γpre-event), so a negative COD is obtained. A positive COD (or coherence loss) indicates surface changes between the co-event scenes spanning the events, such as major damage to a building significantly that increases the interferometric phase variance, causing a decrease in coherence. Hence, the loss of coherence is most effective for detecting damage in built-up areas caused by earthquake. However, COD is generally less effective and less reliable in vegetated areas where coherence changes may be random. Therefore, this study focused on a DPM in urban, built-up areas as changes can be detected easily in SAR imagery. A greater loss in coherence generally correlates with greater severity of the change, such as a fully collapsed building, causing more significant coherence loss than partial collapse [21].

The threshold for significant coherence loss can be chosen by comparing observed coherence changes with reported damage and areas in which it is known that no damage occurred. Yun et al. (2015) compared a DPM with a National Geospatial-Intelligence Agency (NGA) analysis and the United Nations Operational Satellite Application Programme (UNOSAT) damage assessment map. Tay et al. (2020) used high-resolution aerial imagery from the Geospatial Information Authority of Japan (GSI). Here, coherence loss thresholds for DPMs were chosen by considering the damaged area from the Korea Meteorological Administration (KMA) report of the Pohang earthquake [34]. The threshold for significant coherence loss can be chosen by comparing observed coherence changes with reported damage and areas in which it is known that no damage occurred. Yun et al. (2015) compared a DPM with a National Geospatial-Intelligence Agency (NGA) analysis and the United Nations Operational Satellite Application Programme (UNOSAT) damage assessment map. Tay et al. (2020) used high-resolution aerial imagery from the Geospatial Information Authority of Japan (GSI). Here, coherence loss thresholds for DPMs were chosen by considering the damaged area from the Korea Meteorological Administration (KMA) report of the Pohang earthquake [34].

event), so a negative COD is obtained. A positive COD (or coherence loss) indicates surface changes between the co-event scenes spanning the events, such as major damage to a building significantly that increases the interferometric phase variance, causing a decrease in coherence. Hence, the loss of coherence is most effective for detecting damage in builtup areas caused by earthquake. However, COD is generally less effective and less reliable in vegetated areas where coherence changes may be random. Therefore, this study focused on a DPM in urban, built-up areas as changes can be detected easily in SAR imagery. A greater loss in coherence generally correlates with greater severity of the change, such as a fully collapsed building, causing more significant coherence loss than partial collapse

#### *2.4. Selection of Seismic-Related Factors 2.4. Selection of Seismic-Related Factors*

[21].

*Remote Sens.* **2021**, *13*, x FOR PEER REVIEW 7 of 26

We selected factors affecting seismic vulnerability based on the result of previous studies [7,14,15,17]. The factors affecting the seismic vulnerability were prepared based on five main indicators, which were geotechnical, physical, structural, social, and capacity; we selected a total of 15 factors corresponding to these categories. Geotechnical factors included slope and altitude; physical factors included peak ground acceleration (PGA), epicenter distance, and fault distance; structural factors included land use, construction materials, building density, and building height; social factors included elderly population (≥65 years), child population (<15 years), and population density; and capacity factors included distances from hospitals, fire stations, and police stations. The factors were organized into raster-based spatial databases (30 m spatial resolution) and were reclassified using the quantile method to identify and analyze the effect of each class. The data used in this present study are shown in Figure 4. We selected factors affecting seismic vulnerability based on the result of previous studies [7,14,15,17]. The factors affecting the seismic vulnerability were prepared based on five main indicators, which were geotechnical, physical, structural, social, and capacity; we selected a total of 15 factors corresponding to these categories. Geotechnical factors included slope and altitude; physical factors included peak ground acceleration (PGA), epicenter distance, and fault distance; structural factors included land use, construction materials, building density, and building height; social factors included elderly population (≥65 years), child population (<15 years), and population density; and capacity factors included distances from hospitals, fire stations, and police stations. The factors were organized into raster-based spatial databases (30 m spatial resolution) and were reclassified using the quantile method to identify and analyze the effect of each class. The data used in this present study are shown in Figure 4.

**Figure 4.** *Cont*.

**Figure 4.** *Cont*.

**Figure 4.** Factors related to seismic vulnerability: (**a**) slope, (**b**) elevation, (**c**) peak ground acceleration (PGA), (**d**) distance from epicenters, (**e**) distance from faults, (**f**) land use, (**g**) construction materials, (**h**) density of buildings, (**i**) building height, (**j**) elderly population, (**k**) child population, (**l**) population density, (**m**) distance from hospitals, (**n**) distance from fire stations, and (**o**) distance from police stations. **Figure 4.** Factors related to seismic vulnerability: (**a**) slope, (**b**) elevation, (**c**) peak ground acceleration (PGA), (**d**) distance from epicenters, (**e**) distance from faults, (**f**) land use, (**g**) construction materials, (**h**) density of buildings, (**i**) building height, (**j**) elderly population, (**k**) child population, (**l**) population density, (**m**) distance from hospitals, (**n**) distance from fire stations, and (**o**) distance from police stations.

Slope and elevation were extracted from a digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) using basic terrain analysis tools. Slope and elevation are factors affecting the vulnerability of urban environments to earthquakes [14,35]. Degradation in terrain with steep topography, especially at the top of hills and peaks, is greatly enhanced. According to construction standards, a slope of 5 to 9% is suitable for urbanization [16]. The highs and lows of elevation in each area are highly correlated with landslide susceptibility in each area [36]. Therefore, because of the amount of erosion and its relation to human activity, the higher the altitude of an area, the greater the seismic vulnerability. Figure 4a,b shows the slope and elevation maps of southeast Korea, respec-Slope and elevation were extracted from a digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) using basic terrain analysis tools. Slope and elevation are factors affecting the vulnerability of urban environments to earthquakes [14,35]. Degradation in terrain with steep topography, especially at the top of hills and peaks, is greatly enhanced. According to construction standards, a slope of 5 to 9% is suitable for urbanization [16]. The highs and lows of elevation in each area are highly correlated with landslide susceptibility in each area [36]. Therefore, because of the amount of erosion and its relation to human activity, the higher the altitude of an area, the greater the seismic vulnerability. Figure 4a,b shows the slope and elevation maps of southeast Korea, respectively.

tively. Peak ground acceleration (PGA) (Figure 4c) is the degree to which the ground shakes at the Earth's surface and is related to the amount of fault activity [18]. In this study, raw data from the Korea Institute of Geoscience and Mineral Resources (KIGAM) were converted to acceleration data and interpolated throughout southeast Korea [27,37]. The epicenter, or location where an earthquake occurs, is an important factor related to earthquake occurrence; the level of damage is different depending on the ground condition or the structure of the fault plane, on which the greatest damage often occurs at the epicenter. Therefore, we used distance data from earthquake epicenters (Figure 4d) with a magnitude over 4 from 1978 to 2020, including the 2016 Gyeongju earthquake and the 2017 Pohang earthquake. The epicenter locations were acquired from the United States Geological Survey (USGS). The faults are forms of tectonic factors whose presence or absence can be examined in relation to the seismic hazard of different areas. Fault distance (Figure 4e) plays a key role in vulnerability to earthquake hazards, as proximity to the fault causes high seismic risk and damage, and distance from it will reduce the risk and consequently Peak ground acceleration (PGA) (Figure 4c) is the degree to which the ground shakes at the Earth's surface and is related to the amount of fault activity [18]. In this study, raw data from the Korea Institute of Geoscience and Mineral Resources (KIGAM) were converted to acceleration data and interpolated throughout southeast Korea [27,37]. The epicenter, or location where an earthquake occurs, is an important factor related to earthquake occurrence; the level of damage is different depending on the ground condition or the structure of the fault plane, on which the greatest damage often occurs at the epicenter. Therefore, we used distance data from earthquake epicenters (Figure 4d) with a magnitudeover 4 from 1978 to 2020, including the 2016 Gyeongju earthquake and the 2017 Pohang earthquake. The epicenter locations were acquired from the United States Geological Survey (USGS). The faults are forms of tectonic factors whose presence or absence can be examined in relation to the seismic hazard of different areas. Fault distance (Figure 4e) plays a key role in vulnerability to earthquake hazards, as proximity to the fault causes high seismic risk and damage, and distance from it will reduce the risk and consequently provide higher resilience [16].

provide higher resilience [16]. Anti-seismic design in South Korea was introduced in 1988, and it is mandatory only for buildings that are three stories or higher [38]. As of November 2016, 29.9% of residential buildings and 23.7% of non-residential buildings in Seoul were designed to be antiseismic. As there is no guarantee that future earthquakes will not exceed the magnitude of the Gyeongju earthquake, most buildings in South Korea are considered highly vulnerable. To assess their vulnerability, we identified four structural factors of seismic vulnerability: land use, construction materials, building density, and building height, depicted Anti-seismic design in South Korea was introduced in 1988, and it is mandatory only for buildings that are three stories or higher [38]. As of November 2016, 29.9% of residential buildings and 23.7% of non-residential buildings in Seoul were designed to be anti-seismic. As there is no guarantee that future earthquakes will not exceed the magnitude of the Gyeongju earthquake, most buildings in South Korea are considered highly vulnerable. To assess their vulnerability, we identified four structural factors of seismic vulnerability: land use, construction materials, building density, and building height, depicted in Figure 4f–i, respectively [7,17]. The greater the number of floors of a building, despite its quality, the

greater the vulnerability. The number of floors in the building, if not in accordance with safety principles, will definitely increase the damage [39]. Even if the height is treated with due diligence and calculations, it is difficult for the evacuation of buildings, and to search for and rescue people. Obviously, structures of high strength and standard materials have good earthquake safety [40]. Proper deployment of land uses on the basis of urban planning principles, such as proper accessibility, proper distance from the biological hotspots, safety, comfort, and utility can substantially reduce the amount of vulnerability, injury, and economic damage [41].

Increasing population growth, population density, and poor distribution of services and infrastructure pose risks to society [42]. In recent earthquakes around the world, it can be said that most of the damage is to humans and with the increase in population, it is predicted that in the future, the mortality rate will be higher. The earthquake hazard coefficient in urban centers is also more complex and riskier due to urbanization without planning and development [43]. In events such as earthquakes, everyone in the community is vulnerable, but older people and children are the most vulnerable groups in a community and more attention is needed to minimize pain and injury [44]. Children do not tolerate disruption well and older people are psychologically fragile because of their disrupted life rhythms. The elderly population, child population, and population density in southeast Korea are presented in Figure 4j–l, respectively.

We identified the locations of social infrastructure facilities that offer aid in the event of an earthquake, and of hazardous facilities that have the potential to cause huge damage. The degree of accessibility following a disaster was analyzed by considering the physical distances to three factors, including social infrastructure facilities (hospital, police station, and fire station). Distance from a hospital (Figure 4m) and access to health services such as hospitals play a key role in controlling post-emergency complications and providing earthquake rescue and hospitalization services. Proper and quick access to medical facilities will increase earthquake resilience [45]. Distance from a fire station (Figure 4n) and access to police stations (Figure 4o) through the communication networks will speed up rescue operations and service the injured. As such, the greater the distance from fire stations and police stations, the greater the vulnerability [46].

#### *2.5. Machine Learning*

To map seismic vulnerability using a machine learning algorithm, several steps must be performed. First, the spatial relationships between the damaged buildings from the DPM and related factors (geotechnical, physical, structural, social, and capacity indicators) are calculated using the frequency ratio (FR) method. In this research, we analyzed the spatial relationship between damaged buildings' locations (1623 cells) and 15 factors related to seismic vulnerability, based on the FR value of each factor. When the ratio is greater than 1, this denotes that the class in each factor has a closer relationship with seismic vulnerability [47]. In the case that each factor has a less close relationship with seismic vulnerability, the ratio is less than 1. The FR for each factor was calculated by Equation (2) [48].

$$\text{FR} = \frac{\% \text{ of class of related factor}}{\% \text{ of total area}} \tag{2}$$

To apply the machine learning algorithm, we used the 1623 cells of damaged buildings generated from the DPM method. Among these cells, 50% (812) were used as a training dataset and 50% (811) were used as a test dataset. We extracted the same number of cells corresponding to undamaged buildings. All cells were randomly sampled and generating models and the accuracy of each model was done based on training (1624) and test datasets (1622). Several seismic-related maps, including geological maps, were produced at a 30 m resolution. Then, all data were classified as categorical or continuous. Continuous variables include the slope, elevation, PGA, distance from epicenters, distance from faults, building density, building height, child population, elderly population, population density, distance in Figure 5.

from hospitals, distance from police stations, and distance from fire stations. Categorical variables include construction materials and land use. faults, building density, building height, child population, elderly population, population density, distance from hospitals, distance from police stations, and distance from fire stations. Categorical variables include construction materials and land use.

Model validation was carried out through ROC curve analysis of the testing dataset (50%). Receiver operating characteristic (ROC) curve analysis, as an index of model performance, is commonly used to assess predictive accuracy [49]. To quantitatively determine the accuracy of the model verification, the area under the curve (AUC) of the ROC curve is calculated for the total area and correct predictive accuracy is obtained. AUC values range between 0.5 and 1; higher values indicate more reliable algorithm performance. The workflow of the seismic vulnerability mapping carried out in this study is provided in Figure 5. Model validation was carried out through ROC curve analysis of the testing dataset (50%). Receiver operating characteristic (ROC) curve analysis, as an index of model performance, is commonly used to assess predictive accuracy [49]. To quantitatively determine the accuracy of the model verification, the area under the curve (AUC) of the ROC curve is calculated for the total area and correct predictive accuracy is obtained. AUC values range between 0.5 and 1; higher values indicate more reliable algorithm performance. The workflow of the seismic vulnerability mapping carried out in this study is provided

**Figure 5.** Workflow of seismic vulnerability mapping. Abbreviations: DPM, damage proxy map; InSAR, interferometric synthetic aperture radar; ROC, receiver operating characteristic. **Figure 5.** Workflow of seismic vulnerability mapping. Abbreviations: DPM, damage proxy map; InSAR, interferometric synthetic aperture radar; ROC, receiver operating characteristic.

ୈ

୧ୀଵ

#### 2.5.1. LogitBoost 2.5.1. LogitBoost

P ൬<sup>C</sup> x

LogitBoost is a boosting algorithm developed by Friedman et al. [50] to reduce bias and variance. The LogitBoost algorithm was modified from AdaBoost, which was the commonly used boosting method for handling noisy data that executes additive logistic regression with least-square fits for individual classes. LogitBoost reduces training errors and enhances classification accuracy by using additive logistic regression for classification with a base-learning regression scheme and an ability to perform multiclass classification. The damaged building inventory map was divided into two classes: damaged buildings and undamaged buildings, using Equation (3): LogitBoost is a boosting algorithm developed by Friedman et al. [50] to reduce bias and variance. The LogitBoost algorithm was modified from AdaBoost, which was the commonly used boosting method for handling noisy data that executes additive logistic regression with least-square fits for individual classes. LogitBoost reduces training errors and enhances classification accuracy by using additive logistic regression for classification with a base-learning regression scheme and an ability to perform multiclass classification. The damaged building inventory map was divided into two classes: damaged buildings and undamaged buildings, using Equation (3):

$$\mathbf{Lc(c)} = \sum\_{i=1}^{D} \mathfrak{f}\_i \mathbf{x}\_i + \mathfrak{f}\_0 \tag{3}$$

where D is the number of building damage-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 with Equation (4): where D is the number of building damage-dependent factors and β<sup>i</sup> is the coefficient of the i-th component within input vector x. Probabilities were constructed using the linear logistic regression method with Equation (4):

$$\mathbb{P}\left(\frac{\mathbb{C}}{\mathbf{x}}\right) = \exp(\mathrm{Lc}(\mathbf{x})) / \sum\_{\mathbb{C}'=1}^{\mathbb{C}} \exp\left(\mathrm{Lc}'(\mathbf{x})\right) \tag{4}$$

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

#### 2.5.2. Logistic Model Tree (LMT)

The logistic model tree combines the C4.5 algorithm [51] and logistic regression (LR) functions. The information gain ratio (IGR) technique is applied to split the tree into nodes and leaves, and the LogitBoost algorithm [52] is used to fit the logistic regression functions at a tree node. The C4.5 algorithm uses the entropy technique for feature selection because it is the fastest method for providing reliable classification accuracy [53]. The over-fitting problem, which is an important challenge in LMT modeling, is overcome using the CART algorithm, which prunes the tree for modeling the training dataset [54]. The IGR can be formulated using Equation (5):

$$\text{Gain ratio (A)} = \frac{\text{gain (A)}}{\text{split info (A)}} \tag{5}$$

where gain (A) is the information after attribute A is selected as a test for classification of the training samples and split info (A) is the information generated when x training samples are categorized into n subsets [51]. In the next step, the LogitBoost algorithm performs additive logistic regression with least-squares fit for each class Ci (damaged or undamaged building) according to Equation (6) [55]:

$$\mathbf{L}\_{\mathbf{C}}(\mathbf{x}) = \sum\_{\mathbf{i}=1}^{\mathbf{C}\mathbf{F}} \alpha\_{\mathbf{i}} \mathbf{x}\_{\mathbf{i}} + \alpha\_{\mathbf{0}} \tag{6}$$

where Lc(x) is the least-squares fit and CF and α<sup>i</sup> are, respectively, the number of seismicrelated factors and the coefficient of the i-th element of vector x. The a posteriori probabilities in the leaves of the LMT are calculated using the linear logistic regression model with Equation (7) [52]:

$$\mathbf{p}(\mathbf{c}|\mathbf{x}) = \frac{\exp(\mathbf{L}\_{\mathbf{c}}(\mathbf{x}))}{\sum\_{\mathbf{c}'}^{\mathbf{c}} \exp(\mathbf{L}\_{\mathbf{c}'}(\mathbf{x}))} \tag{7}$$

where c is the number of building damage classes and Lc(x), the least-squares fit, is transformed in such a way that ∑ c c 0=1 Lc(x) = 0.

#### 2.5.3. Logistic Regression (LR)

The logistic regression (LR) model, developed by McFadden (1973) [56], is a multivariate regression analysis model that describes the relationship between a bivariate dependent parameter and several independent parameters [57] through the estimation of an optimal model. The addition of a link function suitable for a general linear regression model allows the parameter type to be continuous, discrete, or mixed, thus obviating the requirement for a normal distribution [58,59]. Some studies have shown that the LR model is more accurate than other types of models constructed for the same purpose [60–62]. The LR model based on a general linear model can be derived from Equations (8) and (9):

$$\mathbf{y} = \mathbf{b}\_0 + \mathbf{b}\_1 \mathbf{x}\_1 + \mathbf{b}\_2 \mathbf{x}\_2 + \mathbf{b}\_3 \mathbf{x}\_3 + \dots + \mathbf{b}\_n \mathbf{x}\_n \tag{8}$$

$$\mathbf{P} = \frac{\mathbf{e}^{\mathbf{y}}}{1 + \mathbf{e}^{\mathbf{y}}} \tag{9}$$

where y is the linear logistic model, b<sup>0</sup> is the y-intercept, b<sup>n</sup> is the logistic coefficient of each factor, n is the number of factors controlling a seismic event, x is the earthquake conditioning factor, and P is the probability of damage (ranging from 0 to 1) in the event of an earthquake [60].

#### **3. Results** throughout Gyeongju and Pohang, as seen in Figure 6a. According to the land cover map derived from the Korea Institute of Geoscience and Mineral Resources (KIGAM), Figure

earthquake [60].

**3. Results** 

#### *3.1. Building Damage Inventory Map* 6a reveals the DPM after the Gyeongju earthquake over residential, commercial, agricul-

*3.1. Building Damage Inventory Map* 

*Remote Sens.* **2021**, *13*, x FOR PEER REVIEW 13 of 26

The DPMs were generated from the Sentinel-1 dataset and geocoded to the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM), 1 arcsecond. Figure 6 reveals the map view of Sentinel-1 DPMs over Pohang. The assessment technique is most sensitive to the destruction of the built environment. Pixels are set to be relatively transparent where corresponding to areas where decorrelation did not significantly change during the time spanning the earthquake, suggesting little to no destruction. Increased opacity of the radar image pixels reflects increasing ground and building change or potential damage. The color range from yellow to red indicates an increasingly significant coherence change in the area covered by the pixel. Each pixel in the DPMs was registered to the SRTM DEM and had a corresponding dimension of about 30 m. tural, and vegetated areas. Figure 6b,c show the widespread COD in Seonggeon-dong and Gwangmyeong-dong, respectively, corresponding to collapsed houses (Figure 6d [63]),. These areas consist of residential and commercial areas. Figure 6a shows DPM that indicates COD over residential, commercial, industrial, agricultural, and vegetated areas affected by the Pohang earthquake. Some of the spatially large and strong coherence loss after the Pohang earthquake corresponds to Handong Global University (Figure 6e), associated with collapsed walls (Figure 6f), while Figure 6g shows Songdo-dong district, consisting of residential areas. The DPM of the Pohang earthquake was then used as a building damage inventory map to produce a seismic vulnerability map using machine learning.

to the SRTM DEM and had a corresponding dimension of about 30 m.

where y is the linear logistic model, b0 is the y-intercept, bn is the logistic coefficient of each factor, n is the number of factors controlling a seismic event, x is the earthquake conditioning factor, and P is the probability of damage (ranging from 0 to 1) in the event of an

The DPMs were generated from the Sentinel-1 dataset and geocoded to the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM), 1 arcsecond. Figure 6 reveals the map view of Sentinel-1 DPMs over Pohang. The assessment technique is most sensitive to the destruction of the built environment. Pixels are set to be relatively transparent where corresponding to areas where decorrelation did not significantly change during the time spanning the earthquake, suggesting little to no destruction. Increased opacity of the radar image pixels reflects increasing ground and building change or potential damage. The color range from yellow to red indicates an increasingly significant coherence change in the area covered by the pixel. Each pixel in the DPMs was registered

Here, we compared the distribution of coherence loss areas after the earthquakes

**Figure 6.** (**a**) Damage proxy maps (DPMs) after the Gyeongju and the Pohang earthquakes. DPM of (**b**) Seonggeon-dong and (**c**) Gwangmyeong-dong corresponding to (**d**) collapsed houses. DPM **Figure 6.** (**a**) Damage proxy maps (DPMs) after the Gyeongju and the Pohang earthquakes. DPM of (**b**) Seonggeon-dong and (**c**) Gwangmyeong-dong corresponding to (**d**) collapsed houses. DPM of (**e**) Handong Global University corresponding to (**f**) collapsed walls and in (**g**) Songdo-dong. Yellow to red pixels indicate increasingly more significant potential damage in seismic vulnerability mapping. Red and yellow stars indicate the epicenter of Pohang and Gyeongju earthquakes, respectively.

Here, we compared the distribution of coherence loss areas after the earthquakes throughout Gyeongju and Pohang, as seen in Figure 6a. According to the land cover map derived from the Korea Institute of Geoscience and Mineral Resources (KIGAM), Figure 6a reveals the DPM after the Gyeongju earthquake over residential, commercial, agricultural, and vegetated areas. Figure 6b,c show the widespread COD in Seonggeon-dong and Gwangmyeong-dong, respectively, corresponding to collapsed houses (Figure 6d [63]). These areas consist of residential and commercial areas. Figure 6a shows DPM that indicates COD over residential, commercial, industrial, agricultural, and vegetated areas affected by the Pohang earthquake. Some of the spatially large and strong coherence loss after the Pohang earthquake corresponds to Handong Global University (Figure 6e), associated with collapsed walls (Figure 6f), while Figure 6g shows Songdo-dong district, consisting of residential areas. The DPM of the Pohang earthquake was then used as a building damage inventory map to produce a seismic vulnerability map using machine learning.

#### *3.2. Relationship between Damaged Buildings and Related Factors*

FR values can provide information on the relationship between seismic vulnerability and related factors (geotechnical, physical, structural, social, and capacity). A ratio greater than 1 denotes that the class in the related factor has more impact on seismic vulnerability. The FR values calculated in this study are shown in Table 1. Building damage occurred in areas with elevations of 0–74 m. The classes with a strong impact on seismic vulnerability were: slope of 0◦–2.44◦ (FR = 3.50), elevation of 0–74 m (4.09), PGA of 3.37–28.56 gal (5.31), 0–26.90 km distance from epicenter (5.12), 7.36–39.08 km distance from fault (2.65), building constructed from concrete (2.63), building density of 1256–40,065 (5.13), building height of 7.88–422 m (2.31), commercial areas (3.99), child population of 228–10,747 (3.65), elderly population of 2252–31,218 (3.56), population density of 10,913–146,455 (5.17), 0–6.92 km distance from hospital (3.95), 0–3.39 km distance from police station (3.95), and 0–5.24 km distance from fire station (4.18). Areas in these classes are predicted to experience the highest degree of damage due to earthquakes.


**Table 1.** Frequency ratios of seismic-related factors.


**Table 1.** *Cont.*

#### *3.3. Seismic Vulnerability Map*

Seismic vulnerability maps were made using the training dataset compiled using the building damage inventory map from the DPM of the Pohang earthquake and applying machine learning algorithms, as discussed above. A combination of 15 seismic-related factors served as the dependent variables, and can mainly be classified as geotechnical, physical, structural, social, and capacity indicators. LogitBoost (Figure 7a), LMT (Figure 7b), and LR (Figure 7b) machine learning algorithms were applied to produce the seismic vulnerability maps. Each pixel in the study area was assigned a specific building damage value using the natural breaks method [7]. The seismic vulnerability maps were classified as safe, low to moderate, high, and very high vulnerability classes.

**Figure 7.** Seismic vulnerability map generated using three algorithms: (**a**) LogitBoost, (**b**) LMT, and (**c**) LR. **Figure 7.** Seismic vulnerability map generated using three algorithms: (**a**) LogitBoost, (**b**) LMT, and (**c**) LR.

Figure 8 shows the distribution of pixels in each seismic vulnerability map generated by LogitBoost, LMT, and LR models. In the LogitBoost model, 29.53% were classified as safe, 17.59% as low risk, 13.18% as moderate risk, 22.95% as high risk, and 16.74% as very Figure 8 shows the distribution of pixels in each seismic vulnerability map generated by LogitBoost, LMT, and LR models. In the LogitBoost model, 29.53% were classified as safe, 17.59% as low risk, 13.18% as moderate risk, 22.95% as high risk, and 16.74% as very

and low risk. For the LMT model, 19.74% were classified as safe, 18.54% as low risk, 8.25% as moderate risk, 18.27% as high risk, 35.20% as very high risk. Gyeongju and Pohang were found to be the most vulnerable to earthquake damage. Daegu was classified as a safe and low-risk city, while Busan and Ulsan were high and very high risk. In the LR model, 19.91% were classified as safe, 20.08% as low risk, 20.71% as moderate risk, 19.46% as high risk, and 19.85% as very high risk. The distribution of pixels in the low- and highrisk classes in Figure 8 shows similar results for each algorithm, although LMT shows larger numbers of pixels in very high-risk areas. The most vulnerable areas were Gyeongju and Pohang, whereas low and moderate risk areas were Busan and Ulsan. Seismic vul-

nerability classes were evenly distributed in Daegu.

els.

high risk. Gyeongju and Pohang were found to be the most vulnerable to earthquake damage. Among big cities in the study area, Busan, Daegu, and Ulsan were classified as safe and low risk. For the LMT model, 19.74% were classified as safe, 18.54% as low risk, 8.25% as moderate risk, 18.27% as high risk, 35.20% as very high risk. Gyeongju and Pohang were found to be the most vulnerable to earthquake damage. Daegu was classified as a safe and low-risk city, while Busan and Ulsan were high and very high risk. In the LR model, 19.91% were classified as safe, 20.08% as low risk, 20.71% as moderate risk, 19.46% as high risk, and 19.85% as very high risk. The distribution of pixels in the low- and high-risk classes in Figure 8 shows similar results for each algorithm, although LMT shows larger numbers of pixels in very high-risk areas. The most vulnerable areas were Gyeongju and Pohang, whereas low and moderate risk areas were Busan and Ulsan. Seismic vulnerability classes were evenly distributed in Daegu. *Remote Sens.* **2021**, *13*, x FOR PEER REVIEW 18 of 26

**Figure 8.** Distribution of pixels in seismic vulnerability classes for LogitBoost, LMT, and LR mod-**Figure 8.** Distribution of pixels in seismic vulnerability classes for LogitBoost, LMT, and LR models.

#### *3.4. Model Validation*

*3.4. Model Validation*  A validation step was conducted to assess the reliability of the seismic vulnerability map from each algorithm. ROC curve analysis is a standard way of validating the probability models used to generate seismic vulnerability maps, according to the area under the curve (AUC) [7,14,15]. Higher values indicate more accurate and reliable models. If the AUC, which ranges from 0 to 1, is lower than 0.5, the model is considered unacceptably inaccurate [49]. The accuracy of the seismic vulnerability maps generated using the three algorithms was then evaluated based on ROC curve analysis of the testing dataset (50% of all data). As seen in Figure 9, the AUC values were 0.769, 0.851, and 0.749 for LogitBoost, LMT, and logistic regression, respectively. Thus, the LMT model generated the best seis-A validation step was conducted to assess the reliability of the seismic vulnerability map from each algorithm. ROC curve analysis is a standard way of validating the probability models used to generate seismic vulnerability maps, according to the area under the curve (AUC) [7,14,15]. Higher values indicate more accurate and reliable models. If the AUC, which ranges from 0 to 1, is lower than 0.5, the model is considered unacceptably inaccurate [49]. The accuracy of the seismic vulnerability maps generated using the three algorithms was then evaluated based on ROC curve analysis of the testing dataset (50% of all data). As seen in Figure 9, the AUC values were 0.769, 0.851, and 0.749 for LogitBoost, LMT, and logistic regression, respectively. Thus, the LMT model generated the best seismic vulnerability map in this study. The results indicate that the algorithms are useful to map seismic vulnerability in the southeastern Korean Peninsula. Since all of the AUC values were higher than 0.5, the seismic vulnerability maps produced by all algorithms used in this study are acceptable for predicting vulnerable buildings in southeast Korea [49].

mic vulnerability map in this study. The results indicate that the algorithms are useful to map seismic vulnerability in the southeastern Korean Peninsula. Since all of the AUC values were higher than 0.5, the seismic vulnerability maps produced by all algorithms used in this study are acceptable for predicting vulnerable buildings in southeast Korea [49].

**Figure 9.** Receiver operating characteristic (ROC) curves associated with seismic vulnerability maps generated by LogitBoost, LMT, and LR algorithms. **Figure 9.** Receiver operating characteristic (ROC) curves associated with seismic vulnerability maps generated by LogitBoost, LMT, and LR algorithms.

#### **4. Discussion**

#### **4. Discussion**  *4.1. Building Damage Inventory Map*

*4.1. Building Damage Inventory Map*  The building damage inventory map was generated through the DPM method after the Pohang earthquake using Sentinel-1 SAR data. We found a good correlation between the DPM result and the released map from the KMA report about the Pohang earthquake in our qualitative validation. Therefore, the building damage inventory map is reliable. The KMA map was derived from field survey and damage survey data from local governments. The map shows the distribution of the Pohang earthquake's magnitude using the modified Mercalli intensity (MMI) scale [34]. The results were similar to the KMA report that these areas suffered MMI VII to VIII. KMA defines MMI VII to VIII as damage to major structural parts such as pillars, walls, and roofs, even in well-designed and well-The building damage inventory map was generated through the DPM method after the Pohang earthquake using Sentinel-1 SAR data. We found a good correlation between the DPM result and the released map from the KMA report about the Pohang earthquake in our qualitative validation. Therefore, the building damage inventory map is reliable. The KMA map was derived from field survey and damage survey data from local governments. The map shows the distribution of the Pohang earthquake's magnitude using the modified Mercalli intensity (MMI) scale [34]. The results were similar to the KMA report that these areas suffered MMI VII to VIII. KMA defines MMI VII to VIII as damage to major structural parts such as pillars, walls, and roofs, even in well-designed and well-built buildings. KMA map also shows some areas that suffered MMI V to VI; however, the DPMs did not show any CODs in these areas. The scales showed that the damage was inside buildings, such as minor cracks in walls and damage caused by dropping objects or tiles; therefore, SAR cannot detect a significant coherence change.

built buildings. KMA map also shows some areas that suffered MMI V to VI; however, the DPMs did not show any CODs in these areas. The scales showed that the damage was inside buildings, such as minor cracks in walls and damage caused by dropping objects or tiles; therefore, SAR cannot detect a significant coherence change. Here, the total damaged areas from DPMs were calculated by multiplying the total number of pixels by the pixel cell size, which for the Gyeongju earthquake yielded an area of 1.09 km2 and the Pohang earthquake yielded an area of 1.32 km2. The KMA reported a total of 69 damaged buildings by the Gyeongju earthquake, and 504 damaged buildings by the Pohang earthquake (MMI VII and VIII). The Pohang earthquake caused more dam-Here, the total damaged areas from DPMs were calculated by multiplying the total number of pixels by the pixel cell size, which for the Gyeongju earthquake yielded an area of 1.09 km<sup>2</sup> and the Pohang earthquake yielded an area of 1.32 km<sup>2</sup> . The KMA reported a total of 69 damaged buildings by the Gyeongju earthquake, and 504 damaged buildings by the Pohang earthquake (MMI VII and VIII). The Pohang earthquake caused more damage than the Gyeongju earthquake, some of which was due to the depth of the epicenter. The shallower the epicenter of an earthquake, the more damage it causes. The epicenter depth of the Pohang earthquake was 7 km, while the epicenter depth of the Gyeongju earthquake was 15 km. Additionally, the surface deformation in urban areas with high buildings caused by the Pohang earthquake affected the occurrence of damage in Pohang.

#### age than the Gyeongju earthquake, some of which was due to the depth of the epicenter. *4.2. Seismic Vulnerability Map*

*4.2. Seismic Vulnerability Map* 

The shallower the epicenter of an earthquake, the more damage it causes. The epicenter depth of the Pohang earthquake was 7 km, while the epicenter depth of the Gyeongju earthquake was 15 km. Additionally, the surface deformation in urban areas with high Estimating the seismic vulnerability of an area is vital for environmental management and land use planning, among other applications [19]. Although many methods and techniques have been developed to assess earthquake hazards around the world to date,

Estimating the seismic vulnerability of an area is vital for environmental management and land use planning, among other applications [19]. Although many methods and techniques have been developed to assess earthquake hazards around the world to date, the goals of all these studies are to reduce the economic losses and resulting losses. The the goals of all these studies are to reduce the economic losses and resulting losses. The method used to create seismic vulnerability maps affects the quality of the mapping. Machine learning techniques are effective. In particular, the method used to generate training and testing data is important. Accurate damage building inventory maps can be obtained using the DPM method; we combined DPM and GIS spatial data to produce an accurate seismic vulnerability map.

The seismic vulnerability maps of all three algorithms used in this study revealed that Pohang and Gyeongju are the most vulnerable areas to earthquake damage. We analyzed seismic-related factors by comparing general patterns of damaged buildings with factor maps. Pohang and Gyeongju are the cities where the two largest recent earthquakes happened. Several factors in these areas had higher values, such as the value of PGA and the distance from the epicenter. Therefore, Gyeongju and Pohang were areas with a high risk of earthquake vulnerability. Buildings within 0–36.90 km of an epicenter corresponded to damaged buildings. This result confirms that most buildings close to epicenters were damaged. The results revealed that areas with a high risk had high population and building density, including Gyeongju and Pohang. Similar to Gyeongju and Pohang, Busan, Daegu, and Ulsan are cities with a high density of buildings and populations. These cities consist of areas of low to moderate risk of seismic vulnerability. One of the main causes of a high risk of seismic damage in countryside areas is wood as a building construction material, while buildings constructed with steel reinforced with concrete have a lower level of seismic vulnerability. Therefore, attention should be paid to reconstruction, retrofitting, or renovation of the buildings in these areas.

Seismic vulnerability maps were validated based on ROC curves and AUC values to assess the accuracy of the maps using testing data, which comprised 50% of the total dataset. The AUC data showed that the LMT algorithm had the highest accuracy of 85.10%, which was 8.2% higher than the LogitBoost algorithm and 10.2% higher than the logistic regression algorithm. Therefore, the LMT algorithm is better to produce seismic vulnerability maps than LogitBoost and logistic regression algorithms.

A survey was conducted with 1256 Korean citizens using a questionnaire of citizen earthquake literacy (CEL) during spring 2020. The survey was based on three dimensions, including citizen knowledge, awareness, and management. We developed 15 questions associated with three dimensions. Table A1 in Appendix A presents the questionnaire that consists of five questions for each dimension, and all questions are configured to be answered on a 5-point Likert scale (strongly agree, agree, neutral, disagree, and strongly disagree). Principal component analysis as an exploratory factor analysis approach was used to examine if the items corresponded to each other and explore the construct validity of the instrument. Cronbach's alpha when an item was deleted was also calculated to see whether the items in the construct reliably measured the same latent variable. Cronbach's alpha values for each dimension were 0.847, 0.822, and 0.849, respectively. Furthermore, Cronbach's alpha when an item was deleted in this study was found to be between 0.77 and 0.84. This result indicated that the survey data used in this study were strongly reliable and consistent [64].

To characterize the profile of participants with higher (or lower) levels of seismic literacy, difference in means analyses using a *t*-test, ANOVA test, and Tukey's honestly significant difference (HSD) post hoc test were carried out. Table 2 shows the average values associated with earthquake literacy, including seismic knowledge, awareness, and management, broken down by sociodemographic characteristics of the participants. The *t*-test result shows that male citizens have higher earthquake literacy but not significantly more than female citizens in all three aspects. The ANOVA test was conducted to determine whether there were differences in the subcategories of seismic literacy depending on age, risk awareness, and final educational background. The results show that participants 20 years of age and below, and 60 years of age and above, declared a higher level of seismic literacy than participants between 30 and 50 years of age. Participants in their 30s declared the lowest level of seismic awareness and participants in their 40s declared the lowest level

of seismic knowledge and management among the different age groups. However, there was no significant difference in earthquake literacy in all aspects among all ranges of ages (under 20s to over 60s). Only participants who were aware that where they live is absolutely safe from earthquakes had a statistically significantly higher level of seismic management behavior. They also declared a higher level of seismic knowledge. Participants who were aware that where they live is not safe at all from earthquakes declared a higher level of seismic awareness but declared the lowest level of seismic knowledge and management. This result is an important warning sign for local and regulatory authorities to raise the citizen seismic knowledge and management. In the latter case, participants in graduate school and graduate school graduates declared a higher level of seismic literacy in all three aspects. Participants with an education level of less than high school declared the lowest level of seismic knowledge and high school graduates declared the lowest level of seismic awareness and management.


**Table 2.** Mean values for seismic literacy.

Further analysis was conducted to find significant differences through Tukey's HSD post hoc test and showed that citizens who were aware that where they live is absolutely safe from earthquakes had a statistically significantly higher level of seismic management behavior. Participants with a highly educated background had the highest level of seismic knowledge.

The participants of 60 years old and above declared the highest level of seismic literacy. Some authors posit that this could be explained because adults in this stage of life have acquired greater experience and care responsibilities (either for others or their assets), which may give rise to increased interest in involving themselves in preparedness measures [65]. Participants' educational background also influenced their seismic literacy, especially in seismic knowledge, in line with previous studies [66,67]. Groups with less seismic literacy should be a target of intervention in order to raise earthquake risk awareness and motivate them to adopt preparedness actions.

Nevertheless, in general, the participants declared relatively low earthquake literacy, including seismic knowledge, awareness, and preparedness. This means that they are not properly prepared for earthquake hazards. For the areas that were classified as highly vulnerable to earthquakes, the focus should be on building reconstruction and retrofitting [16], whereas for those in moderate-risk areas, the focus should be on building renovation to reduce their seismic vulnerability. Areas with high building and population densities, such as Daegu, Busan, and Ulsan, should develop education programs to improve earthquake risk awareness and seismic literacy. The government could increase earthquake literacy, starting with a conscious approach to Korean citizens, especially in high vulnerability areas.

Finally, the institution responsible for developing local disaster risk reduction plans and programs should appropriately characterize their target audiences and areas if they expect to obtain more effective and efficient results. We expect that the results reported in this study will be useful input to achieve this. However, this study has certain limitations. First, the questionnaire could provide more information of the participants profile such as marital status, children in the house, annual household income, and household type and could be analyzed to provide specific information about participants' preparedness. Future studies could investigate the interaction between these variables to find specific patterns of seismic literacy.

#### **5. Conclusions**

We conducted a survey of 1256 participants to investigate the seismic literacy of Korean citizens, including seismic knowledge, awareness, and management, using a questionnaire of citizen earthquake literacy (CEL) following the 2017 Pohang earthquake. The results declared that the citizens had low literacy, which means that they are not properly prepared for earthquake hazards. To develop an earthquake risk reduction plan and program efficiently and effectively, not only must one appropriately characterize the target audience, but also indicate high potential earthquake zones and potential earthquake damage. Therefore, this study mapped and analyzed the seismic vulnerability in southeast Korea using LogitBoost, logistic model tree (LMT), and logistic regression (LR) machine learning algorithms based on a building damage inventory map. The building damage locations were generated after the 2017 Pohang earthquake using the damage proxy map (DPM) method from the Sentinel-1 synthetic aperture radar (SAR) data. The DPMs manage to detect coherence loss, which indicates damaged buildings in residential and commercial areas due to the Pohang earthquake and show a good correlation with the Korea Meteorological Administration (KMA) report with modified Mercalli intensity (MMI) scale values of more than VII (seven). The damage locations were randomly divided into two datasets: 50% for training the vulnerability models and 50% for validating the models in terms of accuracy and reliability. Fifteen seismic-related factors were used to construct a model for each algorithm. Model validation based on the area under the receiver operating curve (AUC) was used to determine model accuracy. The AUC values of seismic vulnerability maps using the LogitBoost, LMT, and LR algorithms were 0.769, 0.851, and 0.749, respectively. We suggest that earthquake preparedness efforts should focus on reconstruction, retrofitting, renovation, and seismic education in areas with high seismic vulnerability in South Korea. The results of this study are expected to be beneficial for engineers and

policymakers aiming at developing disaster risk reduction plans, policies, and programs due to future seismic activity in South Korea.

**Author Contributions:** Conceptualization, C.-W.L. and K.-Y.L.; methodology, M.S. and M.H.; software, M.S. and M.H.; validation. K.-Y.L., M.S. and J.H.; formal analysis, C.-W.L.; investigation, J.H. and M.H.; resources, M.S.; data curation, M.S., M.H. and J.H.; writing—original draft preparation, A.S.N.; writing—review and editing, M.S., A.S.N. and J.H.; visualization, M.S.; supervision, K.-Y.L.; project administration, J.H.; funding acquisition, K.-Y.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 provided by the government of Korea (No. 2019R1A2C1085686) and also supported by a 2020 Research Grant from Kangwon National University.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The dataset of seismic-related factor can be found on: https://www. bigdata-environment.kr/user/main.do (accessed on 31 March 2021). Sentinel-1 A/B SAR images are freely available by European Space Agency and distributed and archived by Alaska Satellite Facility (https://search.asf.alaska.edu/#/ (accessed on 31 March 2021)). The DEM SRTM data and epicenter location used in this study were obtained from the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/ (accessed on 31 March 2021) and https://earthquake.usgs.gov/ earthquakes/map (accessed on 31 March 2021)).

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

#### **Appendix A**


**Table A1.** Survey questions asked by a researcher to prompt discussion.

1. Yes, absolutely; 2. Yes; 3. It's normal; 4. No; 5. Absolutely not.

#### **References**


*Article*

## **Land Subsidence Susceptibility Mapping in Jakarta Using Functional and Meta-Ensemble Machine Learning Algorithm Based on Time-Series InSAR Data**

#### **Wahyu Luqmanul Hakim , Arief Rizqiyanto Achmad and Chang-Wook Lee \***

Division of Science Education, Kangwon National University, Gangwon-do, Chuncheon-si 24341, Korea; wahyulhakim@kangwon.ac.kr (W.L.H.); ariefrizqiyanto@kangwon.ac.kr (A.R.A.)

**\*** Correspondence: cwlee@kangwon.ac.kr

Received: 10 September 2020; Accepted: 3 November 2020; Published: 4 November 2020

**Abstract:** Areas at risk of land subsidence in Jakarta can be identified using a land subsidence susceptibility map. This study evaluates the quality of a susceptibility map made using functional (logistic regression and multilayer perceptron) and meta-ensemble (AdaBoost and LogitBoost) machine learning algorithms based on a land subsidence inventory map generated using the Sentinel-1 synthetic aperture radar (SAR) dataset from 2017 to 2020. The land subsidence locations were assessed using the time-series interferometry synthetic aperture radar (InSAR) method based on the Stanford Method for Persistent Scatterers (StaMPS) algorithm. The mean vertical deformation maps from ascending and descending tracks were compared and showed a good correlation between displacement patterns. Persistent scatterer points with mean vertical deformation value were randomly divided into two datasets: 50% for training the susceptibility model and 50% for validating the model in terms of accuracy and reliability. Additionally, 14 land subsidence conditioning factors correlated with subsidence occurrence were used to generate land subsidence susceptibility maps from the four algorithms. The receiver operating characteristic (ROC) curve analysis showed that the AdaBoost algorithm has higher subsidence susceptibility prediction accuracy (81.1%) than the multilayer perceptron (80%), logistic regression (79.4%), and LogitBoost (79.1%) algorithms. The land subsidence susceptibility map can be used to mitigate disasters caused by land subsidence in Jakarta, and our method can be applied to other study areas.

**Keywords:** Jakarta; land subsidence susceptibility mapping; time-series InSAR; StaMPS processing; machine learning

#### **1. Introduction**

Several cities in Indonesia suffer from degradation at the ground level of buildings, known as land subsidence [1,2]. In Jakarta, this process has had a severe impact on urban infrastructure, leading to cracks in buildings, roads and damage to drainage systems [3]. These conditions are problematic because land subsidence may expand coastal flood areas due to sea-level rise [4]. Heavy monsoon rainfall [5] has caused frequent river flooding; if this occurs again, Jakarta could be submerged entirely underwater [1,4].

Recent studies of land subsidence in Jakarta have used various geodetic measurement methods, such as leveling surveys [3] and a global positioning system (GPS) surveys [6,7]. These studies indicated that excessive groundwater extraction is the leading cause of land subsidence and compaction to the vulnerable aquifer system [8]. This compaction may be exacerbated by natural consolidation since Jakarta's landform mostly comprises young alluvial soils that cannot support the weight of

human-made structures [9,10]. Therefore, it is essential to monitor land subsidence in Jakarta to predict further possible occurrences and mitigate damage [11,12].

Over the last decade, land subsidence susceptibility maps have been generated using geological, geomorphological, topographical, and hydrological data; these are considered the main factors influencing land subsidence [11,13,14]. Various methods are used to generate land subsidence susceptibility maps, including frequency ratios (FR) [12,15], weight of evidence (WOE) [16,17], logistic regression (LR) [18], evidential belief functions [11], analytical hierarchy processes (AHP) [19], support vector machines (SVM) [14,20], decision trees [21,22], fuzzy logic [12,23], adaptive neuro-fuzzy inference systems (ANFIS) [24,25], and artificial neural networks (ANN) [26,27]. In general, using a single modeling method leads to lower predictive accuracy than an ensemble method that uses a combination of models and machine learning algorithms [13,28]. Machine learning algorithms have the advantage of finding unpredictable relationships in datasets at multiple scales and have, thus, been recommended to obtain accurate land subsidence susceptibility maps [22].

Another challenge in generating land subsidence susceptibility maps is the low availability of land subsidence inventory maps. In this study, a land subsidence inventory map was generated via time-series interferometry synthetic aperture radar (InSAR) analysis. This technique can be applied to measure the displacement of the earth's surface with an accuracy of up to millimeters per year by improving the selection of coherent pixels and reducing atmospheric noise [29–33]. It has been widely used to monitor land subsidence and generate land subsidence maps, for example, in large cities in Mexico [33,34], Kurdistan, and Iran [35], in Yuncheng Basin, China [36], and in coastal cities and areas such as Venice, Italy [37], New Orleans, United States [38], and Shanghai, China [39]. The InSAR algorithm has been successfully applied to earthquakes [40], volcanic activities [41,42], crustal deformation [31], landslides [43], manmade deformations [44], excessive groundwater extraction [3], mining activities [45], and natural consolidation of young alluvial soil [9].

The recent studies of monitoring land subsidence in Jakarta using interferometry synthetic aperture radar (InSAR) techniques was conducted using the Small Baseline Subset (SBAS) algorithm from 2007 to 2009 [1] and using Geodesy and Earth Observing Systems-Persistent Scatterer Interferometry (GEOS-PSI) algorithm from 2007 to 2010 [46]. Both studies utilized Advanced Land Observing Satellite (ALOS) phased array type-L synthetic aperture radar (PALSAR) data to produce a land subsidence map, and both compared their results with GPS survey data. However, Jakarta's land subsidence map requires updating as it has remained unchanged for the past 10 years, and research on land subsidence susceptibility maps was not found in Jakarta.

Therefore, this study's objective was to generate the updated land subsidence map in Jakarta using the Stanford Method for Persistent Scatterers (StaMPS) with the Sentinel-1 synthetic aperture radar (SAR) dataset from March 2017 to May 2020 in ascending track and March 2017 to April 2020 in a descending track. The mean vertical deformation maps in ascending and descending tracks were compared to validate Jakarta's land subsidence location. After that, the land subsidence map obtained with this method was used as an inventory map to generate a land subsidence susceptibility map in Jakarta that predicted the areas at risk of land subsidence in the future. Two meta-ensemble machine learning algorithms (adaptive boosting (AdaBoost) and Logit Boost), and two functional machine learning algorithms (logistic regression and multilayer perceptron) were used. The result of the land subsidence susceptibility map produced by these algorithms was evaluated to compare all models' accuracy and reliability using receiver operating characteristic (ROC) curve analysis.

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

Jakarta is Indonesia's capital city (Figure 1a), located at 6◦17' south (S) and 106◦82' east (E), on the northern coast of western Java. It is considered a lowland area, with an average altitude of ±7 m above sea level [3,47]. In 2019, its population was 10.5 million, with a population growth rate of about 1.19% per year. The population density was 15,900 people per km<sup>2</sup> , within a total land area of 662.33 km<sup>2</sup> [47]. Historically, the population of Jakarta gradually migrated to the other districts and municipalities that are part of the Jakarta Metropolitan Region (JMR), which covers a total area of 5897 km<sup>2</sup> and includes Bogor, Depok, Tangerang, and Bekasi [48] (Figure 1b). This migration to Jakarta's outskirts led to increased urban development that could cause land subsidence in these areas [3]. Our study area covers the JMR; however, for simplicity, we refer to it as Jakarta. The three administrative areas were chosen as study area due to land subsidence reported on those areas from the recent study using GPS surveys, leveling surveys, and InSAR techniques. Training and testing data in Figure 1b used in this study to generate a land subsidence susceptibility map were assessed from the land subsidence inventory map.

**Figure 1.** (**a**) Study area location in Indonesia and (**b**) training and testing data sets from land subsidence location within Jakarta Metropolitan Region (JMR) depicted from Sentinel-2 satellite imagery taken on 28 August 2020, divided into three administrative areas: Jakarta, Tangerang, and Bekasi.

The geological and geomorphological area in Jakarta lithologically was dominated by alluvium landform (50.20%) and alluvium fans (19.66%), followed by Serpong form (7.03%) that dominated by fragmented pumice sandstones, some limestone, and andesite. A volcanic pyroclastic flow formed tuff Banten (6.77%) and upper Banten tuff (4.88%), Sandstone unit (4.05%), swamp deposit (1.75%), Cihoe form (1.61%), beach ridge deposit (1.53%), and old alluvium (0.74%). Subang form (0.64%) was dominated by layered claystone lithology with limestone and marl found locally, marine deposit (0.52%), and young volcanic rocks (0.22%); Bojongmanik form (0.22%) was dominated by alternating sandstones and claystone inserted by limestones and coastal deposit (0.13%); Parigi form (0.05%) was dominated by medium limestone, lake (0.01%), and sandstone tuff (0.01%). The domination of alluvium landform has a risk of land subsidence due to the compaction of natural consolidation worsened by a human-made structure [3,9] covering the Jakarta area's land use. Land use in Jakarta was dominated by settlement area with 45.72% and followed by 39.79% of rice field, 5.94% of dryland

agriculture, 4.75% of fish pond, 1.98% of shrub-mixed dryland farms, 0.84% of airport, 0.51% of swamp, 0.35% of estate crop plantation, 0.09% of barren land, 0.03% of secondary mangrove forest, and 0.004% of swamp shrub.

Research on land subsidence in Jakarta has been conducted since 1982 using leveling surveys. There were two main monitoring periods: (1) 1982–1991 and (2) 1991–1997. Land subsidence was found in the period 1982–1991 in three regions with the highest accumulated subsidence value compared to other regions, namely, two regions in the northwestern part (Cengkareng and Kalideres districts) and the third in the northeastern part of Jakarta (Kemayoran-Sunter district); accumulated subsidence was found in Kalideres district with up to 68.5 cm between 1982 and 1991 with an annual rate of subsidence around 8 cm/year. Accumulated land subsidence was found in the Cengkareng district, with up to 60 cm between 1982 and 1991, with an annual subsidence rate of around 7 cm/year. Accumulated land subsidence was found in the Kemayoran-Sunter district, with up to 70 cm with an annual subsidence rate of around 6 cm/year. The second period of monitoring land subsidence using a leveling survey between 1991 and 1997 highlighted one region with the highest accumulated subsidence value than other regions in the Kalideres district, with up to 154.1 cm with an annual subsidence rate of around 23 cm/year [7].

Following this research, land subsidence in Jakarta was monitored using GPS surveys from 1997 to 2005, and two regions affected by land subsidence were reported in two stations: (1) Kwitang district and (2) Pantai Mutiara district. The accumulated land subsidence was found in the Kwitang district, with up to 48 cm with an annual rate of subsidence rate of 5 cm/year; accumulated land subsidence in Pantai Mutiara district was around 50 cm with an annual rate of subsidence around 4.6 cm/year. [6]. Land subsidence in Jakarta measured using leveling and GPS surveys positively correlated with excessive groundwater extraction and sea-level rise [6,7,49].

Land subsidence was also reported using InSAR techniques based on ALOS PALSAR satellite data in 2007–2009 using the SBAS method, with land subsidence found in Pluit district with an annual subsidence rate of 21.6 cm/year. Land subsidence was found in Cengkareng district with an annual subsidence rate around 21.8 cm/year, in Bekasi district with an annual rate of subsidence of around 10.6 cm/year, and in Karawang district with an annual rate of around 16.4 cm/year [1]. Research was also conducted with ALOS PALSAR satellite data within the period 2007–2010 using the GEOS-PSI method, finding land subsidence in the coastal area and lowland area in northwestern Jakarta within Penjaringan and Cengkareng districts with an annual subsidence rate of up to 26 cm/year, with accumulated subsidence rates of up to 86.5 cm between 2007 and 2010. Land subsidence was observed in the Bekasi district with an annual subsidence rate of up to 11.5 cm/year [46].

#### **3. Material and Methods**

#### *3.1. SAR Datasets*

A land subsidence inventory map for generating the land subsidence susceptibility map in Jakarta was generated using Sentinel-1 SAR C-band data (5.5 cm wavelength) provided by the European Space Agency (ESA). The SAR data were acquired from March 2017 to May 2020 (91 datasets in the ascending track with path number 98 and frame number 1160, with vertical–vertical (VV) polarization) and in the period of March 2017 to April 2020 (89 datasets in the descending track with path number 47 and frame number 614, with vertical–vertical (VV) polarization). The ascending and descending datasets are listed in Tables 1 and 2, and the reference dates with zero delta day and zero perpendicular baselines from the ascending track are shown on 15 October 2018, whereas those from the descending track are shown on 16 November 2018; both reference dates are shown in bold text in table number 45. A perpendicular baseline graph (generated from Tables 1 and 2 and shown in Figure 2) was used to visualize the temporal baseline from the reference date.


**Table 1.** Acquisition dates (format ddmmyyyy) of the Sentinel-1 synthetic aperture radar (SAR) datasets in ascending track. Delta days = number of days between each acquisition date. B⊥ = perpendicular baseline. The reference dates in ascending tracks are shown in bold text.

**Table 2.** Acquisition dates (format ddmmyyy) of the Sentinel-1 Sentinel-1 synthetic aperture radar SAR datasets in descending track. Delta days = number of days between each acquisition date. B⊥ = perpendicular baseline. The reference dates in descending tracks are shown in bold text.


**Figure 2.** Perpendicular baseline graph from (**a**) ascending track and (**b**) descending track.

#### *3.2. Land Subsidence Conditioning Factors*

In total, 14 land subsidence conditioning factors consisting of geological, geomorphological, topographical, and hydrological data were chosen as the conditioning factors (Table 3) in this study as they are considered the main factors influencing land subsidence [11,13,14]. Each factor was classified using the quantile method, and the factors with different cell size resolutions were resampled into raster datasets with 30 m cell size from each conditioning factor to standardize each factor's resolution.


**Table 3.** Land subsidence conditioning factors in the study area. DEM, digital elevation model; SRTM, Shuttle Radar Topography Mission.

A recent study in Jakarta found that the leading causes of land subsidence in Jakarta are groundwater extraction, the load of buildings and construction, natural consolidation of alluvium soil, and tectonic activities [3]. In this study, groundwater drawdown data were collected from the Groundwater Conservation Center, The Ministry of Energy and Mineral Resources of Indonesia. The observation of groundwater drawdown data is carried out periodically through an automatic water level record (AWLR) system. The annual change in groundwater level from 2019 to 2020 was calculated from 15 monitoring wells, and the map was constructed (Figure 3a) using the inverse

distance weighting (IDW) interpolation method from the annual change in groundwater level data. The obtained groundwater drawdown data in this study were insufficient compared to the study area due to the limited monitoring wells over the study area. Although the available data are few, the use of groundwater drawdown data is essential to determine the relationship between land subsidence and groundwater extraction.

**Figure 3.** Land subsidence conditioning factors: (**a**) groundwater drawdown; (**b**) rainfall intensity; (**c**) distance to roads; (**d**) distance to rivers; (**e**) distance to faults; (**f**) drainage density; (**g**) land use; (**h**) lithology; (**i**) elevation; (**j**) slope; (**k**) aspect; (**l**) profile curvature; (**m**) plan curvature; (**n**) topographic wetness index.

The groundwater level can increase due to conditional factors indirectly associated with land subsidence, such as rainfall, distance from a river, and river density, which can recharge the groundwater level [13,50]. Four years of daily rainfall data (from 2017 to 2020) from seven weather stations were acquired from the Indonesian Meteorology, Climatology, and Geophysical Agency of Indonesia, and the annual rainfall intensity was calculated and interpolated using the IDW tool in geographic information system (GIS) software (ArcGIS; ESRI, Redlands, CA, USA) (Figure 3b). However, the availability of rainfall data was limited compared to the study area due to the cost and constraints of acquiring data. Nevertheless, the utilization of the rainfall data in this study is important to correlate the topographical factor related to the infiltration flow affecting the soil's strength.

Road, river, and fault networks were acquired from the Geospatial Information Agency of Indonesia from the main road map, main river map, and fault map at a scale of 1:250,000 of polygonal-shaped data using the Atlas map of Indonesia. The information on the location of roads, rivers, and faults was used to construct a distance map (Figure 3c, Figure 3d, and Figure 3e, respectively) using the Euclidean distance tool, and the maps were classified using the quantile method to provide suitable classes within a 30 m cell size. The drainage density or river map density was estimated using the kernel density tool (Figure 3f). Distance to the road and land use are related to urban development in Jakarta, which can affect land subsidence [3]. The land-use map (Figure 3g) was acquired from the Ministry of Environment and Forestry of Indonesia that used Landsat data to generate a land cover map. The geological parameters describe the spatial correlation of lithological landform and land subsidence in Jakarta as being caused by the compaction of the alluvium soil landform [3]. The lithology map (Figure 3h) was acquired from the Geological Atlas map of Indonesia from the Ministry of Energy and Mineral Resources, and the polygonal-shaped map was converted into a raster map with 30 m cell size using GIS tools.

The topographical map, which includes elevation, slope, aspect, profile curvature, plan curvature, and topographic wetness index (TWI) (Figure 3i, Figure 3j, Figure 3k, Figure 3l, Figure 3m, and Figure 3n, respectively) data extracted from the digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) [51], was constructed using the basic terrain analysis tools. Elevation influences the hydrological properties and soil moisture, whereby a lower-elevation area possibly gains more precipitation than a higher-elevation area [13,14]. The slope is associated with land subsidence because it affects the infiltration of rainfall (a steeper surface slope decreases infiltration) [13,21], and the aspect is the second derivative of the slope that has a relationship with land subsidence because the slope aspect affects the strength of the soil due to the moisture preservation and the amount of vegetation [11]. Profile curvature is associated with flow speed, sediment, and erosion quantity, while plan curvature is perpendicular to the slope and indirectly affects land subsidence by influencing convergence and divergence of flow across the surface [13]. The TWI defines the degree of deposition of water at a specific site [52]. The topographical factors extracted from the SRTM DEM are widely used as conditioning factors in land subsidence susceptibility mapping [11–13,28].

The relationship of land subsidence occurrence with the conditioning factors in each class is described in Table 4. A ratio greater than 1 denotes that the class in the conditioning factor is more correlated with the land subsidence occurrence [11]. The calculation of frequency ratio shows that the land subsidence occurred in areas with groundwater level data between 20.27 and 21.00 and between 23.31 and 34.55 m below ground level, while land subsidence also occurred in areas with more annual rainfall intensity due to the recharge of groundwater level and the usage of groundwater. Areas between 0 and 126 m with roads were more correlated with land subsidence occurrence. The land subsidence areas correlated with fault distance were between 15,944 and 68,975 m from the fault location. There were three drainage density classes correlated with land subsidence occurrence and settlement areas correlated with land subsidence. In terms of lithological factors, alluvium, alluvium fans, beach ridge deposits, and sandstone landforms were considered more correlated with the land subsidence occurrences. There were three classes in the elevation map, four classes in the aspect map, two classes in the slope map, one class in the plan curvature, one class in the profile curvature, and three classes in the topographic wetness index map correlated with land subsidence occurrences.


**Table 4.** Relationship between land subsidence occurrence and conditioning factors using frequency ratio (FR) model.


**Table 4.** *Cont.*


**Table 4.** *Cont.*

#### *3.3. Illustration of Methodology*

The methods to generate land subsidence susceptibility maps using functional and meta-ensemble algorithms are described below and illustrated in Figure 4.


quantile methods in GIS tools with a similar environment of 30 m cell size for each factor; then, the number of subsidence occurrences in each class was calculated using the cross-tabulation tool in GIS. Next, we calculated the ratio between the percentage of pixels of each conditioning factor class and the percentage of subsidence occurrence pixels to obtain the FR value as follows:

% of class conditioning factor

% of land subsidence occurrence (1)

Frequency Ratio =


#### *3.4. StaMPS Processing*

StaMPS (Stanford Method for Persistent Scatterers) is an analysis method used to facilitate the generation of time-series deformation images of all terrains, including nonurban areas. The StaMPS algorithm uses the spatial correlation of phase measurements rather than a functional temporal model to identify PS pixels. StaMPS processing does not use a model to describe how the displacement rate varies with time. To identify PS pixels in a single master stack of interferograms, StaMPS uses the phase characteristics from the dominant point scatterer in each area and creates interferograms from SAR images. It also reduces decorrelation [54]. Thus, the StaMPS algorithm can identify and extract the deformation signal from stable pixels in all terrains.

The slave images in the acquired SAR datasets were resampled to perfectly match the master image through a co-registration process before generating an interferogram. The co-registration process was applied to two different images to produce refined SAR images, which were then cropped to focus only on the area of interest. Next, differential InSAR (DInSAR) images were generated by subtracting the topographic InSAR images generated using the interferograms from the topographic phases of the SRTM DEM [55,56]. We used the PS method to measure the displacement of the earth's surface [57] to generate the time-series deformation map. The main processes generated a single master stack of interferograms and topographic phase removal [33].

The StaMPS algorithm is shown in Figure 5. Multiple images were co-registered to generate a single master image for the ascending track on 15 October 2018 and the descending track on 16 November 2018; co-registered images with the topographic phases removed were subjected to amplitude and phase noise analysis to derive a subset comprising all PS pixels, with weeding performed using a threshold value of 0.4. The wrapped phase of the selected pixels was corrected for the spatially uncorrelated look-angle error. The phase was then unwrapped, and PS outputs were generated. The deformation map from the line of sight (LOS) displacement could be converted into vertical deformation data [58,59] by assuming the horizontal deformation as very small compared to the vertical deformation caused by land subsidence [60–62]. Recent studies using GPS and leveling surveys reported that the land subsidence in Jakarta shows a vertical deformation [7], and a vertical deformation pattern was also found in research using InSAR to monitor land subsidence in Jakarta [1,46]; hence, the deformation from the line of sight (LOS) in this study could be assumed as negligible and could be converted directly into vertical deformation value using Equation (2) by dividing the displacement or deformation from the line of sight (dLOS) by the cosine of incident angle (θ) from the radar signal. The results of vertical deformation were assigned a negative value from the initial ground-level observation point, indicating that the land subsidence occurred vertically at that point [9].

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

**Figure 5.** Flowchart of Stanford Method for Persistent Scatterers (StaMPS) Processing.

#### *3.5. AdaBoost*

AdaBoost is a machine learning algorithm introduced by Freund and Schapire (1997) [63]. AdaBoost's classifier uses an adaptive resampling technique that produces a series of individual classifiers to classify training samples accurately. The frequency of variables selected by a weak learner was examined, and the relative importance of the variables could be determined [64]. AdaBoost combines multiple weak learners to derive a single strong learner by repeatedly calling a weak classifier and adjusting the attributed weight to the sample. The new classifier focuses more on the misclassified

sample as the misclassified sample is increased compared to the correct sample. The final weight is obtained by adding or subtracting the updated weight in each iteration. The final model is obtained by dividing each final weight by the total adjusted weight [65]. The general form of the AdaBoost algorithm is as follows [66]:

	- a. Fit the classifier fm(x) ∈ {−1, 1} using weights w<sup>i</sup> with the training data;
	- b. Compute errm= E<sup>w</sup> h <sup>1</sup>(y,fm(x))<sup>i</sup> , cm<sup>=</sup> log 1 − err<sup>m</sup> errm
	- c. Set w<sup>i</sup> <sup>←</sup> <sup>w</sup><sup>i</sup> exp<sup>h</sup> <sup>c</sup>m1(y,fm(x))<sup>i</sup> , i = 1, 2, . . . , N, and renormalize so that P <sup>i</sup> wi= 1;

;

3. Output the classifier: signhP<sup>M</sup> m=1 cmfm(x) i .

#### *3.6. LogitBoost*

LogitBoost is a modified version of the AdaBoost algorithm, introduced by Friedman, Hastie, and Tibshirani [66], where the exponential loss function is replaced with the log-likelihood loss function. This method reduces classification error and is less sensitive to noise [28,67]. LogitBoost can handle multiple class problems and uses a regression scheme as the base learner for classification [64]. The general form of the LogitBoost algorithm is as follows [66]:

	- a. Compute the working response and weights:

$$\mathbf{l}\_{\mathbf{i}} = \frac{\mathbf{y}\_{\mathbf{i}}^{\*} - \mathbf{p}(\mathbf{x}\_{\mathbf{i}})}{\mathbf{p}(\mathbf{x}\_{\mathbf{i}}) \left(1 - \mathbf{p}(\mathbf{x}\_{\mathbf{i}})\right)},\tag{3}$$

$$\mathbf{w\_{i}} = \mathbf{p(x\_{i})} \left(\mathbf{1} - \mathbf{p(x\_{i})}\right);\tag{4}$$


$$\mathbf{f}(\mathbf{x}) \leftarrow \mathbf{f}(\mathbf{x}) + \frac{1}{2} \mathbf{f}\_{\mathbf{m}}(\mathbf{x}),\tag{5}$$

$$\mathbf{p}(\mathbf{x}) \leftarrow \frac{\mathbf{e}^{\mathbf{f}(\mathbf{x})}}{\mathbf{e}^{\mathbf{f}(\mathbf{x})} + \mathbf{e}^{-\mathbf{f}(\mathbf{x})}};\tag{6}$$

3. Output the classifier: sign[F(x)] = signhP<sup>M</sup> m=1 fm(x) i .

#### *3.7. Logistic Regression*

Logistic regression is a statistical method used to find the best model to describe the correlation between a dependent variable and several independent variables. This method's advantage is that the variables do not need to be normally distributed [68]. It also offers several ways of selecting the best predictor for use in the P probability model [69,70]. The equations describing the logistic regression are as follows [28,69–71]:

$$\mathbf{f}(\mathbf{x}) = \text{logit}(\mathbf{P}) = \ln\left[\frac{\mathbf{P}}{\mathbf{1} - \mathbf{P}}\right] = \mathbf{c}\_0 + \mathbf{c}\_1 \mathbf{x}\_1 + \dots + \mathbf{c}\_n \mathbf{x}\_n \tag{7}$$

$$\mathbf{P} = \frac{1}{\mathbf{1} + \mathbf{e}^{-\mathbf{f}(\mathbf{x})}} = \frac{1}{\mathbf{1} + \mathbf{e}^{-(\mathbf{c}\_0 + \mathbf{c}\_1 \mathbf{x}\_1 + ... + \mathbf{c}\_n \mathbf{x}\_n)}},\tag{8}$$

where f(x) is a linear combination function called logit(P), P is the probability of subsidence occurrence, 1 − P is the probability that subsidence will not occur, x1, x<sup>2</sup> , . . . , x<sup>n</sup> are input variables, c<sup>0</sup> is the model intercept, and c1, . . . , c<sup>n</sup> are the approximated coefficients of regression.

#### *3.8. Multilayer Perceptron*

The multilayer perceptron is a machine learning algorithm based on the ANN technique that consists of input layers, hidden layers, and output layers [72]. The advantages of the multilayer perceptron algorithm are as follows: the distribution of the training data does not require any assumptions, most inputs are selected during the training process based on the weight adjustment, and the relative importance of the various input measures does not need to be determined [73]. The equation representing the multilayer perceptron for land subsidence classification is as follows [28,72,73]:

$$\mathbf{m} = \mathbf{f}(\mathbf{c}),\tag{9}$$

where f(c) is a hidden function that is optimized during the training process for a given network architecture via the adjustable weights, and c = c<sup>i</sup> for i = 1, . . . , 14, which is a vector containing 14 land subsidence conditioning factors.

#### **4. Results**

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

The mean vertical deformation maps for Jakarta shown in Figure 6a,b were derived by time-series InSAR analysis based on the StaMPS algorithm. We overlaid a true color red/green/blue (RGB) composite image from Sentinel-2 taken on 28 August 2020. Six areas from both ascending and descending tracks showed high deformation; thus, we plotted a deformation trend for these areas.

The vertical deformation rate in the ascending track at Point P1 (Figure 6c) represents the Kosambi area, which subsided 189.48 mm from 2017 to 2020. Point P2 (Figure 6c) represents the Cengkareng area, which subsided 184.02 mm from 2017 to 2020. Point P3 (Figure 6e) represents the Ciledug area, which subsided up to 155.22 mm from 2017 to 2020. Point P4 (Figure 6e) represents the Penjaringan area, which subsided 148.36 mm from 2017 to 2020. Point P5 (Figure 6g) represents the Bekasi area, which subsided 128.17 mm from 2017 to 2020. Point P6 (Figure 6g) represents the Cikarang area, which subsided 271.84 mm from 2017 to 2020.

The vertical deformation rate in the descending track at Point P1 (Figure 6d) represents the Kosambi area, which subsided 210.07 mm from 2017 to 2020. Point P2 (Figure 6d) represents the Cengkareng area, which subsided 216.19 mm from 2017 to 2020. Point P3 (Figure 6f) represents the Ciledug area, which subsided up to 155.92 mm from 2017 to 2020. Point P4 (Figure 6f) represents the Penjaringan area, which subsided 148.31 mm from 2017 to 2020. Point P5 (Figure 6h) represents the Bekasi area, which subsided 107.19 mm from 2017 to 2020. Point P6 (Figure 6h) represents the Cikarang area, which subsided 257.94 mm from 2017 to 2020.

Figure 6g,h show that the deformation pattern was mostly linear. However, Figure 6c–f, representing the Kosambi, Cengkareng, Ciledug, and Penjaringan areas, show quite periodic subsidence with the standard deviation of the vertical deformation rate being higher than the vertical deformation rate in Figure 6g,h. These results occurred due to the seasonal effect of groundwater extraction, and those areas were surrounded by a residential area that mostly used groundwater as the water source. Meanwhile, P6 is one of the most significant industrial areas in Indonesia. The mean vertical deformation maps from ascending and descending tracks were compared to validate the accuracy of the land subsidence inventory map using the StaMPS algorithm, and the result showed a good correlation in Figure 6i with a coefficient of correlation (R<sup>2</sup> ) up to 0.9584 between ascending and descending tracks.

42

**Figure 6.** Mean vertical deformation map for Jakarta depicted using Sentinel-2 image in (**a**) ascending tracks and (**b**) descending tracks; the vertical deformation rate at P1 (Kosambi) and P2 (Cengkareng) in (**c**) ascending and (**d**) descending tracks; the vertical deformation rate at P3 (Ciledug) and P4 (Penjaringan) in ascending (**e**) and descending (**f**) tracks; the vertical deformation at P5 (Bekasi) and P4 (Cikarang) in ascending (**g**) and descending tracks (**h**). (**i**) The comparison of mean vertical deformation between ascending and descending tracks. (**j**) Kriging interpolation of time-series deformation map from mean vertical deformation map of descending track, resulting in the extension of the land subsidence inventory map; the blue area is considered as the nonoccurrence area of land subsidence.

The persistent scatterer density of both tracks in the east and west areas was relatively low due to them being wetland areas and more vegetated than other areas. The SAR dataset used in this study from Sentinel-1 SAR C-band data with 5.5 cm wavelength could not deeply penetrate beneath the trees. Thus, to overcome that limitation, the interpolation of the persistent scatterer points from the descending track was constructed using kriging interpolation in GIS tools to provide land subsidence information over the study area as shown in Figure 6j [74].

#### *4.2. Land Subsidence Susceptibility Map*

The land subsidence susceptibility model's performance depended on the calculated parameters (Table 5) for optimization used in this study.

A land subsidence susceptibility map was generated using 14 land subsidence conditioning factors, training data from our land subsidence inventory map, and four different algorithms: LogitBoost (Figure 7a), AdaBoost (Figure 7b), logistic regression (Figure 7c), and multilayer perceptron (Figure 7d). Land subsidence susceptibility indices were generated for all unique pixels in the study area. The susceptibility indices were reclassified using the quantile method to identify feature pairs in five susceptibility classes: very low, low, moderate, high, and very high [28,64].


**Table 5.** Calculated parameters for the algorithms used in this study.

**Figure 7.** Land subsidence susceptibility map for Jakarta, generated using four different algorithms: (**a**) AdaBoost, (**b**) LogitBoost, (**c**) logistic regression, and (**d**) multilayer perceptron.

The proportion of very high susceptibility land was quite similar for all four algorithms. However, the map generated by AdaBoost differed from those of the other three methods because the AdaBoost model could only classify susceptibility into four classes due to the limits of the probability range. Therefore, the very low susceptibility class was excluded. The susceptibility class proportions (pixel distributions) are shown in Figure 8. For the Adaboost algorithm, the proportions were 0%, 57.88%, 2.54%, 22.02%, and 17.56% for the very low, low, moderate, high, and very high classes, respectively; the respective values for the LogitBoost algorithm were 33.33%, 27.09%, 7.03%, 16.77%, and 15.78%, and those for the logistic regression algorithm were 32.62%, 13.86%, 16.44%, 18.51%, and 18.57%. Finally, the multilayer perceptron algorithm's respective values were 40.64%, 18.70%, 14.02%, 13.42%, and 13.23%.

The distribution of pixels in the very high class and high class in Figure 8 showed similar results for each algorithm, resulting in the maps of these classes being quite similar. The very high class was considered to be the land subsidence areas shown in the mean vertical deformation map in Figure 6a,b, with similarity seen because of the training data used in this study being acquired from the land subsidence inventory map with a large spatial resolution, which allowed more effectively defining the land subsidence area in the susceptibility map. The moderate and high classes were considered land subsidence areas in the future, and the very low and low classes were areas with the lowest probability of land subsidence in the future. The difference between the moderate class from AdaBoost and that from other algorithms is that AdaBoost did not consider other areas far from the land subsidence occurrences. Other algorithms showed that residential areas located in alluvium landforms had a reasonable possibility of land subsidence.

**Figure 8.** Proportions of susceptibility classes for land subsidence susceptibility maps generated using four different machine learning algorithms.

#### *4.3. Model Validation*

The accuracy of all four algorithms in this study was evaluated by ROC curve analysis [12,17]. ROC curve analysis is a standard way of validating the probability models used to generate land subsidence susceptibility maps, according to the area under the curve (AUC) [22,28]. Higher values indicate more accurate and reliable models. If the AUC, which ranges from 0 to 1, is lower than 0.5, the model is considered unacceptably inaccurate [75].

Land subsidence susceptibility maps produced using functional (AdaBoost, LogitBoost) and meta-ensemble (logistic regression and multilayer perceptron) algorithms were compared. The ROC curves for the four algorithms are shown in Figure 9. The largest AUC of 0.811 was from the AdaBoost algorithm (blue line in Figure 9). The multilayer perceptron algorithm had the next largest AUC (0.800; purple line in Figure 9), followed by the logistic regression (0.794; green line in Figure 9) and LogitBoost algorithms (0.791; red line in Figure 9).

**Figure 9.** Receiver operating characteristic (ROC) curves for the land subsidence susceptibility maps produced by functional (logistic regression and multilayer perceptron) and meta-ensemble (AdaBoost and LogitBoost) algorithms.

Variables having at least one tie between the positive actual state group and the negative actual state group can bias the results. Since all of the AUC values were higher than 0.5, the land subsidence susceptibility maps produced by all algorithms used in this study are acceptable for predicting high-risk subsidence areas in Jakarta [75]. However, the AdaBoost algorithm map had the best performance.

#### **5. Discussion**

#### *5.1. Land Subsidence Inventory Map*

A land subsidence inventory map was successfully created through time-series InSAR analysis of Sentinel-1 datasets, using the StaMPS algorithm. The map is displayed as a vertical deformation map in Figure 6a,b for descending and ascending tracks. The comparison was made to validate the mean vertical deformation's accuracy used as a land subsidence inventory map. A coefficient of correlation of 0.9584 was found. As stated above, there were six high-deformation areas (Figure 6c–h) within ascending and descending tracks, which were equally affected (in terms of subsidence) by increased groundwater usage and existing buildings [3].

Figure 6g,h show areas with largely linear deformation rates, the standard deviations of which were smaller than those of the area shown in Figure 6c–f. The periodic subsidence shown was affected due to the variable climate in Jakarta, which is generally characterized by high rainfall but a dry summer. In the rainy season, the groundwater level of all aquifers beneath Jakarta increases, which manifests as deformation, i.e., an uplift in the ground level due to a change in the underlying material's thickness [34,76–79].

As shown by the land subsidence occurrence (Figure 6a,b) and groundwater drawdown (Figure 3a) data, land subsidence for P1–P5, which are between 20.58 and 34.55 m below ground level, correlated with the groundwater level. The correlation between land subsidence and groundwater level might have been stronger if groundwater level data were available for the whole study area. The groundwater level is deepest for P1 and P2, ranging from 22.47 to 34.55 m. P6 has more significant deformation because it contains the largest industrial estate in Indonesia.

Jakarta is mostly situated on young alluvial soil (Figure 3h), which cannot tolerate the maximum compression force of many buildings [9]. Thus, compaction of the unconsolidated alluvial soil occurs and is exacerbated by groundwater extraction (because pore pressure is reduced, leading to further clay compaction) [80,81].

#### *5.2. Land Subsidence Susceptibility Map*

The method used to create land subsidence susceptibility maps strongly affects the quality of the mapping. Machine learning techniques are effective [12,13,28,64]. In particular, the method used to generate training and testing data is important. Accurate land subsidence inventory maps can be obtained using InSAR; we combined InSAR and GIS spatial data to produce an accurate land subsidence susceptibility map [12]. Nevertheless, the distribution data of groundwater drawdown and rainfall intensity data from the interpolation process might not represent the whole study area. They could have affected the spatial distribution of the raster map. Access to the monitoring wells and rainfall station outside the Jakarta area could provide more accurate analysis to determine the relationship between those factors and land subsidence.

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; fulkifadhillah@kangwon.ac.kr (M.F.F.); ariefrizqiyanto@kangwon.ac.kr (A.R.A.)

**\*** Correspondence: cwlee@kangwon.ac.kr

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 km<sup>2</sup> [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/km<sup>2</sup> , 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**<sup>⊥</sup>

ሺሻ**<sup>1</sup> No.** 

**Date (yyyymmdd)**  **Days B**<sup>⊥</sup> ሺሻ **<sup>1</sup>**

**Date (yyyymmdd)** 

#### **No. Date (yyyymmdd) Days B**<sup>⊥</sup> *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

ሺሻ **<sup>1</sup> No.** 


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

<sup>1</sup> <sup>B</sup>⊥: 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. The SAR images were then resampled such that the slave images matched the master image. 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 β*<sup>i</sup>* 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 m<sup>3</sup> /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 <sup>1</sup> , Fabiola D. Yépez 1,\* , Diego F. Lozano-García 2 , Víctor H. Guerra Cobián 1 , Adrián L. Ferriño Fierro <sup>1</sup> , Héctor de León Gómez <sup>1</sup> , Ricardo A. Cavazos González <sup>1</sup> and Adriana Vargas-Martínez <sup>2</sup>**


**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 km<sup>2</sup> 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.

<sup>1</sup> 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 km<sup>2</sup> . 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 km<sup>2</sup> , 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 <sup>3</sup> , Yunjee Kim <sup>2</sup> , Namwook Cho <sup>2</sup> 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 km<sup>2</sup> (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.

<sup>a</sup> Aerial photograph before and after Mt.Umyeon landslides from Kakaomap [47]. <sup>b</sup> Topographical factors were extracted from digital topographic map by National Geographic Information Institute. <sup>c</sup> The detailed soil map produced by Rural Development Administration. <sup>d</sup> Field investigation data produced by Korean Geotechnical Society [38]. <sup>e</sup> 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

**Figure 6.** Map showing input parameters constructed using kriging interpolation based on the geotechnical properties of the sampling location (**a**) cohesion; (**b**) friction angle; (**c**) hydraulic conductivity; (**d**) unit weight. **Figure 6.** Map showing input parameters constructed using kriging interpolation based on the geotechnical properties of the sampling location (**a**) cohesion; (**b**) friction angle; (**c**) hydraulic conductivity; (**d**) unit weight.

In sequence, automatic weather system (AWS) data were used to construct rainfall intensity data. The Korea Meteorological Administration (KMA) operates over 510 AWS sites to monitor the atmospheric conditions near the ground in real time. To construct rainfall data in the study area, rainfall data was measured at 1-h intervals at seven AWS stations near Mt. Umyeon (Table 3). The rainfall intensity over the 16 h from 19:00 26 July to 09:00 27 July, 2011 was obtained to calculate the rainfall intensity at each AWS. Kriging interpolation was also applied to construct rainfall intensity data for the entire study area (Figure 7). **Table 3.** Automatic weather system (AWS) observatory information. In sequence, automatic weather system (AWS) data were used to construct rainfall intensity data. The Korea Meteorological Administration (KMA) operates over 510 AWS sites to monitor the atmospheric conditions near the ground in real time. To construct rainfall data in the study area, rainfall data was measured at 1-h intervals at seven AWS stations near Mt. Umyeon (Table 3). The rainfall intensity over the 16 h from 19:00 26 July to 09:00 27 July, 2011 was obtained to calculate the rainfall intensity at each AWS. Kriging interpolation was also applied to construct rainfall intensity data for the entire study area (Figure 7).


**AWS Observatory Name Number of Stations Latitude Longitude Height above Sea Level (m) Table 3.** Automatic weather system (AWS) observatory information.

Yongsan 415 37.52038 126.97611 31.73

Gwacheon 590 37.44028 127.00249 47

**Figure 7.** Map showing rainfall intensity using kriging interpolation based on automatic weather system observatory locations. **Figure 7.** Map showing rainfall intensity using kriging interpolation based on automatic weather system observatory locations.

Finally, a 2011 landslide occurrence map was prepared to validate the analysis. A total of 103

landslide locations were extracted by superimposing a 1: 5000 digital topography onto an aerial photograph with a spatial resolution of 50 cm before and after the landslide (Figure 8). Finally, a 2011 landslide occurrence map was prepared to validate the analysis. A total of 103 landslide locations were extracted by superimposing a 1: 5000 digital topography onto an aerial photograph with a spatial resolution of 50 cm before and after the landslide (Figure 8).

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

**Figure 8.** Landslide locations.

**4. Methodology**

*4.1. Physically Based Model* 

to shear strength. The shear stress is the stress that acts on the fracture surface, taking into account the weight and influence of the slope material, whereas the shear strength is based on the action of the pore water pressure on the vertical stress acting perpendicularly to the fracture surface. Therefore, an FS value of less than 1 indicates a destructive state; and is expressed as [53]:

∅tan ߙଶݏܿ (௪ݖ௪ߛ − ܦߛ) + ܿ

 ߙ cos ߙ sin ܦߛ where *c* is the cohesion [kN/m3], *γ* is the unit weight of the slope material [kN/m2], *D* is the depth of the slope [m], *γw* is the unit weight of water [kN/m2], *zw* is the groundwater level [m], *α* is a gradient of the slope [degree], and *φ* is a frictional angle [degree]. The SHALlow landsliding STABility model (SHALSTAB) was used, which is based on the distributed hydrological model [36] and wet index of the infinite slope stability model, to calculate the groundwater level (*zw*). SHALSTAB is a hydraulic model that considers only the ground flow while ignoring the outflow of the surface in

(1)

FS =

We used a physical slope model, an infinite slope, in which the depth of the ground is shorter

#### **4. Methodology**

#### *4.1. Physically Based Model*

We used a physical slope model, an infinite slope, in which the depth of the ground is shorter than the length of the landslide. The infinite slope model is the most suitable model for GIS-based landslide analysis and is widely used to analyze landslides caused by rainfall [49–52]. The stability of the infinite slope model is expressed by the factor of safety (FS), which is the ratio of shear stress to shear strength. The shear stress is the stress that acts on the fracture surface, taking into account the weight and influence of the slope material, whereas the shear strength is based on the action of the pore water pressure on the vertical stress acting perpendicularly to the fracture surface. Therefore, an FS value of less than 1 indicates a destructive state; and is expressed as [53]:

$$\text{FS} = \frac{c + (\gamma D - \gamma\_w z\_w) \cos^2 \alpha \tan \mathcal{Q}}{\gamma D \sin \alpha \cos \alpha} \tag{1}$$

where *c* is the cohesion [kN/m<sup>3</sup> ], γ is the unit weight of the slope material [kN/m<sup>2</sup> ], *D* is the depth of the slope [m], γ*<sup>w</sup>* is the unit weight of water [kN/m<sup>2</sup> ], *z<sup>w</sup>* is the groundwater level [m], α is a gradient of the slope [degree], and ϕ is a frictional angle [degree]. The SHALlow landsliding STABility model (SHALSTAB) was used, which is based on the distributed hydrological model [36] and wet index of the infinite slope stability model, to calculate the groundwater level (*zw*). SHALSTAB is a hydraulic model that considers only the ground flow while ignoring the outflow of the surface in the hydraulic model proposed by [54], which considers the flow and surface runoff in shallow ground. It is expressed as follows:

$$\frac{q \times a}{T \times \sin \alpha \times b} = \frac{z\_w}{D} = w \tag{2}$$

where *q* is rainfall [m/day], *a* is the watershed area [m<sup>2</sup> ], *b* is the width of the contour line [m], and *T* is transmissivity [m<sup>2</sup> /day]. Furthermore, *w* is the wetness index, which is the ratio of the groundwater level to the depth of the slope (*zw*/*D*), the relative depth of the actual groundwater level with respect to the slope. The wetness index *w* is between 0 and 1; its maximum value is 1 because the groundwater level does not exceed the depth of the slope. The groundwater level is expressed as:

$$\mathbf{w} = \text{Min}\left(\frac{q \times a}{T \times \sin \alpha \times b}, 1\right) \tag{3}$$

#### *4.2. Monte Carlo Simulation*

In order to calculate the safety factor through a mathematical model, input data such as shear strength (cohesion force and friction angle) of the ground or groundwater level are required. However, these data are usually obtained through limited field surveys or indoor experiments, which is absolutely insufficient compared to the size of a wide area. Uncertainty intervenes in these data [55]. In this study, the probabilistic analysis technique was applied to quantitatively consider uncertainty and reflect it in the analysis. Such methods can be used to estimate the probability that the value of a state function satisfies a threshold by assuming that the input variable is randomly selected from a specific distribution. The probabilistic analysis method replaces the safety factor to evaluate the risk of slope using the probability of failure and is evaluated as the most effective method to quantify uncertainty among the various techniques proposed so far.

Probabilistic methods include MC simulation, first-order second moment (FOSM), the point estimate method (PEM), etc. FOSM and PEM are approximate methods, so their accuracy is relatively low compared to MC simulation [31,56]. MC simulation can represent the variability of input variables by randomly generating input variables, and is suitable for analyzing one or more random variables [57]. In this study, analysis was performed using Monte Carlo simulation among the probabilistic analysis techniques.

The MC simulation method first determines the model of the state function. The state function used in landslide susceptibility analysis depends on FS; and we calculate the probability of failure, which is the probability that FS is less than 1. Second, the probability distribution of the input variables is calculated. Because information about the probability distribution of input variables is not available, this is considered to be a normal distribution with an appropriate average and standard deviation [28,32,33,58]. Third, N random values were extracted from the distribution of input variables, and N values of FS were calculated by substituting the values of the randomly generated input variables. Finally, the probability of failure, i.e., the proportion of the N FS values that are less than 1, was calculated as follows:

$$P\_f = p(FS < 1) \tag{4}$$

#### *4.3. Applications and Validation*

In our MC simulation, the cohesion and frictional angle were considered as random variables. According to previous studies, the probability distributions of geotechnical characteristics follow a normal or lognormal distribution, so the probability distribution of the input variables was assumed to be normal in this study [59–61]. Since the mean and standard deviation are required to define the distribution of input variables, the average was calculated based on the constructed data, and the standard deviation was calculated from the coefficient of variation. In previous studies, the ranges of the coefficient of variation of the internal friction angle and cohesion were 10–20% and 25–30%, respectively [56,62,63]. Because the entire study area was composed of gneiss and thus the geological characteristics were relatively similar, the minimum value of each coefficient of variation was assumed to be 10% for the internal friction angle and 25% for the cohesive coefficient, and the number of repetitions was set to 100,000.

The probability of failure is established by probabilistic methods, in contrast to the deterministic analysis in which we interpret an area to be unstable when FS is less than 1 [64]. Therefore, based on the results of previous studies, landslide-susceptible areas were classified based on failure probabilities of 1%, 5%, and 10% [64]. We then calculated the landslide prediction rate, which is the ratio of the number of landslides in a landslide-susceptible area. Furthermore, we carried out a deterministic analysis of the same dataset and calculated the FS by taking a simple average of the data. Finally, we compared the results from the MC simulation to those obtained using deterministic techniques.

Finally, the area under the curve (AUC) of the receiver operating characteristic (ROC) was calculated for validation purposes. The *x*-axis of the ROC graph is the ratio of the expected landslide area where it shows high susceptibility of landslides, and the *y*-axis is the probability of landslides. The AUC is calculated from the area under the ROC graph and has values between 0 and 1. The closer it is to 1 (100%), the higher the accuracy [45].

#### **5. Result**

The results of the MC simulation are summarized in Table 4 and shown in Figure 9. When a 1% probability of failure was considered to indicate a susceptible area, 57.75% and 42.25% of the study area was found to be unstable and stable, respectively, and the landslide prediction rate was 95.15%. When the susceptible area was set based on a probability of failure of 5%, the unstable area was 54.23%, the stable area was 45.77%, and the landslide prediction rate was 91.26%. When we defined the susceptible areas based on a probability of failure of 10%, the unstable area was 52.02%, the stable area was 47.98%, and the landslide prediction rate was 90.29%. By contrast, according to the deterministic analysis method, the unstable area was 42.72%, the stable area was 57.28%, and the prediction rate, which in this case indicated the proportion of the predicted landslides that occurred, was 81.55%. Additionally, the AUC calculation of the ROC graph of the MC simulation was 0.7332 (73.32%) (Figure 10).

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

**Figure 9.** Map showing (**a**) the probability of failure evaluated using Monte Carlo simulation and (**b**) the factor of safety evaluated using the deterministic analysis. **Figure 9.** Map showing (**a**) the probability of failure evaluated using Monte Carlo simulation and (**b**) the factor of safety evaluated using the deterministic analysis. <sup>107</sup>

**Figure 10.** Receiver operating curve of the results from the Monte Carlo simulation.



simulation more than 5% 45.77 54.23 91.26 Monte Carlo simulation more than 10% 47.98 52.02 90.29 Deterministic analysis less than 1 57.28 42.72 81.55 Based on the results of both analytical methods, we conclude that the stable area occupies more than 40% of the study area. According to the results of the MC simulation, as the reference Based on the results of both analytical methods, we conclude that the stable area occupies more than 40% of the study area. According to the results of the MC simulation, as the reference probability of failure decreases, the proportion of unstable regions increases with increasing landslide prediction rate. Our MC simulation predicted relatively high proportions of unstable areas and landslide prediction rates compared to the deterministic analysis. An AUC value of 0.7 or more can be interpreted as indicating good predictive performance [65,66]. The AUC obtained from the MC simulation (0.7332) was greater than 0.7, which could explain the high landslide prediction rate. The good predictive performance of the MC simulation is attributed to the variability of the input variables.

#### probability of failure decreases, the proportion of unstable regions increases with increasing landslide prediction rate. Our MC simulation predicted relatively high proportions of unstable **6. Discussion**

Monte Carlo

areas and landslide prediction rates compared to the deterministic analysis. An AUC value of 0.7 or more can be interpreted as indicating good predictive performance [65,66]. The AUC obtained from the MC simulation (0.7332) was greater than 0.7, which could explain the high landslide prediction rate. The good predictive performance of the MC simulation is attributed to the variability of the input variables. In this study, we conducted a landslide susceptibility analysis of the Mt. Umyeon area, a representative example of urban landslide damage. An interpolation method was used to construct data for the entire study area because we could only obtain subsurface data and geotechnical characteristics for a few points. In addition, it was infeasible to obtain sufficient data to infer the probability distribution due to the cost and constraints of acquiring the actual property values, so we

assumed the probability distribution to be normal. Since it was difficult to obtain a sufficient amount of field survey data compared to the scope of the study area, but it was difficult to apply the unsaturated soil theory, which requires many experimental values and well reflects the behavior of the groundwater level of the actual ground; saturated soil theory with relatively simple interpretation was applied in this study. In addition, several hydraulic parameters were assumed based on existing studies. Therefore, a more accurate groundwater model can be constructed with sufficient data on hydraulic parameters obtained through experiments. Future studies should increase the number of data acquisition points and experiments at each point to compensate for these limitations.

The results of the landslide susceptibility analysis of the Mt. Umyeon area can be summarized as follows. The MC simulation results showed that the landslide prediction rate (90.25, 91.26, and 95.15%) was significantly higher than the landslide prediction rate (81.55%) of the deterministic method. Comparing the results of the deterministic analysis method and the MC simulation, which is a probabilistic analysis method, by using the coincidence ratio with the location coordinate of the landslide, the landslide prediction rate of 8.7% to 13.6% is higher than the deterministic analysis method. Through these results, the probability analysis method can be applied considering the variability of the rainfall intensity, so accuracy can be improved by considering the uncertainty inherent in the rainfall intensity. MC simulation result also yielded a high AUC value (0.7332). We attributed this improvement in results to the fact that MC simulation can include the variability and uncertainty of the input variables. Moreover, the proportion of unstable areas around Mt. Umyeon exceeded 40%. This is thought to be due to the weathered gneiss distributed throughout the area and the high rainfall intensity.

#### **7. Conclusions**

Following the Mt. Umyeon landslide, a representative landslide case in the Seoul metropolitan city, South Korea, it has become necessary to analyze landslide susceptibility. The purpose of this study is to analyze the susceptibility of landslide disasters considering uncertainty before the occurrence of widespread landslides. Analysis based on physical slope models can be used to evaluate landslide susceptibility in areas without prior landslide occurrence. A spatial database of Mt. Umyeon, the study area, was constructed and analyzed by applying it to the infinite slope model that is similar to the characteristics of landslide occurrence in the study area. In addition, landslide susceptibility was applied and analyzed based on a physically based model and MC simulation, a probabilistic analysis technique considering uncertainty. An infinite slope model was used as the physically based model. In the GIS environment, it was possible to analyze landslide susceptibility in the study area, and by using an infinite slope model and using probabilistic techniques, it was possible to evaluate the landslide susceptibility as a quantitative indicator of probability of failure. Furthermore, to evaluate the accuracy of our landslide predictions, we applied a deterministic method and compared the results. Finally, the accuracy of the landslide prediction was calculated based on the AUC.

This study confirms that it is possible to evaluate high-accuracy landslide susceptibility without prior information of landslide occurrences by combining a physical slope model with probabilistic method. By varying the reference probability of failure from 1% to 10% in the MC simulation, it was possible to adjust the safety level as needed. This means that the reference failure probability can be varied according to the purpose of analysis. For example, a susceptibility map with a high standard probability of failure should be used to determine the locations of disaster prevention structures to minimize costs. Conversely, if the danger zone is temporarily set to minimize human damage, the susceptibility map should be based on the minimum probability of failure. In this way, the same dataset and probabilistic technique can be used for different purposes.

To prevent a repeat of the damage incurred by the Mt. Umyeon landslide, it is necessary to carry out landslide susceptibility studies of areas where landslides have not occurred. In particular, prior landslide susceptibility analysis should be carried out in areas with high population densities to minimize large-scale damage. The methodology presented herein can be used to prepare measures to reduce the damage caused by landslides by analyzing landslide susceptibilities in areas without landslide occurrence data, where landslides have not occurred previously. Furthermore, this methodology can be applied to various regions by extracting input factors by setting an infinite slope model that reflects regional characteristics in consideration of landslide characteristics.

**Author Contributions:** Conceptualization, S.L. and M.-J.L.; data curation, S.L. and J.J.; formal analysis, Y.K. and N.C.; funding acquisition, M.-J.L.; investigation, J.J. and S.L.; methodology, J.J., Y.K. and N.C.; resources, S.L. and M.-J.L.; software, J.J.; validation, J.J., Y.K. and N.C.; supervision, M.-J.L.; visualization, S.L. and J.J.; writing—original draft, J.J.; writing—review and editing, S.L. and M.-J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was conducted at Korea Environment Institute (KEI) with support from Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2018R1D1A1B07041203). This research was also supported by a grant (20010017) of R&D Program of Responding Technology to Climate Disaster such as Heat wave, funded by Ministry of Interior and Safety (MOIS, Korea).

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

#### **References**


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

*Letter*

## **Intelligent WSN System for Water Quality Analysis Using Machine Learning Algorithms: A Case Study (Tahuando River from Ecuador)**

**Paul D. Rosero-Montalvo 1,2,\* , Vivian F. López-Batista <sup>1</sup> , Jaime A. Riascos 3,4 and Diego H. Peluffo-Ordóñez 4,5,6**


Received: 24 April 2020; Accepted: 12 June 2020; Published: 20 June 2020

**Abstract:** This work presents a wireless sensor network (WSN) system able to determine the water quality of rivers. Particularly, we consider the Tahuando River from Ibarra, Ecuador, as a case study. The main goal of this research is to determine the river's status throughout its route, by generating data reports into an interactive user interface. To this end, we use an array of sensors collecting several measures such as: turbidity, temperature, water quality, pH, and temperature. Subsequently, from the information collected on an Internet-of-Things (IoT) server, we develop a data analysis scheme with both data representation and supervised classification. As an important result, our system outputs a map that shows the contamination levels of the river at different regions. Furthermore, in terms of data analysis performance, the proposed system reduces the data matrix by 97% from its original size, while it reaches a classification performance over 90%. Furthermore, as an additional remarkable result, we here introduce the so-called quantitative metric of balance (QMB), which measures the balance or ratio between performance and power consumption.

**Keywords:** prototype selection; river pollution; supervised classification; WSN

#### **1. Introduction**

Rivers are natural watercourses that commonly come from both precipitation (surface runoff), and snowpacks (e.g., water stored in glaciers). Regularly, they flow towards lakes, sea, oceans, or another river. Urban rivers are responsible for providing water resources to crops and human beings as well as navigation purposes. Certainly, this natural resource may not be everlasting. As a matter of fact, there is currently a great deficit of water reserves due to deforestation, inappropriate and excessive use of fertilizers and pesticides, causing environmental issues [1,2]. Likewise, the urbanization and industries have had collateral adverse impact directly on the water quality of river ecosystems worldwide [3]. Besides, the population growth produces enormous wastewater that enters into the rivers without any environmental control. United Nations (UN) settled that 90% of such waste is not correctly treated, and 70% of the industries discharge contaminant content without any adequate standards or rigorous inspections [4,5]. Water pollution contains high levels of biochemical oxygen demand (BOD), nitrogen, and phosphorus. So it is necessary to develop systems that support the

detection and measuring of the contamination levels in rivers to maintain an optimal ecological balance, limiting environmental damage and preventing diseases spread [6]. Consequently, city governments have stated environmental policies intended to create urban regeneration initiatives around the care of their rivers [7]. In this connection, Ecuador, as our case study, has no any short or long-term plan to improve either the urban or rural river conditions [8].

Traditionally, water quality monitoring uses collected samples for laboratory testing, enabling then a wide range of analyses. Notwithstanding, it results impractical to manually measure water pollution at different points along the river. Moreover, this sort of tests may take a few days, and probably not reaching as a good precision as that of in-situ sampling [9]. Nowadays, the use of sensors for monitoring environmental conditions has received significant attention due to the real-time data collection, flexibility, and portability [4,10]. Following from this, the creation of a wireless sensor network (WSN) that combines several sensors with a data processing system and wireless communication can allow for an adequate measure of the water quality, where each sensor becomes a node that shares information among them as well as to a central server [11,12]. Thus, these data are greatly useful for further robust analyzes of water pollution in rivers. However, the large amount of data demands the implementation of machine learning algorithms to create systems that automatically can detect high levels of water pollution and make proper decisions. For that purpose, historical data (training data) become valuable to turn WSN nodes into intelligent systems [13,14].

Consequently, this work presents a novel system composed of three WSN nodes for monitoring in real-time the water pollution present in the Tahuando River (located in Ibarra, Ecuador) using machine learning algorithms. To do so, we establish different measurement points wherein each WSN node acquires the river's conditions data to be later processed internally by the system. In this sense, we consider water-quality variables, namely pH, turbidity, temperature and dissolved solids. Additionally, we carry out a sensor integration and calibration stage for eliminating reading errors. Finally, we sent these data to a cloud server, using a mobile network, where we visualize the node's information with its proper geo-location. As relevant results, a reduction of the required training set of 97% is accomplished by using is the condensed nearest neighbor (CNN) method as a prototype selection approach, as well as the classification stage—with k-NN—reaches 90.6% of performance. Then, our work is an exploratory study on different methods for both prototype selection and data classification applied to water treatment. Therefore, we have no gold standard result or benchmark method. Instead, an exhaustive comparison of representative methods is presented.

The fact that the data analysis process is implemented directly into the WSN represents a novelty itself for the development of both intelligent embedded systems, and data analysis platforms under low-computational resources. The rationale of creating an intelligent system including in-situ data analysis tasls (e.g., data classification) lies in the fact that an embedded systems can perform automatic decision-making processes with no requiring an external server. As well, it enables the possibility that even non-expert operators can readily interact with the system. In addition, it represents a solution to one of the main open issues of WSNs design, namely: information redundancy, which constraints the battery life-time, and often requires the incorporation of an external server for decision-making procedures. Additionally, to display a report of the current river's status, we implement an interactive user interface.

The rest of the manuscript is organized as follows: Section 2 gathers some remarkable related works. Section 3 describes both the system design and the data analysis proposed for implementing the machine learning algorithms. Section 4 presents the tests and results. Finally, Section 5 gathers the concluding remarks.

#### **2. Related Works**

Some works [5,6,9,15] have extensively worked on the estimation of water pollution, presenting different solutions for determining pollution state and its levels along several rivers located

in China. Other works [16,17] analyze river status using satellite photographs. Meanwhile, in [4,10] WSN are instead preferred for data acquisition.

The work presented in [18] develops a WSN to determine the water quality level for human consumption through GPRS-generated data analysis, which is carried out on an external-to-WSN server holding a communication module. Similarly, another work [19] uses a high-performance external server. Specifically, it presents a system able to measure the quality of the water stored in tanks or reservoirs. In this connection, other works have proposed alternatives to improve the data processing aimed at reaching an admissible performance while involving a lower computational burden. An approach to do so is by minimizing the communication load, as done in [20] wherein an additional data compression stage is incorporated—particularly, the principal component analysis (PCA) algorithm is used. By compressing (or reducing the dimensionality of) data, the sending-packets process through WSN is enhanced in terms of performance and processing time. Similarly, the work presented in [21] performs a data analysis including temperature, pH, electrical conductivity (EC) and dissolved oxygen (DO) sensors, whose data are processed on a server and its result is sent back to the proposed WSN for decision-making. Another approach, which is becoming a new embedded systems paradigm is the design of intelligent systems performing an in-situ data analysis. For instance, in [22] the redundancy is minimized following a data fusion criterion to better manage the WSN computational resources, and bring an adequate energy consumption. Under this new paradigm where data analysis is carried out into the same system handling the data acquisition, the design of a system related to water quality monitoring results not only novel but proper. Indeed, on doing so, there would be enabled an affordable, large-coverage and easy-to-use WSN system, which along with right sensors will help environmental or health-related agencies or bodies to effectively make decisions regarding the quality of natural water from a specific source. Following from this, the work [22] involves stages for data acquisition, processing, and visualization.

Nonetheless, no one of these solutions presents an in-situ data analysis. From the reviewed literature, only [23] presents an analysis of rivers in Ecuador. All of the aforementioned works presented appealing solutions to determine the water conditions of different rivers. However, in spite of all these efforts, there are still many open issues, such as: real-time data analysis, sensor calibration, and sending information to storage servers located far away from the acquisition point, among others.

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

Broadly, the proposed system consists of the following stages: initial conditions of the study region (Section 3.1), WSN design for accurate data acquisition (Section 3.2), and the data analysis with both the criteria for prototype selection, and supervised classification (Section 3.3).

#### *3.1. Initial Conditions of the Study Region*

The city of Ibarra (Ecuador) is the capital of the province of Imbabura with a dry-temperate climate of 18 ◦C on average. The urban population is 109 thousand and a rural population of approximately 45 thousand inhabitants. Its main commercial activity is the production of wooden articles and services to medium-scale industries. Regarding its water supply, 90% is carried out through the public distribution network, while the rest is for the use of river and vertier water [24]. Tahuando river is an important water resource in the Imbabura province, being part of the natural system of Ecuador. Due to its ability to transport and the flowing of its waters, it can withstand a large number of pollutants. However, there are several modifications at the ecological level, such as the loss of aquatic species, foul-smelling, and watercolor changes, among others. In Ecuador, only 10% of wastewater is treated. In Ibarra, around 600 liters per second of these waters are discharged into the Tahuando River, causing that no urban regeneration based on the increase in tourism can be carried out [25]. The Tahuando River is located at 0.4◦ latitude and 78.13◦ longitude. It encompasses an extension of 12 km from the community of Pesillo towards Salinas, in the Ibarra city. Figure 1 depicts the geographical location and basin.

**Figure 1.** Geographical description of Tahuando river. Zoomed view of the river route highlighting remarkable surrounding communities in Ibarra city's urban sector (**right**), and a widespread view regarding the rural section of Imbabura province (**left**).

#### *3.2. Wireless-Sensor-Network Design*

The design of our WSN approach is followed from the considered water-quality-related variables: pH, turbidity, temperature and dissolved solids. The considered sensor network is as follows: Firstly, we measure the turbidity and identify what kind of pollutants can be found in the river, such as: wastewater, chemicals, among others. Secondly, we use a pH sensor to determine if the water composition is acidic or basic as well as a quality sensor (total dissolved solids, TDS) to assess the level of dissolved oxygen in the water (cleanliness) [9]. Thirdly, we incorporate a temperature sensor to determine the water's changes and its relationship with the rest of the variables. To suitably develop the WSN network, we consider several operational requirements in the selection of sensors, such as reliability, precision, availability, ease-of-use, and scalability. Furthermore, in the selection of the WSN network processor system, we consider the number of pins and sensor libraries, as done in a previous work [11]. Specifically, the considered sensors are: SKU: SENO189 (turbidity), SKU: PH-7BNC (pH), Ds18b20 (temperature), RB-Dfr-797 (TDS). As well, the Arduino Uno is selected as processing system. Additionally, we use both global position module (GPS) and mobile communications (GSM) Sim808 to send data. Finally, there is a Lipo rider battery manager for power supply with a solar charging system. Figure 2 presents the considered sensors along with the processor system (Arduino Uno).

Likewise, we calibrate each sensor as follows: sensor SKU:PH-7BNC (pH) has a linear response, so its tuning is based on measuring the voltage of several pH solutions. Particularly, we use two solutions, the first one was pH = 4.01, getting a voltage of 2.98 volts; meantime, the second one was pH = 6.86, obtaining a voltage of 2.53 v. Thus, the equation to obtain the estimated pH is:

$$\text{pH} = -5.65 \ast (v\_1) + 21.15,\tag{1}$$

with *v*<sup>1</sup> as the voltage obtained by the sensor SKU:PH-7BNC. Likewise, the turbidity sensor SKU: SENO189 gives a reading ranging between 2.5 to 4.3 volts with values between 3000 and 0 turbidities (NTU), respectively. According to its datasheets, we can write the following equation:

$$\text{NTU} = -1120.4 \ast v\_2^2 + 5742.3 \ast v\_2 - 4352.9 \tag{2}$$

where *v*<sup>2</sup> is the voltage registered by the sensor SKU:SENO189. On the other hand, the datasheet of the temperature sensor Ds18b20 indicates that each Celsius degree can be transformed using the equality 10 mv = 1 ◦C; thus, the equation is:

$$Temp = \frac{v\_3 \ast 5}{1023 \ast 0.01} \,\text{\AA}\tag{3}$$

with *v*<sup>3</sup> as the voltage obtained by the sensor Ds18b20.

Finally, the TDS RB-Dfr-797 sensor provides a flexible calibration protocol, with a reset button, we can return to the initial conditions, that is, a TDS value of 23 mv. Consequently, we refresh the Arduino program and use the next equation:

$$\text{TDS} = \frac{(30 \ast 5 \ast 1000) - (75 \ast v\_4) \ast 5 \ast (1000/1024)}{75 - 0.23},\tag{4}$$

where *v*<sup>4</sup> is the voltage obtained by the sensor RB-Dfr-797.

**Figure 2.** Demonstrative diagram of the proposed WSN system. Considered sensors (SKU: SENO189 (turbidity), SKU: PH-7BNC (pH), Ds18b20 (temperature), RB-Dfr-797 (TDS)), and the processor (Arduino Uno).

Upon sensor configuration, each *v<sup>i</sup>* value will correspond to a digital-analog converter (DAC) with a resolution of 10 bits, already in the microprocessor Arduino Uno. Furthermore, we implement the moving average recursive filter to reduce the acquisition errors and smoothing the signal from each DAC. This filter takes a subset (window) of *N* samples, and calculate its arithmetic average to estimate a filtered sample as [26]. This filter is implemented in each DAC separately through the following equation:

$$y\_n = (2n+1)^{-1} \sum\_{i=n-d}^{n+d} x\_{i\prime} \tag{5}$$

where **x** = (*x*1, . . . , *xL<sup>x</sup>* ) is the input signal, **y** = (*y*1, . . . , *yL<sup>y</sup>* ) is the filtered signal, *d* is the window size, and *L<sup>x</sup>* and *L<sup>y</sup>* are respectively the input and filtered signal lengths. To accounting for a reduction of the computational resources usage, we experimentally define *d* = 11.

With the aim of verifying the data obtained by each sensor and validating the reliability thereof, samples obtained from the river are taken to the Environment Services Laboratory of the Technical

University of the North (Universidad Técnica del Norte (Universidad Técnica del Norte official web site: https://www.utn.edu.ec/web/uniportal/) from Ibarra-Ecuador, as they count on the technology and reagents to make comparison against the data obtained by the WSN. In this sense, following reliability criteria for each sensor, some recommended performance measures are considered, such as: (i) Accuracy: ability to provide the same reading by repeatedly performing the same experiment (standard deviation), (ii) Reproducibility: ability to reproduce the same results when modifying initial conditions of the experiment, and (iii) Stability: ability to produce the same output value in a long time. Overall obtained results are gathered in Table 1, which correspond to 10 tests over controlled environments to assess the data stability. As can be appreciated, the collected data from the sensors exhibit an error average of 5% in contrast to the those generated at the laboratory—such an error is acceptable enough for implementation purposes.


#### *3.3. Data Analysis Paradigm*

For a proper and wide data acquisition, we establish three node points in different locations, based on the population density of Ibarra, as follows: (i) *La Rinconada*, with low population and located at the river's beginning; (ii) *El Tejar*, with middle population rate and some wastewater discharged into the river; and (iii) *La Victoria*, with a larger population density and more discharge of pollutants from the city. Figure 3 shows the geographic locations of the nodes. Furthermore, we label each data from the nodes with a localization tag. For the data acquisition procedure, we design a collection protocol as follows: A schedule consisting in four collecting times is set, namely: in the morning, afternoon, night and early morning. Such a schedule is timed with Timer2, which is an Arduino internal clock. So, the system is timed for alerts at 08:00, 12:00, 17:00 and 00:00. On those times, the system records the sensor readings every 10 min for two hours (amounting to 6 samples per hour). Finally, these captured data are sent to the remote server through the GSM/GPS sensor. This collection protocol was performed during 3 months, generating an enough amount of information to be used in the subsequent data analysis stage.

Once the data are stored in an external server, a two-stages data analysis process is carried out: The first stage is the training set size reduction—via prototype selection—involving the least or no affectation to the intrinsic knowledge they hold. The second one is the classification task, in which the the algorithm that best fits the first stage while keeping a high accuracy is sought. Both stages are set and performed under low-computational cost criteria (given the device conditions). This process is carried out in order to be compiled within each WSN node (including both prototype selection, and classification). Then, system is able to make their own decisions based upon the reduced, stored dataset as well as the implemented classification algorithm. Therefore, on the one hand, the adaptability criterion required by an intelligent system is met, by making it able to be used anywhere on the river. On the other hand, the resulting system requires no re-run the data analysis process and thus it can be readily used by any system operator whom is not required to hold an expertise on embedded systems or data analysis, but only knowledge in water treatment itself.

**Figure 3.** Geographic location of the WSN nodes. Spots strategically selected to acquire data from, in order to encompass representative zones, as well as different types and levels of river pollution.

#### 3.3.1. Proposed Quality Measure: Quantitative Metric of Balance (QMB)

Algorithm analysis is an important part of designing thereof. Traditionally, the analysis of programming code or algorithms lies in applying theoretical and mathematical procedures. Indeed, when selecting supervised classification algorithms, efficient programs must be ensured to be created, as this translates into better power consumption and therefore battery life time usage. In this sense, the here-introduced Quantitative metric of balance (QMB) is aimed at quantifying how proper is the ratio between classifier performance and the data size reduction by the prototype selection stage. In this connection, the closer its value is to 100%, the better the ratio. As these three individual measures have an increasing nature, we multiply them to state a single value, namely, the rate of removed instances (*RI*) times the classification performance (*CP*), and divided by the response time of the classification algorithms (*RT*), as follows:

$$\text{QMB} = \frac{(RI \ast CP)}{RT} \ast 100\%. \tag{6}$$

Certainly, some classification criteria make use of mathematical functions or recursive functions of model adjustment that, when coded in a low-level language (assembler), generate response time delays, memory saturation and an excessive battery consumption. In this sense, the proposed QMB is aimed at penalizing the excessive computational cost in order to make it more feasible the implementation of data analysis algorithms into an embedded system. Besides, since it takes into consideration the number of removed training set instances to quantify the overall performance, this metric rewards the classification algorithm if it requires the least memory capacity when performing the decision-making procedures. When operating under real conditions, the system acquires the data from the sensors, filter the acquisition errors, make the decision through its compiled classification algorithm, and use the selection of prototypes to determine if this new reading improves the prediction ability of the system. If so, it is added into the training matrix otherwise it is only sent to the external server for visualization purposes.

Figure 4 shows the proposed data analysis scheme.

**Figure 4.** Data analysis scheme including prototype selection and classification stages.

#### 3.3.2. Prototype Selection

Since WSN systems have limited computational resources, its battery consumption is directly related to the amount of data to be processed, and therefore the implementation of machine learning algorithms into thereof is limited. In this connection, the prototype selection (PS) techniques may take place by reducing the training matrix size, while utmost maintaining as good classification performance as that obtained when considering the original size. Regarding PS algorithm designing, technical literature reports at least three main methods (namely, compensation-based, edition-based, and hybrid) [27]. As have been mentioned throughout this paper, the whole process is carried out in such manner that the prototype selection results (reduced data matrices) can be stored directly into every WSN node.

In this work, in order to account for an enough coverage, we have chosen three representative algorithms of each method, as follows:


#### 3.3.3. Classification Algorithms

Classification algorithms can learn based on different criteria, having each of them representative algorithms [27]. Herein, we consider four criteria and their respective representative algorithm, namely:


Given that the four aforementioned criteria are essentially different, a comparison of individual performances is necessary to identify the one(s) best fitting the nature of data and classification task. As well, it is of crucial interest to measuring the computational cost that each algorithm involves to be further implemented within the WSN node.

The database—obtained according to the pollution level—has been divided regarding the information acquired by the WSN nodes into 3 types (being our training labels): high, medium and low contamination. Therefore, if the system is located at different spots along the river, it can generate a map of the pollution status and estimate the river's course. Alternatively, if it is located statically, the system can determine, in hours, how the level of contamination varies with respect to the time of day.

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

In order to evaluate the behavior of each stage, we firstly discuss the data reduction in the training matrix. Subsequently, we show the outcome of our proposed analysis scheme, namely, the performance analysis using our defined metric (QMB) for determining the ideal algorithms for its implementation in the WSN nodes. Finally, we present the results of the final implementation of the system and the tests in real environments.

#### *4.1. Data Reduction*

The sensors were acquiring data during the months of July, August and September on random days. As a result, we obtained the data matrix called **<sup>Y</sup>** <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* , where *m* is the number of instances, and *<sup>n</sup>* the number of measured variables (sensors). While, **<sup>L</sup>** <sup>∈</sup> <sup>R</sup>*m*×<sup>1</sup> is the tag vector. Thus, we have that *m* = 507, and *n* = 4. With these data, we implemented the PS algorithms in order to reduce the training matrix and processing time. In addition, to validate the classification criteria, we retained 20% of the **<sup>Y</sup>** matrix for performance testing. In succession, the matrix for the data scheme is **<sup>X</sup>** <sup>∈</sup> <sup>R</sup>*p*×*<sup>n</sup>* , where *p* = 405. Table 2 shows the summary of the PS algorithms results and find a new reduced data matrix **Z**.

Accordingly, we have selected the CNN, DROP1 and DROP3 algorithms as they reach the highest percentages of reduction in the database. Figure 5 shows scatter plots of the initial data set and the reduced versions generated by CNN, DROP1 and DROP3.

**Figure 5.** 2D scatter plots of resulting data matrices *Z* of the chosen prototype selection algorithms. (**a**) Data matrix *X*, (**b**) CNN, (**c**) DROP1, (**d**) DROP3.


**Table 2.** Analysis of PS algorithms in relation to optimization embedded computational resources.

#### *4.2. Classification Performance*

With the reduced data sets, we compared the classification performance using the aforementioned algorithms. Table 3 summarizes the results of the classifiers with cross-validation with ten random folds.


**Table 3.** Classifier's metrics.

To graphically appreciate the results of the whole data processing scheme, just as done in previous works [11,14], we use the principal component analysis conventional algorithm as a dimensionality reduction approach to represent the original data over a lower-dimensional domain. Figure 6 presents scatter plots regarding the two first principal components to depict the decision borders generated by every considered classifier. This process is carried out for demonstration purposes in order to know the algorithms' ability to differentiate each label in an understandable way for the human being perception (visual-type in this case).

**Figure 6.** Decision borders for each classifier. Original data are embedded into a bi-dimensional space using PCA to graphically depict the classification ability of the considered algorithms. (**a**) k-NN, (**b**) Bayesian classifier, (**c**) SVM (Sigmoid kernel), (**d**) SVM (Polynomial kernel).

Numerical results of the joint performance of the prototype selection and data classification are summarized in Table 4.


**Table 4.** QMB analysis for every classifier along with the previously identified PS algorithms.

**Discussion on performance measures:** As can be seen in the Table 3, VSM reaches the best classification performance based on the considered metrics (100%). Nonetheless, its algorithm involves mathematical functions (known as kernel functions), which are not able to readily processed in a WSN. In this connection, the proposed QBM allows for warning about this computational cost in relation to the amount of data used to train the classification algorithm and the system response time when assigning the corresponding label to a new data from the sensors. This can be appreciated from the fact that by reducing the training matrix its performance decreases significantly. The same occurs

for all the considered algorithms excepting for k-NN, whose distance-based nature is non-expensive in terms of computational cost. Furthermore, by using a reduced data matrix, k-NN considerably maintains its performance. Furthermore, it is clearly noted that DROP1 is the best-suited algorithm for prototype selection although its computational cost is very high. Hence, given the design settings and the embedded systems conditions, CNN is preferred and therefore selected as the algorithm for prototype selection, while k-NN is considered as the selected classification algorithm reaching a performance of 90.6% and a QBM value of 72.85%.

#### *4.3. Implementation and Testing*

Figure 7 depicts the functional architecture of the nodes using the proper, selected prototype selection algorithms, which are to be compiled within thereof. As can be appreciated, each node holds the data-acquisition sensor set. The data analysis and processing is as follows: The raw data is first filtered by using the Moving average filter, which, in this case, is enough to remove the components (artifacts) related to reading errors and noise. Subsequently, data are classified by the algorithm k-NN, which assigns a label and decides about the predicted level of water contamination according to the training database and following a distance-based, majority-vote-driven approach. Then, data undergo an additional processing via CNN to determine whether the training database can be improved by removing instances exhibiting negligible relevance regarding either the subsequent classification task or the intrinsic knowledge they may hold. Finally, the output information is converted into a character string together with its label to be sent by the GSM network to the external server and display the data obtained from each sensor and the decision made. It is worth highlighting that the node to be monitored can be selected through the interface.

**Figure 7.** WSN node functional architecture incorporating the workflow of the in-situ data analysis and processing and mainly consisting in filtering, prototype selection and classification.

In the overall work-flow of our approach, the need for using an external sever lies in the fact that optimizing resource consumption at the in-situ analysis (directly on WSN Nodes) entails performing offline data processing tasks, mainly, at three specific points. The first one is when collecting data from each WSN node, being its main function the storing of such information (which—at this extent—corresponds to the outcomes of reading-errors-filtering stage produced by the moving average filter). The second one is the offline, exhaustive running, and comparison of classification algorithms to identify the ones reaching a good compromise between accuracy and computational cost, and therefore, being adequate to be directly implemented into the WSN nodes. Finally, as the third point, the server is used for information visualization purposes (displaying numerically and graphically the acquired data, the decision (classification) made by each node and the river pollution historical). This information is also stored in the server. Of course, those algorithms identified as adequate ones at the second point are the ones that are finally incorporated into the WSN nodes.

Once performed the data analysis procedures, we integrate all sensors into a PCB board incorporating an Arduino Uno as a processor unit. A view of the developed WSN node can be seen in Figure 8.

**Figure 8.** View of the WSN node including the four considered sensors and the processor.

The developed WSN has a considerably high operating consumption for a LiPo-type battery. To increase the life time of both the system and the battery, energy saving modes are used inside the Arduino board that handles the sensor activation. To enable such modes, we consider the use of timers, which work as an internal clock determining the data-acquisition-and-sending timing, and therefore limit current consumption. Hence the power consumption of every single sensor and the processor should be considered. In normal operation conditions, the total electric current consumption (considering all the sensors) amounts to 110 mA, while the GPS-GSM module and the Arduino require 40 mA and 45 mA, respectively. Meanwhile, when the battery saving system is enabled, the sensors and the GPS module are not is used, and thus only the Arduino works and is fed with 15mA. As stated in [28], the following equation relates the battery life time with the total power consumption (*P*):

$$P = \frac{\left(T\_{on} \ast I\_{on}\right) + \left(T\_{sleep} \ast I\_{sleep}\right)}{T\_{on} + T\_{sleep}}\,\prime\tag{7}$$

where *Ton* , *Tsleep*, *Ion*, and *Isleep* stand respectively for Normal Consumption Time, Sleep Consumption Time, Current Consumption at Normal Conditions, and Current Intensity Sleeping Consumption.

As explained in Section 3.3, the system is on during 10 min and then remains in battery saving mode. As a result, the system consumes 78.45 mA per hour. If the used battery is 5 volts at 1000 mA, the system can work continuously for 12.73 h. However, the system is activated only four times per day (early morning, mid-morning, afternoon and night), that is, it only works for 4 h a day. As a result, the system can remain for at least 3 days with no requiring battery manager support. As an advantageous aspect of our system we may say that, when implemented with a solar panel powering the battery, there is experimental evidence that it can work up to 4 months with no discharging or critical battery issues.

Subsequently, over the implemented system, we store the training dataset obtained after running the CNN algorithm, which is to denoted **<sup>Z</sup>** <sup>∈</sup> <sup>R</sup>*s*×*<sup>n</sup>* , by setting the number of prototypes as *s* = 11. At this extent, CNN algorithm is considered as an recommendable approach, since its execution time is the least while its ability to reduce the dataset instances is proper enough. Consequently, if the system requires to be reconfigured to train the classification algorithm model, the CNN algorithm can be compiled readily on the WSN network with no entailing extra battery consumption or diminishing the system performance. Then, we implemented the Bayesian classifier so that it can make system decisions concerning the tag assigned by location. Thus, we can determine the contamination levels (high, medium, low) using the nodes along the river. Since the system is intended to be waterproof, we use a river buoy to keep the system afloat. At its upper part, we install the solar panel and the GPS-GSM communication antenna. Furthermore, the nodes are anchored using an ironwork attached to the river stones, as shown in Figure 9.

**Figure 9.** Anchored node acquiring and sending data to interface. (**a**) Simulation. (**b**) Real conditions.

Besides, for displaying purposes, we develop a monitoring interface in Processing using a local server that downloads and visualizes the information from the server. In this interface, we show the status of each sensor, the node location, and the level of contamination of the river. Figure 10 summarizes both the sensor testing and the visual interface with the decision taken.

(**a**) (**b**)

**Figure 10.** System testing and visual interface. (**a**) Testing embedded system developed in the rural sector. (**b**) Testing embedded system developed in urban sector. (**c**) Visual interface showing low level of contamination. (**d**) Visual interface showing high level of contamination.

For a more extensive analysis, we move the nodes throughout the river to assign a color label, based on the contamination level, as follows: red refers to high contamination, yellow to medium, and green to null or low pollution. Accordingly, Figure 11 shoes the contamination levels along the case-study river. As a relevant result, we identify that at the Campiña church zone there is already a high level of pollution.

**Figure 11.** Tahuando river conditions along its stream bed.

Finally, with all nodes running, we daily capture data to observe the maximum values, in order to detect the hours of the day with highest contamination, which are in line with the human's work schedules. Figure 12 shows the pH, Temperature, and NTU values registered by the sensors during a whole day.

It is worth mentioning that our system may exhibit failures regarding the loss of signal from the GPS-GSM module when restarting it to carry out the data acquisition. To overcome this drawback, we follow a heuristic sensor calibration procedure as follows: On one hand, when activated, the system first turns on the GPS-GSM module so that there would be enough time to re-link to the GSM network and send back a status indicator signal. On the other hand, the length of the cables connected to the sensors was initially very long. This caused that when the volume of water decreased, cables descended to the bottom of the river and got brushed against stones. Consequently, since the length of the system-incorporated sensor is between 2 and 5 cm, an excessive wear on the sensors is induced. To cope with this issue, we search for and identify points where the river depth is the least possible varying, and is not prone to water stagnation.

**Figure 12.** Sensor-generated data acquired per hour during a day.

#### **5. Final Remarks**

In this work, we present the complete design and validation of an intelligent wireless sensor network (WSN) system to measure the contamination levels of a river. Particularly, the Tahuando River is of interest. Broadly speaking, the proposed system involves two stages: electronic device implementation, and data analysis.

For the electronic design, since the case-study river may have high levels of pollution, as well as it may occur significant variations depending on the hours of the day, and zones of its route, we implement several WSN nodes for acquiring the river's conditions information by covering a meaningful zone and within a wide enough range of time. In this sense, we both calibrate and tune the sensors for a correct data collection. Additionally, we experimentally demonstrate that our data reading schedules were adequate for detecting higher pollution hours. Furthermore, we highlight that the river buoys is a key element to meet the node's permeability requirements as well as to enable the proper functioning of each WSN node.

Regarding the proposed data analysis scheme, we demonstrate that a classifier together with a prototype selection is suitable for a WSN-based water-quality monitoring system. It is reached a good trade-off between the computational resource usage (as the training matrix size is reduced to meet the system operation conditions), and the classification performance at detecting the pollution levels along the river. In addition, given the network coverage, the proposed system is able to send information from the WSN node to the server. Therefore, the filtered data can be visualized in an interface, and an in-situ analysis becomes possible. It is important to mention that the server is only for data visualization purposes and does not have the implementation of machine learning algorithms.

As a future work, the battery life is to be more carefully considered by exploring both different methods of extending its duration and alternatives sources of energy to supply the nodes (i.e., using the water flow to generate energy). A large number of nodes and wider coverage (located at different water resources around the province of Imbabura, Ecuador) is highly desirable for further In addition, we are intended to a seek for alternatives to mitigate system affectations due to disturbances caused by the presence of unexpected individuals (either people or animals), as so far our readily solution has been to locating the system in a hardly visible and difficult-to-access spot.

**Author Contributions:** P.D.R.-M.: Conceptualization, methodology, software, formal analysis, investigation, writing—original draft preparation, visualization, resources, V.F.L.-B.: investigation, supervision, project administration, J.A.R.: formal analysis, writing—review, visualization, D.H.P.-O.: validation, writing—review, methodology supervision, project administration. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This work is supported by the Smart Data Analysis Systems Group—SDAS Research Group (http://sdas-group.com).

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

#### **References**


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

## *Article* **Groundwater Potential Mapping Using Remote Sensing and GIS-Based Machine Learning Techniques**

#### **Sunmin Lee 1,2 , Yunjung Hyun <sup>3</sup> , Saro Lee 4,5 and Moung-Jin Lee 6,\***


Received: 17 February 2020; Accepted: 7 April 2020; Published: 8 April 2020

**Abstract:** Adequate groundwater development for the rural population is essential because groundwater is an important source of drinking water and agricultural water. In this study, ensemble models of decision tree-based machine learning algorithms were used with geographic information system (GIS) to map and test groundwater yield potential in Yangpyeong-gun, South Korea. Groundwater control factors derived from remote sensing data were used for mapping, including nine topographic factors, two hydrological factors, forest type, soil material, land use, and two geological factors. A total of 53 well locations with both specific capacity (SPC) data and transmissivity (T) data were selected and randomly divided into two classes for model training (70%) and testing (30%). First, the frequency ratio (FR) was calculated for SPC and T, and then the boosted classification tree (BCT) method of the machine learning model was applied. In addition, an ensemble model, FR-BCT, was applied to generate and compare groundwater potential maps. Model performance was evaluated using the receiver operating characteristic (ROC) method. To test the model, the area under the ROC curve was calculated; the curve for the predicted dataset of SPC showed values of 80.48% and 87.75% for the BCT and FR-BCT models, respectively. The accuracy rates from T were 72.27% and 81.49% for the BCT and FR-BCT models, respectively. Both the BCT and FR-BCT models measured the contributions of individual groundwater control factors, which showed that soil was the most influential factor. The machine learning techniques used in this study showed effective modeling of groundwater potential in areas where data are relatively scarce. The results of this study may be used for sustainable development of groundwater resources by identifying areas of high groundwater potential.

**Keywords:** groundwater potential; specific capacity; machine learning; boosted tree; ensemble models

#### **1. Introduction**

Because groundwater has less exposure to pollution than surface water, it is considered a valuable natural resource for agriculture in many communities [1]. Especially during the drought season, a continuous supply of groundwater is important in agricultural areas. The study area in this investigation, Gyeonggi-do, has recently suffered from damage to agricultural land due to increasing drought. In 2018, widespread damage to crops due to heat waves and drought continued throughout the year, and the average storage rate in 339 reservoirs in Gyeonggi-do was 59% of capacity, which was only 76% of the normal level [2].

Groundwater is a good water resource because it can stably supply the required amount of high-quality water; thus, appropriate water conservation plans are essential for the sustainable use of groundwater [3]. In many areas, the main causes of groundwater depletion are excessive groundwater extraction and unsuitable aquifer recharge [4]. Therefore, accurate estimation and prediction of groundwater recharge should be carried out to support efficient use and systematic management of groundwater resources. From this perspective, groundwater potential mapping using yield data is important. Yield data include extraction volume and the velocity of groundwater at various measurement points. Groundwater yield depends on geological, topographic, and anthropogenic factors specific to the area, and is also related to groundwater potential [5].

In practical terms, groundwater is less accessible than surface water. Groundwater can be presumed by detecting gravity anomalies such as Gravity Recovery and Climate Experiment (GRACE) [6–8]; however, a local groundwater potential map is essential for regional management of groundwater. Thus, studies on the distribution and prediction of groundwater resources have been limited to local scales based on data obtained from point measurements (e.g., meteorological stations, flow measurement points, and groundwater level monitors) [9,10]. In recent years, areal distribution analysis data obtained through remote sensing have been used for global prediction of the water resource distribution in combination with various machine learning techniques, albeit with high uncertainty. To overcome the limitation of groundwater resource surveys based on local information, these data can be converted into global distribution data using satellite imagery. Remote sensing generally produces data in the form of grids or regions, which can be converted into distribution patterns through various processing methods such as machine learning algorithms. By applying the characteristics of remote sensing data to groundwater resources, point-based groundwater hydrological modeling can be extended to the global scale. Therefore, using existing groundwater yield data, it is possible to make regional and local predictions with remote sensing-based methods.

For groundwater potential mapping, a variety of techniques have been applied, including direct drilling for hydrological testing and geophysical models [11,12]. Such methods are suitable for identifying the hydrological characteristics of groundwater, but have high costs in time and money [13,14]. In recent years, studies related to groundwater potential have been conducted using machine learning models with available historical data on groundwater wells with geographic information systems (GIS) [15,16]. GIS technologies have been used for quantitative analysis of spatial distributions in environmental, geological, and hydrological studies [17–19]. One limitation of data-based analysis of groundwater is insufficient availability of data for analysis [20]; groundwater yield varies with hydrological conditions and recharge sources, which have been measured in a limited number of groundwater wells [21]. Therefore, using various models to predict groundwater yield accurately and identifying the optimal model for water resource evaluation in a given region are essential to effective water resource management.

For this reason, studies related to groundwater potential mapping with various data models have become increasingly common [22–24]. Numerous factors that affect groundwater potential have been proposed based on various data modeling methodologies, including statistical models, probabilistic models, machine learning models, and data mining models; yield and spring or well location data are also widely used as groundwater potential indicators. Due to the characteristics of remote sensing and groundwater, groundwater could be indirectly monitored by using remote sensing; much research has been conducted through thematic maps related to groundwater based on remote sensing data and groundwater potential was estimated by reducing the uncertainties [25–27].

The frequency ratio (FR) model is a representative statistical model applied to groundwater potential mapping [26,28,29]. The relationship between groundwater conditioning factors and groundwater potential could be analyzed using basic statistical and probabilistic models, including FR, weight of evidence [30], evidential belief function [31], and logistic regression [32] models. Furthermore, the recent exponential increase in available data has led to identification of data types and data processing techniques that can support decision-making. Several studies in this area have applied machine learning methods such as machine learning models, while artificial neural networks [33] and support vector machines [34] have been widely applied to groundwater potential mapping. Some studies have also used analytical hierarchy methods, which are expertise-based methods requiring a deep understanding of the study area [35,36]. Recently, hybrid and ensemble models that combine or develop existing methodologies have been applied for groundwater potential mapping [37–39]. This paper also uses a hybrid methodology in this respect.

When performing groundwater potential mapping through modeling, the results show poor generalizability without proper training samples. In such cases, the accuracy for training data is high but the testing results show significantly lower accuracy. To overcome the lack of data, robust models built upon basic models have recently been developed and compared [40]. Typically, an ensemble model using learner sequences is developed; voting, bagging, and adaptive boosting are representative ensemble methods that can be applied to various base learners [41]. In this way, unlabeled cases are identified via self-learning by combining information from labeled cases so that the labeled training set is magnified in each iteration until the entire dataset is labeled. This method, which was applied in the present study, could be effective for data-scarce areas because it allows modeling using less data than other approaches.

Previous studies conducted on groundwater recharge and yield have used enough field survey data targeted at adjacent areas. However, these studies are subordinate to field surveys and are not intended to reduce spatial uncertainty on groundwater. Therefore, the purpose of this study was to map and test groundwater yield potential in Yangpyeong-gun, South Korea, using spatial data analysis in a GIS environment. This study processed and analyzed officially published groundwater yield data using remote sensing and GIS to reduce the uncertainty of the data itself. In addition, one of the latest machine learning models, boosted tree method, was applied to predict large areas of low uncertainty using pumping test data from 53 wells; groundwater yield potential is the major issue of this study. The results of this study could provide a scientific basis for efficient use and systematic management of groundwater resources.

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

South Korea consists of eight administrative districts, labeled '-do', which are made up of local administrative districts, labeled with '-si', '-gun', and '-gu'. The study area, Yangpyeong-gun, is located about 50 km from Seoul, in the northeastern part of Gyeonggi-do (Figure 1). Yangpyeong-gun is surrounded by Hongcheon-gun in Gangwon-do to the northeast, Hoengseong-gun in Gangwon-do to the east, Wonju-si in Gangwon-do to the southeast, and Gapyeong-gun to the north. Yangpyeong-gun contains rugged mountainous areas such as Yongmunsan (1157 m), Bongmyun (856 m), and Baekunbong (940 m), and the Namhan River flows from the south to the northwest of the district. About 90% of the total area of Yangpyeong-gun is a green zone covering the protected headwater area of the Han River; this area has a well-preserved and clean natural environment due to legal and institutional regulations [42].

Yangpyeong-gun covers approximately 878 km<sup>2</sup> , and the amount of groundwater used in this area is 41,503,946 m<sup>3</sup> /year. The groundwater use per unit area is 47,258 m<sup>3</sup> /km<sup>2</sup> annually and 129 m<sup>3</sup> /km<sup>2</sup> daily [43]. Groundwater in Gyeonggi-do is used primarily for agricultural purposes in numerous agricultural areas, including Anseong-si, Yangpyeong-gun, Icheon-si, and Yeoju-si. Among all districts in South Korea, Yangpyeong-gun (10,725) has the second highest number of groundwater facilities for agricultural use after Anseong-si [43].

In Yangpyeong-gun, a preliminary survey of available groundwater resources was conducted from December 2017 to June 2018 for drought response and to prevent unplanned development. Among 35 districts prone to drought, 25 were selected based on the feasibility of surveying the target district and the response rate of residents. A resistivity survey (vertical and dipole survey) was conducted to select locations for large-scale groundwater storage. from December 2017 to June 2018 for drought response and to prevent unplanned development. Among 35 districts prone to drought, 25 were selected based on the feasibility of surveying the target district and the response rate of residents. A resistivity survey (vertical and dipole survey) was conducted to select locations for large-scale groundwater storage. In Yangpyeong-gun, Gyeonggi-do, Kyonggi massif metamorphic rocks of Precambrian age and

In Yangpyeong-gun, a preliminary survey of available groundwater resources was conducted

In Yangpyeong-gun, Gyeonggi-do, Kyonggi massif metamorphic rocks of Precambrian age and an intrusive body of Mesozoic Triassic gabbro and syenite are found. Precambrian Kyonggi massif metamorphic rocks consist of the Paleozoic sequence of Yongmunsan and unconformity of Jang-Rak. The main constituent rocks are banded gneiss, migmatitic gneiss, augen gneiss, mica schist, and quartzite. These rocks underwent metamorphism in the Paleozoic and Mesozoic Triassic, when the landmasses of North China and South China collided. an intrusive body of Mesozoic Triassic gabbro and syenite are found. Precambrian Kyonggi massif metamorphic rocks consist of the Paleozoic sequence of Yongmunsan and unconformity of Jang-Rak. The main constituent rocks are banded gneiss, migmatitic gneiss, augen gneiss, mica schist, and quartzite. These rocks underwent metamorphism in the Paleozoic and Mesozoic Triassic, when the landmasses of North China and South China collided. Groundwater development requires continuous management for sustainable supply of water

Groundwater development requires continuous management for sustainable supply of water rather than short-term measures at the time of drought. Specifically, preliminary investigation is needed in drought-prone areas and areas of high importance for agricultural water usage in Gyeonggi-do. To mount an effective response to agricultural drought, a groundwater management plan that ensures sustainable use of agricultural groundwater prior to drought is needed [44]. In this study, continuous groundwater potential data in the study area were used as primary data for a groundwater abundance survey, and could further be used to establish a groundwater development plan. rather than short-term measures at the time of drought. Specifically, preliminary investigation is needed in drought-prone areas and areas of high importance for agricultural water usage in Gyeonggi-do. To mount an effective response to agricultural drought, a groundwater management plan that ensures sustainable use of agricultural groundwater prior to drought is needed [44]. In this study, continuous groundwater potential data in the study area were used as primary data for a groundwater abundance survey, and could further be used to establish a groundwater development plan.

**Figure 1.** Study area. **Figure 1.** Study area.

#### **3. Data 3. Data**

#### *3.1. Groundwater Potential Analysis Based on Remote Sensing Data*

*3.1. Groundwater Potential Analysis Based on Remote Sensing Data*  Various thematic maps constructed using remote sensing source data were applied to machine learning techniques in this study. Recently, high-resolution aerial photographs were used to produce thematic maps of spatial data. Topographic maps were produced through numerical Various thematic maps constructed using remote sensing source data were applied to machine learning techniques in this study. Recently, high-resolution aerial photographs were used to produce thematic maps of spatial data. Topographic maps were produced through numerical mapping using aerial photographs taken in 2006, with corrections and supplemental data collected through field surveys. Forest and soil maps were also constructed using spatial data generated through field surveys along with aerial photography. For land use maps, aerial photographs taken in 2012 were classified using image classification techniques, and their quality was verified using additional high-resolution satellite images from KOMPSAT-2 and KOMPSAT-3 as well as digital topographic maps. Meanwhile, geological maps were produced from field surveys and historical records using base maps generated from aerial photographs. Groundwater yield is a measure of groundwater pumping capacity, which could be stored in aquifers. In this study, groundwater yield potential modeling using machine learning was performed with spatial data generated via remote sensing and GIS such as soil, land cover, and geological maps, as described above.

#### *3.2. Groundwater Well Data from in Situ Sampling*

Groundwater pumped from wells in the study area is used mainly for agricultural purposes and domestic drinking water. Groundwater well data were collected for specific capacity (SPC) (53 wells) and transmissivity (T) (53 wells) from the basic survey report of Yangpyeong-gun [45]. The main use of the groundwater in this area is agricultural, so groundwater surveys are conducted between spring and summer, and our data was obtained between June and August. In the training and testing subsets, yield values above 3.8 and 3.42 (30 m<sup>3</sup> /h) above the median value were considered for yields based on the dependent variables of SPC and T, respectively, which are two different indexes measured in different ways. Groundwater pumping test data used in this study were generated and published from the national groundwater observation and survey data by local governments conducted by Korea Water Resources Corporation (K-water).

SPC data include geographic location coordinates of individual wells and groundwater yield derived from pumping tests. SPC often indicates well performance, because it refers to the amount of water that a well can produce per unit of drawdown. SPC is calculated by dividing the pumping discharge by the drawdown, in units of liters per minute (LPM) per meter, as follows:

$$\text{SPC} = \frac{\text{Q}}{\text{S}} \tag{1}$$

where Q is discharge (unit: LPM) and S is drawdown (unit: m). A low SPC value indicates that more energy is required for pumping. During a drawdown test to determine SPC, pumping should be maintained at a constant speed for a certain period of time, at least 24 h, with little change in drawdown. SPC data acquired during the pumping test can be used to estimate T and identify potential aquifer issues.

T represents the flow rate under a unit hydraulic gradient through a unit width of aquifer of a certain thickness [46]. Hydraulic conductivity (K) is a measure of the water transmission capacity of an aquifer. T of an aquifer is equal to the hydraulic conductivity multiplied by the thickness of the aquifer.

$$\mathbf{K} \mathsf{V}(\mathbf{x}, \mathbf{y}) = \frac{1}{\mathsf{b}} \int\_{0}^{\mathsf{b}} \mathsf{K}(\mathsf{x}, \mathbf{y}, \mathbf{z}) d\mathbf{z} \tag{2}$$

$$\mathbf{T} = \mathbf{K}\mathbf{b}\_{\prime} \tag{3}$$

where T is transmissivity, K is hydraulic conductivity, and b is aquifer thickness. Less drawdown and a thicker aquifer lead to higher T values. It is possible to estimate the amount of water flowing through the unit thickness of the aquifer by combining Equation (3) with Darcy's law.

SPC and T data were separately applied to the FR, boosted tree (BT), and ensemble models in this study; both SPC and T are used in this study in order to consider various aspects of groundwater. The locations of groundwater wells in the study area are shown in Figure 2. Yield data were randomly divided into a training data subset (70%) and a testing data subset (30%), as is the usual division in machine learning methodologies [16,47]. In the training data subset, 37 wells each were represented in SPC and T data, respectively; 16 wells were used to test the models.

each were represented in SPC and T data, respectively; 16 wells were used to test the models.

SPC and T data were separately applied to the FR, boosted tree (BT), and ensemble models in this study; both SPC and T are used in this study in order to consider various aspects of groundwater. The locations of groundwater wells in the study area are shown in Figure 2. Yield data were randomly divided into a training data subset (70%) and a testing data subset (30%), as is

**Figure 2.** Locations of groundwater wells sampled for: (**a**) SPC and (**b**) T data. **Figure 2.** Locations of groundwater wells sampled for: (**a**) SPC and (**b**) T data.

#### *3.3. Groundwater Conditioning Factors 3.3. Groundwater Conditioning Factors*

Various groundwater conditioning factors were used for groundwater potential modeling in this study (Table 1). Topographical, geological, hydrological, and land cover factors are commonly applied to predict groundwater yield potential. Conditioning factors should be considered depending on regional characteristics. For this reason, the correlation between the factors and groundwater potential were analyzed preferentially through the frequency ratio model and the factors were selected; groundwater potential was estimated using 16 factors in this study. The 16 conditioning factors were constructed into a groundwater inventory, including nine topographic factors (convergence index, convexity, mass balance index (MBI), slope angle, slope height, topographic texture, topographic position index (TPI), topographic ruggedness index (TRI), and valley depth), two hydrological factors (flow path length, and slope length and steepness (LS)), forest type, soil material, land use, and two geological factors (lithology and distance from fault) (Figures 3,4). The conditioning factors were calculated and prepared using ArcGIS 10.3 software (ESRI, Redlands, CA, USA). Each dataset was converted into a grid format with 30-m spatial resolution for use in the groundwater inventory of the study area. Various groundwater conditioning factors were used for groundwater potential modeling in this study (Table 1). Topographical, geological, hydrological, and land cover factors are commonly applied to predict groundwater yield potential. Conditioning factors should be considered depending on regional characteristics. For this reason, the correlation between the factors and groundwater potential were analyzed preferentially through the frequency ratio model and the factors were selected; groundwater potential was estimated using 16 factors in this study. The 16 conditioning factors were constructed into a groundwater inventory, including nine topographic factors (convergence index, convexity, mass balance index (MBI), slope angle, slope height, topographic texture, topographic position index (TPI), topographic ruggedness index (TRI), and valley depth), two hydrological factors (flow path length, and slope length and steepness (LS)), forest type, soil material, land use, and two geological factors (lithology and distance from fault) (Figures 3 and 4). The conditioning factors were calculated and prepared using ArcGIS 10.3 software (ESRI, Redlands, CA, USA). Each dataset was converted into a grid format with 30-m spatial resolution for use in the groundwater inventory of the study area.

Topographic factors were calculated from a 1:5000 scale topographic map provided by the Korean National Geographic Information Institute. Spatial data, such as location and topography, were structured using ground control point measurements taken from digital aerial photographs and ground surveys. Aerial photographs were analyzed through numerical mapping, and further calibration was carried out through field surveys to create the topographic map. A digital elevation model (DEM) was first generated from the topographic map and then used to derive topographic factors, including convergence index, convexity, MBI, slope angle, slope height, topographic texture, TPI, TRI, and valley depth. Slope factor impacts groundwater recharge, with gentle slope areas having relatively high percolation and low surface runoff rates and steep areas having high surface runoff [48]. Soil moisture content is also related to slope, which affects precipitation direction [49]. Slope angle is strongly related to groundwater potential; therefore, groundwater-related topographic factors derived from DEM data with SAGA-GIS software [50] were used for modeling. Acceleration and deceleration, as well as flow convergence and divergence of flow, are mainly affected by the curvature of the area [51]. The hydrological factors flow path and LS factor were considered conditioning factors for hydrological features. Topographic factors were calculated from a 1:5000 scale topographic map provided by the Korean National Geographic Information Institute. Spatial data, such as location and topography, were structured using ground control point measurements taken from digital aerial photographs and ground surveys. Aerial photographs were analyzed through numerical mapping, and further calibration was carried out through field surveys to create the topographic map. A digital elevation model (DEM) was first generated from the topographic map and then used to derive topographic factors, including convergence index, convexity, MBI, slope angle, slope height, topographic texture, TPI, TRI, and valley depth. Slope factor impacts groundwater recharge, with gentle slope areas having relatively high percolation and low surface runoff rates and steep areas having high surface runoff [48]. Soil moisture content is also related to slope, which affects precipitation direction [49]. Slope angle is strongly related to groundwater potential; therefore, groundwater-related topographic factors derived from DEM data with SAGA-GIS software [50] were used for modeling. Acceleration and deceleration, as well as flow convergence and divergence of flow, are mainly affected by the curvature of the area [51]. The hydrological factors flow path and LS factor were considered conditioning factors for hydrological features.

A forest map was also used, which was generated from field investigations and interpretation of aerial photographs. To construct the forest map, the near-infrared band was used for image analysis, in addition to the red-green-blue image. Moreover, soil material characteristics can impact the rate of surface water penetration into aquifers, which drives groundwater potential [52]. The soil material factor was extracted from a soil map published by the National Institute of Agricultural Sciences at 1:25,000 scale. Similarly, land cover has an impact on soil conditions such that storage and movement

of groundwater change when land cover changes; the land use factor was extracted from a digital land cover map provided by the Korea Ministry of Environment at 1:25,000 scale. Land use maps were classified into 22 medium-level categories through application of automatic image classification to aerial photographs, and the accuracy was enhanced using additional high-resolution satellite images from KOMPSAT-2 and 3. The land cover map was reclassified into seven land cover categories: urban, farmland, forest, grassland, wetland, bare land, and water.

Geological factors, including lithology and distance from a fault, were also considered in relation to groundwater characteristics. The lithology factor was extracted from a digital geological map produced by the Korea Institute of Geoscience and Mineral Resources at 1:50,000 scale. The study area was composed of 22 lithological units differing in lithology type and geological age. Distance from a fault was also calculated based on the geological map.


#### **Table 1.** Data layers describing groundwater potential.

*Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 8 of 24

**Figure 3.** Groundwater conditioning factors I: (**a**) Convergence Index, (**b**) Convexity, (**c**) Mass balance index (MBI), (**d**) Slope angle, (**e**) Slope height, (**f**) Texture, (**g**) Topography position index (TPI) and (**h**) Topography ruggedness index (TRI).

**Figure 4.** Groundwater conditioning factors II: (**a**) Valley depth, (**b**) Flow path, (**c**) Slope length and steepness (LS) factor, (**d**) Forest type, (**e**) Soil, (**f**) Land cover, (**g**) Geology and (**h**) Distance from fault.

1

#### **4. Methodology**

To be more specific, the purpose of this study was to map and test groundwater yield potential in Yangpyeong-gun, South Korea, using spatial data analysis in a GIS environment. This was performed by four main steps: First, groundwater yield data of specific capacity (SPC) and transmissivity (T) collected from 53 well locations were used. For the training data, 70% of each groundwater yield dataset was selected randomly, and FR and boosted tree (BT) models with classification were applied to the groundwater inventory using Statistica software (Dell Software, Aliso Viejo, CA, USA). Second, the inventory was constructed from nine topographic factors, two hydrological factors, forest type, soil material, land use, and two geological factors. All factors used in this study were generated and processed from remote sensing-based data, such as aerial photographs or imagery from KOMPSAT-2 and -3. Third, this study involved probabilistic analysis of FR, and two machine learning models: the boosted classification tree (BCT) and FR-BCT ensemble models, which were applied to groundwater yield data. Comparative analysis was conducted to compare the models used in this study. Finally, to quantitatively evaluate the performance of the models, the receiver operating characteristics (ROC) and area under the curve (AUC) were used. The study was conducted, as shown in Figure 5.

**Figure 5.** Data flow of this study.

#### *4.1. Frequency Ratio (FR) Model*

FR is an effective stochastic method for evaluating the effects of various factors on the occurrence of a particular event [53]. Thus, the FR value represents the ratio of occurrence of a particular event to the area ratio for each class [54]. A larger FR value represents a stronger relationship between the probability of occurrence and the specific variable [55,56]. This method allows for the clear and simple analysis of the relationship of each factor to the event [57].

To carry out spatial FR analysis, factors related to groundwater potential were classified into ten classes. Among numerous available classification techniques, factors in this study were classified using the quantile technique, which divides classes into equal areas. FR values were calculated using training data for each factor. Each class of each modulator was weighted. Higher FR values represent a stronger relationship between the class of each factor and groundwater potential, whereas for lower FR values, the effect of the class of each factor on groundwater potential is small. If FR is greater than 1, the effect is significant; if FR is less than 1, the effect is not significant [56]. To construct a groundwater

potential map using FR to represent the relative magnitude of the groundwater potential, the FR values calculated for each factor were determined as follows:

$$FR = \frac{P\_{tm}}{P\_{total}}\tag{4}$$

where *Ptrn* is the ratio of the number of SPC data points above a certain level and *Ptotal* indicates the ratio of the number of pixels in a certain class to the total number of pixels in the study area. A greater FR value for potential indicates higher groundwater potential; a lower value indicates a lower groundwater potential. In this study, FR values for each conditioning factor were used to weight the ensemble FR-BCT model.

#### *4.2. Boosted Classification Tree*

In recent years, decision tree models have been used in various fields as a machine learning method [58], including for groundwater potential mapping [52]. Decision tree models perform attribute tests on non-terminal nodes to represent the results on the terminal node, using a tree-like hierarchy that constructs a classification tree of a simple structure [59]. One of the benefits of this method is that the classification process can be graphically represented. However, the results cannot be formed into multiple outputs and the performance of the model depends on the type of data. Many algorithms have been developed from decision trees: classification and regression tree [60], chi-square automatic interaction detector decision tree [61], Iterative Dichotomiser 3 [62], and J48 (C4.5 decision tree) [63]. In addition, ensemble models using sequences of classifiers have been widely developed. Representative ensemble methods such as voting, bagging (sub-sampling), and boosting have also been applied to the decision tree method, including BT algorithms. Therefore, in this study, representative decision tree algorithms of BT models were used to compare the performance of each model's groundwater potential modeling and prediction accuracy.

The BT model is a tree-based machine learning model using the stochastic gradient boosting method. In the last few years, this algorithm has become one of the most powerful machine learning techniques used for prediction. In the BT algorithm, continuous or categorical input factors can be used for classification and regression problems [64].

The BT algorithm is implemented by applying a boosting method to the regression tree. The basic method involves calculating a simple tree sequence in which each successive tree is built against the prediction residual of the preceding tree. This method creates two trees of data for two samples at each split node. Even if the relationship between predictive and dependent variables is nonlinear, the weighting of such trees can support high accuracy of the predicted value. Thus, the gradient boosting method for weighted expansion of simple trees is one of the most common and powerful machine learning algorithms.

All machine learning algorithms are prone to overfitting, which involves a good fit for learning data but a lack of improvement in the predictability of each model. In other words, this is a common problem that applies to most algorithms used for predictive machine learning. A common solution to this problem is to evaluate the quality of the model fit by predicting observations from test samples of "used" data before evaluating each model [65,66]. The accuracy of each solution can be measured in this way to determine when the overflow occurred.

To overcome this difficulty, which is a major problem facing most machine learning algorithms used in predictive models, a specific approach was selected for the BT models. A continuous simple tree is generated using only subsamples selected randomly from the entire dataset. That is, each successive tree is created for the predicted residuals of an independently extracted random sample. Randomness can be added to any degree to protect against overfitting and can provide good predictability. Continuous boosting calculations for independently sampled input samples are known as probabilistic gradient boosting techniques.

#### *4.3. Ensemble Modelling*

Using the two methodologies described above, ensemble methods of FR and BCT were applied in this study. The probabilistic method FR was used to assess the impact of all types of regulatory factors and assign appropriate weights to each class according to their impact on groundwater yield. Using the FR method, individual weights were derived for each factor. Each conditioning coefficient was then reclassified using the derived weight values, and the reclassified dataset was analyzed using the BCT tree-based machine learning models. Finally, a groundwater potential map was constructed using the BCT and FR-BCT ensemble techniques for comparative analysis.

#### *4.4. Assessment on Model Performance*

The performance of groundwater potential classification was assessed using two statistical indicators: sensitivity and specificity. Sensitivity is the percentage of correctly classified pixels in areas with high groundwater potential; specificity is the percentage of pixels classified as having a low groundwater potential. Sensitivity and specificity are calculated as follows [67]:

$$Sensitivity = \frac{TP}{TP + FN} \tag{5}$$

$$Specificity = \frac{TN}{FP + TN'} \tag{6}$$

The numbers of correctly classified pixels are denoted as true positives (TP) and true negatives (TN). Conversely, the numbers of misclassified pixels are expressed as false positives (FP) and false negatives (FN).

In this study, ROC curves were used to evaluate the overall performance of the groundwater potential model. The ROC curve has been applied in various fields as a standard method for evaluating the general performance of a model [68]. This curve is plotted using sensitivity as the x-axis and 100 − specificity as the y-axis. The general performance of the model can be quantitatively assessed based on the AUC value, representing the area under the ROC curve. AUC values range from 0.5 to 1. A value of 0.5 represents a model with very low accuracy. In contrast, 1 represents a perfect model with the highest possible accuracy, and an AUC close to 1 indicates good performance. Generally, when the AUC value is greater than 0.8, the model shows adequate performance [69].

#### **5. Results**

#### *5.1. Results from the Frequency Ratio Model*

Table A1 presents the correlations of FR values between groundwater data (SPC or T) and groundwater conditioning factors derived from the FR model. The FR is a representative value of the statistical proportional position of well locations with SPC values above a specific level. Correlation between groundwater well data and each factor could be shown from the distribution of values biased according to each class. Areas with high FR values are of great importance for groundwater management because they have high groundwater potential. The characteristics of land cover in the area of this study are high in forest area and agricultural area, and relatively low in urban area. Although there are many groundwater wells in urban areas, the urban area is mixed with rural areas, so it requires a different approach from metropolis.

The topographic factor convexity showed a strong correlation with groundwater potential in the 1.1–43.19 class for FR values of over 1.89 and 2.63 for SPC and T, respectively. Similarly, MBI showed a high correlation with SPC (2.16) and T (1.84) in the -0.33 to 0.1 class. The highest FR values of 4.32 for SPC and 4.21 for T were observed when the slope angle was greater than 0 m and less than 0.05 m, indicating that this factor is strongly correlated with groundwater potential. FR values tended to decrease with increasing slope angle and slope height. For topographic texture, the 0.04–29.08 class

on the FR-BCT models.

**Factor** 

exhibited the highest FR values with SPC (2.97) and T (3.95). Low flow path values also led to FR values over 1, indicating that this factor was correlated with groundwater potential.

Among land cover types, urban area showed the strongest relationship with groundwater potential (SPC: 6.66; T: 7.92), followed by wetlands. These results could also be interpreted as showing that the use frequency of wells in urban areas is high. Meanwhile, distance from a fault had FR values of 2.16 for SPC and 3.16 for T in the 0–530.75 class. Among geological factors, alluvium showed a strong correlation with the groundwater data (SPC: 2.93; T: 3.80), followed by granite porphyry (SPC: 1.45; T: 1.01).

#### *5.2. Construction of Groundwater Potential Maps*

The groundwater potential map was modeled using training datasets of SPC and T. The performance of a groundwater potential model depends on the selection of factors. The groundwater potential map was constructed by training the groundwater potential model. First, a groundwater potential value was generated for each pixel in Yangpyeong-gun. Each pixel was indexed by its predicted groundwater potential value. The results of groundwater potential were reclassified using the 1.0 standard deviation method, which is based on the distribution of individual values in the results for each model. In the groundwater potential map, areas with high (low) groundwater potential are shaded red (blue) (Figure 6). All models showed similar distributions of groundwater potential, and the north, southwest, and southeast areas surrounding the central valley region of the study area all showed low potential. *Remote Sens.* **2019**, *11*, x FOR PEER REVIEW 14 of 24

**Figure 6.** Groundwater potential maps based on (**a**) boosted classification tree (BCT) and (**b**) frequency ratio (FR)-BCT models with specific capacity (SPC) data, and (**c**) BCT and (**d**) FR-BCT models with transmissivity (T) data. **Figure 6.** Groundwater potential maps based on (**a**) boosted classification tree (BCT) and (**b**) frequency ratio (FR)-BCT models with specific capacity (SPC) data, and (**c**) BCT and (**d**) FR-BCT models with transmissivity (T) data.

Furthermore, the predictor importance values of each factor were calculated from the BCT modeling results by summing the decreases in node-impurity values (Table 2). All predictor Furthermore, the predictor importance values of each factor were calculated from the BCT modeling results by summing the decreases in node-impurity values (Table 2). All predictor importance values

texture was the second most important factor in the BCT models, with values of 0.3101 and 0.4206, for SPC and T data, respectively. Meanwhile, FR-BCT models showed that forest type and land cover were the second strongest predictors, with importance values of 0.1704 and 0.2295 for SPC and T data, respectively. The importance of TPI, MBI, and valley depth were low in all FR models; convergence index, valley depth, and distance from a fault fell into the third lowest positions based

**Table 2.** Predictor importance values of each factor for the BCT and FR-BCT models.

Topography Convergence index 0.1689 0.0400 0.1518 0.0494 convexity 0.1285 0.1223 0.1913 0.1698 Mass balance index (MBI) 0.0566 0.1345 0.0703 0.1680 Slope angle 0.1909 0.1443 0.2734 0.1850 Slope height 0.1245 0.1393 0.1711 0.1750 Topographic texture 0.3101 0.1196 0.4206 0.1975 Topographic position index (TPI) 0.0387 0.0990 0.0565 0.1246 Topographic ruggedness index (TRI) 0.1967 0.1658 0.2917 0.2003

**Predictor Importance Values SPC T BCT FR-BCT BCT FR-BCT**  were scaled to a maximum of 1.0, as the value assigned to the largest sum among all factors, indicating the most strongly related factor, relatively. For both SPC and T, soil showed the highest predictor importance values in all models, with a value of 1.0. Topographic texture was the second most important factor in the BCT models, with values of 0.3101 and 0.4206, for SPC and T data, respectively. Meanwhile, FR-BCT models showed that forest type and land cover were the second strongest predictors, with importance values of 0.1704 and 0.2295 for SPC and T data, respectively. The importance of TPI, MBI, and valley depth were low in all FR models; convergence index, valley depth, and distance from a fault fell into the third lowest positions based on the FR-BCT models.


**Table 2.** Predictor importance values of each factor for the BCT and FR-BCT models.

#### *5.3. Model Performance Evaluation*

In this study, the groundwater potential model was evaluated based on statistical indices; AUC was used to quantitatively assess the mapping accuracy. As aforementioned, testing was performed based on the 30% of the groundwater well data collected by field investigation; and since groundwater has less seasonal change than surface water, this study did not consider seasonal change for groundwater. Figure 7 presents the model accuracy rate for the SPC (BCT model: 80.48%; FR-BCT model: 87.75%) and T (BCT model: 72.27%; FR-BCT model: 81.49%) well data. In general, all groundwater potential mapping results and modeling of groundwater potential showed good performance; however, the ensemble models showed improved accuracy by approximately 6%. Figure 7 also shows the performance of the groundwater potential models using the ROC curve method. All groundwater potential models performed well in terms of groundwater potential evaluation results (AUC > 0.7). The testing results of the BCT ensemble model show that 20% of the groundwater potential area includes approximately 80% of the valid groundwater wells for SPC, whereas the testing results of the ensemble model for T show that 30% of the groundwater area includes over 80% of the valid groundwater wells. Compared to groundwater potential mapping with the single machine learning model, BCT, all groundwater potential models using the ensemble method with both FR and BCT showed better performance, with 7.27% and 9.22% higher accuracy, respectively, than the BCT model alone. The difference in AUC results showed that the ensemble model provided better results than the individual modeling process. *5.3. Model Performance Evaluation* 

 Valley depth 0.0887 0.0696 0.0929 0.0513 Hydrology Flow path 0.1123 0.0917 0.1819 0.1407 Slope length and steepness (LS) factor 0.1835 0.1466 0.2747 0.1935 Forest Forest type 0.1821 0.1704 0.2084 0.1694 Soil Soil 1.0000 1.0000 1.0000 1.0000 Landcover Landcover 0.1946 0.1572 0.3285 0.2295 Geology geology 0.1002 0.0996 0.1619 0.1386 Distance from fault 0.1271 0.0497 0.3105 0.1061

In this study, the groundwater potential model was evaluated based on statistical indices; AUC was used to quantitatively assess the mapping accuracy. As aforementioned, testing was performed based on the 30% of the groundwater well data collected by field investigation; and since groundwater has less seasonal change than surface water, this study did not consider seasonal change for groundwater. Figure 7 presents the model accuracy rate for the SPC (BCT model: 80.48%; FR-BCT model: 87.75%) and T (BCT model: 72.27%; FR-BCT model: 81.49%) well data. In general, all groundwater potential mapping results and modeling of groundwater potential showed good performance; however, the ensemble models showed improved accuracy by approximately 6%. Figure 7 also shows the performance of the groundwater potential models using the ROC curve method. All groundwater potential models performed well in terms of groundwater potential evaluation results (AUC > 0.7). The testing results of the BCT ensemble model show that 20% of the groundwater potential area includes approximately 80% of the valid groundwater wells for SPC, whereas the testing results of the ensemble model for T show that 30% of the groundwater area includes over 80% of the valid groundwater wells. Compared to groundwater potential mapping with the single machine learning model, BCT, all groundwater potential models using the ensemble method with both FR and BCT showed better performance, with 7.27% and 9.22% higher accuracy, respectively, than the BCT model alone. The difference in AUC results showed that the ensemble

**Figure 7.** Testing results of the BCT and FR-BCT models for (**a**) SPC and (**b**) T groundwater data. **Figure 7.** Testing results of the BCT and FR-BCT models for (**a**) SPC and (**b**) T groundwater data.

#### **6. Discussion**

**6. Discussions**  In this paper, the relationship between conditioning factors and groundwater was first analyzed through the stochastic method of FR. By applying the ensemble technique to the BCT In this paper, the relationship between conditioning factors and groundwater was first analyzed through the stochastic method of FR. By applying the ensemble technique to the BCT model based on the stochastic weighting, it showed effectiveness in the study of groundwater with high uncertainty. In terms of data, this study was based on data created by governments and public institutions and released to the public; at the same time, it is bound by limitations in data collection. Since the importance of data used for training in data-based learning is very high, model accuracy will be improved if more well data is used in future studies.

Few case studies have applied ensemble models from machine learning algorithms in South Korea. The results of this study confirm that the performance of a groundwater potential model can be improved using an existing probability model and machine learning ensemble. Model performance was evaluated based on the ROC, and the prediction rate of the BCT model showed an improvement of 6.1% with FR-BCT for SPC and 6.0% for T compared to the single machine learning model, BCT, indicating that the ensemble method greatly improved model performance. This improvement occurred because the ensemble model could reduce bias using the BT model and improve its predictive ability by avoiding the overfitting problem of basic classification [70]. This finding is consistent with other studies that concluded that the predictive performance of models was improved with a machine learning ensemble model [71].

Remote sensing is a powerful data source that is widely used for monitoring environmental issues; however, since groundwater does not exist on the surface, groundwater can only be indirectly estimated by using remote sensing. Heretofore, many studies have attempted to reduce the uncertainty of groundwater spatially. As a result of applying the proposed FR-BCT model with existing probability models and the machine learning method of the BCT model, the accuracy was relatively improved or similar to previous studies [3,25,34,68]. In addition, by showing accuracy improvements in single and composite models, it has shown potential for reducing the uncertainty of groundwater potential mapping.

#### **7. Conclusions**

The modern global water shortage requires effective water management and planning. Indiscreet use of water resources and inadequate water management can disrupt the continuous and reliable supply of water. The first step in properly planning water resource usage is to accurately predict and respond to the current status of critical resources. Groundwater represents an excellent water source, especially in water-scarce regions. However, the uncertainty of groundwater availability is high; therefore, estimation of groundwater potential is essential. Mapping of groundwater potential is an essential challenge facing effective groundwater resource management and conservation planning.

Various methods of groundwater potential mapping have been proposed. Improvement of the groundwater potential model is one method for estimating the uncertainty of a groundwater model. Although new machine learning technologies are continually improving in predictive performance, not all methods can be effectively applied in areas where data are scarce, because it may not be possible to generalize from a small labeled dataset. Therefore, FR analysis and the BCT model were applied along with the proposed FR-BCT model, which is an ensemble model of these two machine learning models. For this purpose, 16 groundwater control factors based on remote-sensing data were applied to the models: nine topographic factors, two hydrological factors, forest type, soil material, land use, and two geological factors. The model was trained and tested using groundwater well data; 53 wells were separated into training (70%) and testing (30%) datasets. The proposed FR-BCT model was compared with existing probability models and the machine learning method of the BCT model.

These results are useful for supporting comprehensive management of groundwater exploration and groundwater recharge. The method used in this study can be applied to other areas reliant on groundwater use. Managers and policymakers can effectively analyze groundwater potential modeling results to maximize the benefits of management. However, further testing is required in other research areas to determine how reliably the proposed ensemble model reflects groundwater potential.

**Author Contributions:** Conceptualization, S.L. and M.-J.L.; Data curation, S.L. and S.L.; Formal analysis, Y.H.; Investigation, Y.H. and S.L.; Methodology, S.L. and Y.H.; Project administration, M.-J.L.; Resources, S.L.; Software, S.L. and S.L.; Supervision, M.-J.L.; Validation, S.L.; Visualization, S.L.; Writing—riginal draft, S.L.; Writing—review and editing, M.-J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was conducted at the Korea Environment Institute (KEI) with support from the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2018R1D1A1B07041203). This research was conducted by the Basic Research Project of the Korea Institute of Geoscience and Mineral Resources (KIGAM) funded by the Ministry of Science and ICT. This research was also conducted with support from 'Future Forecasting Government Management Support: A Study on the Modeling of Climate and Environment for the Transition to Data Science' funded by the National Research Council for Economics, Humanities and Social Sciences.

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

#### **Appendix A**


**Table A1.** Results of the frequency ratio model.

**Table A1.** *Cont.*


**Table A1.** *Cont.*



**Table A1.** *Cont.*

#### **References**


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

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Remote Sensing* Editorial Office E-mail: remotesensing@mdpi.com www.mdpi.com/journal/remotesensing

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-1603-5