**Novel Approaches in Landslide Monitoring and Data Analysis**

Editors

**Jan Blah ˚ut Michel Jaboyedoff Benni Thiebes**

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

*Editors* Jan Blah ˚ut Czech Academy of Sciences Czech Republic

Michel Jaboyedoff University of Lausanne Switzerland

Benni Thiebes German Committee for Disaster Reduction (DKKV) Germany

*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 *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special\_ issues/landslide\_monitoring\_data).

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-3787-0 (Hbk) ISBN 978-3-0365-3788-7 (PDF)**

Cover image courtesy of Jan Blah ˚ut

© 2022 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**



### **Preface to "Novel Approaches in Landslide Monitoring and Data Analysis"**

Recent progress in landslide science is shown not only by the constant increase in the number of published papers but also by the novel methods, approaches or best practices applied. This Special Issue, which is focused on the novel methods in landslide monitoring, modelling and data analysis, brings together scientists from around the world and shows the state-of-the-art research in this fast-evolving branch of landslide science. We hope it will help to manage the significant landslide hazards around the world.

#### **Jan Blah ˚ut, Michel Jaboyedoff, and Benni Thiebes**

*Editors*

### *Editorial* **"Novel Approaches in Landslide Monitoring and Data Analysis" Special Issue: Trends and Challenges**

**Jan Blah ˚ut 1,\*, Michel Jaboyedoff <sup>2</sup> and Benni Thiebes <sup>3</sup>**


**Keywords:** landslide; monitoring; modelling; susceptibility; InSAR

#### **1. Introduction**

The purpose of this Special Issue is to bring together recent studies related in particular to landslide monitoring and data analysis. In engineering geology, geotechnical engineering and geomorphology, landslide monitoring using standard techniques is quite common. However, the rapid development of both hardware and software solutions, including miniaturization or remote sensing techniques, brings new possibilities for increasing monitoring accuracy, real-time or near-real-time data analysis and early warning.

#### **2. Summary of the Special Issue Contents**

The Special Issue topics can be sub-divided into three groups according to the main topic covered by the articles. The majority of the articles (seven) are focused on landslide monitoring, monitoring data analysis and surveying, while a further two papers are focused on slope stability modelling using large-scale analog models and the remaining four papers deal with landslide susceptibility and detection.

#### *2.1. Landslide Monitoring, Monitoring Data Analysis and Surveying*

Thiery et al. [1] performed airborne electromagnetic measurements for rapid surveyance of the volcanic tropical environment of La Martinique, an island in the Caribbean. They combined their findings with a physical-based model to obtain improved and integrated information about the internal structure of landslides, founding a better understanding of landslides' initiation conditions.

Gili et al. [2] have monitored the Vallcebre landslide in the Pyrenees in NE Spain since 1987. A range of classical and novel methods have been used to that end (e.g., triangulation, photogrammetry, wire extensometers, GNSS-GPS, satellite DInSAR and terrestrial GBSAR). They conclude that while some methods give higher-precision results than others, all systems play valuable roles in landslide movement interpretation, and provide meaningful monitoring results at different study stages.

Fang et al. [3] present a monitoring system installed on the Pingding landslide in Taiwan. Their system consists of a GPS array combined with inclinometers, extensometers and rainfall data. The system is emergency response-centered and provides a basis for local early-warning indices. The paper tackles the important issue of multilateral cooperation among different subjects and disciplines involved in landslide disaster management.

Qiao et al. [4] focus on early-warning methods for largely abandoned rockfill slopes. These pose a significant threat in areas of large construction works. The authors used

**Citation:** Blah ˚ut, J.; Jaboyedoff, M.; Thiebes, B. "Novel Approaches in Landslide Monitoring and Data Analysis" Special Issue: Trends and Challenges. *Appl. Sci.* **2021**, *11*, 10453. https://doi.org/10.3390/ app112110453

Received: 3 November 2021 Accepted: 4 November 2021 Published: 7 November 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/).

ground-based InSAR to monitor slope deformations and verified the method on five landslides in the area of Huangdao, China.

Dabiri et al. [5] used object-based image analysis to map geomorphological features, and assessed the applicability of Sentinel-1 data to the fast creation of post-landslide digital elevation models. Their findings revealed that—without further post-processing—the automatically derived results need to be interpreted with care, as the automatic generation of a digital elevation model is influenced by several factors.

Blah ˚ut et al. [6] propose a methodology for analyzing time-series monitoring data from a large, slow-moving San Andrés landslide on El Hierro, Canaries, Spain. They used precise 3D dilatometric data and compared them with possible landslide-triggering factors (e.g., seismic, rainfall) to allow for fully automatic processing, thus decreasing the subjectivity of the analysis.

Troiani et al. [7] applied different surface analysis and monitoring methods to decipher the structural controls of rock slope stability in coastal areas. They worked on the Adriatic coast of the Conero promontory in Central Italy, and their results stress the need to analyze slope stability over a long timescale, to understand the current processes.

#### *2.2. Slope Stability Modelling Using Large-Scale Analog Models*

Feng et al. [8] determined a soil-water characteristic curve for landslide seepage under varying hydrodynamic conditions. They used large-scale experiments combined with finite element modelling. Consequently, they evaluated the uncertainties in the modelling using the Bayesian approach.

Tang et al. [9] assessed the influence of an intermediate coarse layer on slope stability during heavy rainfall. They found that the unsaturated hydraulic conductivity in the coarse layer was much lower than that of a fine layer, which led the capillary barrier to work at a lower water content. They also revealed that the coarser layers may have negative effects on slope stability.

#### *2.3. Landslide Susceptibility and Detection*

Fabbri and Patera [10] searched for uncertainties associated with the prediction patterns of landslide susceptibility maps. They conclude that the properties of prediction patterns are mostly unknown, but nevertheless, are critical for interpreting and justifying prediction results.

Lai [11] performed an automated data-mining procedure to differentiate the landslide sources and runout zones of landslides triggered by Typhoon Morakot in Taiwan. The author's models revealed that the detection of landslide sources provided accurate results, while the extraction of the runout areas achieved excellent accuracies.

Li et al. [12] explored the influence of multitemporal digital elevation models on the generation of susceptibility maps in the southern Sichuan Province in China. They conclude that the susceptibility assessment level of an area with historical landslides decreases in the short-term and that the usage of multitemporal digital elevation models has a serious impact on susceptibility results.

Li et al. [13] prepared a spatial, proximity-based, geographically weighted regression susceptibility model for the Qingchuan area in China. Their results suggest that the newly developed model shows higher predictive accuracy than five other commonly used models.

#### **3. Bibliometric Analysis of Current Trends**

A simple bibliometric analysis was performed in the ISI Web of Knowledge: "Web of Science Core Collection" (1900–Present) to capture the main trends in the Special Issuerelated topics. The statistics were gathered on 2 November 2021, so only the first ten months of the year 2021 were included in the search. We followed a similar approach to Jaboyedoff et al. [14] and selected keywords related to the topic of this Special Issue. A number of papers published in every year have been analyzed for the query "TOPIC" (Table 1).


**Table 1.** List of queries used in simple bibliometric analysis. Timespan: 1900–Present. Indices: SCI-EXPANDED, SSCI, A&HCI, CPCI-S, CPCI-SSH, BKCI-S, BKCI-SSH, ESCI, CCR-EXPANDED, IC.

\*: Search keywords.

Before 1990, only a small share of the total number of papers had been published. For the topic "landslides", only 585 papers (1.65%) were published before 1990 out of 35,440 (till 10/2021); for the "landslide + monitoring" topic, only nine papers (0.17%) were published before 1990 out of 5210 (till 10/2021); and for the topic "landslide + susceptibility, modelling or assessment", only 13 papers (0.11%) were published before 1990 out of 12,248 (till 10/2021). In the case of the "landslide + InSAR" topic, the first publication appeared in 1998, and up to the end of October this year, 1079 papers had been published. Figure 1 shows the yearly growth in the number of papers in the respective topics since 1990. The presented graphs illustrate the fast-growing exponential trend in all the analyzed topics. However, the fastest growth in recent years can be seen in the "landslide + InSAR" topic.

One of the basic variables for assessing the trends in a time-series is the growth rate [15]:

$$k\_i = \frac{y\_i - y\_{i-1}}{y\_{i-1}} \cdot 100 \tag{1}$$

where *ki* is the growth rate in %, *yi* is the value at the time *i* and *yi*−<sup>1</sup> is the value at time *i* − 1.

Positive values mean a growth in percentage terms compared to the previous year while negative values mean a decrease compared to the previous year in percentage terms. As can be seen from Figure 2, rapid growth in topics "landslides", "monitoring" and "susceptibility" in the 1990s has stabilized since 2006; studies on these topics are now increasing at a relatively stable rate of around 10–15% per year. The "InSAR" topic is a special case as the first paper appeared in 1998, which consequently resulted in very high growth rates. However, during the last five years, this trend seems to have stabilized, and the topic now seems to follow the trend of the other three topics with a yearly growth rate of between 15 and 20%.

**Figure 1.** Number of published papers per year on the selected topics since 1990.

**Figure 2.** Growth rate of the number of published papers per year on the selected topics since 1990. Lines show the five-year moving average.

#### **4. Future Challenges**

Future trends are always hard to predict. It can be expected, however, that the current trend of yearly growth in the number of papers published on the Special Issue topics, of between 110% and 120%, will continue during the next five to ten years. This could result in more than 10,000 papers published per year on the "landslides" topic by 2026, or 20,000 papers published by 2030. It is questionable whether this enormous number of papers could be published as to do so would place enormous pressure on all the persons involved in the publication process, especially the editors and reviewers. Even now, it is hard to ensure a rigorous peer-review process as scientists are often overloaded with review requests. This situation is unsustainable in the long-term and indicates that we might expect important changes in the publication process of landslide-related scientific papers.

One possible step toward solving this situation is to better classify the type of paper published. Papers should be distinguished between novel, innovative papers and applications of existing techniques and methodologies (case studies). While the innovative papers are potentially highly citable, the case studies bring new examples. Thus, both types are important, but do nothing to decrease the overload faced by editors and reviewers.

**Author Contributions:** J.B. prepared the manuscript, M.J. and B.T. revised and edited the text. J.B., M.J. and B.T. were Guest Editors of this Special Issue and, therefore, reviewed all the submissions. All authors have read and agreed to the published version of the manuscript.

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

**Data Availability Statement:** The data analyzed in this paper have been downloaded from "Web of Knowledge Core Collection" (https://www.webofscience.com/wos/woscc/advanced-search, on 2 November 2021).

**Acknowledgments:** The Guest Editors would like to acknowledge the external reviewers that assured the quality standards of the Applied Sciences journal were reached through their much-appreciated iterative review process. The authors, assistant editors and academic editors are also gratefully acknowledged for their interest in contributing to this Special Issue.

**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**


### *Article* **Airborne Electromagnetics to Improve Landslide Knowledge in Tropical Volcanic Environments**

**Yannick Thiery 1,\*, Pierre-Alexandre Reninger <sup>2</sup> and Aude Nachbaur <sup>3</sup>**


**Abstract:** Caribbean areas are particular volcanic territories in tropical environments. These territories juxtapose several landslide-prone areas with different predisposing factors (poorly consolidated volcanic materials, superimposition of healthy materials on highly weathered materials, high heterogeneity of thicknesses, etc.). In these environments, where rapid development of slopes and land use changes are noticeable, it is necessary to better characterize these unstable phenomena that cause damage to infrastructure and people. This characterization has to be carried out on the materials as well as on the initiation conditions of the phenomena and requires complementary investigations. This study, focusing on La Martinique, proposes a landslide analysis methodology that combines new information about landslide-prone materials acquired by an airborne electromagnetics survey with a physical-based model. Once the data are interpreted and compared with field observations and previous data, a geological model is produced and introduced into the physical model to test different instability scenarios. The results show that geophysical investigations (i) improve the knowledge of the internal structure of landslides and surficial formations, (ii) specify the spatial limits of the materials that are sensitive to landslides, and (iii) give a better understanding of landslide initiation conditions, particularly hydrogeological triggering conditions.

**Keywords:** airborne electromagnetics; landslide; physical-based modeling; tropical volcanic environment; La Martinique

#### **1. Introduction**

Landslides are ubiquitous phenomena in the Caribbean [1–6], particularly in La Martinique [7–9]. With more than 600 events [10,11], this territory is the most affected area in the French Caribbean islands [8]. Phenomena can be shallow, deep, rotational, translational, or complex. The many landslides in the Caribbean are mainly due to the following reasons:


Landslides regularly strike the coasts and the hinterlands, and because the island has much built-up land, resulting in anarchic development of the slopes [10,18], landslides can generate damage to the population and infrastructure, creating high rehabilitation costs [8]. Among the most remarkable recent events was (i) the Bellefontaine collapse (vol. = 15 × <sup>10</sup><sup>4</sup> <sup>m</sup>3) in 1991 [18], which required €7 million in work to rehabilitate the slope after the event, and (ii) the Morne Callebasse landslide (vol. = 2 × 105 <sup>m</sup>3) in 2011, which destroyed more than 20 buildings and the road 'RD 48- , bringing 75 expulsions and more than €17.1 million in works [15,18]. Therefore, anticipating landslides and improving their prevention in this French overseas territory has become a major challenge [8,15].

**Citation:** Thiery, Y.; Reninger, P.-A.; Nachbaur, A. Airborne Electromagnetics to Improve Landslide Knowledge in Tropical Volcanic Environments. *Appl. Sci.* **2021**, *11*, 3390. https://doi.org/ 10.3390/app11083390

Academic Editor: Jan Blahut

Received: 28 February 2021 Accepted: 6 April 2021 Published: 9 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/).

The first step in achieving these actions consists of assessing hazards [19–21]. Several methods that are more or less complex and range from qualitative to quantitative can be used. A large review of their uses can be found in the scientific literature [19–23], and some of them are particularly suited to the French regulatory context [24]. Among them, methodologies using physical-based models (PBMs) at the site scale (<1:5000) are the only methods designed for works [25] or that modified and revised existing regulatory hazard maps, taking into account triggering factors [24]. At this scale of work, it is possible to reach a good understanding of the mechanisms and probabilities of failure or to analyze the runout phenomena by taking into account physical processes. However, in Martinique, the parametrization of PBMs can be difficult because of a lack of information about the nature and depth of materials. Indeed, the island presents a very complex geology with different types of lavas deposited from the Oligocene to the present and exhibits extreme spatial variability and heterogeneity with a high degree of weathering [26–30].

Geophysical investigations can be an alternative to classical field investigations (i.e., geotechnical surveys; [31,32]) to obtain quick information about grounds in large areas. These investigations presuppose correlations between measured properties and physical or geotechnical characteristics [33,34]. In the case of landslides, it is possible to delineate the body of the moving mass of stable grounds due to the changing geophysical parameters [34,35]. However, with heterogeneous grounds, the indirect information that is provided can be biased or difficult to interpret. Therefore, calibrating them with direct observations is indispensable [31], thereby limiting investigations in time and space for very large areas [32].

Airborne electromagnetics (AEM) provides information on lithology and regolith over large surfaces and/or hard-to-access areas [36–40]. It is able to provide plentiful pseudo-3D information about geological structures reaching a few hundred meters [41,42] and has been successfully used in many environmental studies [30,43–47] particularly in volcanic settings. AEM is therefore useful to analyze the internal structure of large landslides [48–51].

In 2013, an AEM survey was conducted over Martinique Island [52]. The SkyTEM system was designed for mapping geological structures and for hydrogeological and environmental investigations [53]. AEM data, acquired along flight lines, provide information to a depth of 150 m and allow continuous and homogeneous imagery of resistivity variations [30]. If the subsurface resistivity has an indirect relationship with the soil characteristics [31], it can also provide relevant information on both formation thicknesses and their spatialization [32]. To be fully relevant, correlation with direct observations or independent geological datasets (outcrop, borehole geological log, etc.; [31,32]) remains mandatory. Recently, in the framework of future development [15], there has been demonstrated interest in AEM to delineate landslide-prone areas and to improve landslide hazard assessment. New information brought by AEM has allowed a landslide hazard map to be modified by combining a 3D geological model derived from a joint analysis of field observations and AEM results and a PBM. This study, for regulatory purposes, focused on only one type of island environment and on shallow (failure < 2 m) and moderately deepseated landslides (failure from 2 m to 10 m). Therefore, deeper landslides (failure > 10 m) and associated formations inducing recurrent and very costly damage were not assessed. Considering these results and despite an exhaustive inventory and landslide hazard maps available for the entire island, there is still a lack of knowledge about (i) the nature of certain landslide-prone materials and (ii) the failure mechanisms and triggering conditions for landslides with different failure depths. This work suggests a landslide analysis methodology that combines information about regolith and bedrock derived from AEM results with a spatialized PBM (SPBM) adapted for different types and depths of failures.

Two different areas that are well known for different recurrent slope instabilities, are typical of the island and have caused significant damage for more than 20 years, were chosen to challenge the methodology. The latter is divided into 3 main steps:


The test sites have benefited from previous field, geophysical and geotechnical studies [54–65]. These studies enable (i) a criticism of the AEM results and their interpretation and (ii) the development of realistic failure scenarios based on historical observations and measurements.

#### **2. Martinique and Study Sites**

#### *2.1. Generalities*

Martinique is part of the French West Indies or the Lesser Antilles (Figure 1a) and results from the westward subduction of the Atlantic plate under the Caribbean plate [26,27]. As the largest island of the archipelago (i.e., 1080 km2), it can be divided into two parts: that with mountainous relief in the north (with Pelée Mountain, a.s.l. 1397 m) and that with a gentler slope in the south (Figure 1a). Orographic effects control rainfall; for instance, the average annual precipitation in the northern part ranges from 5000 to 6500 mm.yr−<sup>1</sup> at the highest elevations and from 1200 to 1500 mm.yr−<sup>1</sup> in the southern part (Figure 1c). Because the climate is characteristic of a humid tropical climate, a humid season from July to November and a dry season from January to April can be delineated. The annual temperature varies between 18 ◦C and 32 ◦C at Fort de France.

**Figure 1.** Location of La Martinique (the names of main cities are given in italics); (**a**) slope map; (**b**) simplified geological map; (**c**) precipitation map (mean annual rainfall per year); (**d**) landslides: type and locations (from French national landslide database: BD-MVT, https://www.georisques.gouv.fr/ accessed on the 18 May 2020).

#### 2.1.1. Geology

In terms of the geological history, this island is unique because it is located between the old volcanic arc to the east and the recent volcanic arc to the west. Due to this specific position, it is possible to study the chronology of the volcanic front in detail [29]. Starting more than 25 Ma, the volcanic history is complex with alternating eruptive and marine sedimentation phases [26–29]. These successions allowed the edification of volcanic complexes, which were weathered and dismantled by erosion [30]. Figure 1b gives an overview of the old and recent volcanic deposits.

#### 2.1.2. Landslides

More than 600 landslides have been inventoried during the last 20 years (Figure 2a). They are reported in the French National database (i.e., BD-MVT; https://www.georisques. gouv.fr/ accessed on the 18 May 2020). Each phenomenon is recorded with a minimum set of information, such as the date, the location of the event (centroid in the local geodesic system), the type of formation involved (lithology and/or regolith), and the associated damage, if any. Three types of landslides can be depicted: (i) landslides (i.e., debris slides, rotational and translational slides), (ii) mudflows, and (iii) rockfalls. Figure 1d gives an overview of the location of each phenomenon. Among the different landslides, seven are the subject of special attention because they have generated some damage, and they may continue to generate damage despite some engineering works [8,10].

**Figure 2.** La Médaille landslide; (**a**) location of test site; (**b**) picture of the early warning system to block road traffic in case of landslide activity; (**c**) panorama of compartment A of the landslide; (**d**) structural scheme of the landslide (map is produced with the hillshade from Helimap DTM, 2013).

#### *2.2. Study Sites: Presentation and Previous Works*

The two selected sites, which are located in the municipalities of Fort-de-France and La Trinité (Figures 2a and 3), are characteristic of the observed landslides over 20 years and still regularly generate damage despite works and monitoring systems. The two sites have benefited from various geomorphological, geophysical and geotechnical studies that can help guide the various interpretations resulting from airborne electromagnetic surveys.

**Figure 3.** Morne-Figue area; (**a**) picture of the limit north of the Morne-Figue landslide (the bumps on the road are produced by the shear cracks of the landslide); (**b**) main scarp of the Morne-Figue landslide in 1989; (**c**) panorama of the Morne-Figue area; (**d**) main morphological landslide features (the map was produced with an orthophoto from IGN, 2015, and hillshade map from Litto3D DTM, IGN, 2010).

#### 2.2.1. La Médaille Landslide

The La Médaille landslide is located in Morne Balthazar to the north of the Fort-de-France Municipality (Figures 2a and 3a). With a volume of approximately 260,000 m3 and an area of approximately 8 ha (Table 1), the landslide occurred in 1916 [54,55], removing the old village of La Médaille and causing five casualties. The landslide is bounded to the north and south by two water courses marking two parallel faults and by a cliff of approximately 110 m representing a normal fault in the west. It is probably located on a paleolandslide (rock avalanche) represented by dacitic deposits in the small-scale lithological map [26,27]. Since 1916, several peaks of activity have been recorded (i.e., 1958, 1966, and 1993; [54–62]), each time with the rehabilitation or the relocation of road RN3 due to its crossing the moving mass. It is a slow-moving landslide but is likely to experience accelerations due to GWL variations. More recently, a debris slide in weathered andesitic formations occurred in the upper part of the cliff overhanging the landslide (Figure 2b,c).

Since the 1960s, several investigations (i.e., field surveys, geotechnical studies, and monitoring; Table 1) were engaged to improve the knowledge of the landslide [54–63], especially for the upper part (compartment A). Indeed, the topography of the lower part (compartment C) of the landslide is chaotic and complex, and the very dense tropical vegetation prevents the deployment of different investigation devices. The different field surveys allowed delineation of four main compartments functioning more or less independently. The lower compartment (C) seems to have an influence on the upper compartment (A), while the compartments on the edges are considered independent of the main landslide stricto sensu [58]. All drilling campaigns (five campaigns between 1967 and 1996) were located in compartment A of the landslide, which presents a high blockage potential (Figure 2). Thus, this compartment benefited from 14 boreholes reaching a depth of approximately 20 m. These campaigns located the faults bordering the landslide (Figure 2d).


**Table 1.** Investigations and monitoring of the two sites since the 1960s.

The first campaign in 1967 revealed a more or less clayey upper layer, up to 3 m thick, followed by approximately 19 m of dacitic screes upstream and 3 m of dacitic screes downstream of compartment A (Figure S1). Under this layer, a thin clay layer of a thickness of 0.5 to 0.9 m was observed. Finally, the lower layer was interpreted as andesite bedrock that was more or less fractured [54,55]. Hazmoune et al. [57] questioned this statement following new observations based on new boreholes. Indeed, the last layer would be a mixture of weathered andesites and breccias and could be part of the landslide, which would then be composed of two superimposed bodies. This information corroborates the conceptual scheme of landslide development established in [56,57,63], which mentioned the implementation of a rock avalanche in andesites. This layer would have been weathered and would then have been covered by a thick layer of dacite debris (breccias), which would have also been weathered.

The different boreholes defined a groundwater level between 1967 and 1968 in the upper part of the landslide with a maximum piezometric level of approximately 10 m under the topography. The fault system and a nearby spring probably feed this water level [57]. Finally, a series of rainfall and piezometric records attempted to prove the relationship between precipitation, groundwater level (GWL) and landslide activity. Unfortunately, despite the punctual implementation of a measurement network, shortcomings in the precipitation series or in displacement monitoring did not allow this relationship to be clearly proven. Nevertheless, when the GWL is high, the activity of landslides increases [54,55,57].

#### 2.2.2. Morne-Figue Area

The Morne-Figue area is located in the commune of La Trinité on the east coast of the island (Figure 3b). The study site covers an area of 0.36 km<sup>2</sup> and is limited to the north and west by the Morne-Figue (129 m a.s.l.) and the Morne-Congo (232 m a.s.l.) areas, respectively, which are carved in compact andesite, to the south by the Gué stream and to the east by Crosmy Bay. The site is hilly with elevations varying between 10 m and 100 m a.s.l. and both steep and gentle slopes varying between 5◦ and 45◦. The site, which is very anthropized (approximately thirty houses), is bordered to the west by national road RN1 and is crossed by a small road connecting the neighborhood to the center of the municipality (Figure 2c,d). Since 1977, the east-facing slope has been subject to landslides (rotational and shallow translational landslides). The main phenomena is located below RN1. This landslide, which is approximately 2.99 ha and was triggered in 1988 [64], is associated with rotational failure in the upper part and a translational component in the lower part ([64]; Figure 3d). This phenomenon has two very active periods: from 1987 to 1988 and in 2004. Since 2004, new geomorphological features such as cracks (i.e., traction and compression), small scarps and new tension cracks on the road were observed intermittently (Figure 3).

One drilling campaign [64] and geophysical investigations allowed the involved materials and their thicknesses to be defined (Figure S2). Four types of formations from the

topographic surface were defined: (i) an upper layer composed of clays or backfill with a thickness between 1 m and 3 m; (ii) a second layer with a thickness from 8 m to 10 m and represented by weathered very clayey materials; (iii) a third layer with a thickness between 10 m and 12 m composed of weathered basalts (Pré-vert basalts); and (iv) a weak horizon of basalt whose thickness is not known. A piezometer was implemented for the period of 1987–1988. During the humid period (from May to October), the GWL lies between 1.4 m and 0.5 m below the topography from upstream to downstream. In the dry season, the GWL is approximately 8 to 9 m below the topography. The piezometric levels and displacement measurements were implemented for a short period of one year, showing a correlation between rainfall and landslide activity [64]. Consequently, the landslide was more active during the humid period, with a maximum displacement of approximately 1.4 cm for the period of 1987–1988. In 2016, the geological knowledge of the surroundings was refined by field observations along cross-sections. Beyond the acquisition of new geological information, two units within the landslide were delineated with a very active unit upstream and a latent unit downstream. New shallow parallel landslides were observed within the site.

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

The suggested methodology is based on a transdisciplinary approach in three steps. This approach reduced the uncertainties raised by previous studies, particularly concerning the internal structure of the La Médaille landslide and the different material types and thicknesses of the Morne-Figue area. The interconnected steps are (i) the identification and improvement of the knowledge of involved materials; (ii) the production of a geological model for each site, including the different regolith thicknesses and/or internal structures of landslides; and (iii) landslide susceptibility analysis by a spatialized physically based model based on limit equilibrium equations with research of best fitting geotechnical parameters. Step (iii) integrates empirical triggering factors according to field observations.

#### *3.1. AEM Data*

From 29 January to 16 March 2013, SkyTEM ApS conducted a survey over Martinique Island. This survey, supervised by the BRGM (French Geological Survey) and totaling 4233 line-kilometers, was flown along the N–S direction with a 400-m spacing; locally, this spacing was refined to 200 m over areas of interest for the hydrogeology or risk assessment [52]. Along each line, EM measurements were spaced approximately 30 m apart, with an average ground clearance of approximately 64 m due to the sharp topography of the island. Figure 4 shows the location of the flight lines flowing over the two studied sites.

SkyTEM is an AEM system developed by the Hydro Geophysics Group of Aarhus (Denmark) for hydroenvironmental investigations [53]. This system is composed of (1) a transmitter coil exciting the subsurface, (2) a receiver coil to measure the ground response, (3) a generator as a power source, and (4) several navigation instruments, such as GPS, tiltmeters, and laser altimeters, to locate the loop in space. The SkyTEM system operates in dual transmitter mode. The low moment, with a magnetic moment of approximately 2826 Am<sup>2</sup> and time gates from 11 μs to 1 ms for the present survey, provides early time data for shallow imaging, and the high moment, with a magnetic moment reaching 144,440 Am<sup>2</sup> and time gates from 109 μs to 8.9 ms, allows measuring later time data for deeper imaging. Locally, the depth of investigation of the method depends on the emitted magnetic moment, the bandwidth used, the subsurface electrical conductivity, and the signal/noise ratio.

**Figure 4.** Location of the Airborne electromagnetics (AEM) flight lines. (**a**) La Médaille landslide (hillshade map is produced with Helimap DTM, 2013); (**b**) Morne-Figue area (hillshade map is produced with Litto3D DTM, IGN, 2010).

The AEM method allows imaging of the conductivity contrasts of the subsurface. To obtain usable measurements, several processes are applied to remove couplings with man-made installations and ambient noise from the signal. The processing scheme used is described in [39,40] and is based on singular value decomposition [39,40]. Data are then inverted using the spatially constrained inversion algorithm (SCI) [44]: (i) each usable AEM measurement is translated into a 1D (EM sounding) model divided into n layers, each defined by its thickness and resistivity, and the resistivity variations are displayed according to depth; during the inversion, constraints are applied vertically and spatially between nearby soundings (independently of flight lines). (ii) The ground clearance of the AEM system is also inverted, and the depth of investigation (DOI) is assessed as a final step in the inversion [66]. The results were obtained by running a smooth inversion for 25 layers from 0 m to 300 m deep. Each layer has a fixed and equal thickness. Only the resistivity values vary within a layer [39]. This approach is effective for imaging complex geological structures with the lowest dependency on the starting model, but it only displays a smoothed view of the subsurface. At this step, each flight line can be displayed as a resistivity profile composed of all the associated EM soundings. More 2D information can be obtained by interpolating, on raster grids, the resistivity of layers falling into a depth or elevation range. This interpolation is generally repeated to obtain slices over the entire range of investigations. Slices can then be merged to obtain a 3D resistivity model [45], (i) drawing profiles in any direction for confronting geological data (maps, boreholes, and field observations) and (ii) deriving interfaces for each imaged horizon.

#### *3.2. Interpretations and Conceptualization*

#### 3.2.1. Confrontation with Independent Data

To interpret the different imaged horizons for defining interfaces of interest, a comparison with boreholes and field observations is carried out (Figures 2 and 3) by projecting the

geological data on selected resistivity profiles. The distances between the boreholes and AEM lines vary from 2 to 10 m and 3 to 12 m for the La Médaille landslide and Morne-Figue area, respectively.

#### 3.2.2. Geological Modeling

Based on interpreted resistivity profiles, boreholes, geological maps and field observations, a geological model is produced for each site [67,68]. AEM results are useful to constrain the model for different environments [37,38,41]. The principle is based on the first two steps of the approach described in [67] and the different positions of geological limits (contacts between two geological formations or lithologies). They are used to guide the interpolation by the kriging method as regular grids representing the bottom of each type of identified formation.

#### *3.3. Landslide Modeling*

#### 3.3.1. ALICE Presentation

Landside modeling is performed with ALICE® (Assessment of Landslide Induced by Climatic Events) developed by the French Geological Survey (BRGM, [69]). This tool supports landslide susceptibility mapping for areas ranging from local sites (catchments) to large areas (several municipalities; [15,69–71]). Developed in a GIS environment (MAPINFO®), it is a SPBM described in [69] and summarized below (Figure 5).

**Figure 5.** ALICE® concept (adapted from [69]).

The geometry of the studied area is introduced in raster format with (i) the topography and (ii) the geometry of different material layers. Geomechanical characteristics, cohesion (c), friction angle (ϕ) and specific bulk unit weight (γ), are given for each lithology and surficial formation. These parameters can be implemented by a constant value or by probabilistic distributions to take into account environmental variability and uncertainties [69,70]. The tool supports different landslide geometries and failures (i.e., rotational, translational, and complex with different lengths and depths). Two triggering factors can be used: (i) the groundwater level (GWL) or (ii) seismic acceleration. The GWL represents the saturation ratio (m = h/z), where h is the height of the water table and z is the depth of the material(s) taken into account. The GWL can be implemented empirically in one or several formations by increasing the saturation level from 0 (dry conditions) to 1 (saturated conditions) or with the help of a hydrogeological model taking into account the effective rainfalls. Seismic acceleration is represented by the peak ground acceleration (PGA).

The slope stability computation is based on a limit equilibrium method (LEM) and slice theory described in [72] and in Figure S1. The iteration process is based on the concept of reducing the number of iterations about the interslice function and therefore the computation time [73]. The hypothetical failure surface is divided into *n* vertical slices, and each slice *I* is subject to the normal shear interslice forces and to the shear resistance (Figure S3), where:

$$R\_i = [w\_i \cos \alpha\_i - u\_i b\_i \ \sec \alpha\_i] \tan \phi\_i' + c\_i' b\_i \ \sec \alpha\_i$$

and the moving forces:

$$T\_i = w\_{i^{
sim a\_i}}$$

where *Wi*: weight; *αi*: base inclination; *ui*: average water pressure; *bi*: width of the slice; *ϕi*: effective friction angle; *c*'i: cohesion along the base; and *Ri*: sum of the shear resistances, except the normal shear interslice forces. Ti is the component tending to cause instability. The potential failure is expressed by the factor of safety (FoS). If the FoS is below 1, the slope (i.e., the computation cell) is considered instable. The slope stability assessment is performed on regularly spaced 2D profiles automatically produced over the whole area and based on maximum gradient lines from the DTM. Two types of computation are possible, (i) the computation of the FoS or (ii) the computation of the probability, to obtain an FoS below 1 for each cell. In the second case, the tool performs a random selection of each geotechnical value following probabilistic distributions and Monte Carlo simulations [69–71].

#### 3.3.2. Landslide Modeling Protocol

Landslide modeling, split into three steps, consists of defining the best set of representative parameters for each site.


is established with a maximum fixed at −0.5 m from the topography surface for the La Médaille landslide, −0.5 m from the topography surface for the moderately deep rotational landslide and a GWL at the level of topography for the translational shallow landslide for the Morne-Figue area. The different saturation ratios increased iteratively from dry conditions (GWL = 0) to full saturation conditions measured for each site (GWL = 1).

For each step, the failure geometries were fixed following the field observations (Table 2). The validation of simulations is performed by comparison with (i) an expert map derived from field observations and DTM derivatives for the La Médaille landslide [62,76] and (ii) the landslide inventory performed in 2016 for the Morne-Figue area [62]. Different classical statistical tests are computed (relative error and ROC-AUC; [69–71]) to validate the results. These tests were completed by expert (qualitative) verification.

**Table 2.** Landslide characteristics used for modeling.


#### **4. Results**

#### *4.1. Identification of Involved Materials and Geological Models*

For each site, resistivity profiles and grids (at different depths) were combined with boreholes, field observations and geological maps (Figures 6 and 7). Interfaces of interest were then extracted to constrain the geological model.

**Figure 6.** Resistivity grids for the two sites, with landslide areas and main lithological formations from geological maps. (**a**) Resistivity map between a depth of 4 m and 7 m for the La Médaille landslide; (**b**) resistivity map between a depth of 38 m and 47 m for the La Médaille landslide; (**c**) resistivity map between a depth of 2 m and 4 m for the Morne-Figue area; (**d**) resistivity map between a depth of 19 m and 24 m for the Morne-Figue area.

**Figure 7.** Interpreted resistivity profiles. (**a**) La Médaille landslide (hillshade map is produced with Helimap DTM, 2013); (**b**) Morne-Figue area (hillshade map is produced with Litto3D DTM, IGN, 2010).

#### 4.1.1. La Médaille

Figure 6a,b shows two resistivity grids at two different depths, between 4 and 7 m and between 48 and 47 m deep. For each, the geological formations and the landslide were reported. The dacitic formations (breccias and bedrock) show high resistivity values (>65 Ω.m, Figure 6a). It is possible to observe them in the western part of the landslide. At the level of the landslide, the resistivity values are between 20 Ω.m and 35 Ω.m, corresponding to the dacitic screes observed in the boreholes and described in [58,59] and [61]. Deeper (Figure 6b), high resistivity values (>50 Ω.m) are imaged to the west and northwest of the landslide and correspond to the dacitic formations, whereas the landslide area is characterized by a low resistivity, below 15 Ω.m. A resistivity profile is shown in Figure 7a

and is completed in Figure S1. Nearby boreholes are also projected on this profile to identify the different materials (Figure 7a). At the surface, a thin conductive layer (C1) is visible with a resistivity of approximately 15 Ω.m. It would represent the superior clayey layer described in the different boreholes. Under it, a more resistant layer (R2 ~50 Ω.m) can be observed. It is assimilated by the dacitic formations of the old debris avalanche described in [54]. Below R2, a thick layer (R3) with a lower resistivity (from 15 Ω.m to 30 Ω.m) is imaged. This layer is the one observed in Figure 6a and was interpreted as breccias or fossilized screes in [54,55]. The deep conductive thick layer (C2) would correspond to weathered andesite. This conductive layer, present along the landslide (Figure 6b, Figure 7a and Figure S1) and not revealed with the older boreholes, would be 30 m thick and would overcome the bedrock composed of andesite (R4). Table 3 details each formation.


#### 4.1.2. Morne-Figue

Figure 6c,d show two resistivity grids at two different depths, between 2 and 4 m and between 19 and 24 m. The resistivity range is globally low (<50 Ω.m). At the surface, the highest resistivity corresponds to andesites and volcanic-sedimentary formations. Weathered basalts have a lower resistivity (< 8 Ω.m). Landslides occur in the majority of these materials described in [62,63]. In between, it is possible to define the hyaloclastites. Resistivity profiles are displayed in Figure 7b and Figure S2. Nearby boreholes are also shown. Six horizons are delineated: a conductive layer (C1; <5 Ω.m) corresponding to a very weathered clayey layer and backfill materials that is more or less saturated following the groundwater level, and a more resistive layer (R1; from 5 Ω.m to 8 Ω.m) corresponding to less weathered clay materials. Locally, this layer disappears in lieu of a more conductive layer C2 (from 2 Ω m to 5 Ω m).

According to borehole data (Figure 7b), this layer (C2) would correspond to highly weathered lavas, where water was observed [63]. C2 disappears progressively downstream of the slope and is less weathered at depth (R4) in favor of a few resistive layers (R1B and R1). Under R1, C1 and R1B, R2 shows higher resistivity values (between 8 Ω.m and 11 Ω.m). This layer, located 25 m below the topographic surface, would correspond to hyaloclastites. However, it remains difficult to know the exact nature of this layer due to the lack of deep boreholes. Finally, a more resistive layer (R3), with values higher than 20 Ω.m, corresponds to andesite formations overlying the hyaloclastite formations in some places. It should be noted that two sectors, where resistivity values drop sharply, probably correspond to faults delimiting the topographic depression at the center of the study site. Table 4 gives the different characteristics of each formation.


**Table 4.** Characteristics of the different formations observed in the la Morne-Figue area.

#### 4.1.3. Geotechnical Models

From the interpreted resistivity models, different interfaces were derived for each site. Figures 8 and 9 illustrate the different geological models for the two sites.

**Figure 8.** Geological model for the La Médaille landslide. (**a**) A 2D cross-section representing the different materials selected for physical-based modeling; (**b**) interfaces introduced for ALICE® (hillshade map is produced with Helimap DTM, 2013).

**Figure 9.** Geological model of the Morne-Figue area. (**a**) Location of homogeneous areas identified by resistivity data analysis, field observations and boreholes; (**b**) 2D cross-section representing the different materials selected for physicalbased modeling; (**c**) interfaces introduced for ALICE® (hillshade map is produced with Litto3D DTM, IGN, 2010).

For the La Médaille landslide, four layers were defined (Figure 8a): (i) a first layer corresponding to dacitic screes (R2) upstream of the landslide on steep slopes and on compartment A of the landslide. The thickness varies from a few decimeters to approximately 5 m upstream of compartment A. On compartment A of the landslide, the thickness varies from 5 m to 20 m. (ii) A second layer corresponds to andesite screes (R3) with a thickness between 40 m under compartment A and a few meters downstream of compartment C. Therefore, R2 rests on R3 along compartment A and forms compartment C alone. (iii) A third layer corresponds to dacite (R1) located under R2 up to the concave slope failure corresponding to the upstream limit of the landslide. It is limited downstream by the northsouth fault. Its thickness varies from 20 to 90 m. (iv) A fourth layer corresponds to andesite (C2) forming the bedrock. The fault observed in the resistivity profile (Figure 7a) and documented in the geological map [26,27] is not integrated in the geotechnical model. Figure 8a shows the failure hypothesis described in [54] and questioned in [57]. The new information on formations and their thicknesses introduced in the geotechnical model and numerical simulations under ALICE® should help to better understand how destabilization occurs.

For the Morne-Figue landslide, the area was divided into five homogeneous areas, defined by the main formations identified previously, the structure and the geomorphology (Figure 9a). Thus, S1 corresponds to the hyaloclastite formation (R2) resting on andesite formations (C1); this area corresponds to the main relief. S2 is composed of hyaloclastite (R2) marking steeper slopes. S2 is separated from S3 by a fault identified in the resistivity model and identified in [26,27]. S3 marks the beginning of the topographic depression of the site. It is characterized by a succession of conductive and more resistive layers (C1, R1, C2, and R1B), including the C2 layer over almost its entire surface. S4, in the southern part

of the site, differs from S3 by an absence of layer C2 on the entire downstream and eastern part. Therefore, downstream of S4, only C1 and R1 are identified. Area S4 is probably limited to the north by an assumed fault separating it from S3. Finally, S5 represents the downstream of the site with gentle slopes covered by alluvial deposits. Faults in this area are not integrated into the geotechnical model. Numerical simulations with ALICE® constrained by the new information brought by the AEM (nature and extension of the formations and their thicknesses) allowed spatializing landslide hazards at the study site.

#### *4.2. Slope Instability Analyses*

For stability analyses, the best geotechnical parameters to be introduced in the models should be chosen. Moreover, this step tests the contribution of the information derived from the AEM and uses it to build the conceptual models on both the different materials involved and the different destabilization conditions.

#### 4.2.1. Identification of the Best Geotechnical Parameters

The sensitivity analysis along the cross-sections (Figures 8a and 9b) must reduce the range of values of geotechnical parameters and define the consistency as the best fitting intrinsic for the involved materials. A set of 50 model iterations for each crosssection were carried out (i) with independent variation of each parameter within the range of values and (ii) in dry and fully saturated conditions to obtain the best range of combinations. Tables 5 and 6 give the selected formation and the retained boundaries to perform calculations. For the two sites, cohesion appears to be the most influential parameter in the safety factor modeling followed by the internal friction angle, while the weight bulk density appears to be the least influential, which is often noticed in this type of analysis.

**Table 5.** Geotechnical values selected after sensitivity analysis for the La Médaille landslide. IV = initial values; SV = selected values; γ = bulk unit weight; c = cohesion; ϕ = angle of friction; italic values are introduced as triangular probability distributions in ALICE®.


**Table 6.** Geotechnical values selected after sensitivity analysis for the Morne-Figue area. IV = initial values; SV = selected values; γ = bulk unit weight; c = cohesion; ϕ = angle of friction; italic values are introduced as triangular probability distributions in ALICE®.


Therefore, for the La Médaille landslide, the results for rotational failure with a depth of 20 m show that the cohesion must not be below 12 KPa and 6 KPa for layers R2 and R3, respectively, especially if the angle of friction is low for each layer. If the value is below one of these values, the computed FoS remains under 1, corresponding to a recurrent instability that does not correspond to reality.

For the Morne-Figue area, there are two cases: shallow and moderately deep-seated landslides. For shallow landslides occurring in the R1 materials, if the angle of friction is higher than 15◦, then the slope remains stable. For cohesion, when it is greater than 10 KPa, the slope remains stable. On the other hand, when cohesion is close to 5 KPa, the slope switches from stable to unstable according to the saturation scenario. For moderately deep-seated landslides, the situation is more complex, regardless of the conditions, and the slope is close to instability even with high cohesion values and/or high angle of friction values. When the GWL is low and the saturation conditions are null, then the FoS is high, but not higher than 1.5. Therefore, when the angle of friction is less than 10◦, the stability remains low under any condition. If the angle of friction is higher than 20◦, the slope is stable. For cohesion under 5 KPa, the slope is recurrently unstable (FoS < 1) regardless of the angle of friction. If the angle is greater than 10◦, then according to the angle of friction, the stability increases.

Tables 5 and 6 show the initial and retained values to be integrated in the spatial modeling for the two case studies. Values in italics are used in ALICE® simulations. Values in italics and bold are used with a probability distribution necessary to compute failure probabilities with ALICE®. The probability distributions taken into account are triangular (Figure 10). This shape of the probability distribution is classically selected for ALICE®, as explained in [69,70], and offers the best compromise for testing new hypotheses of slope destabilizations.

**Figure 10.** Representation of probability distributions of different materials. (**a**) La Médaille landslide; (**b**) Morne-Figue area.

#### 4.2.2. Identification of the Optimum Cell Size

For this step, computations are performed with 10 m and 5 m cell sizes. The geotechnical dataset is the best previously defined for each site. For the La Médaille site, two computations are carried out with rotational failures with a depth of 20 m. For the Morne-Figue site, four computations are performed (two for shallow landslides with a maximum

depth of 3 m and two for deep rotational landslides with a depth of 10 m). The goal was to define the best cell size given the best results, and it was decided to perform the computations taking into account the failure conditions observed in the field (i.e., with a high GWL corresponding to the saturation of the materials, GWL = 1). The comparison and validation of the results is carried out by two statistical tests and expert verification. The tests are classical: (i) the relative error analysis (performed with the observed failure identified) and (ii) the analysis of the area under an ROC curve (ROC-AUC). For the La Médaille site, failures correspond to the main scarps of the landslide (i.e., the scarp between compartments A and C and the scarps identified downstream). For the Morne-Figue, failures correspond to scarps of phenomena.

Figures 11 and 12 depict the results for computations carried out with the two cell sizes for each site. For the La Médaille test site, the different scarps are well identified by the models with high probability values of failure in the different scarps. The relative errors are very low, with 0.25 and 0.21 for the 10 m cell size and 5 m cell size, respectively. The two computed ROC-AUCs show that the different models have a high degree of fit with values of 0.89 and 0.91 for a 10 m cell size anda5m cell size, respectively. For the Morne-Figue site, for shallow landslides, the different models computed with cell sizes of 10 m and 5 m have a low degree of fit. Indeed, the relative error and ROC-AUC have values of 0.89 and 0.95 and 0.59 and 0.58 for 10 m and 5 m cell sizes, respectively. These results indicate a low degree of model fit and low model representativeness for translational shallow landslides. In contrast, for moderately deep rotational landslides, the results are very good, with low relative errors and high ROC-AUCs of 0.13 and 0.15 and 0.89 and 0.87 for a 10 m or 5 m cell size, respectively. The Morne-Figue failure area is well represented with high probability values from 0.1 and 1, indicating, according to the equations, high failure probabilities when the materials are saturated.

**Figure 11.** Comparison of results with cell sizes of 5 m and 10 m for the La Médaille landslide. (**a**) Map computed with a cell size of 5 m; (**b**) map computed with a cell size of 10 m. (**c**) Expert map of potential failures; (**d**) statistical tests for (**a**,**b**). For each computation for (**a**,**b**), the GWL = 1.

Hillshade maps were produced with the Helimap DTM (2013).

For both cases, the results few differ between computations with a cell size of 5 m or 10 m, either statistically or visually on the computed maps. Thus, for the two sites, a cell size of 10 m appears to be the best compromise for the next phase. For the Morne-Figue site, given the low representativeness of the models, it is necessary to have other information to better adjust the computations. Therefore, the influence of the GWL is not tested for shallow translational landslides.

**Figure 12.** Comparison of results with cell sizes of 5 m and 10 m for the Morne-Figue area. (**a**) Map computed for shallow translational landslides with a cell size of 10 m. (**b**) Map computed for shallow translational landslides with a cell size of 5 m. (**c**) Map computed for moderately deep rotational landslides with a cell size of 10 m. (**d**) Map computed for moderately deep rotational landslides with a cell size of 5 m. (**e**) Statistical tests for the four computed maps. For each map the GWL = 1. Hillshade maps were produced with the Litto3D DTM (IGN, 2010).

#### 4.2.3. Influence of the GWL

The initial hypothesis following antecedent monitoring and observations is that the two case studies are controlled by the rise of the GWL corresponding to the material saturation. Therefore, several scenarios were tested by gradually increasing the GWL. The scenarios take into account the best calculation cell size (i.e., 10 m), and the best set of geotechnical data. For the La Médaille landslide, the goal is to find the most favorable materials prone to failures, confirming the landslide-triggered hypotheses. Thus, several scenarios were tested with computations of failures (i) for R2 in dry and completely saturated situations (i.e., GWL = 0 and GWL = 1), (ii) for R3 in dry and completely saturated situations (i.e., GWL = 0 and GWL = 1), and finally, (iii) for both types of materials in dry and completely saturated situations (i.e., GWL = 0 and GWL = 1). Once the materials to be considered have been defined, the influence of the GWL is tested by gradually increasing its level from 0 to 1 (with a maximum at 0.5 m below the topography). For the Morne-Figue area, considering the stakes, the main goal is to improve the knowledge about the material of the compartment face and the rise of the GWL to best understand future landslide-prone areas. The maximum depth of the retained GWL takes into account the observations (i.e., the GWL reaches a maximum at 0.5 m below the topography for moderately deep rotational landslides).

For the La Médaille area, the results (Figure 13) show that in dry conditions (i.e., with a very low GWL), the landslide is computed as stable with a failure probability under the low class threshold. Only the dacitic screes (R2) on very steep slopes are computed with a high probability of failure. For the same R2 materials located on gentle slopes, by taking into account these materials alone, then the computed failure probabilities are low, even if the GWL = 1 (Figure 13a). Only a very thin area of the upper part and the northern boundary of compartment A are computed with high failure probabilities. The relative error and ROC-AUC are not satisfactory, with values of 0.73 and 0.55, respectively. By taking into account R3 formations, it is possible to better delineate some failure areas (Figure 13b). The tests are better than previous tests (Figure 13a) but not optimum, with a relative error of 0.32 and an ROC-AUC equal to 0.73 (Figure 13b). Nevertheless, the results are not exactly consistent with field observations, particularly the failure area observed upstream of compartment A. Finally, by taking into account both types of materials (R2 + R3), the failure areas observed in the field are well identified (Figure 13c). The relative error and ROC-AUC are acceptable, with values of 0.25 and 0.89, respectively.

**Figure 13.** Influence of the type of materials introduced to compute failure probabilities with ALICE® for the La Médaille landslide. (**a**) Computations with R2; (**b**) computations with R3; (**c**) computations with R2+R3. Hillshade maps were produced with Helimap DTM (2013).

Figure 14 illustrates the sensitivity of materials when the GWL varies inside of those materials. Therefore, water appears as the engine of recurrent landslide destabilization; the more the GWL increases, the more the number of computed cells with high failure

probabilities increase. Nevertheless, the different computations show that under a GWL of 0.9, the landslide remains stable with a high probability of failure located in the foot of the sliding mass. From a GWL level equivalent to 0.9, it is possible to notice that the areas with high probabilities of failure correspond to the main location of scarps (Figure 14) with a relative error and an ROC-AUC of 0.27 and 0.89, respectively.

**Figure 14.** Influence of the GWL on computed failure probabilities with ALICE® for the La Médaille landslide. The materials involved are R2+R3. Hillshade maps were produced with Helimap DTM (2013).

For the Morne-Figue area, Figure 15 illustrates numerical simulations for moderately deep rotational landslides with a GWL variation ranging from 0 to 1 (i.e., reaching 0.5 m from the topography). With a GWL equivalent to 0.5, the slopes remain relatively stable, and the failure probabilities are equivalent to 0, except near the Morne-Figue slope, where few cells are computed with failure probabilities between 0.001 and 0.0001 (i.e., moderate failure probability). When the GWL value exceeds 0.7, the initiation area of the Morne-Figue landslide is computed with high probability values from 0.7 to 0.01, indicating a high hazard

level following the JTC-1 classification of hazards [22]. A significant increase in probability values is noted for a GWL above 0.8. From this level, the model computes high probabilities of failure with a good recognition of the initiation area of the Morne-Figue landslide represented by values from 0.9 to 0.01. The northern hillside of the area, with slopes greater than 25◦, appears to be more susceptible, with probabilities greater than 0.3, indicating a very high probability of material mobilization by landslides. Finally, when the GWL is near the maximum, the Morne-Figue landslide failure area is clearly identified with failure probabilities between 0.9 and 1, indicating a very strong probability of landslide hazards. The northern hillside is also simulated with very high failure probabilities exceeding 0.4. Note that failure is possible along the Gué stream south of the study site.

**Figure 15.** Influence of the type of materials introduced to compute failure probabilities with ALICE® for the Morne-Figue area. (**a**) Map computed with a GWL = 0.5; (**b**) map computed with a GWL = 0.7; (**c**) map computed with a GWL = 0.8; (**d**) map computed with a GWL = 1; (**e**) statistical tests for the four computed maps. Hillshade maps were produced with the Litto3D DTM (IGN, 2010).

#### **5. Discussion**

AEM data have already been used to study large landslides in temperate alpine environments or in Japan [38,48–51]. Few studies have been carried out on multiple types of landslides simultaneously in complex tropical environments with superimposed lavas that are more or less weathered at depth. Despite this lithological and structural complexity, through this study, AEM data appear to be of primary interest to assess these phenomena generating damages and losses in the West Indies [1]. In the following paragraphs, the AEM method, the role of water in slope failures and the differences with regulatory landslide hazard maps are discussed.

Thus, first, the AEM data allowed a better definition of the internal structure and the materials prone to landslides for both the La Médaille and Morne-Figue sites. For the La Médaille site, some hypotheses [57–60] that were difficult to validate due to the lack of deep boreholes were confirmed, such as the existence of two superimposed main bodies (one associated with a very low resistivity <5 Ω.m) composing the landslide body. The extracted information also allows building a specific geotechnical model integrated into ALICE®. Figure 16 gives an overview of the conceptualization of the La Médaille landslide gathering ancient information with the information derived from the AEM data. Without this latter, the geological model should be different (Figure 16), with downstream and upstream compartments (A and C) defined in [54,55] and located on one type of thinner material. Using this model (Figure 16a), the results of the numerical failure simulations would probably show that only the upper part of the landslide is likely to be unstable. Conversely, taking into account the interpretations derived from AEM data, it is possible to show that the shear surface identified in [54,55] and questioned in [56,57] turns out to be true concerning failure probability computations.

For the Morne-Figue site, the data from the AEM allowed us to clarify the geological structure of the site: (i) the hyaloclastite footprint and (ii) the different successive low and high resistive horizons in the basalts with likely different weathered levels. The derived geological model was able to reproduce the unstable behavior of highly saturated materials, especially for moderately deep rotational slides. The results obtained for shallow landslides are less conclusive. However, a previous study conducted in [15] showed that AEM data and modeling under ALICE® could give satisfactory results for shallow landslides. Here, this is not the case, and we decided to stop the modeling of shallow phenomena because of the mediocrity of the results. To improve the results for this type of landslide, two points should be emphasized: (i) the shallow landslide is located on a preferential location of water circulation mentioned in [77,78]. The use of a static GWL in ALICE® did not allow faithful reproduction of the real influence of water at this location. (ii) The landslide was triggered in complex materials mixing road embankments and likely weathered lavas. This type of material was not investigated during the geotechnical laboratory tests, and at this location, the simplification of the model meant that these specific materials were not included in the study.

Second, beyond these considerations, the results of physical numerical modeling show that water (i.e., introduced by the GWL) remains the main factor of slope destabilization in this environment. For both sites, the GWL has to be high to obtain high failure probabilities. This statement confirms the field observations and the punctual measurements and monitoring during the 1960s and 1980s for the La Médaille landslide and Morne-Figue area, respectively. However, the question of the origin of this water remains. For example, for the La Médaille landslide, if the results seem conclusive, there are still uncertainties, such as the role of the probable aquifer of the C2 formations and the probable upwelling of water along the discontinuities in C2. The measured levels during the monitoring of the landslide activity show that the water supply to the landslide from the source is not sufficient to activate it.

**Figure 16.** La Médaille conceptualization. (**a**) Before this study; (**b**) after this study (\* data not used in the computations of failure probabilities with ALICE®).

Careful consideration of the AEM results shows low resistivity bodies (C2) with values of approximately 5 Ω.m under each unstable site for both sites. This low resistivity may come from (i) highly mineralized water; (ii) preferential flow in more or less weathered and fractured materials; and (iii) very weathered materials with likely a high proportion of smectite [30]. Taking into account the recent works in Martinique [12,30], it is possible to retain hypotheses (i) and (ii) with deep water upwelling in fractured materials feeding shallow water tables in landslide-prone materials, as illustrated in Figure 16b. However, to obtain such a resistivity value, this water must be highly mineralized, and currently, it is difficult to know the mineralization processes [30]. Therefore, this hypothesis needs to be confirmed by deep drilling and less punctual monitoring than those previously undertaken at these sites.

Third, the new geological interpretation and physical failure modeling allowed refinement of the landslide susceptibility maps, especially for Morne-Figue. Figure 17 shows the areas identified by the expert approach in 2004 at the 1:20,000 scale of work for all landslide types and the results of numerical failure simulations for an extreme GWL for

only moderately deep rotational landslides. The obtained maps are different, particularly (i) the area between the two landslides at Morne-Figue and the 1977 landslide and (ii) the area along the Gué stream or (iii) the areas around the Mont Morne-Figue carved in the andesite formations. One important point to highlight is that Figure 17a integrates all landslide types, and Figure 17b shows only one type of event. At this stage, if the numerical failure modeling resulting from the geotechnical model can help experts create some failure scenarios and better identify some potentially unstable sectors, it is necessary to further model shallow landslides, especially the triggering factors and the material thickness involved. There is a certain amount of uncertainty in the results of the shallow failure probabilities during the calibration phase. Additional investigations must be envisaged, probably based on the work of [15]. The new results of simulated failure probability maps allow better definition of the boundaries of high-hazard areas and improvement of the future final landslide hazard map. This will be achieved during a project to revise the landslide hazard regulatory maps between 2021 and 2022. The goal is not to reject one approach and its results compared to another. Indeed, they are complementary; the expert approach is based on the subjectivity of the expert who, through his or her experience, can intrinsically bring elements not taken into account in the modeling [19,24]. In addition, the geotechnical model may contain some uncertainties (i.e., interpretation, oversimplification of the GWL, and material not considered) and generate errors in the stability computations under ALICE®. Therefore, this approach based on interdisciplinarity is in accordance with the landslide hazard mapping strategy proposed in [24] for France with local revisions based on numerical spatialization tools able to provide new information on slope stability quickly for decision support. Therefore, this study complements the study conducted in 2017 on a nearby site for shallow landslides.

**Figure 17.** Landslide hazard maps for the Morne-Figue area. (**a**) Expert landslide hazard map obtained by expert approach and carried out at the 1:25,000 scale of work for all landslide types; (**b**) failure probability map obtained with ALICE® for moderately deep rotational landslides; map computed with GWL = 1; hillshade maps were produced with the Litto3D DTM (IGN, 2010).

#### **6. Conclusions**

AEM data are powerful in revealing in-depth resistivity contrasts in a complex tropical environment [30]. Once the results are compared with in situ investigations, such as boreholes and field observations, it is possible to better delineate the different superimposed lithologies and to understand the internal structure of various grounds. This study illustrates an additional possibility to obtain reliable data on the first 50 m by this noninvasive geophysical method. The results are fundamental for landslide studies in these environments with complex accesses limiting the use of classical field investigation techniques. Thus, it would be essential to verify low-resistance bodies located at approximately 5 Ω.m between a depth of 10 and 50 m. The two examples of the La Médaille landslide and the Morne-Figue area show that under landslides, some low resistivity formations could play the role of aquifers concentrating underground flows, punctually recharging the upper layers and influencing the activity of slope instabilities. The study of these bodies over the whole island, in correlation with recent hydrological studies [12,30] using AEM data, must be considered to better define this role on landslides [46]. It may be possible to refine the knowledge of unstable slopes from this information. From an operational point of view, the data acquired over the whole island could thus be exploited to (i) refine the knowledge of large landslides that regularly generate material damage and (ii) improve the spatialization of hazards and thus regulatory maps of landslide risks. Beyond the island of Martinique, the AEM method seems to be one of the most cost-effective methods to obtain geophysical information on shallow and deep formations, particularly in complex environments (topography, vegetation, and complex geological structure), such as the Caribbean. Thus, with the information acquired by this study and in [15], new perspectives to improve landslide hazard analysis and mapping for Martinique and other Caribbean islands can be considered.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/app11083390/s1, Figure S1: EM profiles for the La Médaille landslide. Figure S2: EM profiles for Morne-Figue area. Figure S3: Concept of slope stability computations with a slice method.

**Author Contributions:** Conceptualization, Y.T. and P.-A.R.; AEM data processing, P.-A.R.; methodology, Y.T.; software, Y.T. and P.-A.R.; validation, Y.T. and A.N.; formal analysis, Y.T.; investigation, Y.T., P.-A.R. and A.N.; writing—original draft preparation, Y.T., P.-A.R., A.N.; writing—review and editing, Y.T.; supervision, Y.T.; project administration, Y.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Risk and Prevention General Direction of the French Ministry of Ecological Transition and the French Geological Survey (BRGM), contracts PS14DRP041, PS15DRP017 and AP16DRP017. The heliborne geophysical survey, called "MartEM" program [53], was co-funded by BRGM, the FEDER funds for Martinique, the Regional Office for Environment Planning and Housing (DEAL), the Regional Council and the Water Office of Martinique (ODE).

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

**Data Availability Statement:** Some results of the study are available at http://infoterre.brgm.fr/ rechercher/switch.htm?scope=9 (accessed on 4 April 2021); http://ficheinfoterre.brgm.fr/document/ RP-66605-FR (accessed on 4 April 2021); http://ficheinfoterre.brgm.fr/document/RP-65407-FR (accessed on 4 April 2021).

**Acknowledgments:** We acknowledge the useful comments and suggestions from the three anonymous referees who helped us to enhance the manuscript. We also acknowledge our colleagues from the BRGM, Gilles Grandjean, Rosalie Vandromme and Severine Bernardie, for the constructive discussion about the approach.

**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**


### *Article* **Spatial Uncertainty of Target Patterns Generated by Different Prediction Models of Landslide Susceptibility**

**Andrea G. Fabbri 1,\* and Antonio Patera <sup>2</sup>**


**Abstract:** This contribution exposes the relative uncertainties associated with prediction patterns of landslide susceptibility. The patterns are based on relationships between direct and indirect spatial evidence of landslide occurrences. In a spatial database constructed for the modeling, direct evidence is the presence of landslide trigger areas, while indirect evidence is the presence of corresponding multivariate context in the form of digital maps. Five mathematical modeling functions are applied to capture and integrate evidence, indirect and direct, for separating landslide-presence areas from the areas of landslide assumed absence. Empirical likelihood ratios are used first to represent the spatial relationships. These are then combined by the models into prediction scores, ordered, equalarea ranked, displayed, and synthesized as prediction-rate curves. A critical task is assessing how uncertainty levels vary across the different prediction patterns, i.e., the modeling results visualized as fixed, colored groups of ranks. This is obtained by a strategy of iterative cross validation that uses only part of the direct evidence to model the pattern and the rest to validate it as a predictor. The conducted experiments in a mountainous area in northern Italy point at a research challenge that can now be confronted with relative rank-based statistics and iterative cross-validation processes. The uncertainty properties of prediction patterns are mostly unknown nevertheless they are critical for interpreting and justifying prediction results.

**Keywords:** landslide susceptibility; ranking; cross validation; prediction model; prediction pattern; target pattern; uncertainty pattern

#### **1. Introduction**

Predicting landslide susceptibility of a region or study area has become a critical necessity with the continuing expansion of urbanizations across hazardous landscapes, increasing soil deterioration and the extensive damages inflicted by landslides [1]. Regional approaches to quantitative predictions have developed following the availability of thematic maps in digital form, including that of interpreted aerial photography and remotely sensed images. Examples of themes related to mass movement processes are expressions of various aspects of the rocks, soil, land cover, land use, soil permeability, groundwater table, and topographic surfaces [2,3].

In any application, the themes are hopefully representing the gravity-induced physical context of the process of slope failure. For this reason, the landslides are considered spatially related with the thematic map units or values over the study area of concern. These spatial relationships are established and integrated by mathematical models that transform raw data from specific landslide distribution maps and contextual maps into prediction maps or maps of the likelihood of future landslide occurrence [4].

Obviously, the detail and accuracy of the digital maps used for modeling must be congruous and the statistical relationships significant to obtain convincing likelihood maps of landslide occurrences that, hopefully, are more informative than what is already known. The term "convincing" must reflect the interpretable quality of the integrated likelihood

**Citation:** Fabbri, A.G.; Patera, A. Spatial Uncertainty of Target Patterns Generated by Different Prediction Models of Landslide Susceptibility. *Appl. Sci.* **2021**, *11*, 3341. https:// doi.org/10.3390/app11083341

Academic Editor: Jan Blahut

Received: 29 January 2021 Accepted: 5 April 2021 Published: 8 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/).

maps as spatiotemporal predictions and the certainty—or conversely, the uncertainty associated with that quality.

The applications discussed in this contribution have the purpose of exposing aspects of spatial predictions that are commonly ignored—validation, robustness, and uncertainty of the resulting susceptibility levels expressed as relative ranks. They are the modeled values that make up the prediction patterns.

Analyzed are a database of landslide occurrences and their spatial context in a study area in northern Italy. The database has been the focus of experiments on comparison of models and analysis of modeling structure. Here, we will focus exclusively on the uncertainties affecting the ranks representing relative predicted levels of susceptibility to active landslides.

The following section deals with the predictive methodology, the terminology, and the analytical procedures applied. The next section describes the study area, its database for modeling, and previous research on landslide susceptibility. Experimental results follow on uncertainty identification and its relative measures. They consist of modeling prediction, target and uncertainty patterns, their interpretation, and comparison. Conclusions deem that uncertainty assessments are a desirable and necessary endeavor for making spatial prediction modeling a worthwhile practice.

#### **2. Predictive Methodology, Terminology, and Analytical Procedures**

A unified framework for modeling prediction patterns was proposed by Chung and Fabbri (1993) [5]. It was termed the "favorability function" and is used in the applications that are the focus of this contribution. The term favorability was selected for its generality in the meaning intended to be comprehensive and unified to cover a variety of mathematical interpretations. The following is a summary of the concepts in the modeling, the terminology used, and the processing strategy that enables us to assess spatial uncertainty.

A study area (SA) is first assumed to be selected by experts as sufficiently representative of the processes that are considered hazardous within it. Hopefully then, the statistics from spatial quantitative data available in the SA are suitable for information extraction by modeling. A second assumption is that we can have at least two subareas. First, an occurrence subarea is identified as affected by the hazardous landslide process, with the known presence of landslide trigger areas of a specific and congruous dynamic type. Then, another subarea is selected of non-occurrences in which landslides are known to be absent (or are unknown if present). The occurrence subarea is where the occurrences are located and characterized, and it is considered a training area (TA) for establishing the spatial relationships. Generally, the nonoccurrence subarea is the complement of the occurrence subarea within the SA. In some cases, however, the nonoccurrence subarea is only part of the SA or even outside the SA. This could also be the case with the occurrence area. Commonly, the occurrence subarea is a very small portion of the SA. A third assumption is that future occurrences within the SA will be similar in type and setting to the known ones.

To allow favorability function prediction modeling, a proposition is constructed as follows:

F*i*: that a point *i* in the SA will be affected by a future landslide |the presence of spatial evidence. (1)

The symbol | indicates "given." The proposition is to be supported as true by means of the known occurrence distribution and their setting. The setting is considered spatial evidence. Therefore, the proposition links the occurrence locations (or their neighborhoods) and their setting as multivariate context. To satisfy the proposition as a mathematical statement, both the digital images representing the occurrence locations and their setting have to be transformed into spatial relationships. These new images have been termed the direct supporting pattern (DSP). They are the indexed occurrence locations. Indirect supporting patterns of the proposition (ISPs) are termed the images of their setting. The term "pattern" is used to indicate the results of this functional transformation into spatial relationships.

The landslide occurrences are converted into clusters of adjacent picture elements or pixels with a numerical value sequentially identifying each occurrence. The image pixel values of the settings are transformed into normalized frequencies if derived from categorical maps such as lithology or land use. Continuous field digital images, such as elevation or slope angles, are transformed into density functions. The ratios of the normalized frequencies and the density functions for the TA are divided by those of the nonoccurrence subareas within the SA. These ratios are termed empirical likelihood ratios (ELRs). We could also apply other types of normalization; nevertheless, the end results would have the same purpose and functionality. ELRs, as described, range in value between zero and infinity and provide a measure of support of the proposition in (1).

Many different mathematical models can be used to convert and integrate the ratios into prediction patterns that classify the SA into levels of likelihood of future landslide occurrence. We will be using five different models in our applications but the favorability function modeling and the processing structure are independent of the mathematical models that will be indicated later on.

In essence, the likelihood ratios enable us to separate the presence of occurrence from their supposed absence within the SA. This separation uses another assumption. While the normalized frequencies or the density functions in the TA (the occurrence subarea) are related to observed/mapped landslide occurrences, the ones in the nonoccurrence subarea relate to both the absence of occurrences and the unknown occurrences. For the latter, we assume that their setting, due to being relatively rare, is diluted within the setting of the absence of occurrences. This assumption allows considering the ratios as a form of contrast between areas with presence and areas with a presumed absence of occurrences. Commonly, the TA with respect to the SA is three or four orders of magnitude smaller.

A favorability function at each point of a SA should have two properties for modeling, i.e.: (1) should be able to measure a relative level of likelihood that a point *i* in the SA contains part of a future landslide of the given dynamic type and (2) should be able to provide a measure of uncertainty associated with the function by using only the part of all possible landslide occurrences present in the training area, TA.

Now suppose that we use a given prediction model, still unspecified at this point, that generates integrated scores as predicted values for every pixel of the SA. The scores are in some a-dimensional units or values between a minimum and a maximum. How do we establish the very high, high, intermediate, low, and very low likelihood of future landslide occurrence? This is not a simple question to answer, and it leads to further assumptions depending on the kind of data available for the SA, database of DSP and ISPs, and additional information.

Suppose, for instance, that we have temporal characterization for the many occurrences of the given specific dynamic type in the TA. We could use the occurrences from an older time interval to model a prediction pattern and then overlay it with the locations of the occurrences from the younger time interval. In this way, we obtain some statistics about their proportions within higher likelihood scores. Should some of them fall within the lower likelihood score areas, we would consider them as poorly predicted. Vice versa, well-predicted occurrences should fall on the higher likelihood score areas. We have termed this exercise as cross validation. In practice, this is a natural way to establishing how "good" our prediction pattern is as a "predictor of future occurrences."

Now let us suppose that the information on time partitioning of the occurrences, the DSP, is not available, as it is in most cases. How can we proceed with some other forms of cross validation? We could empirically pretend, for instance, not to know the existence of some of the occurrences, e.g., 25% of them, and use the remaining 75% to obtain a prediction pattern and then cross-validate it with the 25% we pretended not to know and that was not used as DSP for modeling. In this case, the "future" landslide occurrences for cross validation are the "next" 25%. Furthermore, we could operationally devise iterative strategies, depending on the number of occurrences available as DSP, such as (1) sequential exclusion of an arbitrary number of occurrences to be used for cross

validation of the pattern produced with the remaining ones, (2) sequential selection of a number of occurrences for modeling and using the remainder for cross validation, or (3) random selection of a number of occurrences repeated a convenient number of times. All these strategies will provide different ways to predicting the next arbitrarily selected number or proportion of occurrences. Recall that we have not yet selected any particular mathematical model for predicting.

Cross validation is a strategy for assessing the quality of our prediction modeling and also for comparing different prediction patterns, produced either by varying the number of ISPs or the mathematical models. The results of cross validations are tables of prediction scores for numbers or proportions of occurrence pixels in the SA. How do we interpret these dimensionless integrated scores generated for each pixel of a SA and ranging between a minimum and a maximum?

Given the number of transformations from the original map unit names or continuous values used to compute the ISPs and their conversion into relative scores, these are considered here as impossible to interpret directly by recognizing systematic changes or breaks. Instead, the prediction scores are easily converted into equal-area ranks after sequencing them in decreasing order. It was found convenient to obtain 200 ranks each corresponding to 0.5% of the SA. By equal-area ranking the prediction scores, it becomes practical and simple to display the prediction patterns and generate cross-validation tables of proportions of validation occurrences for each rank and cumulative tables of these proportions to be represented as curves.

A consequence of iterative cross validation is that a different prediction pattern is generated each time so that the set of patterns can be used to further characterize the initial prediction pattern obtained using all the occurrences in the DSP. We can then pile up the patterns and, if their numbers are sufficiently large, apply some form of statistics to obtain average and variance for each stuck of pixels in the SA. We have defined "target pattern" as what we wish to have—a representation of all past and future landslide occurrences (as susceptibility scores and associated uncertainty scores). In our case, some form of averaging the ranks of the prediction patterns from the iterations. We have defined "uncertainty patterns" as the expression of deviations from these averages. From practice, we have found that a very robust means for generating a target pattern is selecting the median rank for each pixel from the iteration prediction patterns, and for the uncertainty pattern, the rank of the ranges of deviations from the median ranks.

Furthermore, three assumptions must be reasonably satisfied [6], namely, (a1) the known landslide occurrences, the DSPs, are a "random selection" of all existing ones, known and unknown (allowing to extend the favorability function from the TA to the entire SA); (a2) the ISPs are correlated with the target pattern (allowing to estimate the function using the known part of the target pattern in the TA); and (a3) the process of slope failure is not random and follows a certain rule (allowing to model the favorability function).

Target and uncertainty patterns have opened the way to further characterize prediction patterns. Some of these will be discussed in Section 4. Let us now consider some favorability modeling functions of empirical likelihood ratios commonly applied and used here in our analyses. They are fuzzy set membership function [7], linear and logistic regression functions [8,9], empirical likelihood ratio function [10–12], and Bayesian predictive discriminant function [13]. We will abbreviate them as FZ, LI, LO, LR, and BP, respectively. The modeling functions imply different representation and combination rules of spatial relationships [5,7,14].

We will not discuss theoretical aspects of the modeling functions here because they were amply dealt with in the above-mentioned contributions [5,7,10–14]. The focus of our study is the characterization and interpretation of prediction patterns, independently of the particular prediction models applied. Our concern is that their applications to the same input data generate prediction patterns whose scores are in entirely different units, and these are considered as not interpretable or comparable except in terms of equal-area ranking. Moreover, they require entirely different combination rules. They imply diverse

assumed interpretations of the spatial relationships expressed by the empirical likelihood ratios, and each interpretation combines ISPs with assumptions of conditional independence (between categorical, continuous field ISPs and for integration of the two types) [12]. In geomorphology, and geosciences in general, maps are often spatially correlated so that such independence seldom exists or can be hypothesized. Nevertheless, it has been found that the existence of a correlation between ISPs causes minor alterations to equal-area rankings, mostly to the lower ranks [12,15].

#### **3. Study Area, Database, and Previous Research**

The Tirano South study area, whose location is shown in Figure 1, occupies the southern half of a community area termed "Comunità Montana Valtellina di Tirano" in the Province of Sondrio, Lombardy Region, in northern Italy. It was established in 1971 for socioeconomic and environmental protection (www.cmtirano.so.it; accessed on 7 April 2021). The geomorphology and landslide processes in the area were described by Sangalli (2008) [16], who compiled the available information and constructed a database for an initial landslide susceptibility analysis. Part of the data is being used in this contribution that also covers the same study area.

**Figure 1.** Locations of the Tirano South study area, in the Lombardy Region, northern Italy.

Elevations in this young Alpine geology area range between 300 and 3000 m a.s.l. and rainfall from 700 to 1900 mm, thus providing different climatic conditions. Vegetation consists of broad-leaved forests at lower elevations and coniferous forests at higher elevations. The geomorphology is controlled mostly by glacier activity, but anthropic activity significantly affects land use at valley bottoms with tourism, agriculture, and industries. Glacial and torrent erosion, together with intense rainfalls, are the triggers of slope instabilities. Out of a variety of landslide phenomena in the area, 70 active landslides were identified from published inventories. However, they did not contain sufficient information for separating rotational from translational dynamic types. In addition to the inventories used to generate digital images of landslide trigger area locations, various cartographies related to slope instability were available for compilation and digitization into the database, as shown in Figure 2. The database offers opportunities for experimenting on spatial prediction modeling of landslide susceptibilities, hazards, and risks.

**Figure 2.** The Tirano South database is shown that will be converted into a direct supporting pattern (DSP) the 70 active landslide trigger areas, **70at,** and the three categorical and five continuous field digital maps to be converted into indirect supporting patterns (ISPs), **ulp** and **ACDIS**. The explanation is in the text.

The Tirano South database consists of digital images within a raster of 1090 pixels by 1194 lines with 20 m resolution. The study area proper within the raster covers 646,091 pixels corresponding approximately to 258 km2. The trigger zones of the 70 active landslides, converted into DSP and abbreviated as **70at**, cover 697 pixels, i.e., as a training area (TA) a little less than 0.11% of the study area (SA). Eight cartographies are used to generate, at the same 20 m resolution, the ISPs consisting of maps with 23 land use classes, **u**1–23, with 51 geologic (lithologic) units, **l**1–51, and with eight permeability classes, **p**1–8 (three categorical maps), in addition to aspect (0◦ to 359◦), **A**, topographic curvature (−32 to +29), **C**, topographic digital elevation (350–2906 m), **D**, internal relief (0–111 m for 3 × 3 pixels), **I**, and slope (0◦ to 61◦), **S**, all derived from a 5 × 5 m resolution digital elevation model, DTM, resampled to 20 × 20 m (5 continuous field maps).

Previous research by Poli and Sterlacchini (2007) [17] on landslide susceptibility in the Tirano area used an earlier and simplified database with 28 active complex landslides and five binary and binarized factor maps applying the weights-of-evidence mathematical model and generating arbitrary thresholds of prediction ranks. That and similar applications of spatial prediction modeling were criticized [18] so that in joint contributions that used the same database, kindly shared by those authors, cross-validation procedures were preferred. The purpose was to assess the relative quality of the modeling results with an alternative model and different analytical processes that exposed the uncertainties associated with the prediction patterns [19,20]. A further study with a wealth of controlled information on Alpine landslides and three national–regional mapping initiatives provided new data for the construction of a high-quality database for the Tirano study area. It contained landslide inventories in addition to geological and soil–land use cartographies. For the northern part of the Tirano study area, prediction patterns were thus generated using both active and quiescent translational–rotational landslide scarps, the empirical likelihood ratio function model, and iterative cross-validation strategies based on blind tests [21].

Indeed, the Tirano study area had become the focus of much more research on landslide susceptibility modeling and other types of landslides and associated floods due to the wealth of information available and certified. Blahut et al. (2010a) [22] used two sets of aerial photographic coverages, 20 years apart, to select natural-condition debris flows for the estimation of the spatiotemporal probability of hazard initiation. The resulting distribution map was then used to model runout zone limitations using an empirical GIS-based simulation tool. The experiments indicated that the spatial variability observed needed validation tests for satisfactory interpretation. Furthermore, another study was carried out by those authors [23] of the degree of spatial pattern agreement between different landslide susceptibility maps that showed similar debris-flow predictive power. Their database contained the distribution of 573 scarps, and one half was used for modeling a susceptibility map, while the other half was used to validate the map. The similarity of predictive power, however, did not necessarily produce similar spatial configurations of susceptibility ranks, a fact that pointed at the variability of ranks and related spatial uncertainty.

Two more activities in the Tirano study area make it an exemplary instance of supporting background for natural hazard and risk studies. Blahut et al. (2010b) [24] constructed a reliable operational inventory of landslides from three available official inventories. Their tasks were (1) to prepare debris flow and factor maps inventories for susceptibility modeling, (2) generate indices of accountability/reliability and selection of factor maps, (3) evaluate/validate susceptibility maps, and (4) compare the results of different maps for combining them into an integrated susceptibility representation. Various methods of subdivision of debris flows and factor maps were obtained, including the separation of the Tirano study area into three more congruous subareas—northern, central, and southern. Again, high spatial variability of the susceptibility ranks was observed as related to the different combinations of factor maps. This made it difficult to delimit the ranks into susceptibility classes. Blahut et al. (2012) [25] studied the historical information on landslides and floods in the Tirano area for hazard estimation and definition of tentative risk scenarios. They developed a case study to exemplify the usefulness in their database of 489 records of damaging events (from the years 1600 to 2001), for realistic scenario generation, producing damage classification maps of territorial threats.

Recent analyses of the Tirano South database [16] focused on (1) credibility analysis of a fuzzy set modeled prediction pattern of landslide susceptibility and separation of well predicting from poorly predicting landslide occurrences [26] and (2) a generalized procedural strategy for comparisons of prediction patterns of active and dormant landslides by different models [27]. This contribution wants to expand that procedure and strategy attempting to interpret the uncertainties associated with target patterns and their consequences for understanding the prediction patterns generated by the application of those very same models.

#### **4. Experimental Results**

The Tirano South database, shown in Figure 2, was reanalyzed in order to expose the uncertainties associated with prediction patterns. The stepwise procedure for generating and assessing prediction patterns, proposed by Fabbri et al. (2017) [27], was also followed here, i.e., (1) use models to obtain prediction patterns from likelihood ratios, (2) crossvalidate the patterns, (3) interpret the cross-validation results, and (4) obtain and compare the target and uncertainty patterns via equal-area ranks. Furthermore, step (5) was added, namely, analyze target, uncertainty, and 50% combination pattern relationships.

Table 1 shows the likelihood ratios for the ISPs computed for the Tirano South study area. For the ISP units and value ranges, only the ratios ≥2 or ≈2 are shown in the table. A ratio of 2 represents a normalized frequency in the presence of occurrences that is twice that in their absence. The arbitrary value of 2 was used as an empirical rule of thumb to separate the supportive ISP units or value ranges from the unsupportive ones of the proposition in (1). A ratio of 1 represents that the frequency in the presence of occurrences is the same as the one in their absence, i.e., null support of the proposition. Only the maximum values of the supportive ratios and respective units and value ranges are shown in Table 1: two land-use classes (2.29 and 4.48); eight lithology units (5.19 to 5.52); three permeability classes (1.92, 2.06 and 1.91), (categorical ISPs); one range of aspect angles (maximum ratio 2.07); two curvature ranges (11.07 and 5.62); two elevation ranges (4.24 and 2.31); one internal relief range (5.97); and one slope angle range (4.26). All prediction patterns generated by the five models have been based on the supports provided by all the ratio values of which, for simplicity, only the most supportive are listed in Table 1.

**Table 1.** Subset is listed of categorical and continuous field ISPs in the Tirano South study area and their respective empirical likelihood ratio values. They will be used for predictions using as DSP the distribution of the set of landslides, **70at** with **ulpACDIS** as ISPs. Abbreviations for categorical ISPs are **u**1–23, land-use classes; **l**1–51, lithology units; and **p**1–8, permeability classes. For the continuous field ISPs, A, C, D, I, and S are aspect, curvature, digital elevation, internal relief, and slope, respectively. Values are bold if the ELR ≥ 2.00. In Italics is the corresponding range of values, with the maximum ratio in brackets. In this reduced version, only ratios ≥2 or ≈2 are shown (table modified after Tables 1 and 2 in Fabbri et al., 2017 [27]).


#### *4.1. Use Different Mathematical Models to Obtain Prediction Patterns from Likelihood Ratios*

The empirical likelihood ratios represent the spatial relationships within the database. They were obtained using the trigger zones distribution image of the 70 active landslides, **70at**, as DSP, and the ratio transformed images of the three categorical and the five continuous-field maps, **ulpACDIS**, as ISPs. Each model was applied to generate a different prediction pattern. The models integrated the ratios for each ISP. The patterns, displayed in Figure 3, were generated by the five mathematical models each requiring different assumptions and providing incompatible scores computed from identical inputs of likelihood ratios. To identify the patterns and their inputs in the computations, a sequence of short names was used as follows: MODEL\_DSP\_ISPs, as in **FZ\_70at\_ulpACDIS** to **BP\_70at\_ulpACDIS**. The one-letter ISP abbreviations were used in analyses with subsets of them.

**Figure 3.** Prediction patterns using the different models: FZ, LI, LO, LR, and BP in (**a**–**e**), respectively, all using **70at\_ulpACDIS**. Black patches are the oversized trigger zones of the 70 active landslides, **70at**. Colors in the legend indicate groups of ranks of % of the study area (SA) in ascending order.

A convenient way to interpret the scores resulting from the modeling is by converting them into 200 equal-area ranks after sequencing in descending order. Fixed recognizable groups of ranks are then associated with pseudo-colors, as shown in the legend of Figure 3. The 200 ranks are displayed in wider groups for lower ranks of lesser concern, and in successively narrower groups for higher ranks of greater concern, e.g., 12.5%, 10%, 5%, 2.5%, 1.5%, and top 1%, of the area of the SA.

The illustration shows the prediction patterns from the five models, all overlaid with the distribution as black polygons of the 70 active landslide occurrences. The patterns represent the likelihood of future landslide occurrences in the SA.

Comparing them, we can observe similarities among higher-ranking groups in the northern part of the study area but strong differences are found in the southern part. Wider patches with high values are in Figure 3b,c, from the LI and LO models, respectively. Discontinuous patches are in Figure 3a,c,d, for the FZ, LR, and BP models. Altogether, there is a similarity between Figure 3a,d,e), and some similarity between Figure 3b,c. Note that we have to compare the same classes, fixed groups of equal-area ranks, in zones of concern. They are the zones with relatively higher ranks but located far from the known occurrences. In particular, we can focus on the top 10% ranks (90–99% and top 1%) corresponding, for instance, to higher elevations (1740–2100 and 2480–2630 m a.s.l.), high slope angles (37–57◦), and very low permeability, in addition to the particular types of land use and lithology (see Table 1). How "good predictors" are the prediction patterns? Which one is the best or preferable? For answering these questions, we will have to consider the prediction-rate curves in Figure 4, described in the next step. They will provide a measure of predictability of future landslide occurrences, i.e., in our case, the "next seven."

**Figure 4.** Prediction-rate curves are displayed generated by iterative cross validation using a sequential exclusion strategy of 7 out of 70 active landslides. (**a**) shows the aggregated iteration results for models FZ, LI, LO, LR, and BP, and (**b**) shows the individual curves for the 10 iterations of the FZ model and the aggregated curve as a thick red curve, corresponding to the red curve in (**a**).

#### *4.2. Cross-Validate the Prediction Patterns*

A cross-validation strategy was used to obtain and characterize the patterns as predictors. Iterative cross validation by the sequential exclusion of seven was tentatively selected out of the 70 active-landslide trigger zones in the DSP. It was computed to obtain the corresponding prediction-rate tables. The tables associate cumulative proportions of equal-area ranks of the prediction pattern with the corresponding cumulative proportions of validation occurrences. From the tables, cumulative curves are obtained, as shown in Figure 4. The illustration provides the prediction-rate curves from the iterative process using the five models. The horizontal axis shows the proportion of study area as cumulative equal-area ranks, each of 0.5% of the study area, i.e., ≈3230 pixels. The vertical axis shows the corresponding cumulative proportion of occurrence pixels in the cumulative class. The 697 occurrence pixels are distributed over the **70at** landslides so that the proportion of occurrence corresponds to the pixel numbers in each.

Note that the curves in Figure 4b represent proportions of seven occurrences, while the thick red curve represents the proportion of all 70 occurrences (curve **FZ\_70atm7\_ulpACDIS**, where m indicates minus).

The widespread distribution of the curves is indicative of the uncertainty affecting their aggregated thick red curve. What is predicted here is the likelihood of the "next" seven occurrences. In the iterations, they are assumed to be unknown. In all iterations, 63 occurrences are used to generate a prediction pattern that is then cross-validated as a predictor of the remaining seven occurrences. Note also that the proportions shown on the vertical axis are three orders of magnitude smaller than those on the horizontal axis (proportions of 697 pixels versus proportions of 646,091 pixels). As a reminder, the vertical axis was kept half the length of the horizontal axis.

We can observe that the prediction-rate curves are not particularly good. For instance, the thick red curve in Figure 4a shows that the top 10% ranks (0 to 0.1) contain about 14% of the occurrences for FZ, LI, and LO, and 22% for LR and BP; the top 20% ranks contain around 50% of the occurrences for FZ, LR, and BP, and about 26% for LI and LO; the top 30% ranks contain 96%, 91%, and 67% for BP, LR, and FZ, respectively, and about 57% for LI and LO.

These relative proportions represent the predictive capability of the modeling results: the proportion of predicted occurrences in the corresponding equal-area ranks. The shallow initial part of the curves is a sign of poor congruity of the setting of the landslide trigger areas used for cross validation. A good prediction pattern should provide an initially steep prediction-rate curve through cross validation in which most of the validation occurrences fall on higher ranks. The curves in Figure 4a show that the FZ, BP, and LR patterns of Figure 3 predict better than those for LI and LO. Figure 4b shows, in addition to the FZ prediction-rate curve identified with a thick red line (**FZ\_70atm7\_ulpACDIS**), the curves for each of the 10 iterations of the FZ model.

Note that a shallow initial curve is the one from the fourth iteration. Recall that in the iterations the cumulative proportions of occurrences in the diagram refers to the seven occurrences being cross-validated, while the proportions for the aggregated FZ curve refer to the 70 occurrences, validated into successive 10 groups of 7. What is being predicted here are the "next" 7 occurrences using the "previous" 63. This is the best approach when not having the time of the occurrences.

#### *4.3. Interpret the Cross-Validated Results*

The next step was to assemble sets of 10 prediction patterns to obtain via rank-based statistics the corresponding target and uncertainty patterns. The term "target" refers to what we are looking for as a validated prediction result. The term "uncertainty" refers to the stability of it. Indeed, the prediction pattern is the most informed prediction because it is using all information available, all the 70 active occurrences as DSP. However, we do not know its prediction capability. This is estimated by generating the target pattern via iterative cross validations using whichever strategy we can apply.

Figure 5 shows the uncertainty patterns obtained applying the five models to the database, the **70at** as DSP and **ulpACDIS** as ISPs. To obtain the target patterns, rankbased statistics was used to select for each pixel the corresponding median rank of the 10 prediction patterns generated by the iterative cross-validation process **70atm7**. The display of the target patterns is not provided here due to the extreme similarity with the corresponding prediction patterns in Figure 3 when using the same color legend. More informative are the uncertainty patterns in Figure 5, generated by ranking the ranges of ranks around the median ranks of the target pattern. The wider is the range the higher is the uncertainty of the corresponding target pattern (and consequently also of the prediction pattern).

**Figure 5.** Uncertainty patterns, all using **70atm7\_ulpACDIS** as cross-validation process, from the different models FZ, LI, LO, LR, and BP in (**a**–**e**), respectively. Black patches are the oversized trigger zones of the 70 active landslides. Colors in the legend indicate groups of ranks of % of SA in ascending order.

We can use the same color legend for the uncertainty ranks as for the prediction or target ranks. The uncertainty ranks relate to all target ranks from high to low. For uncertainty, however, the desirable ranks are the lower ones, indicating lower uncertainties. Obviously, if a high target rank corresponds to a high uncertainty rank, it is considered less credible than one corresponding to a low uncertainty rank. Figure 5 shows areas of high uncertainty, higher ranks, to the North in all the five uncertainty patterns. As to those along the valley edges, they are visible only in Figure 5a,b,d, for FZ, LI, and LR, respectively. Other areas of high uncertainty to the south are common to all the five patterns in the illustration.

At this point, it becomes important to study the relationships between target and corresponding uncertainty ranks with the 70 active landslide occurrences, **70at**, from the **70atm7\_ulpACDIS** process of cross validation.

#### *4.4. Obtain and Compare the Target and Uncertainty Patterns via Equal-Area Ranks*

By cross-validating the target and uncertainty patterns with the 70 active landslides, **70atm7**, the relationship between target and uncertainty ranks were visualized, conveniently expressed in \*1000 units (for instance, rank 900 corresponds to the top 10% equalarea rank). Figure 6 shows as an example the relationship between the target ranks, in descending order on the horizontal axis, and the respective uncertainty ranks, in ascending order on the vertical axis. It was obtained from the cross-validation process **FZ\_70m7at\_ulp\_ACDIS.** In the illustration, the diagram shows that the 70 points are distributed so that the higher target ranks on the horizontal axis appear to correspond to relatively lower uncertainty ranks. It is worth noting that there are five encircled points corresponding to areas of relatively high uncertainty ranks and the lowest target ranks. They are occurrences that contribute to the shallow part of the prediction-rate curves in Figure 4. They could be considered outliers amongst the occurrences used as DSP.

**Figure 6.** Target ranks in \*1000 for the FZ model on the horizontal axis versus uncertainty ranks on the vertical axis, for the 70 active landslides, **70at**, in the Tirano South study area. The five encircled points have the lowest target ranks. See text for explanation.

Such a tendency for the patterns from the five models can also be visualized by sequencing the target ranks of the 70 occurrences in decreasing order, constructing histograms with pair of columns of target and corresponding uncertainty ranks. This is accomplished

in Figure 7 that shows in blue the target ranks and in red the corresponding uncertainty ranks for each occurrence. For all models, the histograms express a preferential distribution of higher uncertainty ranks for lower target ranks.

**Figure 7.** Histograms of decreasing target ranks (blue columns) and corresponding uncertainty ranks (red columns) for models FZ, LI, LO, LR, and BP from (**a**–**e**), respectively, using the cross-validation process **70atm7\_ulpACDIS**.

We may wonder whether that tendency is visible also in the target patterns (or in the prediction patterns). This can be observed by generating combination patterns that relate uncertainty and target ranks by tentatively thresholding the uncertainty patterns. We have generated the 50% combination patterns for each target pattern, as shown in Figure 8. A threshold value was arbitrarily set at the lower 50% uncertainty ranks in the study area to select the corresponding target ranks (or alternatively, we could have used prediction ranks from the respective prediction patterns).

**Figure 8.** The 50% combination patterns, **70atm7**, using the different models FZ, LI, LO, LR, and BP in (**a**–**e**), respectively. Black patches are the oversized trigger zones of the 70 active landslides, **70at**. The dark-gray areas represent 50% of the study area with relatively higher uncertainty. The remaining 50% show the target pattern ranks for areas with lower uncertainty. Colors in the legend indicate groups of ranks of % of SA in ascending order.

By comparing the uncertainty pattern in Figure 5 with the prediction patterns in Figure 3, we can infer that the combination patterns in Figure 8 have retained high values in the southern part of the SA that shows low uncertainty (lower 50% uncertainty ranks). In the northeastern part, however, uncertainty is relatively high (upper 50% uncertainty ranks) so that the high target ranks are not visible in the combination patterns. What are then the characteristics of the 50% combination patterns? How is the uncertainty range threshold shaping them?

Let us consider some revealing details of the rank distribution in a small window of the 50% combination patterns, as shown in Figure 9. Seven landslide trigger areas are displayed as black contours in the illustration. Some are predicted as uncertain and fall on gray areas, in Figure 9b,c,e by LI, LO, and BP models. On the contrary, they have low uncertainty in Figure 9a,d by FZ and LR models. Topologically, similar dispersed patches are visible in Figure 9a,d,e, and more continuous ones in Figure 9b,c. The top 1% combination ranks, purple color, show the high variability of these clusters of pixels. These are the main characteristics of the patterns. They indicate weak robustness of the ranking and the uncertainty at the occurrence location.

**Figure 9.** A particular small subarea is shown of 50% combination patterns with overlaid boundaries of the seven active landslides located in it. In (**a**,**b**) are the patterns from models FZ and LI; in (**c**–**e**) from models LO, LR and BP, respectively. In (**f**) is the location of the subarea. Colors in the legend indicate groups of ranks of % of SA in ascending order.

Next, we can derive a more general characterization of an entire 50% combination pattern by looking at low-uncertainty proportions or pixel numbers in the individual 200 ranks.

#### *4.5. Analyze the Target, Uncertainty, and 50% Combination Pattern Relationships*

Increasing uncertainty ranks with decreasing target ranks, as we have seen with the **70at** landslide trigger areas via the **70atm7** process in Figures 6 and 7, are also found in the target patterns. We can visualize their suspected higher uncertainty (hinted by the curves in Figure 4b) for lower prediction ranks. This is carried out by generating plots similar to the one in Figure 6 but with 646,091 points corresponding to the pixels in the SA, using pairs of target and uncertainty patterns. Instead, the following simpler visualization procedure was followed.

We have generated 200 equal-area ranks out of the prediction, target, and uncertainty ranks, each corresponding to 0.5% of the study area (≈3230 pixels). Being concerned mostly with the highest ranks, e.g., the higher 80%, we have displayed the numbers of pixels in each combination pattern that corresponded with the 50% lower uncertainty ranks, leaving out those with higher uncertainty ranks. We observed, therefore, the decrease in pixel number for each rank due to the elimination of higher uncertainty ranking pixels and detected which ranks were consequently more uncertain. These were the intermediate ranks. Figure 10 compares this decrease in the 50% combination patterns from the five models with the respective target ranks, i.e., the straight lines. In all cases, the intermediate ranks show higher uncertainty, i.e., higher "loss" of target pixels in the combination ranks, between rank 160th and rank 80th. The curves in Figure 10 have strongly variable configurations with one or two concavities and different higher pixel numbers in the vicinity of the 200th and the 40th ranks. Figure 10f contrasts the curves from models FZ and LO. Note the absence of a drop in pixel number for the top ranks in Figure 10e. It shows a sharper pick, as to be expected from the yellow prediction-rate curve from model BP in Figure 4a that appears to be less sensitive to the suspected outlier occurrences.

We may wonder whether the curves indicate some form of database or modeling signature. The curves in Figure 10 show how the combination patterns reach high uncertainties at intermediate ranks. The relative uncertainties increase from initial lows to reach one or two maximum values between rank 190th and rank 60th. Since this is observed for all the patterns from the five models, it appears as a database property. Observe the prediction-rate curves in Figure 4a and consider the curves in Figure 10 showing the number of target ranks corresponding to lower uncertainties close to the two extreme ranks displayed. Would we want to select the top 10% combination ranks, from 200th to 180th, as the acceptable part? Or alternatively, would we prefer the top 20% as more significant? The properties of uncertainty and combination patterns are still unknown and remain a research challenge.

#### *4.6. Considerations on Prediction Patterns as Maps*

What resembles a map, such as the display of a prediction pattern, is not necessarily a map but more a representation of information extraction through normalizations, conversions, assumptions, and integrations. At present, with what we know about a study area, we have to be satisfied by the higher part of the prediction-rate curve and by the corresponding combination pattern, perhaps fine-tuning the process. Otherwise, additional information must be used to choose a satisfactory part of the prediction pattern. Regarding the prediction patterns, would we want to discover instances of unmapped landslide trigger zones in the higher 10% ranked areas away from the known ones? Is the 50% combination pattern helpful in mapping more of what we know or that we do not know? Should we wait for the next seven landslides to see where they will appear? Where should we concentrate on more detailed mapping? How much of the study area should we consider unfit for particular land uses or developments? Would we opt for using cost/benefit or using safety criteria? Providing answers for decision making and

subsequent risk analysis is the logical function of the prediction patterns of landslide susceptibility, as we have discussed.

**Figure 10.** Comparison of pixel numbers in equal-area target ranks with the pixel numbers in the corresponding low uncertainty target ranks of the 50% combination pattern. Models FZ, LI, LO, LR, and BP are in (**a**–**e**). In (**f**) is the overlay of curves for models FZ and LO. The target equal-area ranks contain about 3230 pixels. Only the higher 160 ranks are shown here.

#### **5. Concluding Remarks**

Favorability function modeling was applied to the Tirano South database through five different models of spatial relationships. All the prediction patterns obtained were represented as relative ranks, including their derived target and uncertainty patterns. The 50% lower uncertainty ranks were tentatively used to extract the corresponding target ranks as 50% combination patterns. In them, for all models, the proportion of less uncertain pixels in the rank represents the level of confidence in the prediction. Possibly, the top 10% combination pixels are a significant part of the prediction pattern. We are trying to find answers to questions such as the following: What are the consequences of higher uncertainty for intermediate ranks in combination patterns? How are these diagrams providing a measure of the quality of the prediction? Are they characteristic of the modeling because of the model or the data? How much of the prediction pattern is reliable? How to evaluate the higher ranks in the prediction patterns? Should we select the top 10% or 20% combination ranks? What would our cost/benefit considerations and choices be? All we have generated from the modeling are "relative" integrated equal-area ranks that must be

interpreted as susceptibility to land-sliding. The model preference needs to be a function of the interpretation of prediction patterns and their representation. In our case, FZ, BP, and LR are equally satisfactory but less so for LI and LO, as to the pattern of predictive performance. However, this may be just one of the criteria that might be used.

A five-step procedure is proposed for modeling prediction and uncertainty. The uncertainty ranks are considered important to properly select the susceptible part of the study area, i.e., susceptible to the "next seven" future" landslide occurrences. While various other strategies have also been applied for iterative cross validation to obtain slightly different prediction patterns, the sequential exclusion of seven occurrences appears to do justice to the properties of the database. The procedure is proposed as critical in spatial prediction modeling. Independently of the models used, a necessary research issue is in the interpretation of the uncertainty associated with prediction patterns.

We have described our analysis and modeling results to indicate a way to predict and the assumptions implied. Obviously, we suspect that our experiments on this particular database have a more general significance beyond the specific study area or the five mathematical models used. These considerations, we hope, will be useful to researchers and users of susceptibility maps. This contribution does not provide a solution but poses questions whose answers point at possible solutions.

**Author Contributions:** Investigation, A.G.F. and A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This contribution was initially and partly supported by the European Commission Project "Mountain Risks: from Prediction to Management and Governance" (MRTN-CT-2006-035978, 2007– 2010), Mountainrisk (2007).

**Acknowledgments:** The authors are grateful to four anonymous reviewers and to Jan Blahut (The Czech Academy of Sciences) for helping to improve this manuscript.

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

#### **References**


### *Article* **Past, Present and Future Monitoring at the Vallcebre Landslide (Eastern Pyrenees, Spain)**

**Josep A. Gili 1,\*, Jose Moya 1, Jordi Corominas 1, Michele Crosetto <sup>2</sup> and Oriol Monserrat <sup>2</sup>**


**Featured Application: Monitoring of very slow landslides with Classical and Novel Techniques.**

**Abstract:** Works carried out to monitor the displacements of the Vallcebre landslide (Pyrenees range, NE of Spain) since 1987 are presented. The landslide, which extends over an area of about 0.8 km<sup>2</sup> and affects more than 20 <sup>×</sup> 106 <sup>m</sup>3, has experienced displacements of up to one meter per year in some points and periods. It has been periodically monitored since 1987, using a wide range of surface and in-hole techniques: triangulation with theodolite, Terrestrial Photogrammetry, Electronic Distance Measurement, GNSS-GPS, inclinometers, wire extensometers, piezometers, DInSAR (satellite) and GBSAR (terrestrial). The results obtained using new techniques are compared with those obtained with GNSS-GPS and a wire extensometer, and checked against fixed stable points. From this comparison, we conclude that even though wire extensometers and inclinometers may have the highest precision, in practice, all systems play potentially valuable roles in providing meaningful data for monitoring at different study stages. In the near future, we envisage the installation of a Distributed Fiber Optic array to monitor the risk with a certain space and time continuity. After the evaluation of the precision and advantages of the different methods, the complementary use of some of them is strongly recommended.

**Keywords:** monitoring; landslides; photogrammetry; global positioning system; in-hole wire extensometer; DInSAR; GBSAR

#### **1. Introduction and Site Description**

The Vallcebre landslide is located in the Eastern Pyrenees, approximately 125 km north of Barcelona, Spain (Figure 1). Its situation, geological context and a complete geomorphological description can be found in [1,2]. The present paper is an expanded and updated version of a previous work [3] where we described the landslide as a translational slide with a stair-shape profile. As in most landslides, its structure and behavior are not simple. The landslide is 1200 m long and 600 m wide, involving an area of 0.8 km2, which shows superficial cracking and distinct ground displacements (Figure 2). The mobilized material consists of a set of shale, gypsum and claystone layers gliding over a thick limestone bed. A geological cross-section is presented later in Section 2.4. The landslide consists of three main slide units, which show a decreasing thickness towards the landslide toe. Each unit is formed by a gentle slope surface bounded in its uphill edge by a scarp of a few tens of meters high. At the base of each scarp, an extension area develops in the form of a crack system and a graben [1]. This fact indicates that the lower units move more rapidly than upper ones, which has been repeatedly confirmed by an installed monitoring network [1,2].

**Citation:** Gili, J.A.; Moya, J.; Corominas, J.; Crosetto, M.; Monserrat, O. Past, Present and Future Monitoring at the Vallcebre Landslide (Eastern Pyrenees, Spain). *Appl. Sci.* **2021**, *11*, 571. https:// doi.org/10.3390/app11020571

Received: 23 November 2020 Accepted: 4 January 2021 Published: 8 January 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/).

**Figure 1.** Location of Vallcebre valley (red arrow) within the Pyrenees range, NE of Spain.

**Figure 2.** General view of the Vallcebre landslide from W. The unstable area (yellow contour) is about 0.8 km2, has an average slope of only 17◦ (green line), and is sliding as a translational slide. The white arrows are the average slide directions of the upper, middle and lower units, which are separated by scarp zones (dashed yellow lines).

The average slopes of both the topography and the basal shear surface in the middle and lower unit are gentle, i.e., in the range between 6 and 10◦. The basal shear surface depth is 42 m in the middle unit and 15 m in the lower one, being quite planar for most of the sliding surface.

Complete information on the geology and geotechnics of Vallcebre, and a stability analysis, are beyond the scope of this article, but can be found in [1–8].

Figure 3 shows a geomorphological sketch of the landslide and the location of the monitored points and boreholes. The most active area is the lower unit, whose toe is being eroded continuously by the Vallcebre torrent. Consequently, most of the monitored points were located there.

**Figure 3.** Geomorphological scheme of the Vallcebre landslide superimposed over a vertical aerial image. Several GPS benchmarks, corner reflectors and borehole extensometers have been highlighted with symbols.

The measurement of displacements is very often the simplest way to observe the evolution of a landslide and to analyze either the kinematics of the movement, the response to the triggering conditions (e.g., rainfall) or the efficiency of corrective measures.

A variety of measuring techniques have been developed to track the movements of unstable areas (see i.e., [9,10]). This second work, also known as the 'Red Book' of Geotechnical Instrumentation and Monitoring, is due to John Dunnicliff. Several of these methods have been used in Vallcebre since 1987, beginning with "classical" surveying and photogrammetry, and using Global Positioning Systems (GPS/GNSS) from 1995. In 1996, this site was included in the framework of the NEWTECH Project, funded by the European Commission. Between July 1996 and March 1997, 14 boreholes were drilled in the slope and equipped with inclinometers, wire extensometers and open standpipe piezometers (Figure 3).

In 2004–2005, multipoint piezometers were installed in three additional boreholes. Later, the Vallcebre landslide was included in an EU-FP7-funded research project known as SAFELAND. In 2007, seven artificial radar corner reflectors (CR) were installed to test the DInSAR monitoring capabilities. In 2010–2011, still within the SAFELAND project, the lower part of the landslide was monitored using a Ground-Based SAR (GBSAR). In

this sense, the Vallcebre landslide can be considered to be a real-scale laboratory where the performance of different monitoring techniques can be assessed and compared. For instance, Laser Scanning (both Terrestrial or Aerial) have been disregarded because the ground displacements, mainly parallel to the topography, would be difficult to detect; also the vegetation would mask the results.

In the paper, Section 2 describes the measuring systems. In Section 3, some results are presented and compared in terms of precision and repeatability. Finally, some conclusions are outlined regarding the advantages and drawbacks of the systems tested in Vallcebre. The aforementioned comparison was made in terms of cost, ease of use, precision and continuity of the results.

#### **2. Monitoring Methods Used in Vallcebre and Sample Results**

A summary of the monitoring systems used consecutively from 1987 to the present in Vallcebre is given here. More details can be found in [1,4–7].

#### *2.1. Terrestrial Photogrammetry*

The first monitoring network established on the landslide was based on terrestrial (or "close-range") photogrammetry. A total of seven campaigns were performed at the landslide foot between 1987 and 1992, covering only a small area of about 100 × 50 m (Figure 4). Stereopairs were taken with a Wild P32 metric camera (Figure 4). Each campaign included three photograms that produced two photogrammetric models.

**Figure 4. Left**: example of a 6 × 8 cm glass plate photogram of landslide toe taken from the Vallcebre Torrent. **Right**: Wild P32 camera for terrestrial photogrammetry over a theodolite, which was used for the landslide toe monitoring between 1987 and 1992.

The results of each survey were the change in the coordinates of the main points (displacements). The precision ranges between 1 cm in well-defined points (pointing targets for instance) and 10 cm (rock blocks, trees and so on). Maps with contour lines were also produced (Figure 5). The interpretation of this kind of maps is not simple: the variations of the terrain surface are caused by the landslide movement, but also by the soil erosion during rainfalls and the eroding effect of the Vallcebre Torrent at the base of the toe. Although several blocks exhibit displacement of up to 8 m, the variation of the contour lines is very moderate as a whole (Figure 5).

**Figure 5.** Terrestrial photogrammetry sample result: Map of differences between two campaigns (one-meter contour-line interval) at the landslide toe (Figure 4).

At that time, the use of terrestrial photogrammetry in landslides was not straightforward, because of the difficulty of creating a proper setup with an adequate view over the hill slope. Moreover, specialized and costly equipment had to be available for a precise stereocompilation. Although out of the scope of this paper, it is worth to say that the photogrammetry processing has evolved in the last decade towards automation, ease of use and low cost, thanks to the so-called Structure from Motion (SfM) concept [11–13]. In parallel, the use of drones (or UAS) as aerial photogrammetry platforms permits to overcome the point-of-view issues during acquisition [12–17] In this way, photogrammetry is once again a powerful landslide monitoring technique, providing information that is almost continuous in space, but discontinuous in time.

#### *2.2. Triangulation and Electronic Distance Measurement*

From 1987 to 1995, geodetic measurements with theodolite and EDM (Electronic Distance Measurement) were carried out. Between 1987 and 1992, up to 17 points were triangulated in the toe of the lower unit, mainly for orientation of the photogrammetric models (Figure 6a). In the period 1988–1994, three additional points in the middle unit were monitored with single distance variation measurements (Figure 6b). During 1994 and 1995, a "triangulateration" from new base benchmarks E1 and E2 (Figure 6c) was extended to 16 points spread out through the whole landslide. The angle measurements were carried out with a Wild T2 theodolite (Wild, Heerbrugg, Switzerland, 1975), and the EDM with a Wild DIOR 3002S (Wild, Heerbrugg, Switzerland, 1988).

The rate of movement measured during the period 1987–1995 was strongly dependent on the rainfall and the target position within the landslide. Rates up to 4 m per year were observed at certain points near the toe in the rainy years, while almost no displacement occurred during periods of drought (Figure 7). In the middle landslide unit (Figure 3), the rate of displacement was significantly smaller, in the range of 1 to 30 cm/year.

**Figure 6.** Schemes for the triangulation (**a**), single distance variation (**b**) and triangulateration (**c**) methods as used in the Vallcebre landslide from 1987 to 1995.

In terms of the precision of the observations, the EDM measurements proved to be more reliable (typically 1 cm) than angle measurements with theodolites (around 4 cm at typical distances), at least with the setup, equipment and sighting distances in use in Vallcebre. This is due to the fact that the precise determination of angles needed greater experience and better environmental conditions (no mist or smog in the line of sight, proper target illumination, lack of air vibration due to strong insolation and so on), often not completely available in mountain areas at the time when the measurements have to be taken.

**Figure 7.** First tentative correlation between the displacement of a point in the lower unit (mm, black line) and the cumulative rainfall (mm, blue line), during 1994–1995.

#### *2.3. The GPS-GNSS Surveys*

Following the discussion on the precision of the observations, the expected errors for the Global Positioning System (1 to 2 cm, depending on the method in use [5]) were lower when compared to theodolite and EDM errors. Additionally, the GPS precision was more balanced in the three axes [5]. These facts led us to apply GPS techniques to perform systematic monitoring of the Vallcebre landslide. In December 1995, a complete double survey (EDM and GPS) was carried out in order to link the measurements taken with classical methods with the first GPS campaign. The equipment used then was a Trimble 4000 SSi model with two dual frequency receivers (Trimble, Sunnyvale, CA, USA, 1995, Figure 8). Currently, we are using two Topcon HiperPro dual constellation receivers (Topcon Positioning Systems, Livermore, CA, USA, 2005).

**Figure 8.** Examples of Vallcebre monitoring points with different antenna setup: telescopic pole (**a**); pole with mini-tripod (**b**); ground-plane antenna over a tripod, centered with an optical plummet over an inclinometer borehole (**c**).

Most of the old targets were recovered with minor modifications. As new points have been added since 1996, the monitoring network has now around 50 points (Figure 3), consisting on engravings in rock blocks outcropping in the hillside, steel rods, stakes, and the top of the casing of the inclinometric boreholes (Figure 8); the radar Corner Reflectors vertexes were incorporated to the GPS network as of 2006. There are seven additional points on the limestone around the sliding zone. These were the fixed points used to check the GPS accuracy. This network allowed both the measurement of displacements and the comparison with movements obtained with the borehole equipment (inclinometers and wire extensometers) and the Radar measurements.

Although the Fast-Static GPS method (more precise and robust) was partly applied in the first years, RTK (Real Time Kinematic) is currently in use for productivity reasons [5]. Currently, we can measure the entire network within a single day, travelling across the slope by car and on foot.

Fourteen GPS campaigns were carried out from December 1995 to February 1998, one survey every 2 months approximately (Figure 9a). Later on, the campaigns were continued on a yearly basis (Figure 9b), 36 GPS surveys until now. Average displacement rates derived from GPS campaigns show a sustained regime of velocities: the annual displacement depends on the precipitation of the year (total value and peak periods). More details about the GPS and special considerations for its application to landslide monitoring can be found in [5,18,19].

**Figure 9.** Examples of results of the GPS monitoring (**a**) planimetric displacements from the initial position in the period 1995–1998 (**b**) planimetric movement of one point during the yearly campaigns.

#### *2.4. In-Hole (or Geotechnical) Instrumentation*

Between July 1996 and March 1997, 14 boreholes were drilled in the slope for reconnaissance; 3 more were drilled in 2005 (Figure 10). Some of them were equipped with inclinometers, wire extensometers and open standpipe piezometers. Inclinometers and piezometers were standard devices. As an example, Figure 11 shows one inclinometer head and its typical record, where the slip surface is clearly marked. Measurements were

made every 2–3 weeks until the casing deformation prevented the safe and proper sliding of the probe.

**Figure 10.** NE-SW geological cross-section of the Vallcebre landslide according to [1,2]. Some of the 17 boreholes drilled are shown, which were equipped with inclinometers, wire extensometers (for shear displacement) and vibrating wire piezometers.

**Figure 11. Left**: view of the head of a borehole during the measurement with the inclinometer. **Right**: graph with the successive displacement profiles from August 1996 to December 1997. The main slip surface appeared to be 34 m deep. The pretty constant deformation profile corresponds with a translational slide.

On the other hand, the extensometers were wire-type, specially built following a design described in [20]. They consist of a protected steel wire anchored to the bedrock, below the slip surface, inside a piezometric pipe (Figure 12). After some computation, the wire displacement at the pulley can be related with the horizontal displacement of the landslide (Figure 13). We can quote two major advantages of this device. Firstly, with a potentiometer, it allows the continuous recording of the displacement, particularly necessary to collect information during the concentrated rainfall periods, characteristic of the Mediterranean climates (Figure 14). Secondly, the wire extensometer works properly with big landslide displacements, much larger than the 20–30 cm that would typically break the inclinometric pipe; in fact, in 2012, 16 years after their installation, 3 wire extensometers were still working, with a cumulated total wire displacement of up to 6 m. Full details of the wire extensometer can be found in [4].

**Figure 12.** Head of one of the borehole wire extensometers installed at Vallcebre. The cable for the piezometer (black) can also be appreciated.

**Figure 13.** Sketch of the borehole wire extensometer. (**a**) just after installation; (**b**) after some shear displacement of the landslide (idealized). The relationship between the wire displacement at the pulley and the horizontal displacement of the landslide depends on the width of the shear zone and the borehole diameter [4]. The smaller the borehole diameter, the faster the wire response.

#### *2.5. DInSAR Monitoring Using Corner Reflectors*

The Space Borne SAR monitoring of landslides has been reported in [21–23] and elsewhere. When applying DInSAR to landslide monitoring in forested areas, a condition that is difficult to fulfil is that a sufficient number of targets within the area of interest remain coherent during the observation period. One way to overcome this limitation is the deployment of artificial corner reflectors (CRs) that ensure coherent, high-quality DInSAR estimates. Using CRs, however, demands additional resources; it also prevents historical deformation studies based on archive SAR imagery.

After a detailed discussion of the Vallcebre site suitability for DInSAR with CRs, [6] describe the installation, in 2006, of 7 metallic trihedral CRs (Figure 15) to be observed with the ENVISAT (C-band, with a 35-day revisiting period). The Vallcebre landslide movement

is favorably oriented (towards the west) and has a gentle slope. Additionally, the deformation rates are moderate enough in order to prevent, the phase from "rolling over" (problem due to the ambiguous nature of the DInSAR). One corner (CR\_2) was installed directly on the stable rock to be used as reference. The maximum distance between the corners was 300 m to keep the atmospheric bias negligible. The installation was done with utmost care, with the CR symmetry axis oriented towards the radar line of sight [6].

**Figure 14.** Measurements at borehole S2. Above: rate of displacement derived from the wire extensometer. Below: groundwater table fluctuation measured by the piezometer. A perfect synchronism between rate and water level can be appreciated. This correlation means that groundwater has a big influence in the balance of forces controlling the landslide dynamics.

**Figure 15.** Corner Reflectors (CRs) used for DInSAR monitoring at the Vallcebre landslide [6]. Examples of installation on a big rock block (**left**) or directly over the terrain (**right**).

In Figure 16a, the CRs are clearly distinguishable in an amplitude image of the landslide. Processing interferometrically four 2007 descending ENVISAT images, the map of the rate of movements was obtained (Figure 16b). As an example, Figure 17 displays the resulting displacements, which were small between the end of December 2006 and the beginning of March 2007, whereas for CR1 and CR4 there is a substantial increase in March and April associated with an abundant rainfall period.

**Figure 16.** (**a**) the 7 Corner Reflectors appear as bright spots in a 2007 radar amplitude image; (**b**) the concentric circles on the map represent the rate of movement of the CRs during the first half of 2007.

**Figure 17.** Evolution of the displacements (cm) of the CRs during the first half of 2007. The relative Line-Of-Sight (L.O.S.) displacement was corrected to the direction of the maximum local slope. Precipitation (mm) is also presented as vertical bars.

#### *2.6. Noninterferometric GBSAR Monitoring*

The authors of [24,25] presented the application of the Ground-Based SAR (GBSAR) to landslides. In 2010–2011, the lower part of the Vallcebre landslide was monitored using a GBSAR. In this case, we did not use interferometry based on the phase but a new procedure to process the amplitude component of the GBSAR data acquired in discontinuous mode. This methodology intends to overcome several well-known drawbacks of the radar interferometry, especially in a landslide environment: loss of coherence; lapse of time between acquisitions; aliasing effect; atmosphere influence, and so on. The use of geometric features of the amplitude images combined with a matching technique is fully described in [7].

Between February 2010 and September 2011, this technique was applied to the lower unit of the Vallcebre landslide in order to evaluate the performance of the noninterferometric GBSAR approach. In this period, eight measurement campaigns were carried out using the IBIS-L Ku-band GBSAR (Ingegneria Dei Sistemi S.p.A., Pisa, Italy, 2010; Figure 18). In each campaign, 15 small CRs (Figure 19) were deployed in and around the target area (11 inside, 4 outside as reference).

**Figure 18.** Picture of the IDS GBSAR system. The synthetic aperture is obtained through the movement of the radar sensor (yellow box) along the rail.

**Figure 19.** One of the small Corner Reflectors deployed over the landslide lower unit in order to fix up the vegetation and snow problems. The operator is holding a surveying circular prism, which is observed from a reference point with a Total Station, for validation purposes.

In parallel to the GBSAR measurements, five Total Station surveying campaigns were carried out to validate the results obtained with the GBSAR. The analysis of the total measurements obtained during the 19 months permitted the validation of the technique. Displacements up to 80 cm were measured in the lower part of the lower unit (Figure 20); the movement of the points almost *en bloc* corresponds to a translational landslide.

**Figure 20.** L.O.S. displacements measured with noninterferometric GBSAR between February 2010 and September 2011 at eight small CRs. The CR14, in a stable zone, acts as reference CR.

#### **3. Discussion**

As the different systems were applied to certain common points, the resulting displacements could be compared, and practical data, including advantages and drawbacks, could be derived.

For the superficial methods, we observed that the GPS measurements showed better trends and stability than the Total Station measurements, at least with the Vallcebre conditions and equipment. An example is given in Figure 21. Moreover, the GPS surveying is "all-weather" in practice, which is an additional advantage when working in mountainous areas. A theoretical composition of errors, along with independent checks of the GPS-RTK results obtained at fixed stable points, let us establish the following standard deviations [5]: 16 mm in the horizontal plane and 24 mm in elevation. The GPS network in Vallcebre is still in use.

**Figure 21.** Apparent change in the behavior of point G2 when monitored with classical surveying (dots 1 to 12) or with GPS (dots 12 to 18). The insets with the error ellipses for both systems show a more balanced and better precision for the GPS measurements. Blue arrow points towards the Total Station.

The GPS measurements were also used to check inclinometer and wire extensometer displacements at several boreholes [4]. As an example, in Figure 22, a cross-check is presented. In Vallcebre, the inclinometric readings were possible only until 200 mm of total displacement for the top of the boreholes, due to the casing deformation in the failure zone. Nevertheless, the wire extensometers and GPS measurements have been continued without trouble beyond these figures.

**Figure 22.** Comparison of the corrected wire extensometer displacement (Dh) against the GPS and the inclinometric measurements of the Vallcebre borehole S2.

Compared with the standard GPS (campaign by campaign), the extensometric technique has the advantage of being an automatic and continuous measuring system which can be correlated with rainfall and piezometric heads. On the other hand, the GPS yields the direction of the movement and the ΔZ as well as the ΔDistance. It is worth noting that the GPS can also be installed and operated in continuous (see [18] for instance).

Cross-checking was also useful for the radar measurements, validating or improving the remote sensing techniques. Every GPS campaign included the CR and the borehole heads (Figure 23). For instance, in Figure 24 a phase unwrapping problem is recovered. The figure shows two displacement time series of CR\_4. The continuous blue line represents one phase unwrapping solution, which has an accumulated displacement of 6.9 cm and which shows an upward displacement of 1.1 cm between April and June 2008. Considering the kinematics of the landslide, this type of movement can be considered very unlikely. For this reason, a second solution was chosen, depicted with a dashed line and which shows an accumulated displacement of 11.9 cm in 11 months. The same figure shows the displacement measured by the wire extensometer S5 (green curve), which is located relatively close to CR\_4. It is interesting to compare the time series of CR\_4 (dashed line) and S5. They display quite a similar temporal evolution, which includes a strong deformation gradient between June and July, followed by a stationary period until September and a second, less strong deformation gradient. This similarity confirms the one phase shift applied to the CR\_4 curve. The curves do not fully overlap because one point is 84 m from the other.

**Figure 23.** Cross-checking between techniques: every GPS campaign included the CRs (CR\_2 in the picture) and the top of the boreholes where the wire extensometers and/or the inclinometers were installed.

**Figure 24.** Displacement time series estimated with DInSAR over the CR\_4 Corner Reflector, and measured by the wire-extensometer S5, located at a distance of 84 m. A DInSAR ambiguity is shown in red.

The GBSAR amplitude technique presented in Section 2.6 was also validated against the continuous monitoring given by the wire extensometer system. Prior to performing a proper comparison, the CR5 and CR7 Line-Of-Sight (L.O.S.) results (Figure 20) has to be projected onto the total displacement vector (estimated from GPS measurements in the area). The method can be found in [7] and is outlined in Figure 25.

**Figure 25.** The GBSAR total displacement, DT, is computed from DLOS using the 3D average direction derived from GPS (unit vector UT) [7].

After the projection, the CR5 and CR7 GBSAR displacements can be compared with the closer wire extensometer (S11 and S2 respectively, Figure 26). Dashed straight lines have a 1:1 slope, i.e., the perfect fit, which is almost reached.

**Figure 26.** Charts showing GBSAR total displacement, DT, against wire extensometer displacement, DW: (**a**) for CR5 and S11, and (**b**) for CR7 and S2.

Finally, after all this precision considerations and system crosschecks, the SAR and/or cable displacements can be visualized as 3D vectors, the 3D direction is given by the GPS data (Figure 27). This overall presentation of the results for the landslide lower unit, with 15 displacement vectors, show a displacement pattern corresponding to a very orthodox translational landslide, with a fair general trend towards the NW, and a small gradient (increase in rate) when approaching the torrent.

**Figure 27.** Displacement field of the lower unit of the Vallcebre landslide reconstructed using GBSAR CRs and borehole wire extensometers (S). The yellow figures, in cm, are the total displacements observed from February to November 2010 [7].

#### **4. Summary and Conclusions**

A general overview of the methods used in Vallcebre for measuring the displacements of the landslide has been given. More details of the systems presented here, along with the full formulation, validation and conclusions, can be found in [1,4–7]. To conclude, we present here some lessons learned during the progressive monitoring of the Vallcebre landslide since 1987.

A summary of the main characteristics of the methods is given in Table 1. This table is in accordance with [5,9], but the values are valid only for the specific constraints and equipment in Vallcebre. The characteristics of the methods are summarized along with the precision and main results achieved. The 'Precision' figures seek to characterize the 'real' or field performance of the systems in our landslide instead of the theoretical values that can be found in the instrument's technical specifications. The latter were determined in lab conditions and were found to be too optimistic, in general. The 'Complexity' column aims to quantify the refinement of the equipment and personnel involved, and the 'Cost' column includes the amount of instrumentation, installation and work needed to obtain the data.


**Table 1.** Overview of the methods used in Vallcebre landslide monitoring.

In-hole

a: Wire extensometer. After calibration of the wire equation, 2 mm in the surface displacement can be achieved. b: The application of Terrestrial or Low-aerial photogram. is becoming easier because of drones and VSfM based SW. c: The technique, with a slightly different setup, may become continuous (continuous GPS, Robotic Total St., In-Place-Incl. d: GPS, Wire extensometer and Piezometer are green because they are still operational. e: Distributed Fiber Optics is red because it has not yet been applied.

The methods in the table should not be considered as excluding alternatives, but rather, as complementary systems to be applied progressively or in different areas, as has

ous/continuous download

been explained in Section 2. For instance, the photogrammetry was used in the toe of the slide to monitor a fairly small but very active area. EDM and triangulation (and later the GPS) were implemented to cover the whole hillside and capture mean movements as big as 1 m per year with minimum investment.

The surveying with EDM and GPS were carried out with a given frequency (i.e., campaigns each two months or each year), so the results are discontinuous over time. As stated in note 'c' in Table 1, although not implemented in Vallcebre, it is technically possible to automate the procedure for the continuous monitoring of the displacements.

On the other hand, the GPS results helped to calibrate the parameters of the wire extensometer equation [4]. The installation of these in-hole devices was possible only when additional funding was obtained in order to drill the boreholes and instrument them. The inclinometers and the in-hole wire extensometers complement each other very well. As the landslide was fairly active, the inclinometers had a short 'life', but produced high-quality information on landslide displacement profiles, velocities and the position of the shear surface quite immediately after their installation. In contrast, at the early stages of deformation, the wire extensometers may only record negative displacements, which are not easily related to the superficial ones. However, once the inclinometers were lost, the wire extensometers provided continuous recording, even for very large displacements. Since 1996, when they were first installed, the wire extensometers have provided a measurement readout each 20 min; at some points, the cumulative displacement has been as high as 6 m. The accelerations in the rate of displacement can be easily related to rainfall and groundwater rises (Figure 14), especially during critical rainy events, when other systems are not in operation. Most of the piezometers and 3 wire extensometers installed in Vallcebre were still operational in 2012. Around 2010 an automatic meteorological station was installed in the center of the landslide in order to get the significant precipitation, because the rainfall may vary sharply from one point to another.

Redundancy between methods is advisable in case of malfunctions in the devices or as a means to cover gaps in data. As shown in the previous sections, all the results are in fairly good accordance when adequate corrections are applied (Figures 22 and 23). Prior to the comparison of the different techniques, some error considerations must be made, in particular, the filtering of any systematic errors (Figures 24 and 25, for instance). In order to 'automatically' filter some systematic errors, during the monitoring operation (acquisition and data handling) it is convenient to be very 'systematic' (to follow exactly the same field and processing procedures).

Contrary to structural monitoring, when working with natural materials like soils and rock, we have to accept natural variability. To address this variability, it is better to apply continuous techniques, both in the time and spatial domains, i.e., the two axes presented in Figure 28. It must be highlighted that the tentative classification of the different techniques in the figure (time frequency/spatial resolution) corresponds to the Vallcebre landslide implementation as reflected in this paper; other setups could be adopted at Vallcebre and other sites. The actual installation depends not only on the site, but also on the level of risk and the available budget.

The systems that were operated continuously relies on automatic acquisition. We learned that 'automatic' does not means 'free of maintenance'. On the contrary, the continuous systems need continuous surveillance and maintenance [2]. Another recommendation for long-term monitoring networks is that, when feasible, it is better to avoid points or benchmarks that may exhibit their own movements (large boulders, light poles, trees, concrete pillars, old stone walls and so on). The discreet marking of points diminishes natural disturbance and vandalism. Finally yet importantly, it is mandatory to include within the network several fixed points in stable areas around the landslide for verification purposes.

Over the years, some improvements in the electrical supply (solar panels) and data transmission (remote download) have been introduced in the Vallcebre setup. Currently, wire extensometer and the piezometer readings can be downloaded from a remote location by means of the mobile phone network.

**Figure 28.** Tentative disposition of the different techniques as specifically used in Vallcebre: following [26], systems are classified according to their continuity in the Spatial and Temporal domain. Distributed Fiber Optics is in red because it has not yet been applied. The asterisks indicate that the technique, with a slightly different setup, may become continuous over space or time.

All the works described in this paper are capable of qualifying the Vallcebre landslide as a real-scale "in situ" lab or observatory. The relatively slow rate of motion and the sustained long-term mobility make it suitable for the testing and validation of several new systems within the frameworks of European and National Projects. In the near future, we plan to install an experimental array of distributed fiber optics [27] in order to assess its field performance. The projected array includes around 600 m of FO buried in a trench besides a country road, roughly along the 1000 m contour line shown in Figure 27. This area is subjected to shear displacements produced by the landslide lateral limits and the transition from the middle to the lower unit. The wire extensometer and the GPS control points present in the zone will permit the cross-calibration of the FO array results and performance.

Basic research on the performance of the different instruments and methods that have been applied in the Vallcebre landslide should help to choose a monitoring strategy in new sites. In making such choices, we must bear in mind that the monitoring must address some fundamental geotechnical questions [28,29]. Instruments and methods that help to answer these questions should be selected. As John Dunnicliff stated, "If there is no question, there should be no instrumentation".

**Author Contributions:** Geological investigation and methodology, J.A.G., J.M. and J.C.; Radar monitoring concept and processing, M.C. and O.M.; data analysis and validation, J.A.G. and J.M.; supervision and search for funding, J.C. and J.M.; writing, J.A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The three decades work summarized here has been funded mainly by EC and Spanish Research Projects. In chronological order: TESLEC (*Temporal stability and activity of landslides in Europe with respect to climatic change*; EC; contract EV5V-CT94.0454); NEWTECH (*New technologies for landslide hazard assessment and management in Europe*; EC contract ENV4-CT96-0248 and Spanish CICYT contract AMB96-2480-CE); CICYT project (*Cartografía de la susceptibilidad a los deslizamientos del terreno en el área del Alto Llobregat* ... *Vallcebre mediante los SIG*; Spanish CICYT contract AMB97- 1091-C06); MODEVALL (*Modelación hidrológico-mecánica del deslizamiento de Vallcebre*; Spanish Science and Education Ministry, MEC, contract CGL2005-05282); SAFELAND (*Living with landslide risk in Europe: Assessment, effects of global change, and risk management strategies*; EC, 7thFP, grant agreement 226479); PYRMOVE (*Research network* funded by AGAUR, Government of Catalonia; grant number 2014CTP00051); PYRMOVE (*Prevention and cross-border management of risk associated with landslides*; funded by Interreg V-A, POCTEFA, and by the European Regional Development Fund, ERDF; project EFA364/19). The last project is still active.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data was created or analyzed in this study. Data presented here was gathered from previous published work [1–7]. Most of the Vallcebre monitoring data is not publicly and/or directly available due to the highly heterogeneous nature of the results, which in a significant proportion belong to the last century (paper-based), prior to the standardization of public datasets and digital repositories. However, in case of interest for some pieces of data, the authors are open to requests.

**Acknowledgments:** We would like to thank the individuals from the UPC and from the former Institute of Geomatics who helped us with the field work: J. Spagnolo (in memoriam), A. Lloret, A. Ledesma, T. Pérez-Revilla, J. Marín, D. Serral, A. Pérez-Carreras, E. Sanjuán, J. Rius, D.González, G. Luzi, E. Biescas and M. Cuevas. The authors wish to thank ICGC (Institut Cartogràfic i Geològic de Catalunya, formerly ICC) for the geodesic and photogrammetric work, ESA for the radar data, and Topcon, Trimble and Leica (formerly Wild) for their equipment and software assistance along the years. The comments and suggestions given by three anonymous reviewers during the editorial processing are also greatly acknowledged.

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

#### **References**


### *Article* **New Landslide Disaster Monitoring System: Case Study of Pingding Village**

### **Yao Min Fang 1, Tien Yin Chou 1,\*, Thanh Van Hoang 1, Quang Thanh Bui 2, Duy Ba Nguyen <sup>3</sup> and Quoc Huy Nguyen 2,4**


Received: 19 August 2020; Accepted: 18 September 2020; Published: 25 September 2020

**Abstract:** The Linbeiken area is located in the village of Pingding, Taiwan. Since the Mindulle and Aere Typhoons in 2004, and as a result of the landslide triggered by the continuous heavy rainfall on 9 June 2006, there has been a persistent collapse of side slopes in the area. This paper describes the equipment that was installed to collect on-site topographic and hydrological information in the Linbeiken area upstream of the Pingding River and to monitor changes in the landslide area, as well as the measurements that were collected during the 2008 Typhoon Sinlaku. A case study of a landslide in Pingding, Taiwan was used to monitor the accurate coordinate changes in the potential landslide areas during typhoons. The goal of this study was to establish warning indexes, and to strengthen the software and hardware at the local disaster response center in the hope of gaining a full idea of the surface movement in landslide areas in future flood seasons. This is important for boosting the preparedness to adapt to landslide hazards, for improving disaster warnings, and for reporting efficiently to better protect the lives and property of local residents. The results show that the landslide disaster monitoring and warning system in Taiwan, as applied during Typhoon Sinlaku in 2008, is both effective and comprehensive.

**Keywords:** Global Positioning System (GPS); landslide; disaster prevention monitoring; disaster mitigation

#### **1. Introduction**

Landslides cause significant damage to both people and property. Many studies have researched landslides, such as that of Chen et al. [1], which showed the application of the three-dimensional deterministic model to a landslide event in Taiwan. Hsu and Liu integrated the TRIGRS (Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability) and DEBRIS-2D (debris flow-two dimensional) models for landslides in Taiwan [2], while Bunn et al. and Ramos-Bernal et al. [3,4] presented the LIDAR (Light Detection and Ranging) digital elevation model and ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) imagery research on landslides. Liu et al. [5] used a geographically weighted regression model for studying landslides in the QingChuan area of China. Liu et al. [6] used a variety indexes like the C-, X-, and L-band synthetic aperture radar (SAR) datasets for landslide analysis in China. Assilzadeh et al. applied GIS application for landside prevention Penang Island in the Straits of Malacca [7], while Long et al. [8] analyzed the mapping of rainfall-induced landslides in Vietnam. Chou et al. applied unmanned aerial vehicle technology to produce a digital elevation model-based dataset witha5mresolution for disaster monitoring and management operations [9]. In recent decades, many different techniques have been developed that focus on early warning disaster monitoring and management systems. The study by Anita et al. [10] was based on using sensor nodes to analyze the data for disaster monitoring. Kwak analyzed multiple data sources of various satellite for disaster risk reduction [11]. McCarthy et al. created a GIS (Geographic Information System) expert system framework for the detection and monitoring of hazards [12]. Arattano et al. combined analytical sensors for the warning and monitoring of debris flow [13]. Xuan et al. provided vivid examples of and suggestions for improving monitoring systems [14]. Lucas et al. integrated characterization and monitoring [15]. Trocone et al. analyzed slope response by using the material point method [16].

With the Taiwan Strait on the left, connecting the island and Eurasia, and the Pacific Ocean on the right, Taiwan is located in one of the areas around the world that experiences the most monsoons. Prevailing winds change significantly as seasons alter; the rainy season and typhoons bring abundant rainfall and often cause landslides and landslips. Landslides are referred to as heterotaxy and occur instantly, while landslips are slow and are known as long-term heterotaxy. The Global Positioning System (GPS), which has been widely used in many academic subjects and fields, has proved to be a powerful tool for monitoring artificial and natural structure deformations, as well as slope displacements. GPS has a great number of advantages compared with traditional measurement technologies. Concisely, GPS is more precise, effective, and highly automatic, with low labor intensity required.

Wu analyzed a landslide event in Taiwan (during the 2009 Typhoon Morakot) by combining landslide susceptibility data and drawing a map [17]. In recent years, many different types of technologies have been developed, such as in case studies of Korea and Germany [18–20], which were focused on several technologies for disaster mitigation. Estrela et al. [21] used cameras to observe ground movement and reduce associated risks, while Zhao et al. researched SAR images for landslide monitoring [22]. An operational web-based GIS for the early warnings of landslides in Bangladesh, created by Ahmed et al. [23], is capable of providing alerts five days in advance. Other research works have resulted in the flooding model [24], the flooding forecasting system [25], and dynamic modeling for reservoir watersheds [26], which are all widely applied in disaster management.

On 25 June 2003, the 42–45 K section (Wuwanzi) of A Li Mountain Highway was affected by the most serious landslide in its history due to continuous torrential rain. With ceaseless rain-wash, the slumped section enlarged, and an area of nearly 120 m slipped off of the valley. As a result, the A Li Mountain Public Works Section commissioned Feng Chia University to install various transmission devices and equipment in the 31 K section of the Taiwan 18th Line Highway. As an outcome of this installation, it was suggested that reference values for rainfall levels, water levels of barrier lakes, drilled underground water levels of deposition areas, and overflow water levels be established and then used as alarms when exceeding such values.

Abnormal phenomena such as up-warp and crack-openings once again affected the dipped-slope site of the collapse zone in November 2002. To understand whether the Jiufenershan Collapse Zone was continuously sliding, the Soil and Water Conservation Bureau immediately organized the "Jiufenershan Collapse Zone Observation Project" to gain a primitive understanding of land surfaces. Observations of underground changes via geological drilling, land surface changes observation, the interior inclination of drilling holes, and underground observations can be used as references for emergency responses during typhoons or torrential rain [27].

The village of Pingding in the county of Yunlin, Taiwan was chosen as a study site. This village is located in a valley surrounded by steep mountains. The dynamic sediment process and the arrangement of the sensors are supposed to provide alerts six hours in advance of any potential landslides so as to protect this village. The advantage of this study is its demonstration of a real-time monitoring station in Pingding that provides the residents with a landslide warning and management system. The aim of this research was: (1) To receive real-time data and to monitor the entire study area (high landslide potential area); (2) to integrate software and hardware for preparedness to adapt to landslide hazards; (3) to send warnings to residents to avoid serious landslides. The research procedure is indicated in Figure 1.

**Figure 1.** Flowchart of the research procedure.

#### **2. Material and Methods**

#### *2.1. Study Area: The Village of Pingding*

Pingding is situated in the township of Linnei in the county of Yunlin, Taiwan, and its geographical coordinates are 23◦45 36 N and 120◦36 38 E (Figure 1).

Linnei is divided into mountainous areas and plains by the Taiwan 3rd Line Highway. With mountains on its southeastern side and plains in the northwest, Linnei's altitude is 70–320 m above sea level, and its topography is a slope extending roughly from the southeast to the northwest. Pingding is situated on a mesa-type hill, with an altitude of 325 m above sea level. With a tropical and humid climate, the village has a yearly rainfall of 1500–2000 mm, which is concentrated in the summer. The exposed geological formation is the Toukoshan Formation Huoyanshan Conglomerate Section; the gravel components of the conglomerate are of less than 20 cm in size and are unevenly distributed. The sandstone and mudstone, which are sometimes mixed in the conglomerate, are mainly composed of sandstone and quartz sandstone. The exposed parts of the conglomerate often form cliffs or precipices and gorges if there are any valleys. The Pingding Collapse Zone is located within Pingding

in the mountains east of Linnei. Linnei, which is located in the northeast of Yunlin county (one of the counties in Taiwan), adjoins Pingding, Pingding Mountain, and the township of Zhushan, Nantou.

In terms of its geographical location, Pingding is located between Shexi and Gaoxikou in the township of Gaoshu (Figure 1). The study site is flat, tilting north and west, and the west is similar to a valley (Figure 2).

**Figure 2.** Pingding in Taiwan: (**a**) Pingding in the township of Linnei; (**b**) Pingding in Taiwan; (**c**) satellite image of Pingding.

According to the "Central Taiwan Geological Map—From Zhushan to Jiayi" document from the fourth issue of the Central Geological Survey Article Collection and Annual Report of Taiwan Council of Agriculture by the Soil and Water Conservation Bureau, the exposed strata near the Pingding Collapse Zone include the Huben Stratum and the Sanxing Stratum of the Pleistocene Age, as well as the Accretion of the Holocene. Judging from the positions and patterns of these strata, the Pingding Collapse Zone is located in the range of the Huben Stratum and the Sanxing Stratum, while both of the shores of the Pingding sub-stream are in accretion. According to the on-site geological surveys and the six (BH-1–BH-6) (geographical surveying site from BH-1 to BH-6) 12–30 m deep geological drilling holes, the strata can be further divided into reddish brown silty clay, brown sandy silty soil in the colluvial rocks, brown silty soil in the gravel stratum, and brown or grey siltstone.

#### *2.2. Analysis of the Causes of Collapse in Pingding*

The southwestern hills of Pingding belong to the Toukoshan Formation; the poor cementing of the hills' conglomerate and sandstone tissue raises the possibility of collapse due to washouts. As a result of the frequent headward erosion in the Linbeiken area (Pingding River), the platform has shrunk in the southeastern direction, leaving only the Pingding area behind. Enveloped in serious erosion, the entire mountain area collapses easily when rainfall occurs. The proposed cause of this quick downward flow of the seepage water in the above soil stratum is because of the excellent permeability of the red soil and gravel strata on the surface. The seepage water, however, gets stuck here when confronted with gray sandy mud strata, which are poor in permeability and move horizontally. The underground water finally causes a collapse due to the decompression effect following seepage through the soil from the side slope surface. As the water content in the mud stratum increases, the mud stratum's ability decreases its shear resistance strength. These are key factors that contribute to the damage caused to the side of slopes, and can be judged from the locations and types of previous collapses. Generally, large-scale types of mass collapse are the main factors causing damage, while collapse locations are influenced by the load of slope tops, the water content of local soil (rock) strata, and the gradient.

#### *2.3. Historical Disasters*

The 921 Earthquake in 1999: The serious collapse areas that developed during the 921 Earthquake included the northeast slumping side slopes and the northwest Linbeiken Collapse Zone. The Linbeiken Collapse Pit experienced collapses of different scales, and despite many restoration works being done, the Linbeiken Collapse Pit has doubled in size in the past decade. The Yunlin County Government completed construction of the reinforced soil retaining wall in April 2003.

Typhoons Mindulle and Aere in 2004: Landslides once again struck the south bank of Linbeiken in September 2004 (Figure 3). The areas with reinforced soil retaining walls and part of the north-side tea garden collapsed into the Linbeiken Collapse Pit. The collapse created a stage terrain with a 30 m height difference. The collapse area of the first landslide covered approximately 1 ha. Afterward, it was followed by several collapses and slumps of different scales. The damage caused by this collapse included the slump of a 2 m2 tea garden, but there were no casualties. Residents from four dwellings moved for safety purposes. Regarding public facilities, a diversion ditch of approximately 23 m in length and three areas of PVC (polyvinyl chloride) piping collapsed and were pulled apart.

Continuous torrential rain in June 2006: The flood on 9 June 2006 caused the collapse of the gabion slope protection on the north side of the Linbeiken source bank and the wash-away of part of the Yunlin 61st Line Highway roadbed.

22 May 2007: At 6 a.m., with no influence of rainfall, the original collapse area of the Yunlin 61st Line Highway suffered from another topple, damaging two-thirds of the main access road of the nearby villagers, making it impassable. Emergency work such as fence establishment and warning sign installations were carried out. The main reason for this collapse was the gravity effect influenced by saturation of the soft soil caused by the cracking of tap water pipes of the area, which led to the seepage of tap water and drizzle on previous days.

#### *2.4. Data Collection*

This section shows the organization of the research methods used. Figure 4 provides a flowchart of the research methods adopted to create a disaster monitoring and management system. First, for simulating the landslide disaster monitoring system, we visited the site to carry out surveys. Second, based on the results of the survey, we decided on and installed the appropriate equipment for the study site. Third, the thresholds for the sensors were defined. Fourth, based on the results of the analysis, we identified the accuracy of the measurements and validated the threshold value. Finally, a landslide monitoring displaying system in real time was constructed.

**Figure 3.** Collapse areas of Pingding River upstream, Linbeiken, and Pingding across different periods of time [28]. Note: The yellow parts refer to the collapse on June 2002 (caused by Jiji Earthquake and Typhoon Toraji); the black parts are the collapses on 9 September 2004, while the parts in red are the collapses on 11 January 2005. The arrows indicate the sliding directions of each collapse.

**Figure 4.** Research methodology flowchart.

#### 2.4.1. Equipment for the Collapse Monitoring System

The Linbeiken area in Pingding, upstream of the Pingding River, has faced constant side slope collapses during torrential rain seasons since the 921 Earthquake in 1999. To assess the condition, this study monitored the collapse changes of the Pingding area and established a side slope collapse monitoring system to understand the real-time activities of the collapse zone during flooding (Figure 3). The following work was planned: The establishment of relevant monitoring equipment, including camcorders and devices and equipment related to rainfall monitoring and land surface deformation monitoring, as well as recording, storage, transmission, and display systems for monitoring the above equipment. In light of the simulation analysis results of Kun-yi Chen [29], the completed work included an inbuilt GPS-4 and GPS-5, a retractable cable meter (group 3), a retractable cable meter (group 4) on the sliding surfaces, an inbuilt retractable magnetic induction meter (groups 1), and a retractable meter (groups 1 and 2) in two open spaces that contained obvious cracks and may have suffered from possible continuous displacements. Finally, GPS-1, GPS-2, GPS-3, and slope circles were placed on the top of the fixed soil retaining wall. The details of the monitoring device numbers are shown in Table 1; the monitoring station is equipped with the following sensors: (1) Two sets of rain gauges; (2) three sets of CCD (A charge-coupled device) camcorders; (3) two sets of high-pressure illuminators; (4) one set of remote camcorders; (5) four sets of retractable meters; (6) one set of slop circles; (7) one sets of GPSs. The device arrangement is displayed in Table 1, Figure 5, which highlights the slope deformation monitoring system contained in the GPS receiver base station and five GPS receiver mobile stations. The precision of the real-time coordination solution was ±3 cm, while the solution speed was over one piece of information (including one) every 3 s.



GPS, Global Positioning System; 3D, three-dimensional.

**Figure 5.** Arrangement of the devices presented in Pingding. The monitoring station is equipped with a rain gauge, a CCD camcorder, a high-pressure tiltmeter, a remote camcorder, a geophone, and slope circles (dual direction).

#### 2.4.2. Establishment of Management Values

The fact that side slope failures often take place after rainstorms has led to the belief in a certain relationship between side slope failures and rainfall. Indeed, this relationship between side slope disasters and rainfall has been noticed and proven by researchers such as Brand [31], Lumb [32], and Slosson and Larson [33]. Nonetheless, some other prior hydrological conditions are essential for the triggering of side slope failures. Pre-storm rainfall injects water into side slope surfaces and allows water to flow freely in the slope. In other words, the soil surface needs to be saturated to trigger the following rainstorm's side slope failure mechanism. The pre-storm rainfall required by side slope failures depends on the soil surface cover, the water conduction capacity of the soil, the seepage rate, evapotranspiration, and the hydrological condition of the side slope.

Knowing that the influences of pre-storm rainfall on side slope stability have been studied for years, Lumb [32] discovered the impact that pre-storm rainfall has on side slope failures. In particular, he found out that higher pre-storm rainfall leads to worse side slope failures. He categorized side slope failures caused by rainfall into the following four groups:


#### **3. Results**

#### *3.1. Monitoring Warning Systems*

The data of different devices can be displayed in the following ways: Image data displays, cover realtime image displays, and historical image displays, as demonstrated by Figure 6.

**Figure 6.** Real-time image displays of the status of the study site (green dots are normal status).

Two types of real-time image displays are available: Single-station displays and loop displays of data from all of the observation stations. Historical image displays allow users to browse the CCD images of certain chosen time periods through the single-station display mode. The slope displacement condition is shown through the displays and analyses regarding the space transformation condition of the retractable meters (tensile strain meters) and the GPS slope land deformation monitoring systems. Users are allowed to access the displacement condition through direct clicks on the devices shown in the device arrangement figure, which is a vigilance data reference provided for involved staff, as displayed in Figure 7. There are four groups of retractable meters, while the data displays show the conditions of the devices using red, yellow, and green lights, as per the display mode of the device set up procedures. When the displacement quantities of the retractable meters exceed the set threshold of the vigilance values, the system will switch the light to red or yellow to remind the responsible personnel of the meters' current condition. The GPS slope land deformation monitoring systems also display the devices' condition through red, yellow, and green lights, and the lights are switched to red or yellow as a reminder for personnel responsible for the current GPS conditions, as shown in Figure 8. In both situations, users are also permitted to click on the red or yellow light icons to inquire about the detailed displacement condition of a particular piece of equipment. The slop circle data displays analyze and display data, focusing on the condition of the angle changes of the slop circle. Users are allowed to access the displacement condition through direct clicks on the devices shown in the device arrangement figure, which is a vigilance data reference provided for relevant staff, as displayed in Figure 9. The real-time monitoring data of all of the devices and the arrangement figures of the relevant devices are shown in Figure 10.

#### *3.2. Monitoring Data Analysis*

During Typhoon Sinlaku, the emergency response team of the Soil and Water Conservation Bureau was active from 08:40 on 11 September 2008 to 18:32 on 19 September 2008. It is a hyetograph of the total cumulative rainfall of the entire island, while the cumulative rainfall of the Yulin area, which reached over 300 mm, is demonstrated in Figure 11.

**Figure 7.** A detailed data display of the retractable meter.

**Figure 8.** A detailed data display of the GPS.

**Figure 9.** A detailed data display of a slop circle.

**Figure 10.** Integrated monitoring data display.

**Figure 11.** The cumulative rainfall from 12 to 16 September during Typhoon Sinlaku in 2008.

#### 3.2.1. Rainfall Data

Influenced by Typhoon Sinlaku, the total cumulative rainfall from 11 September to 16 September at the observation stations was 349 mm, with the greatest rainfall intensity being 40.5 mm/h. A hyetograph of the cumulative rainfall during the time when the emergency response team was active, showing the rainfall intensity and the cumulative rainfall conditions, is provided in Figure 12. The yellow line in the figure represents the vigilance values, while the red line refers to the action values. The hourly rainfall reached 12 mm, and the daily cumulative rainfall was 86.5 mm at 6 a.m. on 14 September, which reached the threshold of the yellow level of vigilance. The hourly rainfall reached 20.5 mm, and the daily cumulative rainfall was 140.5 mm at 2 p.m. on 14 September, which reached the threshold of the red level of vigilance. No abnormal phenomena were observed in the comparison between the monitoring results of the GPS and retractable meters (Figure 13). The local village heads, the Yunlin County Government, and the related soil and water conservation departments were informed through the disaster report process of this system, as shown in Figure 14. It is assumed that the topographic features of the Pingding area, which increase the possibility of torrential rainfall within short periods of time, lead to higher rainfall intensity and heavier cumulative rainfall. Therefore, more data concerning local typhoons and torrential rainfalls will be collected in the future to determine whether the vigilance and action values of the rainfall in this area should be raised.

**Figure 12.** Real-time rainfall hyetograph during Typhoon Sinlaku.

**Figure 13.** Relationship between Hourly Rainfall and GPS displacements.

#### 3.2.2. GPS Monitoring Data

Figure 15 shows the GPS monitoring data in normal time without typhoons, as well as a comparison between the data and management values. As shown in the figure, all of the GPS points were in a stable condition, and the vibration of the displacement quantity was within 2 cm, indicating that its precision was approximately 2 cm, thus not exceeding the set threshold of the vigilance values. The GPS displacement of all of the action stations did not exceed the set threshold of vigilance values during Typhoon Sinlaku.

#### 3.2.3. Monitoring Data of the Retractable Meters

The monitoring data from the magnetic induction and cable retractable meters during Typhoon Sinlaku, which are shown in Figure 16, were compared with the management values. Except for the micro-displacements caused by the indium steel wires of the cable retractable meters, all of the retractable meters remained in a stable condition and the set threshold of the vigilance values was not exceeded.

#### 3.2.4. Slope Circle Monitoring Data

The slop circle monitoring data from the Typhoon Sinlaku period are displayed in Figure 17. The figure shows two bigger inclinations, which indicate soil retaining wall inclinations caused by the increase in soil pressure resulting from an increase in the water of the soil after rainfall. The data return to the initial position after the water in the soil evaporates or runs off. The data were compared with those of the management values, and all of the slopes were found to remain in a stable condition and the set threshold of the vigilance values was not exceeded. Figure 18 shows a detailed data display during Typhoon Megi—International Number 1013 is the 13th tropical cyclone named in the Pacific typhoon season in 2010, and was also the strongest storm of the year.

**Figure 15.** GPS-2 horizontal displacement quantities.

**Figure 16.** Detailed data display of the retractable meter during Typhoon Sinlaku (17–18 October 2010).

**Figure 17.** Detailed data display of the slop circle during Typhoon Megi (17–25 October 2010).

**Figure 18.** Relationship between rainfall and GPS displacement data.

#### 3.2.5. Comparison between the Monitoring Devices and Rainfall

A synthetic comparison between the monitoring results is advised if the monitoring data surpass the threshold of the management values to determine whether a disaster has struck or whether an abnormal phenomenon of a single device has happened. As Figure 18 shows, a comparison between rainfall and GPS displacement involves determining the relationship between cumulative rainfall and displacement, as well as the relationship between rainfall intensity and displacement speed. Figure 19 shows a comparison between the CCD images of the time periods with the greatest rainfall intensity, aiming to ensure the rainfall monitoring is mistake-free.

**Figure 19.** *Cont.*

**Figure 19.** Relationship between rainfall and CCD data: (**a**) Cumulative rainfall of Pingding; (**b**) CCD camera 1; (**c**) CCD camera 2; (**d**) CCD camera 3.

During a typhoon and heavy rain, we can observe this system to determine the real-time accumulated rainfall intensity. This monitoring process demonstrated that there is no displacement of the ground or the slope of the retaining wall, which proves that our system is workable.

#### **4. Discussion**

According to the obtained rainfall data, Lumb established the scope of events of different gradations and provided an explanation based on 15-day pre-storm rainfall and 24 h rainstorms. Severe events brought over 100 mm of rainfall in 24 h and over 350 mm of pre-storm rainfall in 15 days, while serious events generated more than 100 mm of rainfall in 24 h and 200 mm of pre-storm rainfall. This proof of a relationship between side slope failures and rainfall encouraged many studies focusing on the threshold of rainfall for causing side slope failures. For example, it was determined that the mechanism of side slope failures is successfully triggered when the rainfall in Southern California reaches 140% of its normal rainfall (an average value of the recorded data of over 100 years) [34–36], suggesting that the threshold for Los Angeles is 125% (1993), which was also confirmed by Wieczorek [37]. Auer and Shakoor studied the collapse distribution of Nelson County, central U.S., after Hurricane Camille in 1969. According to their research, the hurricane, which moved from the west to the east, caused a more serious collapse of the side slopes in the west, northwest, and southwest directions. This indicates that there might be a connection between the directions of side slopes in collapse zones and hurricane routes [38]. Hong Kong has no records of collapse generated by earthquakes; however, the Geotechnical Engineering Office compared the risk of earthquakes and rainfall in terms of their generation of artificial slope collapse [39]. The study results showed that earthquakes bring a lower risk of causing artificial slope collapse than rainstorms.

Brand, Permchitt, and Phillipson [40] examined the side slope failures of Hong Kong through the following three steps:


According to the above research, based on the relationship between disaster numbers and rain intensity, the authors [40] concluded that over half of the studied side slope failures in Hong Kong were caused by local brief rain showers, and that the threshold for triggering side slope failures is crossed when the rain intensity reaches 70 mm/h.

Shi-Jie Jian [41], who discussed the influences of different factors on the mechanism of side slope slips through on-site measurement data, focused on the slopes (28 K + 900 to 31 K + 500) along the roads near Wuwanzi in Gongtian village, which is in the township of Fanlu in Jiayi county. As part of the old landslip region, this area has suffered from landslips for a long time—since the opening of the road. Several measurement devices have been installed for two-year periods since the beginning of 2000, carrying out consistent observations of local strata, land surface deformation, the underground water level, and rainfall. Relevant discussions have also been carried out in coordination with basic theory and indoor experiments. According to the observation results, the condition of this particular area is extremely unstable. It has also been discovered that the slip behaviors of this area are highly related to rainfall. The monitoring data since 2000 indicate that the cumulative rainfall needed for accelerated side slope slips is 80–270 mm, which equals a cumulative rainfall of 3–5 days.

Xing-Fu Ye [42] discussed the influence of rainfall seepage recharge on slope land collapse through the concept of watershed-based water balance using the following steps: (1) Estimating the base flow volume through the river flow volume qualification line method and the base flow estimation mode, and related to base flow as the seepage recharge of underground water; (2) taking into consideration the correlation system of rainfall, seepage, run-off, evapotranspiration, and underground water recharge through the unsaturated soil water balance method. These two modes were both shown to have similar results. The next step was to carry out a sensitivity analysis through STEDWIN (STEDwin is the smart editor to simplify working with Purdue or PennDOT\* STABL programs) and to discuss its influence on the side slope stability. This revealed that the internal friction angle variability affects the side slope stability most significantly, followed by sloped variability and the rise of the underground water level, while the coherence and unit weight come next. From the figure of the area's rainfall and safety factor relationship, it can be seen that as the steepness of the slope reaches 35 degrees and the rainfall reaches 400–500 mm, according to the analysis of Janbu's simplified method, the safety factor reaches 1.0, forming an imminent risk of collapse. Regarding the discussion of practical examples, a correlation was found between the locations in the Cingshuei River basin that have suffered from collapse and the estimated rainfall recharge factors.

Yi-Feng Qiu [43] discussed the effects of the relevant factors that have affected the mechanism of side slope slips through on-site measurements and theoretical analyses. Long-term monitoring work regarding information such as land surface and stratum heterotaxies, rainfall, the underground water level, and the underground water flow volume of the surrounding areas of Wuwanzi along the Taiwan 18th Line Highway (29 K + 900 ~ 31 K + 500) has been carried out since the beginning of 2000. Serious failures took place on 26 June 2003, leading to the depletion of 150 m of the road base of the 31 K + 340 section and a 1.5-month highway disruption. The causes of this event's side slope failures were analyzed, discussed, and compared with the on-site monitoring results. With the theoretical analyses and the estimation of the stability of falling residual slopes, a final possible restoration solution was proposed. Ming-ren Gao [44] studied the slope lands near Wuwanzi along the Taiwan 18th Line Highway. An unexpected slip failure struck the 31 K + 340 section on 26 June 2003, causing disruption of the entire highway. A simulation of dry seasons and rainy seasons based on the underground water level was conducted using SLOPES/W software (SLOPE/W is the leading slope stability software for soil and rock slopes, GeoStudio company, head office in Calgary, Alberta, Canada). The load input was

increased during rainy seasons to achieve the side slope's failure conditions. The failure mode was then integrated, analyzed, and compared with the observation statistics of the tiltmeter.

The Japanese Society for Landslip Solution Technologies [45] suggested a classification of the displacement speed effect on unstable side slope activities, as shown in Table 2. Such a classification is the displacement speed obtained through dividing the obtained stratum displacement by the observation time. The basis of classification depends on the level of speed and the trend of displacement. The table shows that unstable side slope activities can be roughly classified into four categories. The first three categories are side slopes with (1) emergency movements, (2) definite movements, and (3) semi-definite movements, which are moving slopes with confirmed problems that require immediate effective solutions. The fourth type of movement suffers from potential hazards in terms of the stability of side slopes if there is a certain trend in the displacement direction, because the displacement speed of this type of potentially moving side slope is extremely low—i.e., 0.5–2.0 mm every month, thus producing a displacement quantity of 6–24 mm every year. Further measurements are needed to confirm the stability of side slopes with displacement speeds of over 0.5 mm/month but without certain cumulative displacement directions (forward displacements and backward displacements have both been observed). The stability of such side slopes can be affected by measurement errors of devices or by the non-compact surrounding stuffing of inclination observation tubes when installed. The relevant precaution solutions and management principles provided by the Japanese Society for Expressways [46] are displayed in Table 3.


**Table 2.** Table of displacement speed and side slope stability judgment and suggestions [45].

**Table 3.** Monitoring management benchmarks of stratum slips [46].


To date, no domestic regulations have been set for relevant management values; in addition, no long-term monitoring results have been gained. Consequently, referring to Japan's recommended values and the device-measurable precision of side slope stability, judgement is inevitable, since Japan already has a certain amount of related experience. The retractable meters' action values are 5–10 times that of the vigilance values. The tentative monitoring management values are shown in Table 4 based on the premise that all of the monitoring devices' action values are five times that of their vigilance values. In the future, the data obtained by the monitoring of the devices will be analyzed, interpreted, and judged, and interactive comparisons will be made between different monitoring devices.


**Table 4.** Management values of the monitoring devices [30].

In this study, by using these sensors, we integrated software and hardware to monitor landslides in Pingding. We aimed to review the benchmarking management values and to further investigate the management values. In addition, improvement of the real-time monitoring system was also expected after amending the management values based on different viewpoints to better determine the effects and efficiency of monitoring and vigilance.

In Pingding, Taiwan, particularly in summer (i.e., the rainy season), rainfall is the main contributing factor to the occurrence of landslides, rockfalls, and mudflows. Rainfall information is mainly used to monitor whether there is rain alert on the spot. Station rain gauges, automatically adopt the rain station's operation method, and provide a self-sufficient method. Such stations play a role in recording data, storage, transmission, and display systems for monitoring landslides. The results of this research not only monitor the landslide in real-time, but also manage the space transformation conditions when disasters occur.

The advantages of such a station are: (1) The station is equipped with advanced sensors and technologies; (2) provision of early warning information in real-time for residents who live in the area, as well as for managers who are in charge of this issue, by sending an SMS message to a mobile device; (3) a successful landslide monitoring station can support the monitoring network in Taiwan.

The feasibility of monitoring collapses and displacements through GPS was proven by this study. In addition, with the precision of different monitoring devices, the judgment of the management benchmark values of collapse monitoring established by domestic and overseas documentations was found to be reasonable, according to the monitoring results. Further data concerning typhoons and torrential rainfall will be collected in the future to confirm or amend the management benchmark values.

Finally, this study plays a vital role in the early warning system of landslides in Pingding, as highly effective monitoring can provide information on the current status that is not only applicable for weather forecasts, but also for disaster mitigation policy.

#### **5. Conclusions**

This study established the Pingding Monitoring System in the surrounding area of the Pingding Collapse Zone. On-site observation equipment was installed, and topographic and hydrological information was also collected as a reference for local disaster prevention and response.

Centered on emergency response, this study mainly focused on the displacement measurement of land surfaces. In light of the analysis results of side slope stability, the main cause of collapse and slips appears to be the water content of the soil. The previous collapse disaster confirmation suggests that inbuilt meters that can detect the soil moisture in mudstones should be used to serve as a reference for the necessity of advanced vigilance or as an alert when the water content in the soil is on the rise.

Regarding modern disaster prevention management thinking, the strategy for tackling disasters is a process that involves the application of the following four steps: (1) Readiness—establishing emergency response measures and a management system to enable a quick response once a disaster strikes; (2) response—taking immediate actions before, during, and after disasters to reduce the number of casualties and the property loss and to accelerate recovery; (3) restoration, including the restoration of basic vital resource supply systems within short periods of time and the long-term responsibility of restoring people's normal lives; (4) disaster reduction—advocating for certain policies and applying particular measures to alleviate the impacts brought by future disasters. The above thinking suggests that solving disaster management problems through single modes is not sufficient; instead, in addition to traditional engineering treatments, it is essential to consider multilateral aspects such as multilateral and cross-field discussions involving the following issues: Climatology, economics, engineering, geography, geology, law, meteorology, planning, psychology, public policy, and sociology.

In the future, we would like to apply this research to further applications of the system and its subsequent monitoring in order to establish an early-warning system (social-economic goal) and to conduct post-event analysis when a landslide has actually occurred (scientific goal).

**Author Contributions:** Conceptualization, T.Y.C. and T.V.H.; data curation, Y.M.F.; formal analysis, Y.M.F. and T.V.H.; funding acquisition, T.Y.C. and Y.M.F.; investigation, Y.M.F. and Q.T.B.; methodology, Y.M.F. and Q.H.N.; project administration, T.Y.C. and Y.M.F.; software, Y.M.F. and Q.H.N.; supervision, T.Y.C. and D.B.N.; validation, T.Y.C. and D.B.N.; visualization, Y.M.F. and T.V.H.; writing—Original draft, T.-Y.C. and T.V.H.; writing—Review and editing, T.Y.C. and T.V.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This article is the result of a state-level project financed by the Soil and Water Conservation Bureau and Geographic Information Systems Research Center, Feng Chia University, Taiwan (grant number MOST20051110).

**Acknowledgments:** The authors thanks the Soil and Water Conservation Bureau and Geographic Information Systems research Center, Feng Chia University, Taiwan for their support.

**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*

### **Separating Landslide Source and Runout Signatures with Topographic Attributes and Data Mining to Increase the Quality of Landslide Inventory**

#### **Jhe-Syuan Lai**

Department of Civil Engineering, Feng Chia University, Taichung 40724, Taiwan; jslai@fcu.edu.tw; Tel.: +886-4-2451725 (ext. 3118)

Received: 4 September 2020; Accepted: 21 September 2020; Published: 23 September 2020

**Abstract:** Landslide sources and runout features of typical natural terrain landslides can be observed from a geotechnical perspective. Landslide sources are the major area of occurrences, whereas runout signatures reveal the subsequent phenomena caused by unstable gravity. Remotely sensed landslide detection generally includes runout areas, unless these results have been excluded manually through detailed comparison with stereo aerial photos and other auxiliary data. Areas detected using remotely sensed landslide detection can be referred to as "landslide-affected" areas. The runout areas should be separated from landslide-affected areas when upgrading landslide detections into a landslide inventory to avoid unreliable results caused by impure samples. A supervised data mining procedure was developed to separate landslide sources and runout areas based on four topographic attributes derived from a 10–m digital elevation model with a random forest algorithm and cost-sensitive analysis. This approach was compared with commonly used methods, namely support vector machine (SVM) and logistic regression (LR). The Typhoon Morakot event in the Laonong River watershed, southern Taiwan, was modeled. The developed models constructed using the limited training data sets could separate landslide source and runout signatures verified using the polygon and area constraint-based datasets. Furthermore, the performance of developed models outperformed SVM and LR algorithms, achieving over 80% overall accuracy, area under the curve of the receiver operating characteristic, user's accuracy, and producer's accuracy in most cases. The agreement of quantitative evaluations between the area sizes of inventory polygons for training and the predicted targets was also observed when applying the supervised modeling strategy.

**Keywords:** data mining; landslide detection; landslide inventory; Typhoon Morakot

#### **1. Introduction**

Natural hazards occur frequently in Taiwan. Among them, landslides can be triggered by heavy rainfall events, especially in mountainous areas. The precipitous terrain, complicated geology, and dense population may increase vulnerability, causing a serious loss of life, property damage, and economic loss. Furthermore, these events can affect water resource supply and cause other livelihood problems. Therefore, landslide analysis and assessment have become critical in natural hazard and disaster mitigation, prevention, and reconstruction in Taiwan. A growing number of studies have investigated landslide-related topics on different scales, such as slope stability analysis for a specific slope site [1,2]; regional landslide detection and mapping [3–7]; characterization of the relationship between landslides and environmental or triggering factors [8–13]; susceptibility, hazard, risk assessment, and management [14–19]; and modeling or estimating physical, environmental, and rainfall-based parameters [20–24]. Landslide risk assessment and management are crucial, systematic, and extensive frameworks in the related works. Dai et al. [14] divided this topic into

several parts. In this framework, generating landslide inventory (inventory map and database are used synonymously in this study) is essential to connect the following assignments.

A basic landslide inventory should record the landslide's location, date (event-based) or period (multitemporal), and types of movement [25]. The definition, assumptions, requirements, production processes, and statistical properties of landslide inventories have been discussed in detail by Guzzetti et al. [25], Harp et al. [26], Malamud et al. [27], and Shao et al. [28]. Furthermore, Galli et al. [29] and Mondini et al. [11] have assessed the quality and completeness of landslide inventory maps by comparing two maps of the same area. Numerous studies have indicated that the continuing improvements in remote sensing and geographic information systems (GISs) have led to cooperation with data mining and machine learning algorithms to produce a regional landslide inventory. In particular, integrating GIS-based models with geo-spatial data [18,30] and generating event-based landslide inventory [11,25,31] have garnered interest.

Three common features of typical natural terrain landslides have been observed from a geotechnical perspective. These features were the landside source, trail, and deposition fan [32]. The surface of the rupture, comprising the main scarp and the scarp floor, is defined as the source area. A landslide trail may also occur predominantly as a result of the transport of the landslide mass. The majority of the landslide mass is deposited (i.e., deposition fan). The term runout is generally used to indicate the landslide trail and deposition fan [33]. Landslide detection performed using automatic and semiautomatic methods with remotely sensed images usually includes runout areas unless these results have been excluded manually through a detailed comparison with stereo aerial photos and other auxiliary data. The areas detected using remotely sensed landslide detection can be referred to as "landslide affected" areas. Mixing landslide source areas with runout regions may reduce the reliability of a landslide analysis [34], such as susceptibility and hazard assessments. These runout areas should be separate from landslide source areas because they have different mechanisms. More precisely, producing landslide inventory requires further processing after the remote sensing of landslide detections.

Data mining and machine learning-based algorithms have been used increasingly for landslide modeling, especially in landslide susceptibility assessments, such as decision trees, deep learning systems, evolution-based algorithms, fuzzy theory, neural networks, random forests (RF), and support vector machines (SVM) [35–40]. Related studies have also considered various composites and compared the effectiveness of strategies based on established methods in achieving a specific purpose [41–43]. However, few studies have employed these approaches for detecting landslides and producing inventories that consider the separation of runout areas from landslide affected regions. The RF method [44] consists of data mining and machine learning algorithms that have displayed excellent performance in the analysis of numerous remote sensing and landslide topics [45,46]. The RF algorithm is based on the ensembles of various decision tree results and exhibit desirable properties, such as high accuracy, robustness against overfitting the training data, and integrated measures of variable importance [47]. Furthermore, Lai and Tsai [48] and Lai et al. [34] have demonstrated that combining the RF algorithm with a cost-sensitive analysis [49] in landslide susceptibility assessments can reduce extreme omission (missing) or commission (false alarm) predictions because it adjusts the decision boundary.

The main purpose of this study was to explore the feasibility of separating landslide source and runout areas based on landslide affected extents extracted from remotely sensed images in order to upgrade landslide detections into a landslide inventory. The significance of topographic data in relation to the flow velocity, geomorphology, runoff rate, soil water content, and differences between landslide source and runout areas was demonstrated [34,50]. The term "signature" used in this study represents the patterns of landslide source and runout areas in a feature space. Feature space means that the dimension is the used factors with the samples based on a scatter plot form. Therefore, the developed RF-based data mining models and topographic attributes, such as aspect, curvature, elevation, and slope, constructed for the Typhoon Morakot event in 2009 were combined to compare the results with those

obtained using SVM and logistic regression (LR) algorithms. SVM and LR algorithms are commonly applied in the related domains. The performance of cost-sensitive analysis was also evaluated to improve the constructed models by adjusting the decision boundary.

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

#### *2.1. Study Site and Data*

The study area was located in the Laonong River watershed in southern Taiwan and covered approximately 117 km2, as illustrated in Figure 1a. The elevation in the study site ranged from 258 to 1666 m above sea level, measured using the digital elevation model (DEM) modified by Chiang et al. [51]. The average elevation was 716 m. The slope range, average slope, and standard deviation were 0–71.2◦, 25.84◦, and 11.98◦, respectively, which indicated that the terrain of the Laonong River watershed is steep. Three geological formations and four soil types were identified in the study area. The geological formations were Lushan, Sanhsia, and Toukoshan. The four soil types were alluvium, colluviums, lithosol, and loam. Lai et al. [34] reported detailed information regarding geological formations and soil types within the study site. The Laonong River watershed is located in the tropical monsoon region, and the annual precipitation is approximately 3400 mm. Therefore, this area is frequently struck by typhoons. For example, Typhoon Morakot caused extensive rainfall in 2009, which resulted in numerous landslides and debris flow in southern Taiwan. An official report [52] indicated that 769 casualties or missing people were directly or indirectly caused by these landslides. Furthermore, Typhoon Morakot caused losses of approximately \$526 million because of the damage to agriculture, forestry, and fishery. In particular, a riverside village called Xiaolin (sometimes spelled Shiaolin, Hsiaolin) was destroyed by the landslides and debris flows from a devastating landslide nearby, which caused approximately 500 fatalities.

**Figure 1.** Location of the study site (**a**) and landslide inventory map (**b**).

A landslide inventory map of Typhoon Morakot revealed deep-seated landslides, as illustrated in Figure 1b. The landslide source, runout, and channel classes were interpreted manually [34] based on stereo aerial photos and auxiliary data. This study assumed that the used landslide inventory comprised the results of the detection of affected landslide areas extracted from the remotely sensed images. The landslide inventory was also used for quantitative verifications. The 10–m DEM was then analyzed to derive topographic attributes, including the aspect, curvature, elevation, and slope, as illustrated in Figure 2. Furthermore, the landslide inventory map was converted into a grid format of 10 m2 to match the topographic factors for constructing data mining-based models.

**Figure 2.** Topographic layers used in this study. (**a**) Aspect, (**b**) curvature, (**c**) elevation, and (**d**) slope.

#### *2.2. Procedure and Data Mining Algortithms*

Figure 3 illustrates the conceptual procedure implemented in this study. The topographic attributes were derived from the 10–m DEM, which connected the samples extracted from the landslide inventory for modeling (Sections 2.2.1 and 2.2.2) and verification (Section 2.3). Two experiments were designed for the demonstration of separating landslide source and runout signatures. First, the limited training datasets were selected (the five largest area sizes of landslide source and runout classes in the landslide inventory polygons) for polygon by polygon modeling. Second, each model was verified using other four polygon-based samples. Third, these models were used to predict the polygon samples with an area of <sup>≥</sup>100,000 and 50,000–100,000 m<sup>2</sup> to identify the optimal model. Previous results are presented in Sections 3.2 and 3.3. Detailed information regarding the samples extracted from the landslide inventory is presented in Table 1. Finally, the spatial patterns of landslide source and runout areas were produced instance by instance based on the output labels of the optimal model, as demonstrated in Section 3.4.

**Figure 3.** The conceptual procedure used in this study to model and verify the separation between landslide source and runout signatures.


**Table 1.** The training and verification samples extracted from the landslide inventory polygons.

\*Areas and number of samples in the table include the training polygon samples, but these samples were eliminated during the verified process.

#### 2.2.1. Developed Models

The random forests (RF) algorithm [44] is a data mining method used for constructing a classification-based model that employs a supervised strategy. RF is an extension of the tree-based algorithms. The RF employs the information gain (IG) measure or Gini index to determine the degree of impurity of factors or variables, applying the bootstrap and vote operators to improve the results derived from the tree-based rules. The bootstrap method "randomly" selects subtraining data to generate several trees (termed "forests") to avoid the overfitting results. The vote component is used to determine the optimal results according to the developed trees, thus improving the RF classifier. A topographic factor with a larger IG or Gini value had to be selected because higher values indicate a higher priority for constructing a tree node, which should be ignored in the next computation. A sequence of tree-based rules was constructed after several iterations. The rules were then used to classify other instances. Nominal (or discrete) and numeric (or continuous) formats are commonly used in the field. Entropy theory was applied to calculate the IG in the nominal case, as displayed in Equations (1)–(3). In these equations, E(A) represents the entropy of all training data; S and R represent the numbers of landslide sources and runout samples, respectively; E'(a) and v are the entropy and number of subsets of a specific topographic factor, respectively; E(aj) is the entropy of the subset in a specific topographic factor, calculated using Equation (1); IG(a) represents the IG of a specific topographic factor. The Gini index was applied to determine the IG in the numeric cases, as displayed in Equations (4) and (5), where C represents a segmented point for a specific topographic factor used to divide continuous data into two parts, and N1 and N2 represent the numbers of a ≤ C and a > C, respectively. The steps were detailed by Guo et al. [53].

$$\mathrm{E(A)}\, =\, -\frac{\mathrm{S}}{\mathrm{S}+\mathrm{R}}\log\_{2}{\frac{\mathrm{S}}{\mathrm{S}+\mathrm{R}}} - \frac{\mathrm{R}}{\mathrm{S}+\mathrm{R}}\log\_{2}{\frac{\mathrm{R}}{\mathrm{S}+\mathrm{R}}}\tag{1}$$

$$\mathbf{E}'(\mathbf{a}) = \sum\_{i=1}^{\mathbf{v}} \frac{\mathbf{S}\_{\mathbf{i}} + \mathbf{R}\_{\mathbf{i}}}{\mathbf{S} + \mathbf{R}} \mathbf{E}(\mathbf{a}\_{\mathbf{i}}) \tag{2}$$

$$\text{IG(a)} = \text{E(A)} - \text{E'(a)}\tag{3}$$

$$\text{Gini}(\mathbf{a} \le \mathbf{C} \text{ or } \mathbf{a} > \mathbf{C}) \; = \; 1 - \sum\_{i=1}^{\mathbf{m}} \frac{\mathbf{n}\_i}{\mathbf{N}} \tag{4}$$

$$\text{IG}(\mathbf{a}, \mathbf{C}) = \frac{\text{N}\_1}{\text{N}} \text{Gini}(\mathbf{a} \le \mathbf{C}) + \frac{\text{N}\_2}{\text{N}} \text{Gini}(\mathbf{a} > \mathbf{C}) \tag{5}$$

The strategy of adjusting the decision boundary used in the study is termed as the cost-sensitive analysis. The cost-sensitive analysis is a postclassification method based on the cost matrix that reclassifies the instances and balances the accuracies of certain classes when missing or false alarm errors are unreliable [49,54,55]. The dimension of the cost matrix is equal to that of the confusion matrix. The confusion matrix used in this study is presented in Table 2. True negative (TN), false negative (FN), false positive (FP), and true positive (TP) were used to represent an agreement between the classified and reference labels in counts tabulated in a confusion matrix. The cost in this study indicated a weighting without the unit. The diagonal costs in the table represent correct results for the TN and TP, and the remaining costs indicate misclassification costs between the FN and FP. The diagonal costs and other elements are usually set to 0 and 1, indicating unadjusted and adjusted conditions of the decision boundary. A large amount of missing or false alarm errors are severe FP or FN errors, respectively. Adjusting the decision boundary by increasing the cost (weighting) of the FP or FN leads to the inclusion of more samples, thereby balancing the classification result of a certain class. Equation (6) displays the optimal prediction (R) of sample x in class i, where P(j|x) is the likelihood of estimating a classification of a sample into all classes (j), and C represents the costs. In binary cases, the optimal prediction is the landslide source label (class 1 or positive) if the expected cost of this prediction is less than or equal to the expected cost of predicting the runout label (class 0 or negative), as displayed in Equation (7). Given p = P(1|x) and CTN = CTP = 0, Equation (7) can be simplified, as displayed in Equation (8). An adjusted threshold (p\*) of the decision boundary, presented in Equation (10), can be

derived from Equation (9) to classify a sample x as the landslide source label when P(1|x) is larger than or equal to the threshold [56]:

$$\mathcal{R}(\mathbf{i}|\mathbf{x}) = \sum\_{\mathbf{j}} \mathcal{P}(\mathbf{j}|\mathbf{x}) \mathbf{C}\_{\mathbf{j}\mathbf{i}} \tag{6}$$

$$\mathbf{P}(0|\mathbf{x})\mathbf{C}\_{\rm FP} + \mathbf{P}(1|\mathbf{x})\mathbf{C}\_{\rm TP} \le \mathbf{P}(0|\mathbf{x})\mathbf{C}\_{\rm TN} + \mathbf{P}(1|\mathbf{x})\mathbf{C}\_{\rm FN} \tag{7}$$

$$(1 - \mathbf{p})\mathbf{C}\_{\text{FP}} \le \mathbf{p} \, \mathbf{C}\_{\text{FN}} \tag{8}$$

$$(1 - \mathbf{p}^\*)\mathbf{C}\_{\text{FP}} \le \mathbf{p}^\*\mathbf{C}\_{\text{FN}} \tag{9}$$

$$\mathbf{p}^\* = \frac{\mathbf{C\_{FP}}}{\mathbf{C\_{FP}} + \mathbf{C\_{FN}}} \tag{10}$$


#### **Table 2.** The confusion matrix used in this study.

#### 2.2.2. Algorithms for Comparison

Two commonly used data mining methods were used for comparison with the developed models. The support vector machines (SVM) method has been wieldy applied in remote sensing studies [57] to classify land cover or use targets, extract biophysical parameters in different spatial resolutions, and perform landslide analyses [40,41]. The core role of the SVM classifier is to determine the largest margin of the classifier based on factor transformations with linear and nonlinear functions. These transformations are known as the kernel function. Equation (11) describes a concept of the SVM-based classification with a "sign" operator for the binary classification, where l, y, and k(xi,xj) indicate the number of used factors (support vectors), the outputted label, and a kernel function, respectively, and a and b represent the coefficients of the margin. The kernel functions used in this study were linear, polynomial, and radial basis functions.

$$\mathbf{f}(\mathbf{x}) = \text{sign}\left[\sum\_{\mathbf{i}=1}^{l} \mathbf{a}\_{\mathbf{i}} \mathbf{y}\_{\mathbf{i}} \mathbf{k} \begin{pmatrix} \mathbf{x}\_{\mathbf{i}} & \mathbf{x}\_{\mathbf{j}} \end{pmatrix} + \mathbf{b} \right] \tag{11}$$

Logistic regression (LR) is a typical and traditional statistical method used to model landslide events [20,58]. The computation describes the relationship between the dependent and independent variables as defined in Equation (12), where p refers to the likelihood or probability, a is a constant, b is the regression coefficient, and x represents independent variables or used factors:

$$\mathbf{p} = \frac{\exp(\mathbf{a} + \mathbf{b}\_1 \mathbf{x}\_1 + \mathbf{b}\_2 \mathbf{x}\_2 + \dots + \mathbf{b}\_n \mathbf{x}\_n)}{1 + \exp(\mathbf{a} + \mathbf{b}\_1 \mathbf{x}\_1 + \mathbf{b}\_2 \mathbf{x}\_2 + \dots + \mathbf{b}\_n \mathbf{x}\_n)} \tag{12}$$

#### *2.3. Accuracy Assessment*

The quantitative results of comparing the output and ground truth labels, identified based on the threshold, were derived from the confusion matrix to verify the constructed models. An instance was categorized into the landslide (positive) class in most cases when the likelihood of occurrence was ≥0.5; otherwise, the instance was assigned the nonoccurrence (negative or non-landslide) label. Similar to previous classifications, the positive and negative labels in this study were changed from landslide and non-landslide to landslide source and runout classes, respectively. The threshold of 0.5 for classification could also be adjusted by the cost-sensitive analysis (p\*).

Four commonly used quantitative indices derived from a confusion matrix, namely overall accuracy (OA), user's accuracy (UA), producer's accuracy (PA), and area under the receiver operating characteristic (ROC) curve (AUC), were used to quantitatively evaluate the developed models. OA represents the ratio of samples correctly classified as TN and TP, as illustrated in Equation (13). UA and PA reflect the errors for each class (landslide source and runout) as presented in Equations (14)–(17), where R and S represent the runout and landslide source classes. The false alarm (commission errors) and missing (omission errors) assignments can be calculated as 1–UA and 1–PA, respectively. The AUC is commonly used to assess data mining-based models. The binary classification threshold is used for reclassifying instances, calculating the rates of FPs and TPs to draw the ROC curve. In general, the AUC ranges from 0.5–1, and results larger than 0.7 are acceptable, and larger than 0.8 are excellent [59].

$$\text{OA} = \frac{\text{TN} + \text{TP}}{\text{TN} + \text{FN} + \text{FP} + \text{TP}} \tag{13}$$

$$\text{UA(R)} = \frac{\text{TN}}{\text{TN} + \text{FN}} \tag{14}$$

$$\text{UA(S)} = \frac{\text{TP}}{\text{FP} + \text{TP}} \tag{15}$$

$$\text{PA}(\text{R}) = \frac{\text{TN}}{\text{TN} + \text{FP}}.\tag{16}$$

$$\text{PA(S)} = \frac{\text{TP}}{\text{FN} + \text{TP}} \tag{17}$$

#### **3. Results**

The data mining modeling and accuracy assessments used in this study were developed using the WEKA environment (http://www.cs.waikato.ac.nz/mL/weka/), which is a free and open-source platform. A four-stage step was designed to explore the separability between landslide sources and runout signatures. The first stage, as reported in Section 3.1, compared the topographic attributes between the landslide source and runout samples extracted from different sizes of landslide inventory polygons that had an area of <sup>≥</sup>100,000 or 50,000–100,000 m2. In the second stage, the RF and cost-sensitive analysis-based data mining algorithms were used to construct the models polygon by polygon, with the limited training datasets (the five largest area sizes in the landslide inventory polygons were considered), as reported in Section 3.2. The major parameter in the RF-based computation was the number of trees (Ntree). Du et al. [60] determined that 10–200 Ntree did not have any effect on the results, but an increase in Ntree increased the computation loading. Therefore, 100 trees were adopted in this study, as suggested by WEKA, to develop the RF-based models. In the third stage (Section 3.3), the constructed models were applied to predict other landslide source and runout samples extracted from the polygons with an area of <sup>≥</sup>100,000 and 50,000–100,000 m2, respectively. The developed models were compared with the commonly used methods, namely SVM and LR. The fourth stage consisted of visualizing the spatial distributions of predicting the landslide source and runout areas to assess the performance of the constructed model, as reported in Section 3.4.

#### *3.1. Signatures of Topographic Attributes*

The constraints on the area size of the inventory polygons were compared in the pairs of the topographic attributes to preliminarily examine differences between landslide source and runout signatures. The samples were extracted based on area sizes of <sup>≥</sup>100,000 and 50,000–100,000 m2. Furthermore, the value domain of the aspect, curvature, elevation, and slope layers was normalized into the range of 0–1 for the visualization of data distribution and separation. The patterns obtained by comparing the topographic datasets and area constraints are illustrated in Figures 4 and 5, respectively. A larger area size enabled a larger separation of landslide source and runout features. The disagreement was not observed in the smaller landslide inventory polygons. This trend is similar to findings reported by Lai et al. [34]. This trend may reflect that the complete mechanisms of landslide source and runout areas appear in the larger landslide sizes. Moreover, the elevation is a critical factor in this study. These data distributions reveal the feasibility of distinguishing the landslide source and runout signatures when the sizes of landslide-affected areas were sufficient to reveal their mechanisms.

**Figure 4.** Patterns for comparing aspect, curvature, elevation, and slope layers with the samples extracted from the inventory polygon with an area of <sup>≥</sup>100,000 m2, where red and blue dots represent the landslide source and runout classes, respectively. (**a**): Aspect vs. curvature.; (**b**): curvature vs. elevation; (**c**): aspect vs. elevation; (**d**): curvature vs. slope; (**e**): aspect vs. slope; (**f**): elevation vs. slope.

#### *3.2. Construction of Random Forests Based Data Mining Models*

The RF-based data mining models were constructed using the top five polygons of area sizes in the landslide inventory polygons to further explore the feasibility of separating the landslide source and runout areas. Each polygon was transformed into a grid format of 10 m2. These samples could be used to extract the corresponding topographic attributes and produce the training and verification datasets. The quantitative evaluations of OA, AUC, PA, and UA are presented in Figure 6. The accuracies exceeded 80% in most cases. Model No. 3 provided the most accurate predictions. However, the results of predicting (a) the No. 3 and No. 5 polygon samples by using the No. 1 and No. 2 models, (b) the No. 5 polygon samples by using the No. 3 and No. 4 models, and (c) the No. 1 and No. 2 polygon samples by using the No. 5 model revealed a disagreement between the training and verification datasets. Possible reasons for this disagreement include (a) imbalanced sample ratios between landslide source and runout classes, (b) dissimilar data distributions between different inventory polygons, (c) lack of other critical factors, and (d) limited training data.

**Figure 5.** Patterns for comparing aspect, curvature, elevation, and slope layers with the samples extracted from the inventory polygon with an area of 50,000–100,000 m2, where red and blue dots represent the landslide source and runout classes, respectively. (**a**): Aspect vs. curvature.; (**b**): curvature vs. elevation; (**c**): aspect vs. elevation; (**d**): curvature vs. slope; (**e**): aspect vs. slope; (**f**): elevation vs. slope.

**Figure 6.** Accuracies of the random forest (RF)-based models using the top five polygon samples to predict each other, where OA, AUC, PA, UA, S, and R represent overall accuracy, the area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively. (**a**) No. 1, (**b**) No. 2, (**c**) No. 3, (**d**) No. 4, and (**e**) No. 5 polygon samples for modeling.

To address this problem, cost-sensitive analysis procedures, similar to the procedure presented by Lai and Tsai [48], were employed to adjust the decision boundary during RF-based modeling. The results of the cost-sensitive analysis for different costs are illustrated in Figure 7. The developed model without the cost-sensitive analysis (i.e., the cost is equal to 1) had a high missing error for detecting landslide source regions, as displayed in Figure 7a–f, and runout areas, as displayed in Figure 7g,h, because of the assignment of higher FN and FP, respectively. The variations and improvements of accuracies could be observed after increasing the costs of FN and FP for the cases presented in Figures 7a–f and 7g,h, respectively. Based on Figures 6 and 7, Table 3 lists the representative models applying the RF algorithm and cost-sensitive analysis with different inventory polygon samples for further verification, as presented in Section 3.3.

**Figure 7.** Accuracies (*Y*-axis) of the RF-based models with different costs (*X*-axis) and polygon samples for predictions, where OA, AUC, PA, UA, S, and R represent the overall accuracy, area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively. (**a**) No. 1 model predicting No. 3 polygon samples, (**b**) No. 1 model predicting No. 5 polygon samples, (**c**) No. 2 model predicting No. 3 polygon samples, (**d**) No. 2 model predicting No. 5 polygon samples, (**e**) No. 3 model predicting No. 5 polygon samples, (**f**) No. 4 model predicting No. 5 polygon samples, (**g**) No. 5 model predicting No. 1 polygon samples, and (**h**) No. 5 model predicting No. 2 polygon samples.


**Table 3.** The representative models using the random forests (RF) algorithm with the costs for verifications based on pervious comparisons. For example, RF\_3000 indicates that the RF algorithm with a cost of 3000 was used.

#### *3.3. Model Performances and Comparisons*

Based on settings presented in Table 3, the constructed models were used to predict the samples extracted from other inventory polygons with an area of <sup>≥</sup>100,000 or 50,000–100,000 m2, namely the larger and smaller inventory polygons, presented in Sections 3.3.1 and 3.3.2, respectively. Furthermore, the performances between the developed models and the approaches of SVM and LR were compared.

#### 3.3.1. The Case of Larger Inventory Polygons

The prediction results of No. 1 models obtained using the RF and cost-sensitive analysis are presented in Figure 8a. The examination of this figure reveals that a cost of 5000 provided more favorable results, with the accuracy reaching over 80% in most indices. However, the evaluations of SVM-based approaches in consideration of linear, the first- to third-order polynomial, and radial basis functions are displayed in Figure 8b. Overall, the second-order polynomial function outperformed others. Compared with previous LR results, Figure 8c indicates that the developed model in Section 3.2 obtained the highest accuracies, although the UA of the landslide source was slightly lower.

**Figure 8.** Performances (*Y*-axis) of No. 1 model, where OA, AUC, PA, UA, S, and R represent the overall accuracy, area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively. (**a**) Using the algorithm with the costs; (**b**) using the support vector machine (SVM) with linear, the first- to third-order polynomial, and radial basis (RBF) functions; (**c**) comparing the optimal results of (**a**) and (**b**) by using logistic (Logis) regression.

Similar procedures to those used to verify the No. 2 model were applied to assess the algorithms of RF with the cost-sensitive analysis, SVM, and LR. Figure 9a further compares the RF with costs, as displayed in the third column in Table 3. Based on this figure, the RF algorithm appeared to outperform the other models without increasing costs. As illustrated in Figure 9b, the accuracies of all the functions of SVM modeling were similar, with higher missing errors for detecting landslide sources. The RF and the second-order polynomial function of SVM were compared with LR, as illustrated in Figure 9c, which revealed that the AUC of the LR was the highest. The accuracies of the three other approaches were similar.

**Figure 9.** Performance (*Y*-axis) of the No. 2 model, where OA, AUC, PA, UA, S, and R represent the overall accuracy, area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively. (**a**) Using the RF algorithm with the costs; (**b**) using SVM with linear, the firstto third-order polynomial, and radial basis (RBF) functions; (**c**) comparing the optimal results of (**a**) and (**b**) by using logistic (Logis) regression.

For the No. 3 model, Table 4 presents the comparison of the RF algorithm and cost-sensitive analysis. The results indicate that the original RF algorithms provided higher accuracy, whereas the RF method with a cost of 1000 had lower PA and UA for the runout and landslide source classes, respectively. In the SVM case, Figure 10a reveals that the RBF obtained the least accurate results, especially PA and UA for the runout and landslide source signatures. The performances of linear and the first-order polynomial functions were similar. A further comparison of these models, illustrated in Figure 10b, revealed that the results produced by the RF method were more favorable in most cases, especially UA for the landslide source class.

**Table 4.** The performances of the No. 3 model by using the RF method with the cost of 1000, where OA, AUC, PA, UA, S, and R represent overall accuracy, the area under ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively.


**Figure 10.** Performances (*Y*-axis) of the No. 3 model, where OA, AUC, PA, UA, S, and R represent the overall accuracy, area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively. (**a**) Using SVM with linear, the first- to third-order polynomial, and radial basis (RBF) functions; (**b**) comparing the optimal results of Table 4 and (**a**) by using logistic (Logis) regression.

The No. 4 models were also examined as illustrated in Figure 11. The RF method without costs also provided the most accurate results, as presented in Figure 11a. Similarly, the performances of the linear and first-order polynomial functions outperformed the RBF method in SVM-based computation, as illustrated in Figure 11b. The comparison based on the aforementioned methods is illustrated in Figure 11c. The findings indicated that the AUC result of SVM was slightly higher than the other indices of the RF, SVM, and LR, which were similar.

**Figure 11.** Performance (*Y*-axis) of the No. 4 model, where OA, AUC, PA, UA, S, and R represent overall accuracy, area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively. (**a**) Using the RF algorithm with the costs; (**b**) using SVM with linear, the firstto third-order polynomial, and radial basis (RBF) functions; (**c**) comparing the optimal results of (**a**) and (**b**) with logistic (Logis) regression.

In the No. 5 model, the results produced by the RF algorithm with a cost of 2000 displayed high accuracy, as illustrated in Table 5. For SVM-based modeling, Figure 12a indicated that the RBF method had unbalanced predictions with a significant missing error for the landslide source and the highest PA for the runout. However, the results of the linear and polynomial functions were similar. The second-order polynomial SVM and LRs, illustrated in Figure 12b, displayed accuracies similar to the RF\_2000.

**Table 5.** Performance of the No. 5 model using the RF method with a cost of 2000, where OA, AUC, PA, UA, S, and R represent the overall accuracy, area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively.

**Figure 12.** Performances (Y- axis) of the No. 5 model, where OA, AUC, PA, UA, S, and R represent the overall accuracy, area under ROC curve, producer's accuracy, user's accuracy, landslide source, and runout-out, respectively. (**a**) Using SVM with linear, the first- to third-order polynomial, and radial basis (RBF) functions, (**b**) comparing the optimal best results of Table 5 and (**a**) with logistic (Logis) regression.

Significant improvements were observed in the No. 1 and No. 3 models. Moreover, the RF algorithm and cost-sensitive analysis provided favorable results in most cases. To identify the optimal model in which the samples extracted from different size areas of the inventory polygons, Figure 13a compares the RF method for the No. 2, 3, and 4 models, with costs of 5000 and 2000 for No. 1 and 5 models, respectively. Figure 13 indicates that the accuracies produced by the No. 1, 3, and 4 models were similar, whereas the No. 2 and 5 models displayed less favorable results for the detection of a landslide source class. Figure 13b indicated that the AUC of the No. 1 and 3 models were the most favorable. These evaluations indicated that the No. 1 and 3 models provided more reliable prediction results.

**Figure 13.** Comparison of the representative results using the RF method for the No. 2, 3, and 4 models, with costs of 5000 and 2000 for No. 1 and 5 models, respectively. (**a**) User's accuracy (UA), producer's accuracy (PA), and overall accuracy (OA); (**b**) area under the ROC curve (AUC) for landslide source (S) and runout (R) detections.

#### 3.3.2. The Case of Smaller Inventory Polygons

The performances of the No. 1 model based on the RF method with a cost of 5000 and the No. 3 model with the RF algorithm were assessed, as presented in Table 6, to predict the samples extracted from the polygon areas equal to or larger than 50,000 and smaller than 100,000 m2. The accuracies were lower than those displayed in Figure 13a, especially for the runout UA and PA of

the No. 1 model. Furthermore, the No. 3 model with the RF method outperformed the No. 1 model, indicating accurate results. The topographic roughness index (TRI) is commonly used for assessing landslide susceptibilities [61]. TRI was further used to improve this case as also listed in Table 6. However, the significant topographic attributes [62–64] and various indices [65] for separating landslide source and runout signatures need to be further explored. These results reveal the uncertainty between the training data and prediction targets. More precisely, the area size of the training inventory polygon affected the capability and applicability of the constructed models. In practice, manually interpreting the larger targets of landslide source and runout from remotely sensed images to produce the training data was easier than selecting smaller objects, which was a clear limitation of using the supervised strategy with four topographic variables in this study. Addressing this problem requires further consideration regarding the tradeoff between the cost of producing training data and the applicability of the developed models. The developed models were suggested for different area sizes of the inventory polygons.

**Table 6.** Performances of the No. 1 model using the RF method with a cost of 5000 and No. 3 model with the RF algorithm, where OA, AUC, PA, UA, S, and R represent the overall accuracy, area under the ROC curve, producer's accuracy, user's accuracy, landslide source, and runout, respectively.


\* With the topographic roughness index.

#### *3.4. Visualization of Landslide Detection*

To visualize the separated results of the case in Section 3.3.1, the predicted patterns of the No. 3 models for distinguishing landslide source and runout targets were generated, as illustrated in Figure 14. The black and red and blue colors in this figure represent the ground truth and prediction results, respectively. As illustrated in Figure 14a, the developed models displayed suitable detection of the landslide sources (high PA(S)) with some false alarm errors (low UA(S)). Figure 14b also illustrated that the model provided reliable results for detecting runout areas with higher UA(R) and PA(R). These patterns accorded with the results of the No. 3 model in Figure 10b.

The detected distributions with the topographic situation in a three-dimensional (3D) space were further examined. Figure 15 displays the patterns of the No. 3 model for predicting the No. 1 and No. 2 landslide source and runout inventory polygons. In this figure, red and yellow represent the correct and incorrect predictions, respectively. Figure 15a,b displays accurate detections for the landslide sources. The results presented in Figure 15c indicate that the lower part of a runout area can be detected, whereas the upper region may be miscategorized as a landslide source. This finding reveals that elevation is a critical factor for the separation of landslide source and runout areas. This finding also accords with the result of the factor analysis presented in Section 3.1. Future studies can determine the elevation-based threshold to improve the capability of the developed models.

**Figure 14.** Prediction results of the No. 3 model using the RF algorithm to separate landslide source and runout signatures with the areas of inventory polygons of <sup>≥</sup>100,000 m2. Cases of (**a**) landslide source and (**b**) runout polygons.

**Figure 15.** Three-dimensional (3D) view of the No. 3 model for predictions of separating landslide source and runout areas. The landslide source areas of the (**a**) No. 1 and (**b**) No. 2 polygons in the landslide inventory (located region A and B in Figure 14), and (**c**) the runout regions of the No. 1 and No. 2 inventory polygons (located region C in Figure 14). Red and yellow represent correct and incorrect results, respectively.

#### **4. Discussion**

The major finding in this study was that detecting landslides is not fully equal to inventory generation. The detected targets require further examination. Lai et al. [34] suggested that the runout area and topological relationship could be considered independent classes in landslide records when producing a landslide inventory. Therefore, this study focused on the separation of landslide source and runout classes caused by different mechanisms. The feasibility of distinguishing between landslide source and runout signatures to improve landslide detections was assessed. The limited training samples extracted from the top five polygons in the landslide inventory map were used to construct supervised data mining models in order to classify other instances. The top five polygon samples were used for training because manually recognizing the larger targets from a regional scale-based remote sensing image is relatively easier than interpreting smaller objects. Based on this scenario, the developed models with RFs and cost analysis and mapping results demonstrated the stability of the proposed procedure. The hierarchical analysis was designed in this study. Section 3.2 explored the feasible models as listed in Table 3. Section 3.3 demonstrated the performances of developed models verified by the samples of larger and smaller inventory polygons, comparing them with the SVM and LR methods as shown in Figures 8–13 and Tables 4–6. The quantitative evaluations of the optimal results as shown in Figure 13 were determined to be higher than 80% in most cases. The results further demonstrated that the developed models outperform the SVM and LR algorithms in most accuracies.

The primary limitation of this study was that the developed models were more favorable for larger area sizes but may be difficult to apply to smaller areas. The capability of the models based on the supervised strategy was related to the characteristics of the training datasets. The tradeoff between the cost of producing the training data and the area sizes of targets for predictions should be considered. Perhaps a multi-scale model trained by different area sizes of the inventory polygons is a possible solution. This method could also be improved by including more topographic attributes, as described by Sallem et al. [61].

Producing a high-quality landslide inventory is time consuming and expensive. Landslide inventories are crucial for landslide analyses, such as for performing tasks in the framework of landslide risk assessment and management [14]. The potential for reducing the degree of manual interpretation is demonstrated in mapping results. Based on the spatial patterns of predictions, the misclassification for detecting runout areas appeared in the top region. Two possible strategies may solve this problem. First, a threshold based on the elevation attribute could be designed. Second, studies could investigate only extracted landslide sources from the affected areas, and the remainder are assigned the runout label. These runout patterns may be adequate to further assess the related parameters, such as runout distance and damage corridor width [14].

#### **5. Conclusions**

A novel data analysis based on the scenario of separating landslide source and runout areas with the topographic attributes was presented in this paper to improve landslide detections. To explore the feasibility of the proposed approach, four hierarchical experiments were designed. First, elevation was a critical variable according to the factor analysis. Second, the candidate models were constructed using the RF and cost-sensitive analysis and the limited training samples extracted from the top five largest areas of the inventory polygons. The polygon-based cross-validation was also used for the verification. Third, the candidate models verified using different polygon-based samples indicated that the No. 1 model using the RF method with a cost of 5000 and the No. 3 model based on the RF were the most favorable for predicting the larger polygon (area <sup>≥</sup> 100,000 m2) samples. For the smaller polygon (area of 50,000–100,000 m2) samples, the No. 3 model provided acceptable results. The accuracies of developed models were above 0.8 in most cases, outperforming the performances of the SVM and LR. Finally, mapping the predicted distributions revealed that the detection of landslide sources provided accurate results with lower missing rates. The extraction for the runout areas revealed excellent user's and producer's accuracies.

The gap in recognizing landslides has existed in terms of remote sensing and geotechnical engineering. More precisely, the landslide detections extracted from remotely sensed images without additional interpretations to eliminate the runout areas should be termed the landslide affected area and should not be directly treated as a landslide inventory or database for further analyses. This study provided a systematic procedure following a supervised data mining–based approach to separate landslide source and runout signatures. The tradeoff between the size areas of the training data and the predicted range of the developed models was emphasized in the results. More research is warranted regarding the effect of the elevation-based threshold and the binary classification based on the detections of landslide sources only to produce fully automatic landslide inventory productions.

**Funding:** This study was supported, in part, by the Ministry of Science and Technology of Taiwan under project No. 109-2121-M-035-001 and 108-2119-M-035-003.

**Acknowledgments:** The author would like to thank S.-H. Chiang in Center for Space and Remote Sensing Research, National Central University, Taiwan, for providing the landslide inventory and the digital elevation model. The author also thanks J.-Y. Huang in Department of Civil Engineering, Feng Chia University, Taiwan, for refining the figures in the article.

**Conflicts of Interest:** The author declares no conflicts of interest.

#### **References**


© 2020 by the author. 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* **Study on the Early Warning Methods of Dynamic Landslides of Large Abandoned Rockfill Slopes**

#### **Nan Qiao 1,\*, Yun-Ling Duan 1, Xiao-Meng Shi 2, Xue-Fei Wei <sup>3</sup> and Jin-Ming Feng <sup>1</sup>**


Received: 22 July 2020; Accepted: 28 August 2020; Published: 2 September 2020

**Abstract:** The excavation of large-scale underground projects produces a large amount of rubble waste material that is temporarily deposited near the project site, which forms a large-scale waste rockfill artificial slope. The slope has a granular structure, thus, during excavation and trans-shipment, surface shallow landslides may frequently occur. Existing contact monitoring methods such as buried sensors and GPS (Global Position System) are difficult to apply to the monitoring of rockfill landslides. Therefore, there are no appropriate early warning methods for waste rockfill slope landslides during dynamic transfer. Here, we used ground-based interferometric synthetic aperture radar to monitor the deformation of a rockfill slope during the excavation and transfer processes as a proposed method for the early warning against landslides on rockfill slopes during dynamic construction based on the radar interference measurement results. Through data cleaning and data interpolation, the line of equal displacement was generated, and the cross-sectional area of the equal displacement bodies of landslides was calculated. In addition, we established a four-level early warning grading standard, with the rate of change of the cross-sectional area of the equal displacement body as the early warning index, and realized real-time dynamic early warning of waste rockfill landslides during excavation and transportation. Finally, five landslide examples were used to verify the proposed warning method. The results show that the warning method can make an early warning 8–14 min before the occurrence of landslide, which can effectively avoid the appearance of catastrophic events.

**Keywords:**rockfill; ground-based interferometric synthetic aperture radar; construction; cross-sectional area of equal displacement body; landslide warning method

#### **1. Introduction**

In recent years, the construction volume of large underground projects in China has grown [1]. During underground engineering excavation, a large amount of rubble waste material is temporarily piled up near the project area to form a large-scale waste rockfill body. With the development and utilization of rubble waste material, the temporary waste material dump needs to be excavated and transported. Excavation disturbance inevitably causes slope failure. Therefore, the safety of waste rockfill slopes during excavation and transportation has become an important factor [2].

Slope deformation monitoring and landslide warning rank high among the most effective safety prevention and control measures [3]. Displacement monitoring forms the basis of landslide deformation monitoring. Slope displacement and its related parameters are important reference indicators for landslide early warning research [4]. In 1968, Saito proposed the "three-segment" change rule of slope displacement–time curves [5]. Since then, scholars have conducted a vast amount of research on slope deformation monitoring and landslide warning from different perspectives and with different means. For instance, Hoek (1970) analyzed the displacement–time curve measured at the Chuquicanata mine

and proposed an extension method to extrapolate the displacement–time curves and to predict the time of landslides according to the curve trend [6–8]. Zhong proposed the Verhulst inverse function model according to the characteristics of displacement–time curves [9]. Qiang et al. analyzed and studied many typical landslide displacement monitoring results and proposed a landslide early warning method using acceleration as the evaluation index [10]. In addition, many scholars have tried to combine various mathematical models with slope displacement–time curves. Such examples include the golden section method for early warning of the time of landslides [11–13], the markov early warning model [14,15], the fuzzy mathematics early warning model [16,17], the poisson cycle early warning model [18–21], and the gradient-sinusoidal early warning model [20,22,23].

For natural landslides, the displacement vectors at various points on the landslide body can generally be kept consistent. Even under special circumstances, the deformation of each point in the initial stage of slope deformation is disordered, but, once entering deformation acceleration, the displacement of each point on the landslide body and the time curve also tend to be the same. Therefore, the displacement data of one or more monitoring points can be used as the basis for slope deformation analysis and landslide warning research [24]. However, the waste rockfill are massive and unconsolidated and feature dynamic excavation and transportation processes. In slope deformation monitoring, it is difficult to apply contact monitoring methods such as embedded sensors and GPS. Because of the granular structure characteristics of a rockfill body, the rubble particles on the slope will rotate during slope excavation, which leads to inconsistent deformation among the rubble. It is irrational to describe the displacement and warn of the deformation of the slope by a single point. Therefore, traditional slope monitoring and landslide warning methods are not applicable for large abandoned rockfill slopes.

To effectively solve the problem of catastrophic events that may occur on the slopes of the rockfills, it is necessary to select appropriate monitoring methods, to describe the deformation process of slopes with scientific quantitative indicators, and to provide early warning for slope landslides. GB-InSAR (ground-based interferometric synthetic aperture radar) can provide full-field deformation information within the monitored area [25], which makes it possible to use the change of deformation region on the slope as the characteristic value to study the slope deformation. Therefore, the GB-InSAR system (IBIS-L system from Innovative Interferometric Radar for Environmental and Civil Engineering Applications, Italy) was used to carry out on-site monitoring tests on the slope of China's first temporary large-scale groundwater-sealed oil storage dumping site in Huangdao, Shandong Province. Through data processing of the full-field displacement results obtained from monitoring (data cleaning, data interpolation, generating equal displacement lines, etc.), we propose to describe the deformation process of the slope by using the curve of the cross-sectional area of the constant displacement body with time. By summarizing the results of the data, a landslide early warning method with an equal displacement cross-sectional area as an early warning index is proposed to achieve the real-time monitoring and early warning of rockfill landslides during construction.

#### **2. Project Overview**

The waste rockfill is located about 500 m southwest of the underground project entrance, and the floor area is approximately 9.6 <sup>×</sup> 104 <sup>m</sup>2. The ground where the waste rockfill body is located is high in the north and low in the south, and the difference of elevation between the south and north is about 20 m. The east and west directions are relatively gentle, and there is a slightly protruding depression in the middle part. The difference of elevation is about 4 m. The waste rockfill is dumped by a dump truck [26], and it is leveled by a bulldozer. The rockfill body slopes can be divided into two stages based on the formation process. The boundary between the two stages are shown in Figure 1 (the yellow dotted line).

**Figure 1.** Overview of the rockfill. The horizontal boundary width of the red curve is about 100 m, and the maximum vertical level difference is about 50 m.

The waste is mainly rubble produced by the blasting and crushing of granite gneiss. Table 1 shows the petrophysical and mechanical properties of granite gneiss in intact rock mass. The rockfill is highly porous and unconsolidated with high permeability of the rubble. The bottom of the rockfill body is composed of large rubbles with an average particle size of 1 m, and the rubble particles are independent of each other. The slope within the monitoring area (the red polygon area in Figure 1) is the residual slope after the occurrence of the previous landslides, and many large boulders are already broken, having an average size of 10–50 cm (Figure 2). Besides, due to the occurrence of landslide, the rock rolling and impact will be broken, resulting in much coarse-grained soil left on the slope.

**Table 1.** Physical and mechanical properties of intact rocks (from geological survey report).


**Figure 2.** Rock fragments left on the slope after the previous landslide.

Due to the stacking height limits and the economic value of rubble, rubble waste is transferred to a dumping site about 1 km south of the site of rubble transfer, as shown in Figure 3.

**Figure 3.** Relative position of project site (adapted from google maps). Where GB-InSAR is ground-based interferometric synthetic aperture radar.

To prevent uncontrollable large-scale slippage caused by an excavation-related disturbance, rubble excavation is carried out in two phases. In the first phase, a "C"-shaped platform is constructed above the platform of the original rockfill body. Excavation is carried out using construction equipment located on the platform (Figure 4) to control the location and slope height of the area affected by the excavation, thereby controlling the scale of the landslide and the direction of the sliding body. After the removal of most of the rubble via the above the platform, the second phase of the operation is to transport the rubble under the platform.

**Figure 4.** Loading and transportation.

#### **3. Slope Monitoring of Large-Scale Abandoned Rockfill Based on Ground-Based Interferometric Synthetic Aperture Radar**

#### *3.1. Monitoring Equipment*

This research used the IBIS-L system [2,27] for monitoring. The IBIS-L system is mainly composed of four parts: radar sensor, linear scanner, computer, and power supply module (Figure 5). When using IBIS for field monitoring tests, in addition to ensuring the normal monitoring of the target, the power supply of the equipment should also be considered. IBIS has a power of about 100 W. Even though a battery can be used as a power source, it does not effectively guarantee the continuous working time of IBIS.

**Figure 5.** IBIS-L system (from Innovative Interferometric Radar for Environmental and Civil Engineering Applications, Italy) composition.

The theoretical monitoring range of the IBIS-L system is 0–4000 m. Due to the limitation of the radar antenna beam angle, the azimuth monitoring range is −50◦ to 50◦. The specific system parameter settings are shown in Table 2. As for the working process of the system, the computer controls the radar sensor to move from the left end to the right end of the linear guide in an "off–on–off" manner. The length of the linear guide is 2 m, "off" every 5 mm, during which electromagnetic waves are emitted and echoes are received a total of 401 times. Then, the received signal is processed by synthetic aperture technology to realize azimuth focusing. Finally, according to the phase information of each pixel obtained by the synthetic aperture, the deformation information *d* of the target can be acquired by the interferometry of the phase difference between the echo signals of the targets at different times:

$$d = \frac{\lambda(\varphi\_2 - \varphi\_1)}{4\pi} \tag{1}$$

where ϕ<sup>1</sup> and ϕ<sup>2</sup> represent the phase information of each pixel obtained by synthetic aperture at different times, and λ is the radar wavelength. Thus, the displacement of each pixel in the whole field can be obtained.


**Table 2.** The parameters of IBIS-L.

#### *3.2. Monitoring Plan*

To make a preliminary judgment on the position of the site excavation and potential landslides, the optimal monitoring angle of GB-InSAR was first ensured, and, at the same time, the safety of the equipment and the power supply requirements were taken into account. The GB-InSAR was located at the project department entrance about 120 m to the northeast of the rockfill, as shown in

Figure 2. The range indicated by the red line in the figure is the area occupied by the piled rockfill, and the range indicated by the yellow line is the azimuth scanning range of the radar.

#### *3.3. Monitoring Results*

The monitoring lasted 69 days, during which time 11,376 observations were obtained and 72 shallow landslides (the sliding surface depth did not exceed 1 m) were recorded, including 41 small-area landslides that were destroyed within a short time by excavator disturbance. However, when the tops of some slopes were removed, it made it very easy to lose control, resulting in large-scale landslides. In total, 31 larger landslides were recorded, as shown in Figure 6, where the top of the slopes was removed. The records of the landslide area of the 72 landslides are shown in Figure 7, and the numbers represent the landslide areas, where the landslide area is the area of each landslide which was recorded during the field monitoring. The source of the final landslide area record was divided into two parts. Firstly, all the landslide area records that occurred during the day were recorded from the perspective of radar through visual estimation. Secondly, almost all the records of the landslide area from night to morning were obtained by the visual estimation of the construction workers standing near the bottom of the slope. Then, the landslide area visible to the radar could be calculated according to the relative position of the radar position, the construction personnel, and the landslide position. However, the landslide area records shown in this paper are estimated values of the landslide area from the perspective of radar. The landslide areas shown in Figure 7 are the final data that could be measured by radar.

**Figure 6.** Position of the sliding surface. The 72 landslides records can be divided into two types according to the position of the sliding surface as A and B. Examples of the two types of landslide surfaces are shown on the right.

**Figure 7.** Landslide areas (72 landslides). The landslide area is the area of each landslide recorded during the field monitoring. These data are not accurate measurements but intuitive estimates. Since this paper is based on IBIS-L measurement results for analysis, the area data in the figure represent the area of the landslide area that can be seen by the corresponding radar equipment.

According to the field observation, the surface morphology of the 31 landslides can be divided into two types. One is the slope surface above the excavation location of the mechanical equipment, which has no apparent large particle size and dense agglomeration zone. The other is the slope with a significant dense zone (Figure 8). Due to rainfall, the rolling of construction vehicles, and the uneven particle sizes, some layers or areas of the rockfill have a higher density than the surrounding dense zone. Therefore, the 31 landslides were classified based on the surface shape of the slopes (Table 3).

**Figure 8.** Slope with a significant dense zone.

**Table 3.** Classification of landslides.


This article focuses on the abovementioned 31 landslides. First, according to the monitoring results of IBIS-L, the echo signals of each pixel were extracted. Then, we dealt with the pixels that could not meet the precision requirement. The displacement values of the remaining pixels in each monitoring period were calculated by Equation (1). According to the displacement values of existing pixel points, the interpolation method could obtain the displacement values of all pixel points in the area around each cleaned point. Finally, the contour of equal displacement could be generated based on the values. Considering the limitation of displacement monitoring accuracy of IBIS-L, the displacement cross-sectional area of −0.5 mm was selected as the index for the landslide early warning research. The −0.5 mm equal displacement cross-sectional areas of 31 landslides at different times were extracted, as shown in Figure 9.

**Figure 9.** The area–time curve of the equal displacement cross section. The numbers in the legend represent landslide numbers recorded in chronological order. The "0" time represents the last scan before the severe landslides.

The two types of different slope forms exhibit two distinct curves under the same excavation construction conditions. The −0.5 mm equal displacement cross-sectional area of a Type B1 landslide lasts about 30–40 min before a landslide, and it is a cubic curve. The reason for this is the dense zone on the slope has a significant supporting effect on the upper rubble during the deformation of the slope. The −0.5 mm equal displacement cross-sectional area of a Type B2 landslide varies with time for about 30–120 min, and most last less than 60 min. The curve conforms to the exponential form with the accelerated deformation stage in the three sections of the slope displacement–time curve.

#### **4. Early Warning Method for Landslides Based on the Rate of Change of the Cross-Sectional Area of Equal Displacement Bodies**

The uniqueness of waste rockfill bodies makes it unsuitable to use single-point displacement as an indicator for slope deformation analysis or empirical judgment. Therefore, in this study, the deformation characteristics of waste rockfill slopes were characterized by the change of the cross-sectional area of the equal displacement body. On this basis, a warning method for waste rockfill slope landslides during excavation is proposed.

#### *4.1. Ideas on Early Warning*

(1A) Based on the results in [2] and the 26 landslides monitored in the previous period (the remaining five landslides were used to verify the early warning method), the surface characteristics of waste rockfill slopes can be identified, and a preliminary judgment can be made about the possible landslide patterns in the excavation process. During the accumulation of waste rockfill, the density of some accumulation layers and areas caused by rainfall, the rolling of construction vehicles, etc. is quite different, or there might be a certain amount of abnormally large-sized rubble on the slope surface and in the shallow layer. Both conditions play a great role in supporting the rubble particles above the dense zone. Therefore, combined with the surface characteristics of the rockfill slope, a possible landslide pattern can be determined through the analysis of the characteristics of the cross-sectional area of different displacement bodies over time.

(1B) For possible landslides being monitored, based on the real-time constant displacement cross-sectional area data obtained from slope deformation monitoring, a time series prediction model such as ARIMA-SVR [28] (Autoregressive Integrated Moving Average model-Support Vector Regression: a hybrid ARIMA and SVR model of time series prediction mode) can be used. Based on the existing data and the prediction of future data, the future trend of constant displacement cross-sectional areas can be predicted.

(2) Real-time monitoring data can be extracted, the real-time changing process of the cross-sectional area of equal displacement can be studied and analyzed, and the mode of slope deformation can be determined according to the change trend of the cross-sectional area of the equal displacement body.

(3) According to the measured (1A) and predicted values of the cross-sectional area of the equal displacement body (1B), the relationship between the cross-sectional area change of the equal displacement body in different modes and landslides of the rockfill slope can be analyzed, and an early warning judgment of such landslides can be determined. Then, the change process of the cross-sectional area of different landslides and equal displacement bodies can be independently analyzed.

(4) On the basis of the determined early warning criteria of the rockfill slope, an early warning grade division standard of the rockfill slope can be established.

The detailed process is shown in Figure 10.

**Figure 10.** Early warning method of landslides.

#### *4.2. Early Warning Indicators and Grading Standard*

Considering the needs of on-site early warning in terms of the number of data calculations and the data extraction speed, under the condition of ensuring the accuracy of radar displacement monitoring, the deformation of −0.5 mm on the slope is easier to identify in time. Therefore, the cross-sectional area of the −0.5 mm equal displacement body was used as an index to validate the early warning method of the rockfill slope during excavation and transportation [2]. The on-site conditions of waste rockfill slope landslides are complicated, and both different and similar types of landslides often occur in the same area. Therefore, the cross-sectional area of the −0.5 mm equal displacement body was normalized, and then the rate of change was used as a warning indicator of landslides. The cross-sectional area of the −0.5 mm constant displacement body at each moment after normalization is

$$S\_{-0.5}^{\*}(t) = \frac{S\_{-0.5}(t)}{S\_{\text{max}}} \tag{2}$$

where *S*−0.5(*t*) is the measured or predicted value of the cross-sectional area of the −0.5 mm constant displacement body at time *t* and Smax is the cumulative value of the cross-sectional area of the −0.5 mm equal displacement body obtained from the last measurement before the slope starts sliding.

The normalized −0.5 mm equal displacement body cross-sectional area change rate is

$$V\_{-0.5}^{\*}(t) = \frac{S\_{-0.5}(t)'}{S\_{\text{max}}} \tag{3}$$

Based on the field measurement results and the analogy between landslide cases, the characteristics of the cross-sectional area of the −0.5 mm equal displacement body of the two types of 29 landslides that were previously monitored in the rockfill were statistically analyzed. Then, early warning standards were developed according to two different types of landslides. Since the data selected were actually measured landslide data, the lower limit of the 95% confidence band of the mean value of the cross-sectional area of the two types of landslides of the −0.5 mm equal displacement body was selected as the standard to formulate the landslide warning level. According to the characteristic values of 50%, 30%, and 10% of the lower limit of the 95% confidence band of the average displacement rate of the cross-sectional area of the −0.5 mm equal displacement body before sliding, the landslide warning was divided into five levels, as shown in Figure 6. Note that the rate of change of the cross-sectional area of the "B1"-type landslide of the −0.5 mm equal displacement body shows a decreasing trend within a certain period of the initial stage of the development of the displacement field. Therefore, the maximum value of *V*∗ −0.5(*t*) before the start of deformation acceleration is defined as the upper limit of the yellow warning for class B landslides. To make the early warning method more convenient and to ensure the compatibility of the two types of landslides, the early warning standard was uniformly divided into four levels based on the statistical results of the early warning thresholds of the two types of landslides, as shown in Figure 11. The grading results and construction regulations were then established, as shown in Table 4.

**Figure 11.** The lower limit of the 95% confidence band for *V*∗ −0.5(*t*)


**Table 4.** Early warning classification and dynamic control method of the waste rockfill slope during excavation and transportation.

#### **5. Method Validation**

To check the feasibility and reliability of the landslide warning method of rockfill slopes, the last five landslides of two different types of landslides recorded during field monitoring were used to verify the landslide warning method.

The waste rockfill slope was taken as the remaining broken surface after the previous landslide, and the slope angle was approximately 35–38◦. The excavation and loading speed during construction was approximately a 10-t load muck truck loaded every 6–9 min. During the five landslide monitoring periods, there was no rainfall, the daytime temperature was about 20–25 ◦C, and the nighttime temperature was about 15–20 ◦C.

The landslides occurred at 4:55 on 5 October 2013, 2:45 on 10 October 2013, 3:00 on 15 October 2013, 21:00 on 15 October 2013, and 9:38 on 16 October 2013. The final landslide areas were 450, 1000, 1100, 1600, and 1100 m2, respectively. The specific records are shown in Table 5.


**Table 5.** Characteristics of the five test landslides.

First, the slope was initially classified according to whether there was a dense belt formed by the rolling of the slope or whether there was abnormally large-sized rubble in the middle of the slope. Landslides 44 and 53 conformed to the B1-type failure mode, while Landslides 67, 71, and 72 conformed to the B2-type failure mode. Then, ARIMA-SVR model was used to predict the variation trend of −0.5 mm equal displacement cross-section area of these five landslides, and the results are shown in Figure 12. The average deviations between the predicted value and the measured value of five test landslides are shown in Table 6.

**Figure 12.** The measured and predicted values of the cross-sectional area of the −0.5 mm equal displacement bodies.

**Table 6.** The average deviations between the predicted value and the measured value of the five test landslides.


According to Equation (2), the prediction results of five −0.5 mm equal displacement bodies of five landslides were normalized. Then, the deformation rates of the cross-sectional area of the −0.5 mm equal displacement bodies were calculated after normalization using Equation (3). The results are shown in Figure 13.

**Figure 13.** The normalized −0.5 mm equal displacement body cross-sectional area change rate.

From the analysis results of the above examples, the landslide early warning method and the early warning classification established in this paper could provide a warning to avoid risks before the occurrence of all five tested landslides.

Compared with the traditional landslide warning method based on slope deformation, it mainly depends on the displacement of some points on the slope. For a loose rockfill slope, it is feasible to express the deformation of the slope by the displacement of a certain point, but it is difficult to find a suitable point. However, as described in Section 3.3, the variation curve of a −0.5 mm equal displacement cross-sectional area over time conforms to the form of a "three-section" curve. Combined with previous research experience, this will make it easier for the equal displacement cross-sectional area to be used as an early warning indicator.

It should be noted that this verification was based on the last instance of landslide of each type on the site, and the actual working conditions were similar to the previous landslides. The establishment of the landslide warning standard also depends on the statistical results of these previous similar landslide cases. Therefore, these four successful landslide warnings are limited by certain working conditions, such as rock lithology, rubble particle gradation, rolling of vehicles during the accumulation of rock mass, rainfall, and so on. These conditions will directly affect the repose angle of the rockfill, and the repose angle is an important factor affecting rockfill landslides. Generally, when the angle of the slope is less than the repose angle, the slope is unstable. Rainfall has a direct impact on the water content of rockfill, which is directly proportional to the angle of internal friction. For loose rockfill, due to the structural characteristics of the large pores, leakage and evaporation are rapid, so water content is usually small. Particle gradation will also directly affect the repose angle of the rockfill body. Small-particle-size rubble plays a particular "lubrication" role in the rockfill body. When the content of small-particle-size rubble is relatively low, large-particle-size rubble relies on friction and the "aggregate interlock capacity of rubble" to keep the slope in a stable state. When the content of small-particle-size rubble is high, due to the "lubricating" effect of small-particle-size rubble particles, the relative movement between large-particle-size rubble particles is more likely to occur, resulting in slope instability. Therefore, for the slopes of other rockfill bodies under similar working conditions, the actual situation should be considered to make appropriate adjustments to the early warning standard.

In general, most landslide warning models are based on the displacement–time curve of a single point on the slope. According to the curve change that conforms to the three stages of "initial deformation–uniform deformation–accelerated deformation", the landslide warning levels are divided into "caution–vigilance–danger". For example, Wang [29], He [30], and Qiang [31], among others, all studied the landslide early warning index and landslide early warning classification standard based on the "three stages" of the single point displacement–time curve. However, which point of displacement data should be selected as the basis in the process of landslide warning is a fundamental question. To solve this problem, Qiang [24] proposed that the deformation of each point on the slope is disordered before the basic penetration of the sliding surface. However, after the basic penetration of the sliding surface, the slope begins to slide as a whole. As mentioned above, the rockfill body is loosely piled, and the deformation of adjacent positions on the slope is disordered before the occurrence of landslide. Moreover, under the condition of excavation disturbance, a shallow landslide is likely to occur in a short time, and the time between the connection of the sliding surface and the occurrence of landslide is generally quick. The landslide early warning method established in this paper describes the deformation process of the slope through the area of equal displacement cross-section, avoiding the problem of selecting appropriate points on the slope, which has significant importance for landslide warning of loose accumulation.

#### **6. Conclusions**

Due to the granular structure of waste rockfills, it is difficult to find a suitable single-point displacement as an early warning index. The advantage of GB-InSAR full-field scanning is that it helps to obtain slope displacement and other displacement information (the electromagnetic wave reflection intensity of the measured target, etc.) through data processing. In this paper, a landslide early warning system based on the rate of change of the cross-sectional area of equal displacement bodies was proposed for the landslides of unconsolidated coarse sedimentary rocks.

By summarizing the previouslandslidemonitoring and data processing results, using the normalized rate of change of the cross-sectional area of a −0.5 mm equal displacement body as the early warning criteria, two types of landslides and four levels of early warning classifications were established as follows: (a) Level 1 (green) warning area, *V*∗ −0.5(*t*) <sup>&</sup>lt; 0.01; (b) Level 2 (yellow) warning area, 0.01 <sup>≤</sup> *<sup>V</sup>*<sup>∗</sup> −0.5(*t*) <sup>&</sup>lt; 0.02; (c) Level 3 (orange) warning area, 0.02 ≤ *V*<sup>∗</sup> −0.5(*t*) <sup>&</sup>lt; 0.03; and (d) Level 4 (red) warning area, *<sup>V</sup>*<sup>∗</sup> −0.5(*t*) <sup>≥</sup> 0.03.

The validity of the early warning method was verified based on five landslides to compare the predicted and real values in this paper. Further research is needed to take into consideration comprehensive factors such as rock mechanics and rainfall for early warning of rockfill landslides.

**Author Contributions:** N.Q. contributed to conceptualization, numerical and experimental investigations, and preparation of the paper; Y.-L.D. contributed to conceptualization, numerical and experimental studies, and preparation of the paper; X.-M.S. contributed to conceptualization, numerical and experimental investigations, and preparation of the paper; X.-F.W. contributed to conceptualization, numerical and experimental investigations, and preparation of the paper; and J.-M.F. contributed to experimental investigations. All authors have read and agreed to the published version of the manuscript.

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

**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*
