Semi-automatic recognition and mapping of rainfall induced shallow landslides using optical satellite images

Semi-automatic recognition and mapping of rainfall induced shallow landslides using optical satellite images

Remote Sensing of Environment 115 (2011) 1743–1757 Contents lists available at ScienceDirect Remote Sensing of Environment j o u r n a l h o m e p a...

4MB Sizes 4 Downloads 342 Views

Remote Sensing of Environment 115 (2011) 1743–1757

Contents lists available at ScienceDirect

Remote Sensing of Environment j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / r s e

Semi-automatic recognition and mapping of rainfall induced shallow landslides using optical satellite images A.C. Mondini a,b,⁎, F. Guzzetti a, P. Reichenbach a, M. Rossi a, M. Cardinali a, F. Ardizzone a a b

CNR IRPI, via Madonna Alta 126, 06128 Perugia, Italy Dipartimento di Scienze della Terra, Università degli Studi di Perugia, piazza dell'Università 1, 06123 Perugia, Italy

a r t i c l e

i n f o

Article history: Received 11 November 2010 Received in revised form 4 March 2011 Accepted 5 March 2011 Available online 5 April 2011 Keywords: Satellite Optical image Landslide Automatic mapping Statistics Terrain classification

a b s t r a c t We present a method for the semi-automatic recognition and mapping of recent rainfall induced shallow landslides. The method exploits VHR panchromatic and HR multispectral satellite images, and was tested in a 9.4 km2 area in Sicily, Italy, where on 1 October 2009 a high intensity rainfall event caused shallow landslides, soil erosion, and inundation. Pre-event and post-event images of the study area taken by the QuickBird satellite, and information on the location and type of landslides obtained in the field and through the interpretation of post-event aerial photographs, were used to construct and validate a set of terrain classification models. The models classify each image element (pixel) based on the probability that the pixel contains (or does not contain) a new rainfall induced landslide. To construct and validate the models, a procedure in five steps was adopted. First, the pre-event and the post-event images were pan-sharpened, ortho-rectified, co-registered, and corrected for atmospheric disturbance. Next, variables describing changes between the pre-event and the post-event images attributed to landslide occurrence were selected. Next, three classification models were calibrated in a training area using different multivariate statistical techniques. The calibrated models were then applied in a validation area using the same set of independent variables, and the same statistical techniques. Lastly, combined terrain classification models were prepared for the training and the validation areas. The performances of the models were evaluated using four-fold plots and receiver operating characteristic curves. The method proved capable of detecting and mapping the new rainfall induced landslides in the study area. We expect the method to be capable of detecting analogous shallow landslides caused by similar (rainfall) or different (e.g. earthquake) triggers, provided that the event slope failures leave discernable features captured by the post-event satellite images, and that the terrain information and satellite images are of adequate quality. The proposed method can facilitate the rapid production of accurate landslide event-inventory maps, and we expect that it will improve our ability to map landslides consistently over large areas. Application of the method will advance our ability to evaluate landslide hazards, and will foster our understanding of the evolution of landscapes shaped by mass-wasting processes. © 2011 Elsevier Inc. All rights reserved.

1. Introduction A landslide inventory map shows the location and extent of landslides that have left discernable signs in an area (Guzzetti et al., 2000). An event-inventory shows slope failures caused by a single landslide trigger, e.g., an earthquake (Harp & Jibson, 1996; Lin et al., 2004), a rainfall event (Bucknam et al., 2001; Cardinali et al., 2006; Tsai et al., 2010), and a rapid snowmelt event (Cardinali et al., 2001). Event-inventories are important to document the full extent of a landslide disaster (Bucknam et al., 2001; Guzzetti et al., 2004; Tsai et al., 2010), for geomorphological and erosion studies (Guzzetti et al.,

⁎ Corresponding author at: CNR IRPI, via Madonna Alta 126, 06128 Perugia, Italy. E-mail address: [email protected] (A.C. Mondini). 0034-4257/$ – see front matter © 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2011.03.006

2009; Fiorucci et al., 2011), and to prepare and validate landslide susceptibility models (Guzzetti et al., 2006a, 2006b; Rossi et al., 2010). Landslide event-inventory maps are commonly prepared through the interpretation of stereoscopic aerial photographs taken shortly after an event (Harp & Jibson, 1996; Bucknam et al., 2001; Cardinali et al., 2001), by visual analysis of high resolution digital elevation models (DEMs) obtained from airborne Lidar sensors after an event (Ardizzone et al., 2007; Corsini et al., 2007), and through reconnaissance field surveys (Cardinali et al., 2006; Santangelo et al., 2010). Often, a combination of these techniques is used. Investigators have attempted to use satellite images to identify and map landslides using pixel based and object oriented approaches (e.g., Lin et al., 2002; Hervás et al., 2003; Cheng et al., 2004; Rosin & Hervás, 2005; Metternicht et al., 2005; Barlow et al., 2006; Lee & Lee, 2006; Weirich & Blesius, 2007; Martha et al., 2010; Tsai et al., 2010). Most commonly, change detection techniques applied to high-resolution

1744

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

(HR) and very-high-resolution (VHR) optical images are used to identify landslide areas (e.g., Nichol & Wong, 2005; Lee & Lee, 2006; Weirich & Blesius, 2007; Tsai et al., 2010). This requires that pre-event and post-event images are available for the affected area. Investigators have also attempted to exploit HR and VHR monoscopic optical images to detect individual landslides or groups of landslides visually. The technique works where the visual evidence of the slope failures on the image (i.e., the landslide “signature”) is sufficiently clear (Marcelino et al., 2009; Fiorucci et al., 2011), and where geomorphologists experienced in the identification of landslides on aerial and satellite imagery are available. Synthetic Aperture Radar (SAR) images have also been used to detect new (Czuchlewski et al., 2003) and old (Singhroy & Molch, 2004) landslides that have changed considerably the land cover or the surface morphology. In this work, we present the result of an experiment aimed at testing the use of pre-event and post-event VHR panchromatic and HR multispectral satellite images to help detect rainfall induced landslides, and to prepare an accurate landslide event-inventory map. The paper is organized as follows: we first introduce the study area (Section 2), and we briefly describe a rainfall event that has resulted in abundant shallow landslides in NE Sicily, Italy (Section 3); we then present the satellite images available for the study (Section 4), and the processing performed on the satellite images (Section 5); this is followed by a description of the construction and validation of multivariate terrain classification models used to detect and map the event landslides in the study area (Section 6). In Section 7, we discuss the advantages and limitations of the proposed method for the recognition and mapping of rainfall induced shallow landslides exploiting optical satellite images. 2. Study area The study area is located in the Messina province, along the Ioanian coast of Sicily, southern Italy, between 38°03′ and 38°06′ N, and between 15°26′ and 15°30′ E (Fig. 1). Small, sub-parallel channels drain the area, with elevations raging from sea level to more than 1000 m along the divide between the Tyrrhenian Sea and the Ionian Sea. The hydrological regime of the catchments is ephemeral. In the area, marine sediments consisting chiefly of terraced sand, gravel, and silt, Pliocene to Pleistocene in age, cover a tectonized metamorphic bedrock, comprising phyllite, marble, and paragneiss, Palaeozoic in

age. Climate is Mediterranean, with hot and dry summers, and precipitation falling mostly in the period from October to January. Landslides, including shallow soil slides and debris flows, deep-seated rotational and translational slides, and complex and compound failures, are abundant in the area, and are caused primarily by rainfall. 3. Rainfall event and associated shallow landslides On 1 October 2009, a high intensity storm hit the Ionian coast of Sicily, SW of Messina. The Briga rain gauge (Fig. 1), located two kilometers E of the Ionian coast at an elevation of 139 m, measured 200 mm of rain in 6 h. The high intensity rainfall triggered more than 1000 shallow landslides in an area of about 60 km2, corresponding to a density exceeding 16 slope failures per square kilometer. The rainfall induced slope failures were predominately shallow soil slides and debris flows (Fig. 2). The high intensity rainfall also resulted in widespread erosion and deposition of debris along ephemeral drainage channels, inundation, and local modification of the coastline (Fig. 3). Damage was particularly severe in the Giampilieri, Scaletta Zanclea, Guidomandri, Pèzzolo, Altolìa, and Itàla villages, and at several sites along the transportation network. The shallow landslides and the inundation resulted in 37 fatalities, including 31 deaths, six missing persons, and numerous injured people. 4. Materials For the study area, pre-event and post-event aerial photography, a high-resolution DEM, VHR panchromatic and HR multispectral satellite images were available. A landslide inventory map showing old and recent (event) landslides in the study area was also available. Pre-event photography included stereoscopic, black and white aerial photographs taken in 1954, 1995 and 2005, at about 1:33,000 scale. Post-event photography was taken on 5 to 7 October 2009, 4 to 6 days after the rainfall event, when the Italian national Department for Civil Protection flew an airborne Lidar sensor and a digital optical camera to obtain a high-resolution (1 m × 1 m) DEM, and a coverage of pseudo-stereoscopic, color photographs, at approximately 1:3500 scale. A second set of post-event digital stereoscopic photographs was obtained about one month after the event, at approximately 1:4500 scale.

Fig. 1. Shaded relief image of the study area, NE Sicily, Italy. Landslides triggered by the 1 October 2009 high-intensity rainfall event in the Messina Province, Sicily, are shown in purple. Yellow line outlines tanning area (T), red line outlines validation area (V). Colored rectangles show location of landslides portrayed in Fig. 2. Star shows location of the Briga rain gauge.

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

VHR panchromatic and HR multispectral satellite data consisted in two images taken by the QuickBird satellite on 2 September 2006 (Pre-Event in Fig. 3) and on 8 October 2009 (Post-Event in Fig. 3), seven days after the landslide event. The two images were Ortho Ready Standard Imagery level 2A products, corrected for sensor characteristics, radiometry, and geometry (Digital Globe Inc., 2007). The characteristics of the two images are listed in Table 1. The landslide inventory map covered the Giampilieri, Divieto and Racinazzi catchments (Fig. 1), collectively extending for 9.4 km2, 16% of the area affected by the landslide event. The map was prepared at 1:10,000 scale through a combination of: (i) field surveys carried out in the period from October to November 2009 to document the ground effects of the intense rainfall, (ii) visual interpretation of the pre-event and the post-event stereoscopic and pseudo-stereoscopic aerial photographs, (iii) visual interpretation and digital analysis of the image acquired by the QuickBird satellite sensor on 8 October 2009, and (iv) analysis of the high-resolution DEM and derivative products, including slope and curvature maps. The field surveys were used to check locally the accuracy of the landslide mapping obtained through the analysis of the aerial and the satellite images. Production of the landslide event-inventory required: (i) about 15 days/weeks of two geomorphologists for the interpretation of the aerial and satellite images, (ii) about three days of field work of a team of three geomorphologists, and (ii) about five days of two GIS experts for the digitization of the landslide information and the production of a geographical database. Overall, we estimate that production of the event-inventory required about five days per person, per square kilometer. The inventory map portrays the main geomorphological features caused by intense rainfall in the studied catchments, chiefly soil slips and debris flows (Fig. 2). The source, transport, and depositional areas typical of debris flows (Crosta et al., 1990) were difficult – and locally impossible – to map separately. In the inventory, a single polygon is

1745

used to portray a single slope failure or a group of coalescing failures (Fig. 2). The inventory also shows the types and pattern of the preexisting (old and very old) landslides in the catchments (not shown in Fig. 2). The old and very old landslides are mainly deep-seated slides and compound mass movements (Cruden & Varnes, 1996). For these landslides the source and the depositional areas were mapped separately in the inventory. In the study area, the inventory shows 821 shallow landslides caused by the intense rainfall event, for a total landslide area of 0.75 km2, 8% of the area. 5. Processing of satellite imagery 5.1. Initial steps Preliminary image processing consisted in: (i) pan-sharpening, (ii) ortho-rectification and co-registration, (iii) atmospheric correction, and (iv) identification and removal of urban (built-up) areas. For pan-sharpening, the Gram–Schmidt spectral sharpening algorithm available in ENVI (release 4.7.3) was used. The technique preserves in the new “sharpened” image, with a higher resolution of 0.6 m, the spectral information of the original multispectral information, at the lower resolution of 2.4 m (Laben & Brower, 2000; Aiazzi et al., 2009). To evaluate the impact of noise introduced on the individual pixels by the pan-sharpening algorithm, we adopted the approach proposed by Mangolini et al.(1992) and Munechika et al. (1993), and formalized by Wald et al. (1997).The evaluation was done pixel by pixel. First, the original panchromatic and multispectral images were filtered using a band pass filter, and re-sampled by selecting one every four pixels to obtain (synthetic) images that parallel images acquired by an hypothetical QuickBird-like sensor with a 2.4 m ground sampling distance (GSD) for panchromatic, and a 10.6 m GSD for multispectral. Next, the pan-sharpening Gram– Schmidt algorithm was applied to the synthetic images, obtaining

Fig. 2. Shallow landslides triggered by the 1 October 2009 rainfall event in the study area. Yellow lines show landslide boundaries. a, landslides along the slopes N of the Giampilieri village. b, panoramic view of landslides along the slopes S of the Giampilieri village. c, landslides along the slopes N of the Molino village. Location of individual photographs is shown in Fig. 1.

1746

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

Fig. 3. Pan-sharpened images of the study area taken by the QuickBird satellite on 2 September 2006 (Pre-Event), and on 8 October 2009 (Post-Event). Yellow line outlines tanning area (T), red line outlines validation area (V). Black rectangles show location of enlargements portraying the outlet of the Giampilieri torrent in the Ionian Sea before (upper image) and after (lower image) the rainfall event. The abundant sediments transported by the torrent moved the coastline seawards.

pan-sharpened multispectral images at 2.4 m GSD directly comparable to the original multispectral images acquired by the (real) QuickBird sensor at the same GSD. The approach is used because multispectral images with a 0.6 m GSD are not available for a direct comparison with the (real) pan-sharpened images (at 0.6 m GSD). The pan-sharpened images were ortho-rectified using the satellite Rationale Polynomial Coefficients (RPC) and the 1 m × 1 m DEM available for the study area. For ortho-rectification, two sets of ground control points (GCP) were selected visually on the pan-sharpened images, including 20 GCP on the pre-event image and 15 GCP on the post-event image. The reduced number of GCP used to ortho-rectify the post-event image is a consequence of the lower viewing angle of the image (2.4° off-nadir), when compared to the pre-event image (12.9° off-nadir) (Table 1). To improve the geographical matching between the two images, the pre-event image was co-registered to the post-event image using eight homologous points, a first order transformation, and a nearest neighbor re-sampling. Overall, the accuracy of the ortho-rectification, determined selecting 15 additional

Table 1 Characteristics of the QuickBird images used for recognizing and mapping the event rainfall induced landslides (Fig. 3).

Satellite operating altitude Panchomatic GSD Multispectral GSD Product type Acquisition date Off-nadir viewing angle

Pre-event image

Post-event image

450 km

450 km

0.6 m 2.4 m Ortho ready standard imagery, level 2A 26/09/2006 10:06:20 12.9°

0.6 m 2.4 m Ortho ready standard imagery, level 2A 08/10/2009 09:53:31 2.4°

GCP not used for the ortho-rectification, was 0.8 pixel, corresponding to a mean error of about 0.5 m. This positional error is significantly smaller than the error commonly accepted for landslide mapping at 1:10,000 scale (Guzzetti et al., 2002; Malamud et al., 2004; Santangelo et al., 2010; Fiorucci et al., 2011). Two further adjustments were made to the pan-sharpened, orthorectified, and co-registered images. First, a relative radiometric correction was performed to reduce the effect of different atmospheric and illumination conditions (El Hajj et al., 2008), a result of images taken in different periods and with different viewing angles (Table 1). Next, urban areas mapped visually on the pan-sharpened, ortho-rectified, and co-registered images (less than 3% of the study area) were removed to reduce noise introduced chiefly by the different viewing angles, whose effect was most severe in built-up areas. To perform the atmospheric correction, a dark object subtraction was performed to minimize the effects of path radiance (Hadjimitsis et al., 2004, 2010), and the brightness of the “slave” (post-event) image was adjusted to the brightness of the “master” (pre-event) image through radiometric normalization (Dermanis & Biagi, 2006). For the purpose, the following normalization criterion was adopted (Hong & Zhang, 2008): N = ðS−μs Þ ×

σM + μM σs

ð1Þ

where, for each pixel in the image, S is the radiance of the “slave” image, N is the radiance of the corrected (new) image normalized to the “master” image, μ and σ are the mean and the variance of the radiance values, with subscripts M and S indicating the “master” and the “slave” images, respectively. Using Eq. (1), the probability distribution of the radiance value of the “master” and of the corrected

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

“slave” images have the same mean (μ) and variance (σ). Selection of this normalization criterion was based on the fact that it does not require a specific operator experience in the selection of pseudoinvariant points (Hong & Zhang, 2008), and it can be easily automated.

study, changes in NDVI (δNDVI in Fig. 4), were detected using the following equation:  δNDVI = NDVIpre −NDVIpost =

5.2. Change detection Four change detection techniques (Jianya et al., 2008) were selected to identify modifications (changes) between the pre-event and the post-event images (Fig. 3) attributed chiefly to the occurrence of the new rainfall induced landslides. The first two techniques are image difference techniques (Richards & Jia, 2006) that consist in (i) measuring changes in the Normalized Differential Vegetation Index, NDVI (Hayes & Sader, 2001), and (ii) computing the images spectral angle, θ (Klinker et al., 1992; Masek & Sun, 2004). NDVI, first developed to measure the health of the vegetation and to estimate green biomasses (Hayes & Sader, 2001), is commonly used to detect new landslides (e.g., Nichol & Wong, 2005; Lee & Lee, 2006; Weirich & Blesius, 2007; Tsai et al., 2010). In our

1747

   NIR−red NIR−red − ð2Þ NIR + red pre NIR + red post

where, NIR is the radiance (ρ) in the near-infrared band (0.760– 0.900 μm), red is the radiance in the red band (0.630–0.690 μm), and the pre and post subscripts identify the pre-event and the post-event images (Fig. 3), respectively. Determination of the image spectral angle (Klinker et al., 1992; Masek & Sun, 2004) is an alternative (and complementary) technique to determine differences between two images, on a pixel-by-pixel base. In the multispectral domain (blue, green, red, and near-infrared), the image spectral angle θ is given by the following general equation:  → → u×v θ = arccos juj⋅jvj

ð3Þ

Fig. 4. Changes between the Pre-Event (2 September 2006) and the Post-Event (8 October 2009) images (see Fig. 3). δNDVI, changes in the Normalized Differential Vegetation Index. θ, images spectral angle. PC-2, second principal component. IC-4, fourth independent component. COMP, color composite of δNDVI, θ, and IC-2. Inventory, map of the rainfall induced landslides in the study area.

1748

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

 → → where, u = (ρblue, ρgreen, ρred, ρNIR)u and v = ρblue ; ρgreen ; ρred ; ρNIR v are vectors representing the radiance values of individual pixels in the same geographical position in the two images (Masek & Sun, 2004). In our case, the spectral angle θ measures the angular difference between vectors representing individual pixels in the same geographical location in the pre-event and the post-event images, a quantitative indication of local changes between the two images (θ in Fig. 4). We hypothesize that major changes in δNDVI and θ are due to significant variations in the vegetation or soil cover, and specifically alterations produced by the rainfall induced landslides that have removed the vegetation cover locally or have altered the soil surface, completely or partly (Fig. 2). We further hypothesize that minor changes in δNDVI are related primarily to periodic plant cycles (phenology). The other two change detection techniques are based on geometrical transformations, and include Principal Component Analysis (PCA, Jolliffe, 1986), and Independent Component Analysis (ICA, Comon, 1994; Stone, 2005; Hyvärinen & Oja, 2000). In image processing, PCA and ICA are similar computational methods used to transform (separate) multivariate information represented by (possibly) correlated variables into a new set of additive variables, called principal (for PCA) or independent (for ICA) components. PCA exploits correlation criteria. In PCA, the first component accounts for as much of the variability in the data as possible, and the following components account for as much of the remaining variability as possible.ICA uses independent criteria. In ICA, the independent components are (transformed) variables that represent the original (measured) variables. ICA requires non-Gaussian distributions, and PCA works better with Gaussian data. However, the two techniques are not mutually exclusive, and can coexist for image analysis. The two different geometrical transformations highlight different properties of the images, which cannot be simply related to physical processes. With this respect, they are different from image difference techniques, including changes in NDVI (δNDVI) and the image spectral angle (θ). Due to conceptual and operational differences, when applied to the same dataset, PCA and ICA result in different sets of principal or independent components (variables). We have exploited this opportunity to obtain two new (transformed) variables (one principal component and one independent component) to help detect and map the rainfall induced landslides. The principal and the independent components were extracted from a single stack of images containing, for the entire study area, the pre-event and the post event red (0.630– 0.690 μm) and NIR (0.760–0.900 μm) bands, for a total of four layers. The blue (0.450–0.520 μm) and the green (0.520–0.600 μm) bands were not included in the stack because they were strongly correlated with the corresponding red bands (correlation coefficient, ε N 0.97).Visual inspection of the first principal component (PC-1) revealed that the component was correlated strongly with the presence (or absence) of shadows, and was of little help for landslide detection. Inspection of the second principal component (PC-2) revealed a good visual match with the presence (or absence) of the recent rainfall induced landslides (PC-2 in Fig. 4). Correlation of PC-2 to δNDVI (ε = −0.69) and θ (ε = 0.45) revealed that the second principal component contained information not given by δNDVI or θ, and potentially useful to detect and map the rainfall induced landslides. Similar criteria were adopted to select the best independent component, from a set of four. Visual inspection of the geographical distribution of the four independent components revealed that the fourth component (IC-4) matched closely the distribution of the rainfall induced landslides (IC-4 in Fig. 4). IC-4 was (anti-)correlated to δNDVI (ε = − 0.63), and marginally correlated to the spectral angle θ, (ε = 0.18). We concluded that IC-4 provided information useful to detect and map the rainfall induced landslides, and different from the information given by δNDVI and θ. Further inspection of the correlation

between PC-2 and IC-4 revealed a weak correlation (ε = −0.29). We conclude that the two variables PC-2 and IC-4 provide complementary information for the detection and mapping of the new rainfall induced landslides, in our study area. 6. Landslide detection and mapping Recognition of the recent rainfall induced landslides was performed using a multivariate terrain classification method that involved three steps. First, the study area was divided in a training area and in a validation area (Figs. 1 and 3). Second, the available satellite and landslide information was used to calibrate three multivariate classification models in the training area. Third, the classification models were applied in the validation area, using the same set of explanatory variables. 6.1. Preliminary steps The study area was partitioned in two sub-areas: a northern (training) area covering 6.86 km2 (72.9%) (T in Figs. 1 and 3), and a southern validation area, covering 2.55 km2 (27.1%) (V in Figs. 1 and 3). The two areas have similar – but not identical – geomorphological and physiographic characteristics (Table 2). In the training area the landslide inventory shows 382 recent landslides, for a total landslide area ALT = 0.51 km2, 7.4% of the training area. In the validation area the inventory portrays 439 recent landslides, for a total landslide area ALT = 0.24 km2, 9.4% of the validation area. Next, a terrain mapping unit for the multivariate classification was selected. A terrain mapping unit is a portion of the land surface that contains a set of ground conditions that differ from the adjacent units across definable boundaries (Burrough & McDonnell, 1998; Guzzetti et al., 1999; Luckman et al., 1999). For this study, “grid cells” were selected as the mapping unit of reference. For simplicity, 0.6 m × 0.6 m grid cells were selected, geographically coherent with the individual pixels of the pan-sharpened, ortho-rectified, and co-registered satellite images. As a result, the training area consisted of 18,398,631 grid cells, including 1,681,969 (9.1%) landslide cells and

Table 2 Main characteristics of the training (T) and the validation (V) areas (Figs. 1 and 3). Statistics for elevation and slope obtained from the 1 m × 1 m Lidar DEM.

Total extent Urban area Elevation Min Max Mean St.dev. Slope Min Max Mean St.dev. Main rock types Landslides Number Total area Per cent area Landslide cells

Training area

Validation area

km % km2 %

6.86 72.9 0.16 2.3

2.55 27.1 0.10 3.9

m m m m

0.6 934.4 386.9 196.3

1.4 790.25 356.4 182.9

° ° ° °

0 83.9 30.7 13.3 Metamorphic, sedimentary

0 85.9 32.2 13.2 Metamorphic, sedimentary

382 0.51 7.4 1,681,969 9.1 16,716,662 90.9

439 0.24 9.4 1,424,814 20.9 5,378,594 79.1

2

# # km2 % # % Non landslide cells # %

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

16,716,662 (90.9%) non-landslide cells, and the validation area consisted of 6,803,408 grid cells, with 1,424,814 (20.9%) landslide cells and 5,378,594 (79.1%) non-landslide cells (Table 2). A grid cell was classified as a “landslide cell” if it contained 50% or more landslide area.

6.2. Training of the terrain classification models To train a terrain classification model using a multivariate statistical technique, a single grouping (dependent) variable must be identified, and a set of explanatory (independent) variables must be selected (Michie et al., 1994). For our experiment, the presence (or absence) of rainfall induced landslides in the grid cells in the training area was selected as the grouping variable, and the four variables obtained by processing the pre-event and the post-event satellite images (i.e., δNDVI, θ, PC-2, and IC-4, Fig. 4) were selected as the explanatory variables. Model calibration was performed using standard multivariate classification techniques (Michie et al., 1994; Brown, 1998; Venables & Ripley, 2002), including: (i) logistic regression (LR), (ii) linear discriminant analysis (LDA), and (iii) quadratic discriminant analysis (QDA). The adopted multivariate classification techniques (LDA, QDA, and LR) assume that the dependent (grouping) variable contains the same (or a similar) number of cases in the two dichotomous groups; in our case: (i) the group of stable (no landslide) cells, and (ii) the group of unstable (landslide) cells. In the training area, the proportion of landslide cells (9.1%) is significantly lower than the proportion of stable cells (90.9%), violating the assumption. To overcome the problem, we selected randomly 500,000 cells (from a population of 18,398,631 i.e., 2.5%), including 250,000 landslide cells and 250,000 landslide-free (stable) cells. To evaluate the sensitivity of the classification results to the random selection, we repeated the random selection multiple times, obtaining very similar classification results. Results of the classification models (LDA, QDA and LR) are summarized in Figs. 5 and 6. Fig. 5 portrays (from left to right): (i) maps showing the predicted terrain classification (i.e., landslide/ no-landslide), in five unequally spaced probability classes – [0.0– 0.2], [0.2–0.45], [0.45–0.55], [0.55–0.8], [0.8–1.0] –, (ii) histograms showing the number of grid cells in the same five probability classes, and (iii) receiver operating characteristic (ROC) plots (Green & Swets, 1966; Mason & Graham, 2002) showing the “hit rate” (y-axis) vs. the “false alarm rate” (x-axis) (Fawcett, 2006), and the value for the area under the ROC curve (A ROC ), a quantitative measure of the model performance (Mason & Graham, 2002; Jollifee & Stephenson, 2003; Fawcett, 2006). Fig. 6 shows, for each model: (i) a four-fold plot giving the number and percentage of true positive (TP), true negative (TN), false positive (FP), and false negative (FN), and (ii) a map showing the geographical distribution of the same four cases (TP, TN, FP, and FN). Inspection of Figs. 5 and 6 allows for the evaluation and comparison of the fitting performance of the three classification models i.e., the ability of the individual models to reproduce the known geographical distribution of the recent rainfall induced landslides and of the stable slopes, in the training area (Figs. 1 and 3). Analysis of the ROC plots reveals that all three models performed well, with values of the area under the ROC curve, AROC N 0.860. The LDA and the LR models obtained the best predictive performance (AROC = 0.873), and were (marginally) superior to the QDA model (AROC = 0.864). The proportion of grid cells correctly classified by the three models (TP + TN) was also high (N83%), and similar for the three models (LDA = 86.3%, QDA = 85.3%, and LR = 83.1%). Further inspection of the four-fold plots (Fig. 6) reveals that the classification models were more efficient in predicting the stable slopes than the unstable (landslide) areas. With this respect, the LR model was the

1749

most efficient, predicting correctly 75.3% of the grid cells classified as “landslide cells” in the event-inventory map (Fig. 1). Adopting the scheme proposed by Rossi et al. (2010), the three terrain classifications (LDA, QDA, and LR) were combined using a regression-based approach to obtain an “optimal” terrain classification (combined model, CM). In this approach, the optimal model response is taken as the dependent variable, and the results of the three single classification models (LR, LDA, and QDA) are the explanatory variables. Logistic regression was used to classify each grid cell in the training area, based on the explanatory variables. Result of the combination of the three single models (Figs. 5j–l; 6g,h) indicates that the combined model is marginally superior to the “best” of the single models, i.e. the LR and LDA models. Specifically, the area under the ROC curve is slightly larger (AROC = 0.877). Inspection of the four-fold plots for the combined model (Fig. 6g) confirms that the CM model predicts significantly better the stable areas than the landslide areas. 6.3. Validation of the terrain classification models The three single terrain classification models (LDA, QDA, and LR) calibrated in the training area (T in Figs. 1 and 3) were applied to the validation area (V in Figs. 1 and 3), without modifications or adjustments to the individual model parameters. In other words, the three regression models calibrated in the training area were applied to the same four explanatory variables (δNDVI, θ, PC-2, and IC-4) in the validation area, to obtain estimates of the geographical distribution of the landslide areas and the stable slopes, in the validation area (V in Fig. 3). The distribution of landslides and stable slopes predicted by the three classification models were then compared to the known distribution of rainfall induced landslides in the validation area, to determine the predictive skills of the three terrain classification models. Results are shown in Figs. 7 and 8. Inspection of Fig. 7 provides insights on the ability of the three models to predict the (unknown to the models) geographical distribution of the rainfall induced landslides and of the stable slopes in the validation area. Analysis of the ROC plots reveals that all three models performed well in the validation area, with values of the area and the ROC curve, AROC N 0.800. This performance is only marginally lower than the performance achieved by the models in the training area. The LDA model exhibited the best performance (AROC = 0.826), and was marginally superior to the LR model (AROC = 0.824) and the QDA model (AROC = 0.806). The proportion of grid cells predicted correctly by the three models (TP + TN) was also high (LR = 83.2%, LDA = 84.5%, and QDA = 83.4%), but slightly lower than the performance obtained in the training area. Further inspection of the fourfold plots (Fig. 8a,c,e,g) confirms that the models were more efficient in predicting the stable slopes than the landslide areas, in the validation area. Finally, an “optimal” terrain classification was produced for the validation area (V in Figs. 1 and 3) using the outputs of the three classifications for the validation area (Fig. 7a,d,g). The result indicates that the combined model for the validation area is marginally superior to the performance of the “best” of the single models, i.e. the LR model. Explicitly, the area under the ROC curve for the combined model CM is AROC = 0.826. Inspection of the four-fold plot for the combined model (Fig. 8 g) confirms that the CM model predicts significantly better the stable slopes than the landslide areas. 7. Discussion Recognition and mapping of landslides triggered by a single natural event (e.g., a rainfall event, an earthquake, and a rapid snowmelt event) is an important but difficult, uncertain, time consuming, and resources demanding task (Guzzetti et al., 2000;

1750

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

Fig. 5. Results of terrain classification models in the training area (T in Figs. 1 and 3). Rows show classifications produced (from top to bottom) through: (a–c) linear discriminant analysis (LDA), (d–f) quadratic discriminant analysis (QDA), (g–i) logistic regression (LR), and (j–l) through a combined (optimal) regression model (CM). Columns portray (from left to right) (a, d, g, and j) maps portraying the terrain classifications, (b, e, h, and k) histograms of the count of grid cells (pixels) in five unequally spaced classes, and (c, f, i, and l) ROC plots showing the “hit rate” (y-axis) vs. the “false alarm rate” (x-axis), and AROC, the area under the red ROC curve.

Malamud et al., 2004). We have tested the possibility of using optical (VHR panchromatic and HR multispectral) satellite images, and derivative products, to help the semi-automatic recognition of rainfall

induced shallow landslides. The method involves multiple steps, all of which have limitations that should be considered when applying the method and when using the results. In this section, we first discuss the

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

1751

Fig. 6. Results of terrain classification models in the training area (T in Figs. 1 and 3). Rows show classifications produced (from top to bottom) through: (a and b) linear discriminant analysis (LDA), (c and d) quadratic discriminant analysis (QDA), (e and f) logistic regression (LR), and (g and h) through a combined (optimal) regression model (CM). Columns portray (from left to right) (a, c, e and g) four-fold plots summarizing the number and percentage of true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN), and (b, d, f and h) maps of the distribution of the four categories shown in the four-fold plots.

problems associated with the different steps in the proposed method. Next, we discuss the results obtained, in view of their specific and general relevance for landslide inventory making. The first phase of the analysis consisted in the initial processing of the two satellite images, including pan-sharpening, ortho-rectification, co-registration, atmospheric correction, and identification and removal of urban areas. All these steps involved judgment and decisions that have affected the final result, including the landslide terrain classification. In the proposed method, great care is necessary when performing the image pre-processing for the optimal quality of the landslide recognition.

To detect the new rainfall induced landslides, we selected to “fuse” (pan-sharp) the panchromatic and the multispectral information to obtain synthetic pre-event and post-event images (Fig. 3). Pansharpening was performed in the attempt to exploit jointly the (complementary) radiometric and geographical information captured by the satellite sensor i.e., the improved geographical resolution of the panchromatic information, and the better radiometric information provided by the multi-spectral image. We maintain that the joint analysis of the panchromatic and multispectral information facilitated the identification of the new landslides, and particularly of the smallest slope failures. Due to the reduced GSD (0.6 m), use of the

1752

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

Fig. 7. Results of terrain classification models in the validation area (V in Figs. 1 and 3). Rows show classifications produced (from top to bottom) through: (a–c) linear discriminant analysis (LDA), (d–f) quadratic discriminant analysis (QDA), (g–i) logistic regression (LR), and (j–l) through a combined (optimal) regression model (CM). Columns portray (from left to right) (a, d, g, and j) maps portraying the terrain classifications, (b, e, h, and k) histograms of the count of grid cells (pixels) in five unequally spaced classes, and (c, f, i, and l) ROC plots showing the “hit rate” (y-axis) vs. the “false alarm rate” (x-axis), and AROC, the area under the red ROC curve.

pan-sharpened images resulted in an improved identification of the boundary of the individual landslides, with an accuracy greater than commonly accepted for landslide inventory maps at 1:10,000 scale

(Malamud et al., 2004; Fiorucci et al., 2011; Santangelo et al., 2010). For pan-sharpening, adoption of a “radiometric invariant” technique was essential. Use of non-radiometrically-correct images would have

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

1753

Fig. 8. Results of terrain classification models in the validation area (V in Figs. 1 and 3). Rows show classifications produced (from top to bottom) through: (a and b) linear discriminant analysis (LDA), (c and d) quadratic discriminant analysis (QDA), (e and f) logistic regression (LR), and (g and h) through a combined (optimal) regression model (CM). Columns portray (from left to right) (a, c, e, and g) four-fold plots summarizing the number and percentage of true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN), and (b, d, f, and h) maps of the distribution of the four categories shown in the four-fold plots.

introduced unnecessary noise in the terrain classification. Also, use of pan-sharpened images has facilitated the (i) accurate identification of the GCPs required for ortho-rectification, (ii) co-registration, and (iii) visual mapping of the urban areas. By processing the pan-sharpened images, four explanatory variables were selected (δNDVI, θ, PC-2, and IC-4, Fig. 4) for constructing the multivariate classification models used to detect the landside locations. Selection of the explanatory variables, and particularly of the principal (PC-2) and the independent (IC-4) components, was driven by the analysis of the (linear) correlations between the variables, and by the visual evaluation of the matching

(or mismatching) between the geographical distribution of each variable and the geographical distribution of the rainfall induced landslides shown in the event-inventory map (Fig. 1). Selection of different relevant explanatory variables, or of a different number of relevant variables, was possible. However, we maintain that this would have resulted in only slightly different results. When applied to a different area, the procedure may select different principal and independent components. For PCA, the first component is likely to be related to the dominant land cover type. If new landslides have occurred in most land cover types, and if the slope failures have left distinct signatures in all land cover types, it is likely that the second

1754

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

principal component (PC-2) will be related to the presence of the new landslides. Similarly, for ICA we expect that one of the independent components will be related to the presence (or absence) of the new landslides. In all cases, the heuristic step of variable selection will be difficult to automate. We acknowledge that the pan-sharpening process has introduced noise in the “fused” images. The problem consists in determining if the “signal” in the pan-sharpened images prevails over the “noise”. To determine this, for each pixel (2.4 m GSD) and for each of the four change detection indices (NDVI, θ, PC-2, and IC-4), we computed the difference between the original and the fused images (i.e., pixel-bypixel, original vs. fused pre-event images, and original vs. fused postevent images). On average, the differences were two orders of magnitude lower (and at least one order of magnitude lower) than the differences between stable and landslide areas, in the original images. We conclude that the noise introduced by pan-sharpening is (at least) one order of magnitude lower (and on average two orders of magnitude lower) than the “signal” in the images. We further assumed that the image spectral response was not altered by changes in geographical resolution (from 2.4 m to 0.6 m GSD) (Aiazzi et al., 2007), and we did not evaluate the pan-sharpened images at full resolution (0.6 m GSD). To perform the multivariate terrain classification only the information obtained from the two available satellite images was exploited. We maintain that this was reasonable because the scope of the experiment was to test the possibility of using optical satellite data for the detection and accurate mapping of the new event landslides. However, improved terrain classifications may be obtained by exploiting additional information on the location of the new shallow landslides derived e.g., by processing the DEM of the study area (Passalacqua et al., 2010; Tarolli et al., 2010). When measuring the performance of the terrain classification models (Figs. 5–8) and, more generally, of the proposed method for the semi-automatic recognition and mapping of the recent landslides based on the treatment of satellite imagery, the assumption is made that the reference event-inventory (Fig. 1) is complete and correct i.e., that the inventory shows all the recent rainfall induced landslides in the area, and that the shown landslides are of the correct size and in the exact location. The assumption cannot be verified, because the “true” distribution of the rainfall induced landslides is unknown (Malamud et al., 2004; Fiorucci et al., 2011). Thus, our measures of the terrain classification performance are relative (Rossi et al., 2010). Local inspection of the inventory in poorly illuminated areas (e.g., north facing slopes) or in areas covered by forest revealed that some of the grid cells classified as landslides by the statistical models, and where no landslides were shown in the inventory, could contain recent landslides. In these areas, due to the poor lighting conditions and the presence of a thick vegetation cover, shallow soil slips and channeled debris flows were probably not recognized in the postevent aerial photographs. These areas are remote and difficult to access, and field surveys were not performed. Visual analysis of the location and geographical distribution of grid cells classified as false positives (FP, i.e., grid cells where new landslides were not shown in the event-inventory and that should contain slope failures given the local values of the four explanatory variables) on the large scale aerial photographs and in the field, revealed that many of the FP grid cells were located in areas affected by soil erosion, a consequence of high intensity rainfall and associated runoff. In these areas, soil was removed by the action of water, but the geomorphologists did not recognize individual landslides characterized by distinct mass transport phenomena. In places, the distinction between surface erosion and shallow landslides is faint and difficult to detect from the visual analysis of the aerial photographs, or in the field. Analysis of the terrain classifications revealed clusters of grid cells classified as stable by the models surrounded by unstable cells. Visual

inspection of these areas revealed patches of terrain located inside landslide areas where the soil cover was not removed during (or following) failure. In these areas, landslides did occur during the rainfall event, but the land cover did not change sufficiently for the statistical models to classify the terrain as a landslide. Further analysis of the terrain classifications revealed a number of individual grid cells, or groups of cells with total area not exceeding 5 m2, classified as FP (stable areas recognized as landslides by the classification models) surrounded by TN (stable areas recognized as stable terrain by the classification models). Similarly, the analysis revealed a number of individual grid cells, or groups of cells with total area not exceeding 5 m2, classified as FN (landslide areas recognized as stable cells by the classification models) surrounded by TP (landslide areas recognized as unstable terrain by the classification models). Given the size of the single grid cells (0.6 m × 0.6 m), we argue that these localized conditions are geomorphologically unrealistic. Using “map algebra” techniques, the errors can be easily localized, and corrected. The method is based on the statistical analysis of variables that capture changes in the land cover conditions caused by the occurrence of the new landslides. The result of the analysis will depend on multiple factors, including the time between the pre-event and the post-event images. Ideally, this time should be short, to minimize effects introduced by land cover changes produced by causes different from the event landslides. In practice, a pre-event image taken immediately before an event is seldom available. Selection of the optimal pre-event image (where more than one is available) depends on: (i) the season or time of the year of the image, which should be the same of the post-event image, (ii) the viewing angle, preferring images with a small off-nadir angle to facilitate ortho-rectification, (iii) the area covered by the image, which should be similar to the area covered by the post-event image, and (iv) the proportion of cloud coverage, which should be minimal. In an operational environment, we expect that the selection of the pre-event image will be the result of a compromise between the listed criteria. The procedure assumes that the new event landslides are the only (or primary) cause for the changes in the land cover conditions. This is particularly relevant for δNDVI and θ. Changes introduced by other causes, including e.g. phenological variations, agricultural practices, constructions, and other landslides, introduce noise and may result in a less accurate terrain classification. Where multiple pre-event images are available, their joint analysis may allow for identifying (and resolving) potential problems in the terrain classification. The analysis of multiple pre-event images may also allow for the detection of landslides of different ages (or periods), facilitating the construction of multi-temporal inventory maps (Fiorucci et al., 2011).The procedure further assumes that the event landslides have left discernable features on the land cover, and that the land cover changes were captured faithfully by the post-event images. In places, landslides leave subtle morphological and land cover changes that are difficult to recognize (Fiorucci et al., 2011). In these areas, the proposed procedure may not be capable of detecting some or all of the new landslides. This is because the signal of the new landslides in the image is weak, or it similar to the signal produced by other features unrelated to the presence of the new landslides. Use of the procedure in the study area (9.4 km2) required five days of work, including two days of an expert in remote sensing for image processing and analysis, one day of a geomorphologist for the multivariate terrain classification, and one day of two geomorphologists for model validation and testing of the results. Overall, we estimate that production of the landslide event-inventory through the semi-automatic analysis of the satellite images required about half a day per person, per square kilometer. This rate is significantly faster than the rate required for preparing the event-inventory using the traditional approach of aerial photograph interpretation aided by field checks. The proposed procedure is demanding computationally, a consequence of using satellite images characterized by a very small GSD,

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

which results in large matrixes that require large computer memory to process and produce large output files. Should the procedure be applied to a larger area, extending for hundreds to thousands of square kilometers, the procedure may require large computer resources. With this respect, availability of VHR images may be a limitation, and use of images with a coarser GSD may be the only practical solution. Availability of images with an improved radiometric resolution may compensate for the coarser GSD. Selection of the most appropriate geographical and radiometric resolution will also depend on the available satellite sensors and the scope of the mapping, including: (i) the size of the smallest landslide that needs to be recognized consistently, (ii) the requested degree of completeness for the inventory, and (iii) the statistics of landslide size (Guzzetti et al., 2002; Malamud et al., 2004; Fiorucci et al., 2011). The method used to recognize and map the rainfall induced landslides did not allow for the separation of the different sections of the shallow slope failures in the study area, including the source, transport and deposition areas (Crosta et al., 1990). This is a limitation of the method, particularly if the event-inventory is prepared for susceptibility, hazard, or risk assessment purposes. However, consistent and accurate mapping of the different geomorphological subdivisions of a shallow landslide, and specifically of a debris flow, is problematic and uncertain even in the field (Fig. 2), or through the visual interpretation of large-scale aerial photographs. With this respect, the event-inventory prepared through the semi-automatic analysis of the satellite images is not different from, or inferior to, the original landslide event-inventory available for the study area (Fig. 1). Further, in an operational environment rapid production of an inventory after an event is often more important than the preparation of an extremely detailed inventory that requires large resources and significant time to be completed. During the initial image processing phase, urban areas were removed from the analysis to reduce the noise introduced by the different viewing angles of the pre-event and the post-event images (Table 1). As a result, rainfall induced shallow landslides were not detected in the urban areas. This is a limitation of the landslide eventinventory. The problem must be considered if the proposed procedure is to be applied in areas where urban zones are present. Depending on the characteristics of the images, and the type, height, and density of the buildings, the local terrain conditions, the extent, size and radiometric signature of the landslides, and the possibility of detecting landslides in urban areas varies. For detecting and mapping the event landslides, the proposed method used an approach based on the spectral properties of the individual pixels. External information and heuristic, statistical or phenomenological criteria were not used to help in detecting the new landslides, or to improve the landslide mapping. Use of external information and adoption of a specific set of relevant criteria may improve the mapping. Alternatively, an object oriented approach could be adopted to detect the new landslides, or to improve the mapping obtained by our method. An object oriented approach could exploit information extracted from the satellite images (Barlow et al., 2006; Martha et al., 2010) and information obtained from a high resolution DEM (Passalacqua et al., 2010; Tarolli et al., 2010). As a limitation, an object oriented approach may be less general than a pixel based approach for landslide mapping, and may require additional time and specific expertise to establish and tune an appropriate set of criteria. Finally, we like to stress that availability of a method to recognize and map recent landslides rapidly and semi-automatically using standard satellite images and image processing tools represents an important aid for the production of accurate and reasonably complete landslide event-inventory maps. Where adequate pre-event and postevent satellite images are available, the procedure can be used in the aftermath of an event to obtain preliminary, but reasonably accurate, estimates of the number, extent, and density of landslides in an area,

1755

and of the total landslide area. Analyzing the event-inventory, expert geomorphologists can infer information on the type of the slope failures (Cruden & Varnes, 1996). Application of empirical relationships to link landslide area and volume (Guzzetti et al., 2009; Korup et al., 2010) allows estimates of the volume of the individual slope failures, and of the total volume of mobilized or eroded material. This is important information for erosional studies (Lavé & Burbank, 2004; Imaizumi & Sidle, 2007; Guzzetti et al., 2008, 2009; Fiorucci et al., 2011) and for vulnerability estimates. The event-inventory maps obtained exploiting the satellite images can always be refined (corrected) through the interpretation of post-event aerial photographs (where available), or field surveys. Existence of a preliminary inventory facilitates the preparation and reduces the cost for the production of a new inventory. 8. Conclusions We have presented a new method for the semi-automatic detection and mapping of rainfall induced shallow landslides. The method exploits information obtained from pre-event and post-event, VHR panchromatic and HR multispectral satellite images, using standard image processing techniques. Processing of the images was aimed at detecting areas where changes attributed to slope failures have occurred in the period between the pre-event and the post-event images, chiefly landslides triggered by a high intensity rainfall event on 1 October 2009. The method proved successful for the identification of areas affected by recent slope failures in the 9.4 km2 study area. The good performance of the terrain classification is measured by the area under the receiver operating characteristic curve, AROC N 0.860 in the training area and AROC N 0.800 in the validation area. The large values of the AROC metrics prove the ability of the classification models and – most relevant – the efficacy of the proposed method. We expect that the method will facilitate the production of landslide event inventories over large areas, where pre-event and post-event satellite imagery of adequate quality is available. The method was tested in Sicily, Italy, to detect rainfall induced shallow landslides. However, the method can be used in other areas to detect landslides caused by similar (rainfall) or different triggers, including earthquakes, provided that the event slope failures have left discernable features captured by the post-event satellite images, and that the available terrain information and satellite images are of adequate quality. The proposed method can facilitate the rapid production of accurate landslide event inventories covering large areas. We expect that this will improve our ability to evaluate landslide hazards (Guzzetti et al., 2005, 2006a, 2006b), will foster our understanding of the evolution of landscapes shaped by mass-wasting processes (Lavé & Burbank, 2004; Imaizumi & Sidle, 2007; Guzzetti et al., 2009; Korup et al., 2010), and will augment the ability of civil protection and land management authorities to respond to landslide disasters. Acknowledgments Work conducted in the framework of the ASI MORFEO project. ACM was supported by a grant of the Agenzia Spaziale Italiana, and MR was supported by a grant of the National Department for Civil Protection. We thank G. Iovine (CNR IRPI), the several colleagues at the University of Florence, and A. Gennari (University of Perugia) for contributing to preparing the landslide event-inventory, and L. Santurri (CNR IFAC) for recommendations on image processing. We thank the four anonymous reviewers for their comments and recommendations. References Aiazzi, B., Baronti, S., & Selva, M. (2007). Improving component substitution pansharpening through multivariate regression of MS+ Pan data. IEEE Transactions on Geoscience and Remote Sensing, 45(10), 3230−3239, doi:10.1109/TGRS.2007.901007.

1756

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757

Aiazzi, B., Baronti, S., Lotti, F., & Selva, M. (2009). A comparison between global and context-adaptive pansharpening of multispectral images. IEEE Geoscience and Remote Sensing Letters, 6(2), 302−306. Ardizzone, F., Cardinali, M., Galli, M., Guzzetti, F., & Reichenbach, P. (2007). Identification and mapping of recent rainfall-induced landslides using elevation data collected by airborne Lidar. Natural Hazards and Earth System Sciences, 7, 637−650. Barlow, J., Franklin, S., & Martin, Y. (2006). High spatial resolution satellite imagery, DEM derivatives, and image segmentation for the detection of mass wasting processes. Photogrammetric Engineering and Remote Sensing, 72(6), 687−692. Brown, C. E. (1998). Applied multiple statistics in geohydrology and related sciences (pp. 248). : Springer–Verlag. Bucknam, R. C., Coe, J. A., Chavarria, M. M., Godt, J. W., Tarr, A. C., Bradley, L. -A., Rafferty, S., Hancock, D., Dart, R. L., & Johnson, M. L. (2001). Landslides triggered by hurricane mitch in Guatemala — Inventory and discussion. U.S. Geological Survey Open File Report (pp. 01−443). Burrough, P. A., & McDonnell, R. A. (1998). Principles of geographical information systems. Oxford: Oxford University Press 333 pp. Cardinali, M., Ardizzone, F., Galli, M., Guzzetti, F., & Reichenbach, P. (2001). Landslides triggered by rapid snow melting: The December 1996–January 1997 event in Central Italy. In P. Claps, & F. Siccardi (Eds.), Proceedings 1st EGS Plinius Conference, Maratea (pp. 439−448). Cosenza: Bios Publisher. Cardinali, M., Galli, M., Guzzetti, F., Ardizzone, F., Reichenbach, P., & Bartoccini, P. (2006). Rainfall induced landslides in December 2004 in south-western Umbria, central Italy: Types, extent, damage and risk assessment. Natural Hazards and Earth System Sciences, 6, 237−260. Cheng, K. S., Wei, C., & Chang, S. C. (2004). Locating landslides using multi-temporal satellite images. Advances in Space Research, 33(3), 96−301. Comon, P. (1994). Independent component analysis, a new concept. Signal Processing, 36, 287−314. Corsini, A., Borgatti, L., Coren, F., & Vellico, M. (2007). Use of multitemporal airborne LiDAR surveys to analyse post-failure behaviour of earthslides. Canadian Journal of Remote Sensing, 33(2), 116−120. Crosta, G., Guzzetti, F., Marchetti, M., & Reichenbach, P. (1990, August 6–10). Morphological classification of debris-flow processes in South-Central Alps (Italy). Proceedings VI International Congress IAEG, Amsterdam, The Netherlands (pp. 1565−1572). Cruden, D. M., & Varnes, D. J. (1996). Landslides types and processes. In A. K. Turner, & R.L. Schuster (Eds.), Landslides: Investigation and mitigation, transportation research board special report, vol. 247. (pp. 36−75)Washington: National Academy Press. Czuchlewski, K. R., Weissel, J. K., & Kim, Y. (2003). Polarimetric synthetic aperture radar study of the Tsaoling landslide generated by the 1999 Chi-Chi earthquake, Taiwan. Journal of Geophysical Research, 108(F1), 7.1−7.11. Dermanis, A., & Biagi, L. (2006). Telerilevamento. Informazione territoriale mediante immagini da satellite.Milano: Casa Editrice Ambrosiana 279 pp. (in Italian). Digital Globe Inc. (2007). QuickBird Imagery Products, Product Guide, revision 4.7.3. El Hajj, M., Bégué, A., Lafrance, B., Hagolle, O., Dedieu, G., & Rumeau, M. (2008). Relative radiometric normalization and atmospheric correction of a SPOT 5 time series. Sensors, 8, 2774-2791 ISSN 1424–8220. Fawcett, T. (2006). An introduction to ROC analysis. Pattern Recognition Letters, 27, 861−874. Fiorucci, F., Cardinali, M., Carlà, R., Rossi, M., Mondini, A. C., Santurri, L., Ardizzone, F., & Guzzetti, F. (2011). Seasonal landslides mapping and estimation of landslide mobilization rates using aerial and satellite images. Geomorphology, doi:10.1016/j. geomorph.2011.01.013. Green, D. M., & Swets, J. M. (1966). Signal detection theory and psychophysics.New York: John Wiley and Sons Inc ISBN 0–471–32420–5, 455 pp. Guzzetti, F., Carrara, A., Cardinali, M., & Reichenbach, P. (1999). Landslide hazard evaluation: A review of current techniques and their application in a multi-scale study, Central Italy. Geomorphology, 31, 181−216. Guzzetti, F., Cardinali, M., Reichenbach, P., & Carrara, A. (2000). Comparing landslide maps: A case study in the Upper Tiber River Basin, Central Italy. Environmental Management, 25(3), 247−263, doi:10.1007/s002679910020. Guzzetti, F., Malamud, B. D., Turcotte, D. L., & Reichenbach, P. (2002). Power-law correlations of landslide areas in central Italy. Earth and Planetary Science Letters, 195, 169−183. Guzzetti, F., Cardinali, M., Reichenbach, P., Cipolla, F., Sebastiani, C., Galli, M., & Salvati, P. (2004). Landslides triggered by the 23 November 2000 rainfall event in the Imperia Province, Western Liguria, Italy. Engineering Geology, 73(2), 229−245. Guzzetti, F., Reichenbach, P., Cardinali, M., Galli, M., & Ardizzone, F. (2005). Probabilistic landslide hazard assessment at the basin scale. Geomorphology, 72, 272−299. Guzzetti, F., Galli, M., Reichenbach, P., Ardizzone, F., & Cardinali, M. (2006a). Landslide hazard assessment in the Collazzone area, Umbria, central Italy. Natural Hazards and Earth System Sciences, 6, 115−131. Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., & Galli, M. (2006b). Estimating the quality of landslide susceptibility models. Geomorphology, 81, 166−184. Guzzetti, F., Ardizzone, F., Cardinali, M., Galli, M., & Reichenbach, P. (2008). Distribution of landslides in the Upper Tiber River basin, central Italy. Geomorphology, 96, 105−122. Guzzetti, F., Ardizzone, F., Cardinali, M., Galli, M., Rossi, M., & Valigi, D. (2009). Landslide volumes and landslide mobilization rates in Umbria, central Italy. Earth and Planetary Science Letters, 279, 222−229, doi:10.1016/j.epsl.2009.01.005. Hadjimitsis, D. G., Clayton, C. R. I., & Retalis, A. (2004). Darkest pixel atmospheric correction algorithm: A revised procedure for environmental applications of satellite remotely sensed imagery. Proceedings 10th International Symposium on Remote Sensing, SPIE, 5239. (pp. 464).

Hadjimitsis, D. G., Papadavid, G., Agapiou, A., Themistocleous, K., Hadjimitsis, M. G., Retalis, A., Michaelides, Chrysoulakis, S. N., Toulios, L., & Clayton, C. R. I. (2010). Atmospheric correction for satellite remotely sensed data intended for agricultural applications: Impact on vegetation indices. Natural Hazards and Earth System Sciences, 10, 89−95. Harp, E. L., & Jibson, R. L. (1996). Landslides triggered by the 1994 Northridge, California earthquake. Seismological Society of America Bulletin, 86, S319−S332. Hayes, D. J., & Sader, S. A. (2001). Change detection techniques for monitoring forest clearing and regrowth in a tropical moist forest. Photogrammetric Engineering and Remote Sensing, 67(9), 1067−1075. Hervás, J., Barredo, J. I., Rosin, P. L., Pasuto, A., Mantovani, F., & Silvano, S. (2003). Monitoring landslides from optical remotely sensed imagery: The case history of Tessina landslide, Italy. Geomorphology, 54, 63−75. Hong, G., & Zhang, Y. (2008). A comparative study on radiometric normalization using high resolution satellite images. International Journal of Remote Sensing, 29(2), 425−430. Hyvärinen, A., & Oja, E. (2000). Independent component analysis: Algorithms and applications. Neural Networks, 13(4–5), 411−430, doi:10.1016/S0893-6080(00) 00026-5. Imaizumi, F., & Sidle, R. C. (2007). Linkage of sediment supply and transport processes in Miyagawa Dam catchment, Japan. Journal Geophysical Research, 112, F03012, doi: 10.1029/2006JF000495. Jianya, G., Haigang, S., Guorui, M., & Qiming, Z. (2008). A review of multi-temporal remote sensing data change detection algorithms. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 37. (pp. 757−762): B7. Jollifee, I. T., & Stephenson, D. B. (2003). Forecast verification.A practitioner's guide in atmospheric science. : John Wiley 240 pp. Jolliffe, I. T. (1986). Principal component analysis.: Springer-Verlag978-0-387-95442-4, doi:10.1007/b98835 487 pp. Klinker, G. J., Shafer, S. A., & Kanade, T. (1992). The measurement of highlights in color images (pp. 309−334). USA: Jones and Bartlett Publishers0-86720-295-5. Korup, O., Densmore, A. L., & Schlunegger, F. (2010). The role of landslides in mountain range evolution. Geomorphology, 120(1–2), 77−90, doi:10.1016/j.geomorph.2009.09.017. Laben, C.A., Brower, B.V. (2000). Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening. US Patent 6,011,875. Lavé, J., & Burbank, D. (2004). Denudation processes and rates in the Transverse Ranges, southern California: Erosional response of a transitional landscape to external and anthropogenic forcing. Journal of Geophysical Research, 109, F01006, doi: 10.1029/2003JF000023, 2004. Lee, S., & Lee, M. -J. (2006). Detecting landslide location using KOMPSAT 1 and its application to landslide-susceptibility mapping at the Gangneung area, Korea. Advances in Space Research, 38, 2261−2271. Lin, P. -S., Lin, J. -Y., Hung, H. -C., & Yang, M. -D. (2002). Assessing debris flow hazard in a watershed in Taiwan. Engineering Geology, 66, 295−313. Lin, C. -W., Shieh, C. -L., Yuan, B. -D., Shieh, Y. -C., Liu, S. -H., & Lee, S. -Y. (2004). Impact of Chi-Chi earthquake on the occurrence of landslides and debris flows: Example from the Chenyulan River watershed, Nantou, Taiwan. Engineering Geology, 71(1–2), 49−61. Luckman, P. G., Gibson, R. D., & Derose, R. C. (1999). Landslide erosion risk to New Zealand pastoral steeplands productivity. Land Degradation & Development, 10(1), 49−65. Malamud, B. D., Turcotte, D. L., Guzzetti, F., & Reichenbach, P. (2004). Landslide inventories and their statistical properties. Earth Surface Processes and Landforms, 29(6), 687−711. Mangolini, M., Rankin, T., Wald, L. (1992). French patent No. 92–13961, 20 November 1992; Methods and Device for Increasing the Spatial Resolution of Images from Other Images of Better Spatial Resolution, USA patent serial 081249,882, 26 May 1994. Marcelino, E. V., Formaggio, A. R., & Maeda, E. E. (2009). Landslide inventory using image fusion techniques in Brazil. International Journal of Applied Earth Observation and Geoinformation, 11, 181−191. Martha, T. R., Kerle, N. V. Jetten, van Westen, C. J., & Vinod Kumar, K. (2010). Characterising spectral, spatial and morphometric properties of landslides for semi-automatic detection using object-oriented methods. Geomorphology, 116(1–2), 24−36. Masek, J. G., & Sun, G. (2004). Technical note: A spectral-angle methodology for mapping net forest cover change in northeastern China. International Journal of Remote Sensing, 25, 5629−5636. Mason, S. J., & Graham, N. E. (2002). Areas beneath the relative operating characteristics (ROC) and relative operating levels (ROL) curves: Statistical significance and interpretation. Quarterly Journal Royal Meteorological Society, 128, 2145−2166. Metternicht, G., Hurni, L., & Gogu, R. (2005). Remote sensing of landslides: An analysis of the potential contribution to geo-spatial systems for hazard assessment in mountainous environments. Remote Sensing of Environment, 98, 284−303. Michie, D., Spiegelhalter, D. J., & Taylor, C. C. (Eds.). (1994). Machine learning, neural and statistical classification. Internet version. http://www.amsta.leeds.ac.uk/~charles/ statlog/. Munechika, C. K., Warnick, J. S., Salvaggio, C., & Schott, J. R. (1993). Resolution enhancement of multispectral image data to improve classification accuracy. Photogrammetric Engineering & Remote Sensing, 59(1), 67−72. Nichol, J., & Wong, M. S. (2005). Satellite remote sensing for detailed landslide inventories using change detection and image fusion. International Journal of Remote Sensing, 26(9), 1913−1926. Passalacqua, P., Tarolli, P., & Foufoula-Georgiou, E. (2010). Testing space-scale methodologies for automatic geomorphic feature extraction from LiDAR in a complex mountainous landscape. Water Resources Research, 46, W11535, doi: 10.1029/2009WR008812.

A.C. Mondini et al. / Remote Sensing of Environment 115 (2011) 1743–1757 Richards, J. A., & Jia, X. (2006). Remote sensing digital image analysis (4th edition). : Springer978-3-540-25128-6 439 pp. Rosin, P. L., & Hervás, J. (2005). Remote sensing image thresholding methods for determining landslide activity. International Journal of Remote Sensing, 26(6), 1075−1092. Rossi, M., Guzzetti, F., Reichenbach, P., Mondini, A. C., & Peruccacci, S. (2010). Optimal landslide susceptibility zonation based on multiple forecasts. Geomorphology, 114, 129−142, doi:10.1016/j.geomorph.2009.06.020. Santangelo, M., Cardinali, M., Rossi, M., Mondini, A. C., & Guzzetti, F. (2010). Mapping landslides with landslides with a laser rangefinder binocular. Natural Hazards and Earth System Sciences, 10, 2539−2546, doi:10.5194/nhess-10-2539-2010. Singhroy, V., & Molch, K. (2004). Characterizing and monitoring rockslides from SAR techniques. Advances in Space Research, 33(3), 290−295. Stone, J. V. (2005). A brief introduction to independent component analysis. In B. S. Everitt, & D. C. Howell (Eds.), Encyclopedia of statistics in behavioral science (pp. 907−912). Chichester: John Wiley & Sons Ltd978-0-470-86080-9.

1757

Tarolli, P., Sofia, G., & Dalla Fontana, G. (2010). Geomorphic features extraction from high resolution topography: Landslide crowns and bank erosion. Natural Hazards, doi:10.1007/s11069-010-9695-2. Tsai, F., Hwang, J. -H., Chen, L. -C., & Lin, T. -H. (2010). Post-disaster assessment of landslides in southern Taiwan after 2009 Typhoon Morakot using remote sensing and spatial analysis. Natural Hazards Earth System Sciences, 10, 2179−2190, doi: 10.5194/nhess-10-2179-2010. Venables, W. N., & Ripley, B. D. (2002). Modern applied statistics with S (pp. 495). New York: Springer0-387-95457-0. Wald, L., Ranchin, T., & Mangolini, M. (1997). Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images. Photogrammetric Engineering & Remote Sensing, 63(6), 691−699. Weirich, F., & Blesius, L. (2007). Comparison of satellite and air photo based landslide susceptibility maps. Geomorphology, 87(4), 352−364.