Instability Index Derived from a Landslide Inventory for Watershed Stability Assessment and Mapping

Abstract

Watersheds represent natural units of social–ecological systems and affect crop productivity. Extreme weather events accelerate the natural erosion process by triggering more landslides in watersheds. To achieve the land degradation neutrality set up by the UN’s Sustainable Development Goals, it is necessary to assess and map spatiotemporal landslides in watersheds. This paper proposes an innovative approach to calculating the instability index by preparing an annual landslide inventory, determining the optimum sub-watershed, compensating for shadow effects on the time series of the landslide area ratio, and classifying the standard deviations to different levels of instability. Taking the Qingquan watershed as an example, the instability index calculated for 22 sub-watersheds makes it possible to identify hot spots that are prone to collapse. This new index can also be used to evaluate the effectiveness of watershed management before and after completion of a specific engineering project, as well as to update the latest upriver situation to evaluate current management practices and develop strategies for future planning. Based on this new approach, the Soil and Water Conservation Bureau of Taiwan assesses the stability of 28 watersheds, and the results are made available on the Big Geospatial Information System.

1. Introduction

Watersheds represent natural units of social–ecological systems [1] that channel rainfall to creeks/streams/rivers and eventually into reservoirs, bays, and the ocean. The significant changes in the hydrology of different watersheds are mainly caused by soil erosion, which also affects the infiltration rate, water-holding capacity, nutrient availability and organic matter content, soil depth and soil biota, and, hence, crop productivity [2]. Extreme weather events worsened by global warming modulate the intensity and frequency of precipitation and furthermore accelerate the erosion process by triggering more landslides in watersheds. Although the mechanism is clear [3], and the threat of ever-increasing floods is real [4,5], the consequences and long-term effects of compound disasters have not been investigated in a systematic way. To achieve land degradation neutrality as set up by the UN’s Sustainable Development Goals, it is necessary to assess and map the spatiotemporal degradation of land and soil, particularly in watersheds.

The most common approach to assessing the stability of watersheds is to compute the landslide susceptibility index (LSI) based on preparatory factors weighted according to importance. Koukis and Ziourkas [6] conducted a comprehensive review and discussed more than 64 preparatory factors. Süzen and Kaya [7] reviewed 145 studies spanning the period from 1986 to 2007 and concluded the five most influential natural factors to be lithology/geology, slope angle, landuse/landcover, drainage network, and slope aspect. Such approaches all rely heavily on the landslide inventory used to derive the weight of each preparatory factor. The factors or parameters used in such approaches, however, are mostly topographically-related and are cell-based local parameters that do not consider the influence of neighboring cells [8]. For regions with frequent landslides, the variations in topography are significant. All topographically-related factors, including the slope angle, slope aspect, drainage network, and even landuse/landcover may change in the future due to possible changes in the underlying physical and mechanical processes. They must be remeasured to reflect the latest status of landslide susceptibility. For this reason, the accuracy of the cell-based LSI approach is limited. One indirect way to assess watershed stability is through estimating the temporal and spatial stability of soil moisture through the use of remote sensing, such as was done in works by Cosh et al. [9] and Hu et al. [10]. The soil moisture, however, is closely related to precipitation, which is more like a triggering factor rather than a preparatory factor related to landslides.

Rowbotham and Dudycha [11] pointed out that the use of geomorphometrically significant terrain units extracted from a digital elevation model (DEM) are an efficient alternative to approaches using regular grid cells. This motivates an assessment of stability in terms of a sub-watershed unit rather than a point, using the spatial and temporal distribution and relationships revealed from a long-term landslide inventory rather than using topographically-related and cell-based local factors. Such an approach requires a long-term, detailed landslide inventory to understand the dynamic changes associated with landslides, yet this is not available for most parts of the world. To serve as an ideal area of study, the region under consideration not only has to experience frequent landslides, but it also must be observed intensively by satellites, and there must be a rapid—yet accurate—approach to distinguishing between landslides and the shadowed areas within images of the region [12].

Together with high mountains and frequent earthquakes, Taiwan has one of the highest erosion rates in the world [13] and is extremely vulnerable to landslides [14]. Taiwan’s government, therefore, has been endeavoring to manage watersheds in order to mitigate the issue of flooding by developing various measures to assess watershed stability. The National Space Organization (NSPO) of Taiwan operated Formosat-2 to acquire high-spatial-resolution (2m) imagery of Taiwan from 2005 to 2016, which likely represents the most intensive observation from space for a country with such frequent occurrences of landslides. The Forest Bureau of Taiwan (FBT), in the meantime, funded three multi-year projects in a span of seven years to generate a long-term, detailed landslide inventory for the entire Taiwan region using annual composite Formosat-2 images acquired between 2005 and 2016 [15]. This inventory provides rich spatial information about the ever-changing surfaces of watersheds reformed by the ever-increasing floods accelerated by global warming.

This paper proposes an innovative approach to calculating the instability index from a landslide inventory time series analysis. The instability index is used to categorize the sub-watershed stability in order to identify hot spots that are prone to collapse. Instead of computing the LSI based on some preparatory factors that are not available (such as detailed geological composition) or not measured on a regular basis (topographically-related factors), future landslides are assumed to be highly correlated to historical landslides. A similar assumption was also made by Di Paola et al. [16] and Stamatopoulos et al. [17]. Since a satellite image time series is usually acquired with different viewing geometries under different atmospheric conditions at different times, the results for the area variations in a single landslide derived from these images are inevitably biased by cloud coverage or shadows. Therefore, we propose an optimum sub-watershed by considering both the hydrological process and geological zoning. In each optimum sub-watershed, the time series of the calculated landslide ratio retains the true variation trends, and its standard deviation provides a new index, called herein the instability index, to categorize the sub-watershed stability. This new index makes it possible to quantify the degree of variation for each governance unit, thus identifying the hot spots that are prone to collapse. It can also be used to evaluate the effectiveness of watershed management before and after completion of a specific engineering project, as well as to update the latest upriver situation to evaluate current management practices and develop strategies for future planning.

2. Materials and Methods

2.1. Landslide Inventory

A landslide inventory is a detailed record of the distribution and characteristics of past landslides [18]. Among various landslide mitigation strategies, a detailed landslide inventory is critically needed in landslide-prone regions [19]. Typhoon Morakot (6–10 August 2009) brought severe debris flows/landslides and a significant loss of lives/properties in an extreme case of torrential, sustained rainfall. After this disaster, the FBT decided to make the best use of Formosat-2 imagery, and three multi-year projects were funded over a span of seven years to generate a long-term, detailed landslide inventory for the entire Taiwan region using annual composite Formosat-2 images acquired between 2005 and 2016 [12,20]. All Formosat-2 imagery was preprocessed with the Formosat-2 automatic image processing system (F-2 AIPS) [21], which is able to digest the raw Gerald format data and apply basic radiometric and geometric corrections, as well as conduct rigorous band-to-band coregistration [22], automatic orthorectification [23], multi-temporal image geometrical registration [24], multi-temporal image radiometric normalization [25], spectral summation intensity modulation pan-sharpening [22], edge enhancement, adaptive contrast enhancement, and the absolute radiometric calibration [26].

Based on the characteristics of areas of destroyed vegetation that are spectrally distinct from the background, various classification approaches have been proposed to map landslide areas directly from remote sensing imagery using fully automatic methods [27,28,29], such as supervised classification [30,31], unsupervised classification [32], and the ratio image [33]. However, cultivated lands, roads, houses, riverbeds, and even clouds are usually misclassified as landslides from high-spatial- but low-spectral-resolution imagery. For applications that require accurate information about landslides, the traditional method of manually delineating landslide regions by visually interpreting images is still more reliable [34,35,36,37,38,39]. Such an approach, however, is time-consuming and laborious. The Expert Landslide and Shaded Area Delineation System (ELSADS) was developed to facilitate with determining rather than delineating landslide and shadow areas in a quick and accurate fashion [15]. A true-color composite of Formosat-2 image is illustrated in Figure 1a, where the region annotated by the red box (Figure 1b) is at full resolution with locally enhanced contrast. Meanwhile, the same region is displayed in Figure 1c as a standard false color composite and draped onto a digital terrain model (DTM) that can be transposed, rotated, and scaled with a mouse. The regional thresholds of the histogram of the panchromatic band (Figure 1d), the normalized difference vegetation index (NDVI) histogram (Figure 1e), and the normalized green red difference indices (NGRDI) histogram (Figure 1f) are determined automatically and annotated as the vertical red lines in each window. The non-vegetated areas (green polygons) and the dark areas (white polygons) are determined by these thresholds and overlaid on (Figure 1b) and (Figure 1c) automatically. By switching the window shown in Figure 1b to the previous images or other geospatial layers, such as a geological map or a land use and land cover map, and viewing the window of Figure 1c from different angles or at different scales, the image interpreter simply selects a polygon (turns red) or deletes a polygon (vanishes), and the results are exported and stored automatically. The determined landslide areas and shadow areas are annotated as yellow and blue polygons, respectively [15]. Scanning scene-by-scene, the ELSADS enables the image interpreter to integrate all the useful geospatial information and determine rather than delineate the landslide and shadow areas in a quick and accurate fashion.

To support the open data policy, the FBT released this long-term, detailed landslide inventory. Together with the annual Formosat-2 composites acquired between 2005 and 2016, the SWCB makes them available in the form of vector and image tiles published on the SWCB Big Geospatial Information System [40], as shown in Figure 2. It should be noted that the FBT also sets up a threshold of 10 ha to categorize severe landslides. The geometric attributes such as area, perimeter, mean slope, mean aspect, mean elevation, center coordinates, etc. of each landslide are calculated and registered for future analysis.

2.2. Quantitative Assessment

2.2.1. Governance of Watershed Management and Flood Mitigation Units

The SWCB is responsible for soil and water conservation on upper slopelands and for flood-controlled slopelands. In compliance with policies and verified budgets, the watersheds in Taiwan are divided into governance units (GU) based on classification levels listed as A, B and C, as shown in

Table 1. Classified level of mountain and flood control of management units.

Classified Level	Principle	Strategy
Maintenance area (level A): Previous management is effective; it is a low disaster risk area.	Initiation zone: landslides have been stabilized Transportation zone: the necking zone has been dredged. Accumulation zone: the riverbed changes have been stabilized.	The inspection and maintenance of structures The monitoring of watershed conditions A disaster prevention and evacuation program
Staging management area (level B): Previous management has not been completed, but the assessment shows the management is effective in terms of reducing disaster risk.	Initiation zone: the landslide rehabilitation has not been completed. Transportation zone: partial necking zone is dredged. Accumulation zone: the sediment deposition potential remains high.	Yearly management project and rolling-wave planning A disaster prevention and evacuation program
Fundamental protection area (level C): Previous management has had no obvious effectiveness; it is a high disaster risk area.	Initiation zone: large landslides are difficult to rehabilitate, so engineering projects cannot be carried out. Transportation zone: there is significant sediment production; there is reoccurrence of necking zones. Accumulation zone: the riverbed accumulation is serious.	Adopting low intensity construction methods before relocation Adopting temporary mitigation, disaster prevention, and fundamental disaster control Enhancing disaster prevention and evacuation

. The annual budget is justified and allocated to each GU according to its classified level, protected targets, potential debris flow torrent, and landslide ratio. At the time of this writing, a total of 1725 potential debris flow torrents had been designated and announced by the SWCB (http://246.swcb.gov.tw/debrisinfo/GoogleMap.aspx; browsed on 22 February 2019). Among all these criteria, the classification level is the most critical, yet it is the only condition that cannot be assessed quantitatively. Basically, the classification level should be associated with watershed stability. Therefore, some general principles are employed to form a few descriptive, qualitative questions that can be used, for example, to determine whether the landslide area has been rehabilitated and has achieved stability, whether the necking zone has been dredged, and whether the riverbed changes in the accumulation zone have been stabilized, or the accumulation problem is still serious, and the sediment deposition potential remains high. According to the answers to such questions and the improvements made every year, two extra sub-levels (+ or -) are given for each classification level, resulting in a total of nine levels (C-, C, C+, B-, B, B+, A-, A, A+). This motivated an effort to look for a quantitative approach to assess watershed stability.

Another issue of the current GU is that it is divided with considerations of not just the drainage system and geological zone, but also the administrative district. The coverage of each GU is somehow too large to identify hot spots with lower stability. As also pointed out by Rowbotham and Dudycha [11], the use of geomorphometrically significant terrain units extracted from a DEM are an efficient alternative to approaches using regular grid cells. This also motivated an effort to devise a new strategy by which to specify an optimum sub-watershed.

2.2.2. Optimum Sub-Watershed

One direct, intuitive way to assess watershed stability is to focus on each landslide and track the variations in its area over time. However, even with the same satellite, different images are usually acquired with different viewing geometries under different atmospheric conditions at different times. The area variations in a single landslide derived from these images are inevitably biased by cloud coverage or the presence of shadows. One example is a case of rapid response mapping of significant landslides triggered by Typhoon Morakot in August of 2009. Most of the available aftermath images were actually acquired in the winter of 2009 or in the spring of 2010. During that period of time, the sun elevation was low in the subtropical zone, when shadows could occupy as much as 30% of an entire image over a mountainous area [41]. As a result, some delineated landslide areas decreased after Typhoon Morakot. To avoid such misleading results, it is necessary to introduce some grouping mechanisms.

From the perspective of hydrological processes, slopeland stability is decreased by the number of landslides occurring in the same sub-watershed, both in the upstream and downstream directions. If a broader sub-watershed is specified, there will be a greater number of landslides. From the perspective of geological zoning, each sub-watershed should be comprised of the same geological composition to keep this preparatory factor the same in each sub-watershed. When a smaller sub-watershed is specified, a closer geological composition is likely to be found. This is equivalent to finding an optimized solution: the sub-watershed must be specified according to some criteria. We started by conducting a simple experiment in which various thresholds of accumulated flow were specified to obtain different river systems using the hydrology module that comes with the ArcGIS® spatial analyst tool. As shown in Figure 3, the landslide polygons are overlaid on three river systems calculated with three values of accumulated flow (300, 600, and 900). The number of sub-watersheds increases as the threshold of accumulated flow decreases, and more landslide polygons across the boundary are divided into smaller polygons. This is not preferred because any landslide larger than 10 ha will be categorized as a severe landslide and will be monitored by the FBT. The river system derived from a smaller threshold of accumulated flow is more fragmented, and the total number of severe landslides is less.

We repeat the iteration by increasing the accumulated flow at an interval of 100 and reach an optimum threshold of 1000, which maximizes the total number of severe landslides in the watershed, shown as the yellow dot in Figure 4. Note that this optimum sub-watershed is purely based on the concept of grouping landslides into a unit. It is easy, consistent, and reproducible and is mainly controlled by the hydrological process. To our surprise, the boundary of the out of optimum sub-watershed matches most of the existing forest compartment boundaries quite well, but reflects more realistic situations because it is the latest DEM considered in the hydrological process, and there is no administrative district issue involved.

2.2.3. Instability Index

With the well-defined optimum sub-watershed, the variations in landslides are no longer monitored or discussed one by one, but rather are added together and represented as a whole for each sub-watershed. To support the preparation of the annual landslide inventory, a cloudless scene can be selected within that year. The bias introduced by shadows, however, is inevitable in mountainous areas [41]. This can be compensated for by assuming that the landslide area ratios are the same both inside and outside of the shadows, as discussed in Liu et al. [42]. Figure 5 gives an example of shadow compensation in the watershed for a potential debris flow torrent Tainan DF045, and Figure 6 shows all of the thumbnail images of the watershed. The shadow area ratios are indeed rather large in some years, where the corrections ranged from 16.7% in 2013 to 110.6% in 2009. The overall trends were approximately the same before and after the shadow compensation.

The same approach can be applied to each optimum sub-watershed to obtain its landslide area ratio time series, for which the standard deviation provides a quantitative indicator of the stability, namely the instability index. The computational flowchart revealing the steps used to calculate the proposed instability index is given in Figure 7.

Generally speaking, a lower instability index value implies that the variations in the landslide area ratio are fewer in number and that the optimum sub-watershed is relatively stable, and vice versa. In other words, the stability of a watershed can be expressed as a set of the instability indexes calculated from the 12-year landslide inventory in every optimum sub-watershed. An objective way to rank the stability of all optimum sub-watersheds is to cluster these instability indices by applying the Jenks natural breaks classification method in order to determine the best arrangement of values divided into five classes of instability levels, as listed in

Table 2. Grading scale chart for instability levels.

Instability Levels	Grading Scale	Explanation
1	Stable area	Landslides never occurred, or only a few landslides occurred in this area with good natural restoration. The area is not easily destroyed.
2	Low active landslide area	Some landslides occur, but the landslide area is small, or the restoration after the landslide is obvious.
3	Unstable landslide area	Some landslides occur with massive landslide areas, or the landslide restoration is less obvious.
4	Severe landslide area	Massive landslide area, or an area where landslides occur frequently
5	Large-scale landslide area	An area where large-scale landslides occur

. Noted that the thresholds of the five instability levels determined using the Jenks natural breaks classification method are not fixed/absolute values, and they are not necessarily valid in other watersheds. A future plan has been made to define absolute thresholds of five instability levels that can be applied to any watershed in Taiwan. The instability index calculated in each optimum sub-watershed makes it possible to assess its stability in a quantitative fashion and to identify hot spots that are prone to collapse. It can also be used to evaluate the effectiveness of watershed management before and after completion of a specific engineering project.

3. Results

This innovative approach has been employed by the SWCB to quantitatively assess the stability of a total of 28 watersheds, and the results have been published on the SWCB Big Geospatial Information System [40]. The Qingquan watershed in Wufeng Township is used as an example to explain and discuss the results. Typhoon Aere devastated this watershed in 2004, triggered a total of four potential debris flow torrents, and caused several large-scale landslides. A total of 22 optimum sub-watersheds are specified in the Qingquan watershed, and the main protected targets include the Tuchang Tribe, the Qingquan Tribe, and the Minduyou Tribe. The annual landslide inventory is overlaid on the 22 optimum sub-watersheds to illustrate the spatial and temporal distribution and its relationship to landslides in Figure 8.

One remarkable feature of the spatial distribution is that most landslides actually occurred and reoccurred in a few sub-watersheds. Another feature is the rapid restoration of landslides in some sub-watersheds, which refer to some major disaster events or specific engineering projects. To clarify the cause–effect relationship, the annual landslide areas of every sub-watershed, as well as the derived instability index and level of instability are listed in

Table 3. Annual landslide area, instability index, and level of instability of 22 optimum sub-watersheds in the Qingquan watershed in Wufeng Township, Hsinchu County.

No	Area (ha)	Annual landslide Area (ha)												Instability Index	Instability Level
		2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	2016
20894	11.23	0.37	0.00	0.26	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.01110	4
21735	324.94	9.36	5.81	9.98	6.36	3.81	4.11	4.75	1.46	2.56	4.05	2.71	1.80	0.00843	4
21897	372.25	18.69	11.93	14.92	11.09	11.71	11.19	10.75	7.65	7.41	10.33	9.08	3.93	0.01003	4
22062	108.72	1.13	1.11	1.52	1.04	2.51	1.14	0.95	0.00	0.00	1.30	0.54	0.70	0.00620	4
22253	298.04	16.24	10.24	29.37	29.58	15.86	15.17	21.38	14.67	12.43	16.27	19.87	13.98	0.02062	5
20610	150.57	0.53	0.42	0.28	0.00	0.00	0.00	0.00	0.00	0.00	0.27	0.16	0.00	0.00128	2
20759	11.29	0.18	0.00	0.29	0.33	0.00	0.00	0.21	0.00	0.14	0.00	0.00	0.00	0.01130	5
20928	163.21	0.47	0.10	0.00	0.00	0.19	0.00	0.00	0.70	0.71	0.36	0.36	0.38	0.00162	3
21045	89.31	0.32	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.26	0.00	0.00	0.00126	2
21084	106.23	0.11	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00029	1
21100	155.90	1.84	1.10	1.42	0.00	0.00	0.39	0.00	0.10	0.20	0.00	0.00	0.00	0.00420	3
21191	270.83	1.16	0.23	0.50	0.00	0.51	0.18	0.70	0.47	0.29	1.27	1.09	0.83	0.00153	2
21414	145.41	2.46	0.00	1.47	0.54	0.00	0.27	0.00	0.00	0.10	0.42	0.70	0.00	0.00520	3
21704	277.64	0.83	0.00	0.32	0.36	0.21	0.00	0.00	0.20	0.00	0.11	0.00	0.00	0.00089	2
22057	76.77	4.40	5.41	5.81	4.12	3.36	5.82	2.34	1.40	2.07	0.23	0.00	0.00	0.02899	5
22061	7.39	0.12	0.12	0.17	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00853	4
22191	117.37	0.06	0.21	0.39	0.00	0.31	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00120	2
21628	73.18	0.00	0.16	0.00	0.00	0.39	0.00	0.00	0.77	0.00	0.00	0.00	0.00	0.00327	3
20713	175.81	0.00	0.00	0.45	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00073	2
21512	58.93	0.00	0.00	0.17	0.00	0.00	0.00	0.12	0.00	0.00	0.00	0.00	0.00	0.00098	2
21793	10.98	0.00	0.00	0.00	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.00026	1
21320	207.05	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.20	0.22	0.00	0.00040	1

for comparison.

Three of the 22 optimum sub-watersheds are selected to represent cases of a low active landslide area (No. 20713), a severe landslide area (No. 21735), and a large-scale landslide area (No. 22253). The landslide area ratio time series of these three cases are shown in Figure 9 for comparison, while the annual thumbnail images overlaid with the corresponding landslide polygons are given in Figure 10 (No. 20713), Figure 11 (No. 21735), and Figure 12 (No. 22253). Instead of comparing the landslide area ratios, the annual thumbnail images illustrate the spatiotemporal variations in severe landslide areas in each optimum sub-watershed.

No. 20713 is basically a stable optimum sub-watershed with an area of 175.81 ha. Because this sub-watershed is located on the leeward side of Maipalai Creek, it is seldom influenced by rainfall. Only a few landslides (0.45 ha) occurred in 2007, and the area quickly recovered in later years, as shown in the annual thumbnail images (Figure 10). Note that the landslide area ratio of case No. 20713 is one order less than the one of the other two cases and is plotted against the right y-axis in Figure 9. The level of instability of the optimum sub-watershed No. 20713 is categorized as 2, which reflects the actual situation pretty well.

The optimum sub-watershed No. 21735 has the largest area (324.94 ha). Several severe landslides were triggered in the tributary of Tuchang River and the potential debris flow torrent of DF044, as shown in the annual thumbnail images (Figure 11). The landslide area ratio time series reveals that this sub-watershed was gradually stabilized from an annual landslide area of 9.36 ha in 2005 to 3.81 ha in 2009, yet it increased again to 4.75 ha in 2011, decreased the following year and increased again to 4.05 in 2014. In other words, this sub-watershed exhibits a repeated pattern of landslides. The level of instability of the optimum sub-watershed No. 20713 is categorized as 4 (severe landslide zone), which reflects the actual situation pretty well.

No. 22253 is also a large optimum sub-watershed with an area of 298.01 ha. Quite a large scale landslide was triggered by Typhoon Aere in 2004, and the total landslide area had just reduced from 16.24 ha in 2005 to 10.24 ha in 2006, when a sudden almost 3-fold expansion (29.37 ha) was found the next year, as shown in the annual thumbnail images (Figure 12). This sub-watershed exhibited a similar pattern as that of a repeated landslide, where the average landslide area was never less than 10.0 ha throughout the years. Note that 10 ha is a threshold by which to define a severe landslide by the FBT. Although it is not a single landslide, but a total landslide area that is larger than 10 ha, it is worth paying attention to this sub-watershed. The level of instability of optimum sub-watershed No. 22253 is categorized as 5 (large-scale landslide zone), which reflects the actual situation as well.

Figure 13a shows the locations and regions of 28 watersheds assessed by the SWCB (blue polygons), where Qingquan watershed is annotated as a star. Figure 13b illustrates the quantitative assessment of Qingquan watershed stability based on the instability index derived from a long-term landslide inventory. Each optimum sub-watershed is annotated with one color to indicate its level of instability. Together with the long-term, detailed landslide inventory and the annual Formosat-2 composites acquired between 2005 and 2016, all 28 watershed stability assessments have been made available on the SWCB Big Geospatial Information System [40].

4. Discussion

An innovative approach is proposed in this research to calculate the instability index from a time series analysis of the landslide inventory for each optimum sub-watershed. The landslide area ratio time series calculated in each optimum sub-watershed, as well as the annual thumbnail images overlaid with the corresponding landslide polygons, make it possible to clarify the cause–effect relationship and identify hot spots that are prone to collapse. Taking three optimum sub-watersheds Nos. 20713 (Figure 10), 21735 (Figure 11), and 22253 (Figure 12) as an example, we are able to clarify that the main cause of the large scale landslide areas was Typhoon Aere in 2004, resulting in three optimum sub-watersheds with levels of instability of 2, 4, and 5. The classified level of instability reflects the actual situation pretty well. Note that the Qingquan watershed is classified at the B⁺ level according to the existing regulations. The spatial and temporal distribution and relationship of the landslides (Figure 8), however, indicates that most of the areas are actually quite stable, whereas some landslide events occur repeatedly in a few optimum sub-watersheds. Such an assessment of stability for each optimum sub-watershed derived from a time series analysis of the landslide inventory elaborates relevant and reliable spatial information by which to support decision-making.

Up to the present time, the most common approach to assessing the stability of watersheds has been to compute the LSI. Now, more factors are being taken into account, and the LSI models are becoming increasingly more complicated. As aforementioned, however, such an approach is inevitably limited by the data quality of such factors since some factors are not available (such as the detailed geological composition), and some factors are too expensive to be measured on a regular basis (topographically-related factors). Based on the assumption that future landslides are closely related to the historical landslides, the statistic-based instability index proposed in this work is promising and will be even more reliable, provided the database is large enough to include all kinds of situations. The key of success of this work is the long-term and detailed landslide inventory derived from Formosat-2 imagery (2005-2016), which is probably the most intensive observation from space for a country with such frequent occurrences of landslides. It is much easier to prepare a long-term, detailed landslide inventory for any place in the world now using the free and open Sentienl-2 or Landsat-8 imagery.

A landslide occurs when conditions related to both preparatory factors and trigger factors are satisfied. Note that the calculation of the LSI is only based on preparatory factors weighted according to importance. There is no input from the information for trigger factors at all. The instability index proposed in this work, therefore, should not be confused with the LSI. In other words, a region with a higher LSI value might still have a lower instability index value, provided the conditions triggering rainfall are never satisfied in that region. Considering the difficulty of taking all possible rainfall conditions into the calculation of the LSI, the statistic-based instability index proposed in this work is essentially an indicator of joint influence from both the preparatory factor and the trigger factor conditions. For the case of Taiwan, an optimum sub-watershed with a high instability index value can be assessed and mapped as an unstable region, but an optimum sub-watershed with a high LSI value might be rather stable and free of any landslide under the rainfall conditions during the period from 2005 to 2016.

Earthquakes are another trigger factor that is also very important, particularly in a seismically active area like Taiwan. A significant earthquake will immediately trigger a lot of landslides, such as the case of the 2018 Hokkaido Eastern Iburi earthquake. The soil loosened by an earthquake will also be flushed out easily during the rainy season that follows. Therefore, after a major earthquake, the calculated landslide susceptibility should be increased accordingly to the level of the earthquake. However, it is impossible to have a landslide inventory that is comprehensive enough to include all levels of earthquakes for all regions. Neither the LSI approach nor the instability index approach can be used directly after a major earthquake. The current landslide inventory was compiled from 2005 to 2016, after the Chi-Chi earthquake (ML7.3) occurred in central Taiwan with a focal depth of 8 km on September 21, 1999. This earthquake caused large scale multiple landslides and was the strongest earthquake recorded during the 20th century in Taiwan. Although there have still been some minor earthquakes in different parts of Taiwan in recent years, the LSI and the instability index derived from the same landslide inventory can be employed under the assumption that the variations in the contribution from the earthquake factor are insignificant. This assumption would be invalid if there were to be another big earthquake in Taiwan.

One important assumption of this statistic-based instability index approach is that future landslides are closely related to historical landslides. Although a similar assumption was also made by Di Paola et al. [16] and Stamatopoulos et al. [17], this assumption can be reexamined and better supported by taking into consideration local landslide triggers, geological and morphological conditions, and so on. The SWCB of Taiwan has funded another multi-year project to develop a nowcast system for landslide hazards [12] by establishing an LSI map with four preparatory factors, including slope angle, slope aspect, total flux [8], and geological composition. As aforementioned, the detailed geological composition in most parts of the mountainous area is not available, and the topographically-related factors are not measured on a regular basis to reflect the latest status of landslide susceptibility. Therefore, the accuracy and applicability of existing LSI maps are somewhat limited in some regions. Nevertheless, this LSI map can be compared to the results of our instability index, with the intention to clarify the reasons behind any deviations and to evaluate the levels of correlation. This has been planned in a future work.

5. Conclusions

This paper proposes an innovative approach to calculate the instability index by preparing the annual landslide inventory, determining the optimum sub-watershed, compensating for shadow effects on the landslide area ratio time series, and classifying the standard deviations for different levels of instability. Considering the difficulty of taking all possible rainfall conditions into consideration during the LSI calculation, the instability index is essentially an indicator of a joint influence from both preparatory factors and trigger factors. Based on the long-term, detailed landslide inventory derived from Formosat-2 imagery (2005-2016), as well as the assumption that future landslides are closely related to historical landslides, the statistic-based instability index proposed in this work is promising and more reliable than the LSI approach. Both the instability index and the LSI derived from the existing landslide inventory can be employed under the assumption that the variations in the contribution from the earthquake factor are insignificant. This assumption would be invalid if there were to be another big earthquake in Taiwan. Taking the Qingquan watershed as an example, the instability index calculated at 22 sub-watersheds makes it possible to clarify the cause–effect relationship and identify hot spots that are prone to collapse. The spatial and temporal distribution and relationship of the landslides indicate that most of the areas are actually quite stable, whereas some landslide events occur repeatedly in a few optimum sub-watersheds. Such an assessment of stability for each optimum sub-watershed derived from a time series analysis of the landslide inventory elaborates relevant and reliable spatial information to support decision-making. Based on this new approach, the SWCB assesses the stability of 28 watersheds, and the results are made available on the Big Geospatial Information System [40]. Together with the innovative approach to calculate the instability index, it is possible to assess and map the spatiotemporal landslides in the watersheds and achieve land degradation neutrality set up by the UN Sustainable Development Goals.