The Geomorphometry of Rainfall-Induced Landslides in Taiwan Obtained by Airborne Lidar and Digital Photography

1. Introduction

Taiwan has a land area of 36000 m². 26.68% of the land areas are covered by plain region, whereas 27.31% are hilly and 46.01% are mountainous. By official definition for the purpose of land conservation management, hilly lands refer to the area under 100m but with a slope more than 5% or the area between 100m and 1000m. Mountainous lands refer to the area with an altitude above 1000m. Therefore, 73.32% of the areas are under conservation management. The complicated landscape of Taiwan is characterized by small drainage basins, highly fractured rock, high relief, and steep stream gradients. Frequent earthquakes due to the collision of Eurasian Plate and Philippine Sea Plate in eastern Taiwan further loosen the top surface of the land. Rock formations are highly fractured and jointed. Therefore the lands are particularly sensitive to episodic events such as typhoons and earthquakes, and various types of anthropogenic disturbance.

In addition, Taiwan is located in tropical and sub-tropical zones, often suffering from heavy rainfalls, especially in the summer seasons with typhoons. The average annual rainfall of Taiwan is 2500 mm which is about three times the world average. Landslides are easily induced by the heavy rainfall come along with typhoons. These physiographic settings make Taiwan a fragile land, especially vulnerable to rainfall-induced landslides. The consequence is the sedimentation of the reservoirs. And the turbidity of the water in reservoirs becomes a major factor impacting the sustainable operation of water supply reservoirs in Taiwan. Landslides have to be recovered and their hazards have to be mitigated. The necessity of landslide survey is obvious.

Aerial photo interpretation has long been adopted for landslide inventory (Liu et al., 2001). This conventional method is based on visual perception of colour tone and geomorphometric features of landslides on the aerial photographs. Both manual interpretation and automatic recognition of satellite images are also used. Most of the recent automatic classification methods of landslides using images are based on spectral features other than topographic features. Therefore, landslides cannot be correctly recognized. A recent study is to establish an interactive approach with a software interface for assisting visual interpretation of landslides (Lau et al., 2006). Both spectral and spatial parameters are employed for the inputs of the software to assist the interpreter/operator to correctly recognize and delineate landslides. Automatic recognition of landslides solely on basis of spectral information of digital images is efficient in terms of time consumption, whereas the results usually can not meet the requirements for taking engineering measures (Parise, 2001). Nevertheless, manual interpretation is too slow to meet the requirements for emergency response. A hybrid approach is to combine the advantages of automatic processes with manual interpretation. The extraction of gemorphometric parameters from airborne LiDAR data is thus considered for integrating in the interactive interface to assist the interpreter.

Airborne LiDAR is the state-of-the-art technology for efficiently taking high density and high resolution elevation data for a wide area. This feature is also suitable for emergency response or quick assessment of landslide disasters. Hsiao et al. (2005& 2006) shows that the integration of multi-temporal airborne LiDAR and aerial photography can give detailed change information of large-scaled deep-seated landslide as demonstrated by the Jiu-fen-er earthquake landslide. For establishing an interactive interface for assisting visual interpretation of landslides, morphometric parameters derived from LiDAR are required for setting the internal defaults (Lau et al., 2006). In this interface, four primary parameters are selected, namely the greenness, the slope angle, the object height model, and surface roughness. Normalized Vegetation Index (NDVI) is taken for denoting the greenness if colour IR digital aerial photography is applied.

For these purposes, surveys were carried out with airborne LiDAR and digital camera to obtain digital terrain models (DTM) and digital surface models (DSM) of 1m grid and colour orthophotos of 50cm grid. DTM, DSM and orthophotos are georeferenced and transformed into the local coordinate system with Taiwan Datum 1997 (TWD97). Subsequently, the geomorphometric features of the landslides are analyzed. In this study, the geomorphometric characteristics of three selected events will be examined and these will be taken as reference values for setting the defaults in the software interface.

2. Conventional API Approach of Landslides and Its Implication

Rainfall-induced landslides are in majority shallow-seated in the high relief terrains of Taiwan. Techniques of stereoscopic airphoto interpretation have been adopted for landslide inventory in Taiwan since 1973 when an aerial survey team was established under Agricultural Council of the government. Though it is labour intensive, it is believed to be reliable. The core spirit of this approach is the synergy of human perception to include both 2D and 3D features of the target and its environment. Any automated attempt should take this into account. Therefore, geomorphometric features of landslides constitute important ingredients in the automation process.

2.1. Rainfall-induced Landslides in Taiwan

For practical applications in the physiographic environments of Taiwan, the classification scheme of landslides developed by Varnes (1978) is simplified into five major categories, namely rock falls, shallow-seated landslides, deep-seated landslides, dip-slope and wedge slides, and debris flows. Thus, types of landslides can be differentiated by their physical appearance. It is especially useful for practical applications using remotely-sensed images.

		Type of Materials
Bed rock	Engineering Soils
Debris	Soils
Falls		Rock falls	Shallow-seated slide
Topples
Slide	Translational	Dip-slope and wedge slide
Rotational	Deep-seated slide
Flows		(not applicable)	Debris flow	(not applicable)

Table 1.

A simplified classification scheme of landslides applied in Taiwan.

There are 270 events of natural disasters in Taiwan in 50 years from 1958 to 2007 including categories of typhoons (71.1%), flooding (15%), earthquakes (8.5%), torrential rainfalls (2.2%), wind-storms (1.5%), mountain flooding (0.7%), and landslides (0.7%) (NFA, 2008). As shown in Figure 1, the frequency of natural disasters is in a trend of increasing. In total, 89% of the events are concerning with rainfall hazards and 97% of them are directly or indirectly concerning with landslides. Rainfall-landslides become a critical issue in managing natural distasters.

Figure 1.

Statistics of natural disasters in Taiwan from 1958 to 2007

Remote sensing has been an important tool for landslide inventory. The physical appearance of landslides is the basis of the recognition of the boundary and the type of a landslide. However, the displaced materials of a rainfall-induced landslide are usually washed away from steep slopes. It remains only the fresh scars of the rupture surface. The fresh landslide scars emplacing at various slope gradients and various slope locations would normally include landslide types such as rock falls, debris slides, channel bank failures, and debris flows. In this study, the landslides concerned will cover all these types except debris flows. The exception is due to the reasoning that debris flows are triggered by a different mechanism with more contributions from flowing-water instead of gravity itself. In other words, debris flows can be treated as a transformation of other shallow-seated landslides when high concentration of rainfalls and liquefaction of displaced materials take place.

2.2. Procedures of Air Photo Interpretation

Air photo interpretation (API) is a process of understanding to associate shapes and pattern and other characteristics on vertical images with real features or phenomena on the ground. Interpretation by aerial photographs has been the most efficient and realistic way for identifying landslide topography in a wide area. Currently, researches in automatic extraction of landslides using images and digital elevation data become important topics (Barlow et al., 2003, Chang & Liu, 2004, Fernandes et al., 2004, Parise, 2001, Liu et al., 2008, Mantovani et al., 1996). However, visual interpretation by well-trained personnel is still believed to be more accurate and reliable than by computers. Interpretation process needs high skill and the results largely depend on the expertise of the interpreter. Sense of perception of a specific feature such as landslide can be acquired by practices and by an interpretation key describing visual signature characteristics of the object, including size, shape, pattern, tone, association, and texture. To minimize subjective factors of individual interpreters, cross checks should be implemented for a case covering a wide study area such as a few hundreds of aerial photographs. And, map making should be performed very carefully with, not only aerial photographs, but also site investigation.

The procedures of the conventional API adopted for a wide area of landslide inventory usually include steps as follows:

1. Acquisition and preparation of aerial photographs of the study area.

2. . Aerial photograph interpretation (identifying landslide topography) – A stereoscope is used to pick up accurate landslide topography from aerial photographs. The scale of the panchromatic aerial photographs taken by the Aerial Survey Office of Forestry is about 1:20,000. Since 1976, about 20000 aerial photographs are taken every year. Photo index can be used for choosing the particular cloud-free photographs. Landslides with more than 50m in length were identified and their scarp, moving mass, internal structure, and moving direction are drawn with coloured pencils on the paper-printed photographs. A standard legend should be established.

3. Tracing the identified features on the topographic map – Tracing the features of landslides onto the topographic map by comparing identical landforms both on the photographs and the map. An original map of landslides is thus created.

4. (Digitization and drawing the final map – The landslide features are then digitized. Subsequently, landslide scarps and lineament structures are compiled and printed with a backdrop of conventional contour map in a GIS environment. These maps were examined and revised by the researchers.

5. Field check and update the attribute table from field records.

6. Ancillary materials for interpretation.

7. Final presentation and backups.

8. The second step of the API procedures is the most critical one where stereoscope is usually used to perceive the sense of 3D features and a well-trained interpreter should be acquainted with interpretation key for the study area.

2.3. Interpretation Key

The perception of landslides from a bird-eye view of aerial photographs is also largely depending on the scale or spatial resolution of the photographs. Landslides can not be mapped properly when they are smaller than a minimum mapping unit such as 5mm on the paper prints. Before 2008, the aerial photographs taken by Aerial Survey Office had been the conventional panchromatic photographs in a scale around 1:20000. Therefore, the minimum mapping unit of the landslides will be larger than 100m in the real ground. In general, four factors affect the quality of the mapping results, namely the scale, the time lag between the landslide event and the aerial photography, the type of film used, and the overall quality of the photographs. Table 2 shows the criteria used for the recognition of landslides on aerial photographs. The general feature of a rainfall-induced landslide is characterized by the fresh landslide scars in elongated shape and located in a relatively steep slope. It takes place in any kind of geology so long as there are some weathered overburdens. Features on aerial photographs include the bright tone, the bare surface, and the features shown in Table 2. Manual interpretation uses both 2D and 3D features of the landslides for recognition. The 2D features include tone, location, and shape. The 3D features include location, direction, slope, and shadow effects. A sound consideration of the automation of landslide recognition should be able to take care of all these aspects.

Feature	Description	Discrimination rule
Tone	Light, grey light	Brightness"/Threshold
Location	Near ridges, cut-off slopes, road-sides	Trigger events and buffer zone of the feature
Shape	Spoon-shaped, elongated-oval, dentritic, rectangular, triangular	Location-specific and topography-specific
Direction	The drop direction of the landslide is the gravitational vector on the ground surface.	Roughly perpendicular to the streams and t opography-specific
Slope	Depend on types of landslides. E . G . S hallow-seated landslides "/ 45% ; Deep-seated landslides ~40% ; Debris flows ~10-20% .	Slope "/ Threshold
Shadow	Depend on whether the l andslides are in shadow -side or sunny-side	Solar azimuth in related to slope aspect

Table 2.

The criteria for the recognition of rainfall-induced landslides.

2.4. Geomorphometry of Landslides

Obviously, geomorphometry has been applied in manual interpretation. Geomorphometry, the science of quantitative land surface analysis is also known as geomorphological analysis, terrain morphometry, terrain analysis, and land surface analysis (Hengl & Reuter, 2009). The aims of geomorphometry are to extract surface parameters and objects using input digital terrain models. Pike (1988) listed a dozen groups of parameters used as terrain descriptors using manually digitized digital terrain models and he used a resulting "geometric signature or topographic signature" to categorize terrain characteristics and suggested the degree of danger from landslides. Topographic signature of life and their processes are deemed to be strongly influenced by biota (Dietrich & Perron, 2006). Guth (2001& 2003) took terrain fabric as measures of a point property of the digital terrain models and the underlying topographic surface. This study is also known as topographic fingerprints (Densmore & Hovius, 2000) for characterizing the location of a landslide on the slope. The state-of-the-art technology of high resolution satellite images, digital aerial photography, and airborne LiDAR opens a new era in the automation of landslide recognition, especially the possibility of applying geomorphometrics. And, the extraction of land surface parameters becomes more and more attractive for both stochastic and process-based modelling, making use all the level of detailed digital terrain models.

It is shown that the topographic-based analyses can be used to objectively delineate landslide features, generate mechanical inferences about landslide behaviour and evaluate relatively the recent activity of slides (McKean & Roering, 2004, Glen et al., 2006). Especially, surface roughness derived from LiDAR DTM allows an objective measurement of landslide topography. Eigenvalues of surface normals can be an effective parameter for differentiating shallow landslides and debris flows (Woodcock, 1977).

For establishing an interactive interpretation software interface to assist the interpreter, it is clear that expert knowledge of the morphometric properties of landslides is required. And, data acquisition with the new sensors of aerial digital camera and LiDAR becomes feasible. Therefore, the general properties of slope angles, OHM and roughness of rainfall-induced landslides are included in this study.

3. The New Interactive Approach and Parameters of Geomorphometry

Figure 2 shows some typical rainfall-induced landslides in Taiwan. Landslides are bare in high relief terrains with densely-vegetated surroundings. Typical modernized aerial survey system nowadays is equipped with a digital camera and a LiDAR sensor. The procedures of landslide inventory are subjected to change to adopt the new types of high resolution digital data. Thus, an interactive system for manual interpretation under a digital environment is required. Standard products generated by the new survey system include orthophoto, DTM and DSM. In addition to the functions for data management and manipulation in the interactive system, algorithms for automatic recognition of landslides are also required to assist or guide the interpreter for improving the efficiency.

Figure 2.

Typical rainfall-induced landslides in Taiwan,

3.1. The interactive system for color orthophoto interpretation

On basis of the experiences in airphoto interpretation and national landslide inventory, a man-machine interface is developed using windows software development tools including Visual Studio.NET, Borland C++ Builder, and OpenGL. Figure 3 is the flowchart of the interactive system which includes three data entries and four parameters. The entries and parameters will be modifies when more standard products are available. Parameters of roughness, OHM and Slope are derived from LiDAR data. Parameter 4 the greenness is derived from color orthophoto. These four parameters are used for highlighting potential areas of landslides by default settings of threshold for the parameters. Another option is to manually define training areas to obtain the threshold from the training sample.

The visualization on the screen shows both 2D and 3D perspectives of the results (Figure 4). Final setting of parameter thresholds can be optimized visually. And finally, the interpreter can further edit the results of automated detection. Or, the interpreter can even carry out all the interpretation discarding the automated results. Finally, ground truth can be imported to compare with the results for accuracy assessment.

For practical reasons, only four major parameters which can be easily derived from the standard aerial products available by a national agency are used for the automatic back-processing in the interactive system (Figure 3). Simple thresholds are used to highlight the potential landslides. For example, roughness < 5m, OHM < 10m, slope > 40 degrees, and greenness < -0.40. Default settings of thresholds are set on basis of geomorphometric analysis of rainfall-induced landslides for the specific area in related to physiographic conditions and the triggering event. Another option is to obtain the thresholds from the training sample. In this system, a landslide seed is located on the screen by the interpreter. The values of 25 pixels extracted from a 5x5 window centred at the assigned seed are used to calculate statistical means and standard deviations. Three times of the standard deviations are taken as the thresholds. Any pixel with a value within three standard deviations of the means will be assigned as a pixel of landslide. Thus, the omission and commission errors of landslide recognition can be minimized. In addition, the thresholds can be tuned interactively to see the correctness of matching between the landslide feature on the colour orthophoto and the highlighted area (Figure 4).

3.2. The parameter derived from orthophoto

Because rainfall-induced landslides of natural slopes are mostly covered by densely-vegetated surroundings, vegetation index will be critical for indicating the areas of bareness. The most popular one is the NDVI (Normalized Vegetation Index).

$$\text{NDVI }=\left(\text{NIR}-\text{R}\right)/\left(\text{NIR}+\text{R}\right)$$ E1

where R stands for the grey value of Red band, and NIR stands for grey value of Near Infrared band. Theoretically, if the image digital values are calibrated to stand for the reflectance of the target, the NDVI can be widely applicable. However, the digital numbers of Red band and NIR band of digital aerial camera are not calibrated for this purpose. Therefore, the NDVI value is a relative indicator of biomass. NDVI can be applied for recent digital aerial cameras which usually includes an NIR band. If the colour aerial photographs include only RGB bands, an alternative of greenness parameter can be used. Greenness is also a relative indicator, of which the radiometric values are not normalized.

$$\text{Greenness}=\left(\text{G}-\text{R}\right)/\left(\text{G}+\text{R}\right)$$ E2

where G is the grey value of Green band, and R is the grey value of Red band. The range of the values of NDVI and Greenness is between -1 and 1. Nevertheless, the range for those of landslides may change with natural weather, terrain conditions and type and settings of the camera sensor. A relative low value implies that the area of the pixel is low-vegetated or bare.

Figure 3.

Flowchart of the interactive system.

Figure 4.

Screen shots of the interactive system. (Left) Parameter settings; (Right) Accuracy assessment by comparing classified result with ground truth.

3.3. The parameters derived from airborne LiDAR

Three parameters are derived from airborne LiDAR DTM and DSM, namely the slope, the object height model (OHM) and the surface roughness. The factors in the mechanism of slope stability usually include slope angle, strength of materials, and pore water pressure (Turner & Schuster, 1996). If the slope gradient is high, the slope can be unstable. Slope is thus selected as the first parameter due to its importance and that it can be easily derived from DTM. There are two surfaces which can be easily defined by LiDAR-derived data. One is the digital terrain model (DTM) standing for the bare ground surface. The other is the digital surface model (DSM) standing for the upper envelope of all the objects above the bare ground surface. For an area of rainfall-induced landslide, the difference between these two well-defined surfaces can be minimal. Therefore, the OHM defined as the difference of these two surfaces can be a good parameter for automatic landslide recognition. It is straightforward that, due to the wash out or sliding, the surface of landslides in nature should be smoother than their surroundings. Surface roughness has been proved to be an objective and useful measurement of landslide topography (McKean & Roering, 2004, Woodcock, 1977, Glen et al., 2006).

(a) Slope

Slope angle of a landslide is the angle between the horizontal and the ground surface of the longitudinal axis of the landslide. Slope angle for each of the landslides can be determined by the slope angles derived from LiDAR DTM. If the surface of the ground is

$$\text{Z}=\text{f}\left(\text{X},\text{ Y}\right)$$ E3

the slope (in radian) can be defined as.

$$\begin{array}{l}slope={\tan }^{-1}\left(\sqrt{f_x^2+f_y^2}\right)\\ where{f}_x=\frac{\partial Z}{\partial X},f_y=\frac{\partial Z}{\partial Y}\end{array}$$ E4

In common practice, the DTM is stored in grid form. The slope of a grid element such as Z5 in Figure 5 is computed by using a 3x3 moving window.

Figure 5.

Slope calculation by a kernel of 3x3 moving window.

If the fluctuation of local height becomes too large due to the nature of LiDAR data or due to the nature of local relief, the resulted slope angles will be subjected to heavy noises with discontinuities of slope angles. It is therefore necessary to introduce an image processing method to resolve the problem. The first order of differentiation is applied for convolution operation with DTM. In the principle of image processing, a 2D (x, y) convolution is equivalent to two passes of 1D convolution of both (x) and (y). This simplification can be implemented more efficiently (Sharpnack & Akin, 1969, Parker, 1997). For example, formula (5) is a 1D Gaussian function and formula (6) is the first order of its derivative. Therefore, the slope formula in (4) can be implemented by convolution operations in both x and y directions with DTM grid.

$$f(x)=e^{-\frac{x^2}{2{\sigma }^2}}$$ E5

$$f_x=f\text{'}(x)=\left(-\frac{x}{{\sigma }^2}\right)e^{-\left(\frac{x^2}{2{\sigma }^2}\right)}$$ E6

(b) OHM

OHM is obtained by simply subtracting DTM from DSM for describing the height of objects above the bare ground. DTM is also referred to nDSM, i.e. normalized DSM, denoting the significance of the surface is tightly related to DSM. DTM is the bare ground surface excluding all objects above the ground. In forestry land, the difference between DSM and DTM can be referred to CHM (Canopy Height Model), denoting the general heights of the trees. The surface objects especially in forests are generally depleted in areas of landslides. Therefore, a minimal value of OHM can be expected in landslide areas.

(c) Surface Roughness

Surface roughness can be described by either the variance of DSM or OHM in a local window area. In this study, roughness is defined as one standard deviation in a 5x5 moving window on OHM for describing the relief variation in the local area. This can partly diminish the effects of landscape undulation. A 5x5 window is used for extraction the variance of the OHM values in the moving window and then the value of one standard deviation is used to stand for the surface roughness of the central pixel. In the areas of rainfall-induced landslides, the roughness will be lower than other areas due to the depletion of surface materials.

3.4. Scale Effects of Digital Terrain Models

Because slope angle, OHM, and roughness are generated from DTM, they are subject to the change of DTM grid-size. This poses a requirement to understand the possible scale effect due to the change of DTM grid-size for landslide areas (Claessens et al., 2005). A contraction of 1m grid is carried out to obtain grids of 5m, 10m, and 40m for comparison. A pixel on the grid will cover a larger area when the scale is smaller. There are two approaches for the contraction, namely pixel thinning and pixel aggregation. With pixel thinning, every nth pixel is kept. With pixel aggregation, the new pixels represent averages of the n pixels specified by the contracting factor. In Taiwan, DTMs of 5m, 10m and 40m grids are created on bases of photogrammetry. Therefore, pixel aggregation approach is used in this study for its comparability to image matching.

4. Test Areas and Materials

The landslides induced by rainfall events in Shimen, Alishan and Ilan of northern, middle and eastern Taiwan are selected for this study. Figure 6 is the location map of the three study areas and the landslides of these areas due to the relevant events. Surveys were carried out with both sensors of airborne LiDAR and digital camera to DTM and DSM of 1m grid and orthophotos of 50cm grid. DTM, DSM and orthophotos are georeferenced, co-registered and transformed to the local coordinate system with Taiwan Datum 1997 (TWD97) for the analyses of the induced landslides.

Figure 6.

The location map of the study areas and the landslides of these areas.

Aerial surveys were conducted after rainfall events as shown in Table 3. Although the maximum rainfall in the period of Typhoon Longwang in Shimen was as small as 208 mm, this event was the one followed three larger events in three months of the same year, i.e. Haitang (504mm) on July 16, Matsa (818mm) on August 3, and Talim (384mm) on September 1. The event in Alishan was just a concentrated torrential rainfall. On 9th June 2006, the cumulative rainfall had reached 811mm in 24 hours and 1200 mm in 48 hours. Enormous amount of debris flows and slides took places. LiDAR data and aerial photographs were taken right after the event on 22nd June of 2006. There had been no records of heavy rainfall events one year prior to this event. The landslides observed with these datasets can be solely attributed to this rainfall event. Typhoon Kalmaegi on July 17 took place nine month after Typhoon Krosa on October 4 of 2007 in Ilan area. The rainfall took place after a dry and hot summer season. The occurrences of the three selected study areas are different.

Name and size of s tudy area	Date of data acquisition	Rainfall event	Date of the event	Maximum rainfall (mm)
Shimen (48 sq. km)	Jun. 17, 2006	Typhoon Longwang	Sept. 30, 2005	208
Alishan (36 sq. km)	Jun. 22, 2006	Torrential rainfall	Jun. 9, 2006	1200
Ilan (4 sq. km)	Nov. 4, 2008	Typhoon Kalmaegi	Jul. 17, 2008	1100

Table 3.

Rainfall events related to the study areas

The orthophotos were then generated by the aerial photographs taken by direct-georeferencing technique and ortho-rectified by LiDAR DSM without using ground control points. Photography and laser scanning are synchronized. Because airborne LiDAR is equipped with GPS and IMU, an event mark is given when photography system triggers a transistor-transistor logic pulse. Thus, the instantaneous GPS and IMU information can be used to resolve the exterior orientation of the photo frame, i.e. x, y, z, ω, ψ, κ. Subsequently, the true-ortho ground surface model, i.e. LiDAR DSM, is used for the ortho-rectification of the central projected photograph.

Leica ALS50 airborne LiDAR system used in this study is consisted of 2 major parts, i.e. a laser scanning assembly and a Position and Orientation System (POS). The former one is for triggering laser pulses, controlling the range, the swath, the FOV, the scan rate and the pulse rate. These parameters decide how fast we can make a complete coverage of the survey area. The second part is critical to the positioning accuracy.

Point density is an important indicator for the spatial resolution of LiDAR DTM and DSM. An understanding of the forest closure and crown density can be obtained by preliminary inspection of the point-density distribution of point clouds (Means et al., 2000; Naesset, 2002). In Alishan study area, the point density in average is around 2.3 points/m² with ground point density of 0.6 points/m². The upper envelope of the point clouds is interpolated to form DSM of 1m grid, whereas the point clouds that hit the bare ground or that are filtered to eliminate off-ground points are interpolated to form DTM. In other words, DTM denotes the bare ground surface. The accuracy of the DTM and DSM can be varied due to the change of land-cover types and density of vegetation. For assuring the accuracy, ground survey with total stations was carried out for 347 selected sample points. The RMSE is 0.82m, and mean error is 0.73m (Table 4). The error actually is a bias verified in the field check because this is due to the dense low bushes underneath the tree-canopies. This over-estimation of DTM is noteworthy especially for tropical and sub-tropical forest. In general, the accuracy of bare grounds is about 0.15m. Similarly, Shimen and Ilan areas were flown with looser point density of 1.5 points/m² with ground point density of 0.45 points/m².

Figure 7 is an example of a blown up of 1 square km of the Alishan study area. It is clearly shown that the landslide area can be enhanced on the OHM image where the landslide areas are with low OHM values.

Locations	Sample size	Average error (m)	RMSE (m)	Standard Error (m)
Tree base	219	0.70	0.77	0.33
Open Ground	128	0.79	0.90	0.43
Total	347	0. 73	0.82	0.37

Table 4.

Accuracy assessment of the DTM in forest lands.

Figure 7.

Blown-up of a 1x1 km area of Alishan study area.

5. Results and Discussion

5.1. The resultant images of the four parameters

Greenness can be extracted from RBG orthophoto where the landslide area exhibits lower value (Figure 8C). Local slope can be calculated using 3rd finite difference algorithm (Figure 8D). OHM is a normalized height of objects above the bare ground surface. Because terrain effect has been removed, OHM exhibits a good appearance of landslides (Figure 8E). The roughness of landslide area is obviously lower than that of the environment (Figure 8F). In other words, the smoothness of landslide area is obviously higher than that of the environment.

It also can be observed that the shaded-relief image of DSM gives a better contrast between landslides and their environments than that of DTM due to the contribution of the shading effect of the trees and other above-ground objects (Figure8A and B). In addition, The DSM-shaded image in nature is a true ortho-image, possessing the advantage of no occlusion of object shading when compared with orthophoto of the same area (Figure 7D). It is costly to process an orthophoto to a true orthophoto which needs to incorporate the correction of objects along with the terrain correction. Therefore, if airborne LiDAR survey is carried out alone without an integrated digital camera, the DSM-shaded image can be a good surrogate of panchromatic photograph for manual interpretation.

5.2. Results of Manual Interpretation of Landslides

Landslides of the study areas (Figure 6) are obtained by manual interpretation of colour orthophotos of 50 cm grid and DSM-shaded images of 1m grid using the criteria of expert knowledge for conventional aerial photo interpretation.

The total number of the rainfall-induced landslides in the 36 km² in Alishan of middle Taiwan is 106 with a total coverage area around 1.29 km². The landslide occurrence rate is around 4%. Statistically, 8% of the landslides have a longitudinal length of less than 30m; 36% between 30~60m; 67% less than 100m; 86% less than 150m. If more than 5 pixels are the minimum mapping unit for visual interpretation, usually more than 36% of the landslides will not be mapped using remotely-sensed images in medium resolution. The total number of landslides in Shimen of northern Taiwan is 200 with landslide coverage of 0.76 km² in 48 km² of study area. The landslide occurrence rate is around 1.4%, which is only one third of the rate in Alishan although the total number of landslides is more than that in Alishan. This implies that smaller consecutive rainfall events in Shimen area trigger more small landslides than that in Alishan area. This assertion can be further supported by the evidence observed in Ilan area of eastern Taiwan. The total number of landslides in Ilan area is 12 in 2 km² of study area with landslide coverage of 0.14 km². The landslide occurrence rate is around 7.0%. The average area of a landslide in Ilan area is also larger than that of Shimen area, yet comparable with that in Alishan area (Table 5).

Figure 8.

Resultant images.

Study area	Total area (km 2 )	Total landslide area (km 2 )	Totla number of landslides	Landslide occurrence rate (%)	Average area of a landslide (m 2 )
Alishan	36	1.29	106	4.0	122
Shihmen	48	0.76	200	1.4	38
Ilan	2	0.14	12	7.0	117
Average	-	-	-	4.1	92

Table 5.

Statistics of the landslide distribution of the study areas

5.3. Results of the morphometric analyses of Landslides

Manually-interpreted landslides are overlaid with DTM/DSM derivatives to extract the selected geomorphometric parameters including slope angle of landslides, object height models, and surface roughness. Statistics of the landslides in Alishan area (Table 6) show that the mean slope angle of the areas covered by landslides is 40.99 degrees with one standard deviation of 14.14 degrees. In contrast, the mean slope of the whole study area is 33.97 degrees with a standard deviation of 15.71 degrees. Generally, average slope angle in landslide areas is higher than that of the whole area. Figure 9 shows that the peak of the curve of slopes of landslide areas is higher and when the slopes are more than 31 degrees the faction of landslide slopes is more than that of the general slopes. This tendency holds true for both Shimen and Ilan areas.

		Slope (deg)		OHM (m)		Roughness (m)
		Whole area	Slide area	Whole area	Slide area	Whole area	Slide area
Alishan	Mean	33.97	40.99	14.31	4.4 0	3.25	2.05
Std. Dev.	15.71	14.14	9.69	6.3 0	2.69	2.56
Shimen	Mean	35.15	43.79	13.23	2.15 0	2.37	1.48
Std. Dev.	14.28	12.95	8.01	4.7 0	1.87	2.11
Ilan	Mean	29.00	40.48	10.20	6.15	2.55	0.40
Std. Dev.	20.14	13.14	10.81	8.32	2.82	1.32
Average of the means		32.71	41.75	12.58	4.23	2.72	1.31

Table 6.

Statistics of the geomorphommetric parameters of the rainfall-induced landslides.

The mean value of OHM of the landslide areas in Alishan is 4.40 m with one standard deviation of 6.3 m; whereas for the whole study area, they are 14.31 m and 9.69 m, respectively. OHM of landslide areas are obviously smaller than that of the surroundings where are vegetated with high forests (Figure 8E). Figure 9 shows that the distribution of OHM for the whole study area is in bi-modal with one additional peak between 10~31 m. The peak in the right side is a forestry peak representing the concentration of trees. The mean OHM of Shimen area is as small as 2.15m denoting a cleaning ground surface of the sliding areas, whereas the mean OHM of Ilan area is 6.15m denoting the landslide areas remain some tree residues above the ground surface.

Figure 9.

Statistics of the three selected parameters of Alishan area.

The mean roughness of the landslide areas is 2.05 m whereas it is 3.25 m for the whole Alishan area. The cumulative curve of roughness shows that 83% of the landslides have a roughness less than 2m and 88% less than 3m. In general, the means of the landslide areas are less than those of the whole areas. This indicates that ground surface of landslide areas are significantly soother than their surroundings, reflecting the truth of Figure 8(F). The mean surface roughnesses for both of Shimen and Ilan areas are smaller than 2.0m which are even smaller than that of the Alishan area.

The significance of these three morphometric parameters can also be perceived from the average of the means in Table 6 that the differences of the parameters of the whole test area are substantially different from that of the landslide areas.

5.4. Analysis of Scale Effects of Digital Terrain Models

Table 7 shows the statistics of slope angles, OHM, and roughness of landslides in DTM grids of 1 m, 5 m, 10 m, and 40 m, respectively. Two features can be observed in the table: (1) statistics of 40 m grid are obviously different from others; (2) the roughness in four different grids gives quite different values. The former one reflects the unreliability of the statistics when grid-size is comparable to the lengths of landslides (see also Figure 10). The later one shows that there is a significant relationship between surface roughness and grid-size. In other words, there is a scale effect for this parameter. The value of the parameter is changed along with the grid-size. These can be further observed from Figure 10. When the dimension of landslides is similar to or less than the dimension of DTM grid-size, the computed slope angles become unstable, maybe too big or too small. The OHM shows similar phenomena that in 40 m grid, the pixels become mixed cells, i.e. trees nearby the landslide give contribution to the OHM. Surface roughness exhibits changes in all different grid-sizes.

It is noteworthy that there is no cell with a roughness of more than 22m for the curve of 40m grid. 22m is about the half of the 40m-gridsize. This shows that the distribution of roughness is scale-dependent. In short, DTM with a grid size smaller than 40m will not be suitable for analyzing the rainfall-induced landslides which are usually with an area smaller than 40x40m² as demonstrated in this study (Table 5). Therefore, it should be carefully treated when applying DTM with different resolution for geomorphometric studies.

		grid	grid	grid	grid
Slope	Mean	40.99	40.69	40.25	37.77
Std. Dev	14.14	13.77	13.44	13.20
OHM	Mean	4.40	4.86	5.01	6.61
Std. Dev	6.30	5.99	5.85	5.90
Roughness	Mean	2.05	7.06	13.37	33.96
Std. Dev	2.56	4.66	7.39	7.38

Table 7.

Statistics of slope angles, OHM, and roughness of landslides in four grids.

Figure 10.

Scale effects of slope, OHM and roughness derived from various grid-sizes

6. Conclusions

Conventional airphoto interpretation has long been adopted as a standard approach for reliable national mapping of landslides and it is still applied for this purpose in many places of the world including Taiwan. For establishing an interactive interpretation interface to assist the interpreter, expert knowledge of morphometric properties of landslides are required for entries to automatic detection algorithm to highlight the potential areas of landslides in the system. In this study, for understanding these properties, aerial surveys were carried out with airborne LiDAR and digital camera to obtain DTM and DSM of 1m grid and orthophotos of 50cm grid. The landslides induced after torrential rainfalls in middle, northern and eastern Taiwan are selected for this study. It is proved that the morphometric parameters of rainfall-induced landslides are useful in the automatic detection of landslides for highlighting the potential areas in the interactive system. However, they have to be defined in related to local conditions and the specific events triggering the landslides. It is also observed that scale effects are obvious for roughness but not for slope and OHM. The scale effect takes place when the DTM grid is comparable to the average size of landslides, i.e. 40m in this study.

Acknowledgments

This study was sponsored by the grant of Council of Agriculture, Taiwan. Project ID is 95COA-12.1.1-S-a1.