CHAPTER 5 LANDSLIDE EXTRACTION WITH A GEOMORPHOLOGICAL
5.2 Study Area and Materials ................................................................................ 9 0
5.2 Study Area and Materials
Materials for this research include high resolution satellite images and airborne LiDAR data.
5.2.1 Physiographic Settings of the Study Area
Hsiaolin village is located in Chiahsien District, Kaohsiung City (Figure 5.1). The study area is covered by 9 map-sheets of 1/5000 national photomaps:
95193025~95193027; 95193035~95193037, and 95193045~95193047. The village is located on a river terrace of Chisan River. The geological map in Figure 5.2 (Song et al., 2000) shows that the area is situated in the Western Foothill Zone of Miocene sedimentary formations including Changchikeng Formation, Tangenshan Sandstone, Yenshuikeng Shale, and Peliao Shale. The area is primarily covered by Tangenshan Sandstone and Yenshuikeng Shale.
Tangenshan Sandstone consists of alternate layers of sandstone and shale, whereas Yenshuikeng Shale consists of alternations of siltstone and shale with occasional lens-type conglomerates. The river terrace materials include recent fluvial and colluvial deposits of sand and gravel.
91
Figure 5.1 SPOT image taken on 2009/08/24 after Typhoon Morakot. The 8-digit numbers are the map numbers of national 1/5000 map series.
92
Figure 5.2 A regional geological map near the Hsiaolin village (Song et al., 2000).
5.2.2 Satellite images
This study uses SPOT images taken at approximately the same season as the
93
first LiDAR survey in 2005 used for comparison. The Formosat-2 image taken after Typhoon Morakot was collected and compared with the second LiDAR survey in 2010. In addition, there are several typhoon events from 2007 to 2009.
Therefore, this study also uses SPOT images acquired from 2005 to 2009 (Figure 5.3) to analyze landslide recurrence rate. The resolution of enhanced-mode SPOT images is 2.5 m, pan-sharpened Formosat-2 image have a resolution of 2.0 m.
Figure 5.3 Satellite images of the study area from 2005 to 2009. Bright grey features on the images are mostly landslide scars. Landslide occurrence
increasingly increases in this period of time, as shown in Figure 5.8
94
5.2.3 Airborne LiDAR data
LiDAR data before and after Typhoon Morakot were collected for this study. The LiDAR feature of multiple returns provides a good means for editing the point clouds and produce DSM, DEM, and CHM (Canopy Height Model) or DBM (Digital Building Model). This in turn enables the analysis of multi-temporal datasets. As Figure 5.4 shows, the DEM and DSM in this study are based on 2005 LiDAR survey. The landscape suffered from dramatic changes after Typhoon Morakot (Figure 5.5). The large landslide near Hsiaolin Village is the most conspicuous example. Figure 5.6 shows the DEM and DSM of the study area acquired in 2009 after Typhoon Morakot. Both of the LiDAR datasets in this study were surveyed using a common guideline (MOI, 2006) and a common datum—TWD97 for geodetic coordinates and TWV2001 for vertical system—to maintain the same level of accuracy. The RMSE (Root mean square error) was 16.7 cm with a standard deviation of 16.3 cm for 2005 LiDAR data.
The RMSE was 20.2 cm with a standard deviation of 18.3 cm for 2009 LiDAR data. RMSE is a measure of the dispersion between the coordinates obtained by Airborne LiDAR and those surveyed in the field. Whereas, standard deviation is a measure for the concentration of the differences between these two datasets.
The accuracy of these two datasets meets the requirement set in the MOI guideline (MOI, 2006).
Figure 5.4 DEM and DSM images before Typhoon Morakot
Figure 5.5 3D perspective views of Hsiaolin Village before and after Typhoon Morakot. Hsiaolin Landslide has a volume of ~25 million cubic meters with a maximum depth of 85 m on top area and a maximum length of 3396 m from top
to the other side of Chisan River. The landslide completely destroyed the
95
DEM and DSM images before Typhoon Morakot
3D perspective views of Hsiaolin Village before and after Typhoon Hsiaolin Landslide has a volume of ~25 million cubic meters with a maximum depth of 85 m on top area and a maximum length of 3396 m from top
to the other side of Chisan River. The landslide completely destroyed the village.
DEM and DSM images before Typhoon Morakot
3D perspective views of Hsiaolin Village before and after Typhoon Hsiaolin Landslide has a volume of ~25 million cubic meters with a maximum depth of 85 m on top area and a maximum length of 3396 m from top
to the other side of Chisan River. The landslide completely destroyed the
96
Figure 5.6 DEM and DSM obtained after Typhoon Morakot. As compared to those of Figure 5.4, dramatic landform change can be found in river valley as
well as mountain slopes, especially the example of Hsiaolin Landslide.
5.3 The Geomorphological Model for Landslide Extraction 5.3.1 Introduction
The proposed model includes both global and local detection procedures, and uses a supervised classification method for global landslide detection. The focus of this paper is on global detection. Because of the diversity of the geologic and topographic environments in which landslides occur, omission and commission errors are unavoidable when using the global approach. Thus, local landslide detection is required to increase the accuracy of the resulting landslide
97
map (Liu et al., 2009). With the attendance of geological expert, the local approach employs several interactive manual editing tools to compile landslide information and minimize commission and omission errors. Therefore, the aim of the global detection is to include as much as possible the areas which are vulnerable to landslides. For error analysis, the user accuracy, producer accuracy, average accuracy, and overall accuracy were calculated from a confusion matrix (Kohavi and Provost, 1998).
Figure 5.7 shows the flowchart of the geomorphometric model established in this study. Landslide areas possess geomorphometric characteristics that can be used to establish a geomorphometric model to describe the topographic feature of landslides. As the first step, global parameters based on landslides extracted from satellite images by classifying bare land and then filtering out commission errors produced by bare agriculture lands and debris flows were obtained. Landslide polygons were then overlaid on parametric maps derived from 2005 LiDAR data. The parametric parameters of the extracted samples were then used as training sample globally. Thresholds of various parameters were derived based on statistics of the training samples of landslides. Threshold values of the six geomorphometric parameters (T1~T6) were defined a priori based on some user-defined training areas, that is, the landslide polygons. The mean and standard deviation values of each index were calculated and the threshold values were set to be the mean ± 3 standard deviations. The proposed method classifies a pixel as a landslide pixel if the following expression is true:
98
(Slope > T1)∩(Roughness < T2)∩(Curvature > T3)∩(OHM <
T4)∩(Greenness < T5)∩(Wetness > T6). Otherwise, it is classified as a non-landslide pixel. Because the global landslide detection algorithm is pixel based, isolated landslide pixels were removed by morphological filtering (e.g., opening and closing). Small landslides were eliminated by setting a minimum mapping unit. Finally, the detected landslide pixels were converted into vector-based polygons. In other words, the pixel conforms to the threshold criterion is designated as 1, otherwise it is designated as 0. The area of the intersecting set of all the parameters was categorized as landslide area. This is a dichotomic multi-criteria evaluation approach.
Figure 5.7 Flowchart of the geomorphometric model
slope curvature OHM Topographic wetness index
Interpretation- Landslides before and after the event
Statistics of landslides and Thresholds of selected geomorphometric parameters
Dichotomy
Model prediction of landslide areas
Accuracy evaluation
99
5.3.2 Geomorphometric Parameters of Landslides
For extracting landslides from high accuracy and high resolution LiDAR data, parameters for establishing the model were selected based on the criteria usually used in manual interpretation of landslides, including the 2D and 3D landslide features detailed previously in Table 3.1. The parameters of the geomorphometric model in this study were derived from LiDAR DEM and DSM. The major parameters in this model include slope, surface curvature, OHM (object height model), OHM roughness, and topographic wetness index.
In addition, NDVI (Normalized Difference Vegetation Index) or greenness is one of the most important indexes for landslide recognition due to that fresh shallow-seated landslides are characterized by bare land without or with little vegetation cover. Therefore, it is also included in the model. A number of vegetation indices, such as the NDVI (Jackson et al., 1983), EVI (Enhanced Vegetation Index) (Liu and Huete, 1995), and LAI (Leaf Area Index) (Chen and Black, 1992) have been used in remote sensing for analyzing vegetation cover.
Of these indices, NDVI is the standard method for comparing relative biomass and vegetation greenness in remotely sensed images. A higher NDVI indicates a higher level of healthy vegetation cover. The greenness index is similar to the NDVI, except that it substitutes a green band for the near-infrared band.
These parameters are also closely related to the factors for landslide susceptibility (Tarolli et al., 2011). The control factors of slope stability usually include slope angle, strength of materials, and pore water pressure (Turner and
100
Schuster, 1996). If the slope gradient is high, the slope can be unstable. Slope angle was thus selected as the first parameter because of its importance, and can be easily derived from DEM. Because DEM represents the bare ground surface and DSM represents the upper envelope of all the objects above the bare ground surface, the difference between these two well-defined surfaces is minimal in the area of rainfall-induced landslide. In this case, the OHM, defined as the difference between these two surfaces, can be a good parameter for automatic landslide recognition. After wash out or sliding, the surface of landslides in nature should be smoother than the surroundings. Surface roughness is an objective and useful measurement of landslide topography (Glenn et al.,2006;
Woodcock, 1977; Mckean and Roering, 2004). Landform curvature is another critical factor controlling the susceptibility of landslide occurrence (Pirotti and Tarolli, 2010).
The definition of the parameters is as follows (Wilson and Gallant, 2000;
Zhou and Liu, 2006):
(1) Slope. The slope angle of a landslide is the angle between the horizontal surface and the ground surface of the longitudinal axis of the landslide. The slope angle for each landslide can be derived from LiDAR DEM data. A variety of methods are available for terrain slope gradient estimation. However, the details of a high-resolution terrain model may introduce high variations in changes of local slope gradients (Sharpnack and Akin, 1969). This study adopts the method proposed by Parker (1997) to overcome this problem, that is
101
the derivatives of the Gaussian function are convoluted with the DEM in the x and y directions, respectively, and then combined to estimate the slope.
(2) OHM. Object height models (i.e., OHMs) are obtained by subtracting DSM from DEM to describe the height of objects above bare ground. The OHM describes the heights of above-ground objects in raster format. Objects close to zero in height may represent the bare soil that characterizes landslides.
(3) OHM Roughness. Roughness is a derivative of OHM, defined as one standard deviation in a 5 x 5 moving window. This measure, which is a function of geological structure and lithology, describes the relief variation in the local area. Because most landslides occur in bare soil areas, the surface is smoother than that of forested areas. Thus, a surface roughness index can be used to detect landslide areas. To account for the high terrain variation in mountainous areas, this study uses object heights rather than surface heights. For simplicity, the standard deviation of object heights within a local window serves as the surface roughness index.
(4) Curvature. Curvature is the second derivative of the surface (Schmidt et al., 2003). Two optional output curvature types are possible: the profile curvature is in the direction of the maximum slope, and the plan curvature is perpendicular to the direction of the maximum slope. The curvature is the slope form and has a significant effect on surface runoff, soil erosion, and deposition processes (Stefano et al., 2000). This study applies a 15 x 15 medium filter to the DEM to suppress any accidental height changes in the high resolution elevation model.
102
The curvature along the slope direction was then calculated with a 5 x 5 mask.
(5) Topographic wetness index (TWI). Wetness is derived from the concentration of a small watershed (Kirkby, 1975; Wilson and Gallant, 2000).
Topography is often one of the major controls of the spatial pattern in saturated areas, which in turn is a key to understanding the variability of hydrological processes. The topographic wetness index has become a widely-used tool to describe wetness conditions. The formula is as follows:
tan )
ln( θ
ω = A
(5.1)where A is the local upslope contributing area and θ is local slope.
(6) NDVI or greenness. This parameter is derived from satellite images or orthophotos acquired at a compatible time as the LiDAR survey. In other words, there are no rainfall events between the time that both the LiDAR data and the images or orthophotos are acquired. Because rainfall-induced landslides of natural slopes are mostly covered by densely-vegetated surroundings, the vegetation index is critical for indicating the areas of bareness. The most popular index is the NDVI:
NDVI = (NIR-R)/(NIR+R) (5.2) where R stands for the grey value of the red band and NIR stands for grey value of the near infrared band. Theoretically, if the image digital values are calibrated to stand for the reflectance of the target, the NDVI can be widely
103
applicable. However, the digital numbers of the red band and NIR band of digital aerial cameras are not calibrated for this purpose. Therefore, the NDVI value is a relative indicator of vegetation cover. NDVI can be applied to modern digital aerial cameras, which usually include an NIR band. If color aerial photographs include only RGB bands, an alternative greenness parameter can be used. Greenness is also a relative indicator with radiometric values that are not normalized:
Greenness=(G-R)/(G+R) (5.3) where G is the grey value of the green band, and R is the grey value of the red band. The values of NDVI and Greenness range from -1 to 1. Nevertheless, the range for these values in landslides may change depending on natural weather, terrain conditions and type, and camera sensor settings. A relatively low value implies that the area of the pixel is low vegetated or bare.
5.4 Results and Discussion
5.4.1 Establishing The Geomorphometric Model of Landslides
Bare land has a relatively low reflectance in the infrared region of the electromagnetic spectrum. This feature can be used in unsupervised classification to obtain a preliminary map of landslides. On an interactive screen, manual editing of the results can filter out commission errors such as bare crop fields and debris flows. Figure 5.8(A)-(D) show the distribution of landslides over four different years. Six typhoons affected Taiwan in 2008: Kalmaegi,
104
Fung-wong, Nuri, Sinlaku, Hagupit, and Jangmi. A comparison of the images in 2007 and 2008 reveals more landslides in 2008 (Figure 5.8E). The number of landslides increased substantially after the torrential rainfall of Typhoon Morakot (Figure 5.8F).
The recurrence rate of landslides, defined as the repetitive occurrence of landslides between two different times, was 65% between 2007 and 2005. The recurrent rate was even as high as 95.9% between 2009 and 2008. 64.1% of the landslides in 2008 reappeared in 2009 after Typhoon Morakot. The high recurrence rate between succeeding years shows that landslides happen in similar environmental conditions. To verify the accuracy of the landslides obtained by satellite images, conventional aerial photo-interpretation was conducted. It is shown that the overall accuracy was 92.4% with omission error of 9.2% and commission error of 16.1%.
105
Figure 5.8 Landslide distribution between 2005 and 2009. Landslides on images are high-lighted with yellow polylines. New landslides are in red polylines when comparing images taken in 2005 and 2008 (E) and those in
2008 and 2009, respectively.
5.4.2 Statistics of Geomorphometric Parameters
Figures 5.4 and 5.6 are the primary data of DEM and DSM obtained in 2005 and 2009, respectively. For further understanding the features of landforms, geomorphometric parameters are extracted from these primary datasets. Figure 5.9 shows the distributions of major LiDAR-derived geomorphometric parameters selected for landslide recognition in this study.
Figure 5.9 The distributions of major LiDAR
parameters selected for landslide recognition in this study. The coordinates of the maps are (209810，
Figure 5.10 shows the frequency distribution of geomorphometric parameters based on 2005 landslide data. Fig
these parameters based on 2009 landslide data. The average slope of landslides in 2005 is 31.2 degrees. The surface roughness is generally below 1.5 m, with a cumulative fraction of 90% below 1.5 m (Fig
derived from the difference of DSM and DE
20% and 30% of all the landslide pixels having a value below 0.5 m and 3.3 m,
106
Figure 5.9 The distributions of major LiDAR-derived geomorphometric parameters selected for landslide recognition in this study. The coordinates of
，2566339) and (217609，2557916) for the lower right and upper left, respectively.
.10 shows the frequency distribution of geomorphometric parameters based on 2005 landslide data. Figure 5.11 shows the frequency distribution of sed on 2009 landslide data. The average slope of landslides in 2005 is 31.2 degrees. The surface roughness is generally below 1.5 m, with a cumulative fraction of 90% below 1.5 m (Figure 5.10B). On basis of the OHM derived from the difference of DSM and DEM, the average OHM is 9.1m with 20% and 30% of all the landslide pixels having a value below 0.5 m and 3.3 m,
derived geomorphometric parameters selected for landslide recognition in this study. The coordinates of
2557916) for the lower right
.10 shows the frequency distribution of geomorphometric parameters .11 shows the frequency distribution of sed on 2009 landslide data. The average slope of landslides in 2005 is 31.2 degrees. The surface roughness is generally below 1.5 m, with a .10B). On basis of the OHM M, the average OHM is 9.1m with 20% and 30% of all the landslide pixels having a value below 0.5 m and 3.3 m,
107
respectively. Figure 5.10D is a frequency distribution of OHM. The major fraction of OHM is distributed between 5 m to 20 m. A cumulative fraction is 37% and 92% for OHM under 5 m and 20 m, respectively. Only 8% of OHM exceeds 20 m, indicating commission errors of trees can be as high as 8%.
Figure 5.11 shows the frequency distribution of geomorphometric parameters based on 2009 landslide data obtained from images after Typhoon Morakot. In other words, the training samples of the geomorphometric parameters are obtained from the LiDAR data taken in 2009. The average slope of the landslide areas is 33.8 degrees, with a major range in 25~50 degrees. A cumulative fraction is 25% and 90% for slope under 25 and 50 degrees, respectively. The average roughness is 1.2 m, with 90% less than 1.5 m. The average curvature is -0.008, showing that most of the slope forms are more concave than convex. The OHM ranges from 5~20 m with an average of 9.1 m. Similarly, there are 30% of the landslide pixels having an OHM less than 3.3 m. The average roughness of OHM is 2.6 m, with a standard deviation of 1.2 m.
The frequency distributions of various parameters derived by landslides in 2005 and 2009 show no obvious differences. In both cases, the average slopes fall within the range of 30~50 degrees, with a roughness of 1.1~1.7 m, curvature of -0.04~-0.02, OHM under 17 m, and OHM roughness of 1.5~3.5 m.
Figure 5.10 Frequency distribution of geomorphologic parameters of landslides
Figure 5.11 Frequency distribution of geomorphologic parameters of landslides
108
Frequency distribution of geomorphologic parameters of landslides in 2005.
Frequency distribution of geomorphologic parameters of landslides in 2009.
Frequency distribution of geomorphologic parameters of landslides
Frequency distribution of geomorphologic parameters of landslides
109
When using the landslides in 2008 for training samples, the slope ranges from 25~55 degrees, with an average of 38.2 degrees. As a comparison, the general average slope for 2009 landslides is 33.8 degrees, with an OHM of less than 20 m, roughness less than 1.5 m, and average curvature of -0.018. More concave slope forms were present in 2008 than in 2009. Before the Morakot landslide event, the average OHM was 7.3 m, and the average roughness was 2.4 m with a standard deviation of 1.2 m.
The average slopes of 2008 landslides are higher than those of 2009 landslides. However, the curvature for 2008 is less than that for 2009. There are no obvious differences in OHM and roughness. In 2008, a total of 60% of the landslides have an area of less than 0.5 hectares, whereas the average area of individual landslides in 2009 become larger, with 73% of them possessing an area of less than 1.0 hectare.
The average slopes of 2008 landslides are higher than those of 2009 landslides. However, the curvature for 2008 is less than that for 2009. There are no obvious differences in OHM and roughness. In 2008, a total of 60% of the landslides have an area of less than 0.5 hectares, whereas the average area of individual landslides in 2009 become larger, with 73% of them possessing an area of less than 1.0 hectare.