Next Article in Journal
Estimation of Herbaceous Fuel Moisture Content Using Vegetation Indices and Land Surface Temperature from MODIS Data
Next Article in Special Issue
Characterization of Landslide Deformations in Three Gorges Area Using Multiple InSAR Data Stacks
Previous Article in Journal / Special Issue
Topographic Correction of Wind-Driven Rainfall for Landslide Analysis in Central Taiwan with Validation from Aerial and Satellite Optical Images
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Supervised Method of Landslide Inventory Using Panchromatic SPOT5 Images and Application to the Earthquake-Triggered Landslides of Pisco (Peru, 2007, Mw8.0)

1
ISTerre, IRD, CNRS, Université de Grenoble, 1381 rue de la piscine, F-38041 Grenoble, France
2
INGEMMET, Av. Canada 1470, San Borja, Lima 41, Peru
3
LEGOS, CNRS, Université de Toulouse, 14 av. Ed. Belin, F-31400 Toulouse, France
*
Author to whom correspondence should be addressed.
Remote Sens. 2013, 5(6), 2590-2616; https://doi.org/10.3390/rs5062590
Submission received: 28 March 2013 / Revised: 9 May 2013 / Accepted: 9 May 2013 / Published: 23 May 2013

Abstract

:
Earthquake is one of the dominant triggering factors of landslides. Given the wide areas covered by mega earthquake-triggered landslides, their inventory requires development of automatic or semi-automatic methods applied to satellite imagery. A detection method is here proposed for this purpose, to fit with simple datasets; SPOT5 panchromatic images of 5 m resolution coupled with a freely and globally available DEM. The method takes advantage of multi-temporal images to detect changes based on radiometric variations after precise coregistration/orthorectification. Removal of false alarms is then undertaken using shape, orientation and radiometric properties of connected pixels defining objects. 80% of the landslides and 93% of the landslide area are detected indicating small omission errors but 50% of false alarms remain. They are removed using expert based analysis of the inventory. The method is applied to realize the first comprehensive inventory of landslides triggered by the Pisco earthquake (Peru, 15/08/2007, Mw 8.0) over an area of 27,000 km2. 866 landslides larger than 100 m2 are detected covering a total area of 1.29 km2. The area/number distribution follows a power-law with an exponent of 1.63, showing a very particular regime of triggering in this arid environment compared to other areas in the world. This specific triggering can be explained by the little soil cover in the coastal and forearc regions of Peru. Analysis of this database finally shows a major control of the topography (both orientation and inclination) on the repartition of the Pisco-triggered landslides.

1. Introduction

Earthquake is one of the main triggering factors of landslide [1] and the most damaging trigger of landslides in terms of social cost. As an example, more than 60% of casualties in landslides are encountered in earthquake-triggered landslides over the 2004–2010 period [2]. In this context, understanding and quantifying the susceptibility of slopes to landslide is of major interest.
Triggering of landslides by earthquakes is a complex problem as it encompasses two main domains: seismology and landslide mechanics. The later includes several fields such as rock mechanics, geomorphology, hydrology and hydrogeology. Numerous factors, related both to the earthquake and to the settings in which the seismic waves propagate, drive the number, size and type of landslides. From the study of landslide inventories, it emerges that the triggering of landslides by earthquake is controlled by two main class of parameters: (i) the earthquake source properties (magnitude, duration, fault rupture dynamics) (e.g., [1,35]); (ii) the site properties (geology, topography, geometry, water content) (e.g., [69]). The source magnitude controls the area affected by landslides [1] relating the ground motion decay with distance to the epicenter or to the fault plane (e.g., [3]). However this generic law is affected by different factors: (1) a significant over-abundance of landslides in the hanging walls of ruptured thrust faults relative to the footwalls rate are reported (e.g., [4,5]), showing the effect of the dynamics of the fault rupture on the landslide distribution; (2) the greater erosion rate due to precipitations observed after the earthquake [8,10] shows that earthquakes are found to interact non-linearly with precipitations; (3) the site effect, or amplification of the seismic waves due to shear-wave velocity contrast between the landslide material and the stable material [11,12] or topographical particularities like ridges [9]. This topographic control on co-seismic landslides has also been pointed out through different parameters: convexity [7], slopes, proximity from ridges and crests [6,9].
The effects of all the parameters listed above are fairly well known qualitatively but not quantitatively due to the very few number of cases studied. Indeed, because of the tedious work of detection, the studies of earthquake-triggered landslides have been based on approximately 40 coseismic landslide inventories worldwide [3,13], sampling only few different situations of magnitude, focal mechanism, climatic settings and land-surfaces. All these results emphasize the need for the generation of new datasets in order to draw robust conclusions on co-seismic landslides. Landslide inventories are, most of the time, based on field inventories and/or visual analysis of aerial or satellite images. These methods often lead to incomplete or biased inventories, due to subjectivity of the operator [5,14]. Therefore, automatic or semi-automatic detection methods have been developed in the past years (e.g., [1523]). These methods are either based on the supervised or unsupervised classification of one satellite image (e.g., [18,24]), or on the detection of new landslides in a pair of images acquired at different dates [20,21]. Detection is either based on the comparison of individual pixels or on the analysis of objects extracted from the images. The interest of object-oriented methods is to enable analysis of the object shape or texture that often characterizes the landslides compared with other objects [19]. For a review of the landslide detection methods, see Guzzetti et al.[25].
Accuracy of the methods is generally reported using both the detection error (omission = real landslides that were not detected) and the false alarm rate (commission = “landslides” that were detected but did not exist). Unsupervised methods display errors of commission and omission of 56% and 24% respectively [18,19]. Supervised methods use training data to estimate the parameters of the image segmentation in object oriented algorithms. It allows to increase this accuracy with errors of commission between 12% and 27%, similar to omission errors [20,21,23]. These 3 latter methods showing better accuracy used Very High Resolution (VHR, ≤1 m) sensors, whose past archives have small coverage on the globe. It therefore restrains the use of such datasets to very specific regions.
The detection has been tested with different types of sensors. Multi-spectral optical images are particularly well suited for object-oriented analysis [16,19]. Moreover, the use of Normalized Difference Vegetation Index (NDVI) coming from Landsat TM, SPOT or Quickbird images is of high interest to detect changes of vegetation and thus landslide in areas covered by vegetation. The high precision of Lidar (e.g., [26]) allows the detection of specific textures characterizing landslides, but the Lidar datasets do not cover extensive areas. On the contrary, radar images present the advantages of a good spatial and temporal coverage even in cloudy regions. In particular, polarimetric radar has been used to detect new landslides [27,28], taking advantages of different scattering properties on the bare soil and forested ones. On bare soil area, the InSAR images also allows detection of slow moving landslides [29,30]. The drawbacks of radar images are (1) the presence of shadows or blind areas on the image that prevent monitoring specific slopes, and (2) their medium resolutions not adapted to detect landslides of small sizes, which are the most common [31]. The recent launch of DLR TerrarSar-X and Tandem-X missions (among others) improves the availablity of high resolution radar images, even if past archives are still limited. Most recent methods use VHR data, either Quickbird (0.6 m) or Ikonos data (0.4 m) (e.g., [21,23]). The compensation of such very high resolution instruments is that their images cover only a limited area. Most studies have therefore been performed at a basin scale (typically 50 km2) for rain-induced landslide inventories (e.g., [20]). In case of co-seismic landslides, the distance affected by landslides can be as much as 10 times the fault length for small earthquakes (Mw < 5) and 3 to 5 times the fault length for larger earthquakes [32]. That is, for an Mw 8.0 earthquake, landslides can be triggered up to 600 km away from the source. These large distances require imagery acquired with wide swath that only high to medium resolution satellite can provide.
Moreover, for reduced areas, the detection can benefit from fine Digital Elevation Model (DEM) acquired by lidar or derived from a couple of optical images (e.g., [19]). However, for larger areas, fine DEM (typical resolution of 10 m) is generally not available and/or pricing can be high. In these cases, the 90-m posting SRTM or 30-m posting ASTER GDEM datasets are free alternatives.
Finally, in most cases (especially when revisiting historical events), panchromatic images are the only optical dataset available. Landslide detection using panchromatic images [15,17,33,34] is less accurate than using multi-spectral images. The use of panchromatic images solely makes it difficult to identify clearly the spectral signature of landslides compared with other objects [34]. Moreover, all the proposed methods have been used on vegetated regions. In arid environment such as the coast of Peru, landslides occur on bare soil. Therefore, changes in texture are small and object detection based on texture analysis is of low efficiency. For example, Nichol and Wong [35] showed that a pixel-based method alone does not enable to detect landslides occurring on bare soil on SPOT images.
Therefore, in this study we adapt a pixel-based detection method to be working with (1) medium resolution (30 m) DEM, (2) panchromatic satellite images of high to medium resolutions, and (3) non-vegetated areas. We use 5 m SPOT panchromatic satellite images to cover extensive areas. We apply this method to pairs of co-seismic images, surrounding the Pisco earthquake (Mw 8.0, Peru, 15 August 2007) over an area of 27,000 km2, where a field inventory (Figure 1) has been realized in the months following the earthquake [36]. Using visual interpretation of satellite images, we show that the field inventory is not consistent with the satellite inventory and thus cannot be used to validate our method. The method is thus validated using a comparison with an expert-based inventory realized on two subsets of the satellite images. A first analysis of the topographical characteristics of the database is finally proposed.

2. Data

2.1. Study Area

The Pisco earthquake (Mw 8.0) struck the coast of Peru on 15 August 2007 (Figure 1). The epicenter was 18 km deep and ruptured the megathrust along the subduction zone of Nazca and southern American plates [38]. It broke a fault segment of 170 km, parallel to the coast, composed of 2 main asperities [37]. The analysis of the accelerogram registered in a radius of 200 km showed the complexity of this rupture [39].
This earthquake triggered landslides, 134 of which were detected during a field survey following the remaining accessible roads, the days after the earthquake [36]. This inventory is mostly composed by rockfalls. Rockfalls mainly occurred along roads that run at the foot-cut slope and were thus easily detected during the field work [36]. The resulting number-volume statistics do not exhibit a classical linear relation as plotted in log-log coordinates [31], showing that this inventory is certainly incomplete. Moreover, the total number of landslides reported in this inventory is far out of classical values for Mw 8.0 earthquake. Indeed Keefer [3] reported average number of landslides being on the order of 105 for an Mw 8.0 earthquake. A more complete inventory is therefore needed for better analysis and quantification of the triggering factors of these co-seismic landslides.

2.2. Satellite Images

For this study, we chose to use SPOT5 images, which are a good compromise between the spatial resolution (2.5 to 5 m pixel size) and the area covered (each scene is 60 × 60 km2). Moreover, the data programming since 2002 enables a collection of archive images of high interest for studying past co-seismic landslides.
The method developed in the following is based on the processing of image pairs before/after a particular event. Pairs of images have been chosen optimizing 3 criterions: (1) small B/H values, where B is the baseline between the optical center of the two images and H is the altitude of the satellite; (2) small cloud coverage; and (3) shortest temporal baseline. The first criterion is key for the orthorectification of the two images since it shows how DEM errors propagate in the difference of coregistration of the 2 images [40]. The last criterion is not a limiting factor for the Pisco case because the arid climate prevents the vegetation from changing rapidly. The areas covered by the images have been chosen to encompass most of the areas visited by the field survey after the earthquake.
Using these criteria, 8 pairs of images have been chosen and 2 of them acquired twice after the earthquake (Figure 1 and Table 1). Characteristics of the images are shown in Table 1. The total area covered by the image pairs is 27,000 km2, imaging 75% of the landslides detected by the field survey. Time separations between images of a pair are ranging from 3 to 6 years, which is short compared with the rate of growth of vegetation in this area. Moreover, triggering by other earthquakes is excluded in this time-lapse, even if numerous aftershocks were recorded with magnitude below Ml 6.1 Sladen et al.[37]. Indeed, these aftershocks are clustered updip of the hypocenter and south of the Paracas Peninsula, at least 60 km far from the main area of landslide (Figure 1). This distance excludes systematic landslide triggering by magnitude 5–6 earthquakes [3].

3. Methodology

The proposed landslide detection method follows a classical scheme: (1) precise image orthorectification; (2) change detection; (3) removal of false alarms. The change detection is adapted from a classical pixel-based method using the comparison of multi-temporal images. New developments from our study focus on the removal of false alarms adapted to the dataset used here.

3.1. Orthorectification

Orthorectification of SPOT5 imagery is performed using the Cosi-Corr software [41]. Each pre-earthquake image is first orthorectified using Ground Control Points (GCP) taken on the shaded DEM. GCPs are mostly located at the crossing between steep valleys or at the summit of the mountains. Around 30 GCPs are chosen for each image orthorectification. The ASTER Global DEM (GDEM-v2) is used, with a spatial posting of 1 arc-second (approximately 30 m in our area) and altitude errors (GDEM) between ±7 and ±15 m [42]. Then, the post-earthquake image is orthorectified on the same grid as the pre-earthquake one, using GCPs based on the picking of recognizable features (mostly man-made structures) on both images. The GCPs are first picked manually before running the Cosi-Corr automatic optimization in order to improve the quality of the co-registration [41].
The B/H values range between 0.0014 and 0.0307, except for pair #7 that exhibits a B/H value of 0.1578 (Table 1). Formal errors of coregistration co using the ASTER GDEM-v2 data are given by co = GDEM · B/H, which is between 2 and 46 cm except for pair #7 that can present errors up to 2.4 m. All these values are however much less than 1 pixel of the SPOT images (5 m in our case). The precision of this process is key as coregistration default could cause false change detection. A good coregistration is also key to detect landslides of small sizes. For consistency between all the pairs of image available, the resolution is set to 5 m for all images.

3.2. Clouds Detection

Clouds are first detected in order to remove these areas from the following analysis. First, we calculate 3 types of radiometric indicators over windows W of 50 × 50 m size, i.e., the mean μ and standard deviation σ of radiometric values for each image, and the coefficients of correlation C between the pre and post-images:
C = cov ( W pre , W post ) σ W pre σ W post
where cov is the covariance between the two windows selected. The clouds are detected by selecting the area of weak correlation (C ≤ 0.1), where either the pre- or post- image displays very high (μ is greater than 95% of the maximum histogram value) and homogeneous (σ below 5 pixels) radiometric values inside the 50 m × 50 m window. This initial area is then iteratively expanded to connected 50 m × 50 m areas that display C below 0.8 and μ greater than 50% of the maximum histogram value. These numbers have been chosen after a process of trial and errors on one image from the pair #1. It allows to detect 96% of the total cloud area, as estimated by a comparison with a manual identification of the clouds on one image. Only small clouds (of less than 50 m size) or the borders of the big clouds are not well detected and are thus removed manually at the end of the whole process.

3.3. Change Detection

In panchromatic images, new areas of landslides are found to exhibit higher brightness values (e.g., [34]). Therefore, we apply a classical image differencing technique [43] to detect the changes between the two multi-temporal images. This technique consists in normalizing the two images impre and impost (Equations (2) and (3)) to differentiate them (Equation (4)) and then selecting the pixels exhibiting large variations. We independently normalize each image since we noticed that very few pixels change their radiometry in between the 2 acquisitions. Indeed, the landscape is dominated by rocky mountains where the vegetation cover does not change rapidly. Also, the viewing angles of the 2 acquisitions are very close (less than 1° for all pairs, except for pair #6 that exhibits difference of angles of 7°) and all images experienced little change in solar incidence due to the low latitudes of the studied area. The normalized images are thus expressed through (Figure 2(a)):
im pre = ( im pre μ pre ) / σ pre
im post = ( im post μ post ) / σ post
The difference image is expressed through (Figure 2(b)):
diff = im post im pre
where μpre and μpost refer to mean values of pre- and post-images, and σpre, σpost to their standard deviations. In the difference image, the new landslides exhibit large positive pixel values. We therefore select the pixels above a threshold T. This threshold can be defined through different methods [33]. Due to the Gaussian shape of the difference image (Figure 2(b)), we define the threshold as an affine function of the standard deviation of the difference image diff (noted σd hereafter):
T = A σ d
where A is a parameter of the method.

3.4. False Alarms Removal

This thresholding method applied to our images leads to many false alarms, which can be classified in 3 classes: (1) river bed variations inside deep valleys; (2) anthropogenic changes of several types such as villages growth, new roads on slopes, agricultural fields on valleys and valley sides close to villages; (3) variations of solar reflections. These false alarms are usually removed using a combination of multispectral information, texture, shape and precise topography (e.g., [19,21]). However, in our case not all of these information are available. The proposed methodology therefore mostly focuses on the development of methods for removing false alarms, adapted to the available datasets (SPOT5 panchromatic images of 5 m resolution, GDEM-v2, non-vegetated area).

3.4.1. Objects Definition

The connected pixels detected by the threshold analysis are first gathered. This step leads to the definition of objects, with typical radiometry, slopes, and shape characteristics. We first eliminate objects containing less than 4 pixels. We then compute mean mo and standard deviation σo of the radiometry for each object o. Then the slope inclination β(x, y) and aspect θ(x, y) (orientation of the slope) at each point of the DEM (H(x, y)) are calculated. The slope of an object βo is defined as the median of the series of β(x, y) obtained for all the pixels of each object. Slope aspect of an object θo is defined by the median of the series of θ(x, y) for all the pixels of each object. We also compute the main direction of each object Do through a linear regression between the North and East coordinates of the pixels of each object. We can notice that for an object aligned in the main slope direction, Do and θo will be very similar.
We finally compute the correlation coefficient C(x, y) between the pre- and post-images (See Equation (1)) over windows of parametrized size w, at each point (x,y) separated by S = w/2. We thus define the correlation Co of an object as the mean correlation of its pixels. The separation between the landslides and other objects is done in the following using these different characteristics.

3.4.2. Change Detection at Large Scales

To detect the changes at large scales, we use the fact that the correlation between two images displaying large changes is low (e.g., [44]). As an example, the correlation C(x, y) is shown for a subset of an image where changes occur (Figure 3(a)). These changes are mostly due to anthropogenic factors during the time lapse between the image acquisitions (Figure 3). We notice that anthropogenic changes mostly occur in or near villages, due to the grouping of the housings and agricultural fields. New roads are treated separately in the Section 3.4.4. Peruvian villages are either settled inside valleys or on the flattest area of the mountains with expansion of the agricultural fields on its flanks. Therefore, villages can be detected by first identifying the flattest area where changes occur and then expanding the area to connected steeper zones where the correlation displays also similar values.
Therefore the detection can be summarized in 3 steps:
  • We select all the objects with βo below a slope value α (Figure 3(b)). Each object has a correlation Co.
  • For one of this object, we search its connected points (distant by less than w pixels) having a correlation lower or equal to Co (Figure 3(c)).
  • Step 2 is iterated until no connected points display values lower than or equal to Co.
The window size w is a parameter of our method whose range must be defined in relation with the expected maximum size of the landslides to avoid their removals. For our study, we identified visually the largest landslides and defined the maximum value of w at 128 pixels, i.e., 640 m.
Other changes at large scales can occur, notably large variations of illuminations. We detect these false alarms by selecting objects that displays Co lower than 0.7. This fairly conservative threshold value has been chosen after a process of trial and error to avoid removing objects that could be landslides.

3.4.3. River Bed Variation

The changes in river bed due to sediment transports or river level variations also cause false alarms. Some of these changes, mostly the ones occurring on slopes lower than α, are detected by the previous process (Section 3.4.2). However, some of the rivers flow in steeper valleys. Moreover, because of the low resolution of our DEM, the low slopes are only partly well defined. Therefore, we also identify the narrow valleys by detecting the zero crossings of the smoothed DEM partial derivatives in both directions. The DEM is smoothed using a Gaussian filter G(x, y; σ) of σ = 2 pixels in order to detect the main valleys with better robustness. The partial derivatives of the topography are then calculated using the convolution of H(x, y) and the gradient of G:
I = H ( x , y ) * G ( x , y ; σ ) = H ( x , y ) G ( x , y , σ )
A valley is detected by a negative change of the sign of the angle of I⃗. A positive change indicates that this point is situated on a crest. An illustration of the valley detection is given in Figure 4. We note that due to the low resolution of the DEM, a misfit can exist between the detected valleys and the real valleys. Therefore, an object is considered inside the valley when its distance to the detected valley is lower than a certain value V. V is chosen depending on the DEM resolution (here 30 m) and the standard deviation σ of the Gaussian function used to filter the DEM (here σ = 2 pixels). We define V as V = 30, σ = 60 m in our case.
We also note that small thalwegs/valleys are not detected by this method, due to the low resolution of the DEM combined with the Gaussian filter applied on the DEM (Section 3.4.1). Detecting small thalwegs was not a main objective of the valley detection because many landslides can flow inside and objects situated inside small thalwegs must therefore not be removed.

3.4.4. Road Detection

A number of false alarms arise from roads on the slopes either because of new roads that exhibit higher radiometry or because of changes of lighting. The sinuous roads in the mountain exhibit a large variability of radiometry. Therefore in our analysis, the same road is often detected in many small objects (Figure 5). We detect the roads based on their aspect, as roads have rarely slopes greater than 10% and are therefore perpendicular to the main slope direction. The object direction Do is compared with the mean direction of the perpendicular to the local slope Ds. The objects that exhibit |DoDs| < B are removed. B is a parameter of the model with possible value between 0° and 90°.
A function is first developed to detect anthropogenic changes and changes of solar reflections that affect large scales (Section 3.4.2). Other functions are also developed to remove the main false alarms of smaller scales: roads and river bed variation in steep and narrow valleys.

4. Validation

4.1. Comparison with the Field Inventory

We first compared the results of the detection with the field inventory by checking if the landslides identified on the field are detected by the automatic method. Out of 134 landslides inventoried in the field, 101 are situated on the image footprints and only 5 (5% of the landslides covered by the images) are detected by the automatic procedure. This low score can be explained by various factors:
  • The field inventory is dominated by landslides of small size. Among the landslides from the field inventory covered by the images, 79% are less than 66 m3. Assuming a relation area versus volume of V = 0.146A1.33[45], we infer a size of 100 m2, or 4 pixels of a SPOT5 image. This means they are removed by our method during the step of object definition (Section 3.4).
  • The field inventory is also constituted by many rockfalls (45% of the database), which are difficult to see on satellite images taken at vertical angle
  • Most of the mass movements detected by the field inventory occurred on the road [36], which were already cleaned at the time of the satellite acquisitions (at least 2 months after the earthquake, see Table 1). These deposits therefore cannot be detected on the satellite images. One exception is a small landslide scar (100 m2) detected in the SPOT5 images, which furthermore matches the location of a mass movement detected on the field.
  • There are uncertainties on the field inventory where no testimony exists, i.e., in arid regions the growth of vegetation is slow, and it is therefore difficult to really give a date to the observed landslides.
Figure 6 is an example of the difficult interpretation of the field work. Indeed, this image presents 2 different landslides distant by only 600 m. The first one has been detected by the field inventory, but the multi-temporal images show that it occurred between 19 June 2003 and 26 July 2007, before the 15 August 2007 earthquake. The other one has not been detected by the field inventory because the local relief prevents to see it well from the road. Similarly, many remote landslides cannot be detected by the field inventory, which has been realized by following the roads in the main valleys.
This comparison shows that the two inventories are not imaging the same types of landslides. The field inventory detects mostly the collapses on or very close to the road, with a better detection of rockfalls and mass movements of small volumes than with satellite images. The remote-sensing inventory detects the larger mass movements outside of the roads.

4.2. Sensitivity Analysis

For a better quantification of the method efficiency, we select 2 specific areas of 11 × 9 km2 (Area A1) and 10 × 6 km2 (Area A2), which encompass both areas of valleys and mountains. Those two sub-regions have topographic characteristics that are similar to our whole study region, with slope distribution showing a mean value of 35°/40° and maximum slope of 65°. An exhaustive inventory of landslides is realized by a visual comparison of the pre- and post-images for the two areas. Respectively 174 and 45 landslides are identified for A1 and A2, corresponding to areas of 64,250 m2 and 24,875 m2 respectively.
We test the performance of the automatic detection method on the test areas and calculate the resulting errors of omission (landslides not detected) Eo,i and commission (false alarms) Ec,i relative either to the area (i = a) or to the numbers of landslides (i = n) (e.g., [46]) :
E o , i = 1 i true i total
E c , i = 1 i false i changes
where itrue corresponds to the area or numbers of actual landslides well detected by the method, ifalse corresponds to the area or numbers of changes detected by the method but are not landslides, itotal is the total number or area of actual landslides, ichanges is the total number or area of changes detected.
The automatic detection and the resulting errors of omission/commission are computed for the set of parameters defined in Table 2. The results show the strong sensitivity to the minimum angle of the landslide slope α and to a lesser extent the radiometric threshold value A chosen (Figure 7). In particular we see that α greater than 15° implies a great increase in the area error. This does not indicate necessarily that the landslides are in low slopes, but rather that the precision of the DEM does not enable to pick up the small fluctuations of the topography. Precision of the DEM is thus a critical parameter of the method. Great improvements of the method are thus foreseen with the release of the high resolution global DEM from the Tandem-X mission [47]. Other parameters W, B are found to have less impacts on the results.
We found the optimum errors by minimizing the sum of commission and omission errors. The minimal errors are found to be Eo,n = 20%, Ec,n = 50%, Eo,a = 7%, Ec,a = 35% obtained for α = 15°, A = 3.5, w = 32, B = 20°. This sensitivity analysis also shows that the area errors are less than the number errors. This indicates that the method is more suitable to detect large landslides than the smaller ones. All landslides larger than 700 m2 (=28 pixels) are detected by the method.

5. Application: Landslides Triggered by the Pisco Earthquake

5.1. Inventory Validation

The methodology is now applied to all SPOT5 image pairs using the parameters obtained in the Section 4.2. High number of false alarms remained after the automatic processing (Section 4.2). A visual screening was thus necessary to remove them. False alarms are principally new roads often associated with landslide (Figure 5), small clouds not detected by the cloud algorithm at large scale (Section 3.2), small agricultural fields that are numerous on the slopes of Peru, and the presence of snow notably on the image pairs #7 and #8. All these false alarms are easily recognized visually. Moreover, using a pixel-based criterion to detect landslides often leads to pepper and salt inventories. Therefore, we gather manually the patches of changes detected that belong to the same landslide. We also get an independent validation of the algorithm over the whole dataset. Indeed, we also realize an expert-based inventory of the largest landslides by screening all the images visually. We estimate our visual detection is complete for landslides of size greater than 40 pixels. We see that all the landslides detected visually are also detected by the automatic procedure.
Over one scene (scene 1), 2 inventories are realized taking advantage of the multi-temporal acquisition (2005, 2008, 2010) (Figure 8). The 2005–2008 inventory is composed of 219 landslides (for an area of 215,000 m2). The 2008–2010 inventory is composed of 10 landslides for an area of 5,000 m2. A decrease by a factor 21 in numbers and 43 in area is therefore observed between the periods 2005–2008 and 2008–2010. The multi-temporal inventories show that not all of the landslides detected with 5 years baselines (2005–2010) are co-seismic landslides, but that the effect of the Pisco earthquake is major compared with other triggering factors in this area, and notably the rain effect. This can be explained by the high aridity of the coastal environment of Peru, with less than 3 cm a year of precipitation on average at the meteorological station of Ica (Figure 1). This multi-temporal inventory shows that large temporal baselines (2.9–6 yr), as used in this study, are not a limitation for studying the co-seismic landslides of Pisco.

5.2. Characteristics of the Landslide Inventory

The mass movements detected are of different types (Figure 9), mostly flows, rockfalls and topples. The area detected for each mass movement contains not only the source area but also sediment-transported area and debris-accumulated area (Figure 6).
The regional inventory realized using all the images is composed of 866 landslides of more than 100 m2 representing a total area of 1.29 km2 (Figure 10). These figures are rather low for a magnitude Mw 8.0 earthquake, as Keefer [3] has predicted around 105 landslides being triggered for such event. We reckon it can be explained by 4 main reasons: (1) the distance between the source and the relief is large (80 km), that is, half the length of the ruptured fault; (2) the earthquake is coming from the subduction, therefore not rupturing the surface and producing lower ground movements; (3) the coast of Peru is arid, therefore the soil cover—the part the most prone to slide—is very thin; (4) the constant seismic activity of the central coast of Peru [48] can regularly purge the slopes and therefore stabilize them.
We compute the normalized probability density function of this distribution (pdf) on the Figure 11. A power-law can be fitted over 3 orders of magnitude and displays a slope exponent of 1.63. Since all landslides with area greater than 750 m2 were also identified visually (see Section 5.1), we can conclude that this exponent is fairly robust. This power-law is similar in a mathematical sense to the second part of the inverse gamma law observed for landslide area [31]. Therefore exponents can be compared. The Pisco case displays exponent lower than previous landslide inventories [31,49] with exponents between 2.1 and 2.5 obtained for the United States, New Zealand, Italy, Guatemala and Taiwan. This means that the Pisco inventory is dominated by larger landslides compared with other seismic or rain events in various seismic and climatic contexts.
Very few landslide inventories exist in this region of the world. To our knowledge only two inventories were realized, after the 1946 and 1970 Ancash earthquakes of magnitude Ml 7.3 and 7.7 respectively. No statistics on the surface or volume of landslides exist for the 1970 earthquake [50]. On the contrary [51] shows that the volumes of the 1946 earthquake-triggered landslides are following a power law. The exponent is found to be 0.51 much lower than for other cases with exponent of 0.8 (e.g., Northridge case). This lower exponent for the volume distribution in Ancash and for the surface distribution in the Pisco case shows a common particularity of the Peruvian landslide inventories.
Both inventories are situated in regions with very thin soils due to both steep slopes and the very little vegetation cover. We suggest that the very thin soil cover in this hyperarid region of Peru explains both the few number of mass movements reported and the particular slope exponent of the distribution. The Pisco earthquake occurred on a region where low volumes of available material were susceptible to be mobilized. Indeed, steep slopes of intensely fractured basement rocks are present all along the dry canyons incising the western Cordillera in Pisco region, but the soil, fine grain or superficial sedimentary cover is very rare. Mean annual precipitations along the coastal plain or western cordillera are only of a few millimeters per year and thus vegetation is absent or focused in the valleys. The strong topographic gradient, the regional desert climate and its geology explain why the slopes are only covered with sparse and thin colluvial deposits.
On the contrary, all other published inventories focus on regions with important soil cover or low consolidated materials. For instance, the Northridge inventory (USA) is dominated by landslides occurring in very weakly cemented sedimentary rocks [52]. The highest erosion rates in Taiwan occur in terrains of weak substrate [53]. We therefore suggest that the slope exponent of the landslide distribution is dependent on the soil cover of each region and shows the prime role of the soil cover in the landslides triggered by earthquakes. This explanation is also corroborated by the different slope exponents of volume-area distributions for shallow soil-based landslides and rockslides [45].

5.3. Topographical Properties of the Landslides

We first notice that the majority of landslides is clustered in a 100 km band, from 1,000 to 4,000 m of altitude (Figure 10), where the steepest topography exists. This suggests a main control of the regional slopes on the observed landslide distribution. We therefore compute the distribution of slope gradient of inventoried landslides on Figure 12. This distribution is given for either the landslide number (Figure 12(a)) or area (Figure 12(b)). To be correctly analyzed, this distribution is normalized by the total distribution of regional slopes (Figure 12(c)), that is, we compute for each slope value by step of 5° the proportion of the total area affected by the landslide. This normalization allows us quantifying the susceptibility of the different slope values to the landslide. In particular, this process shows an increase of the susceptibility of the slopes to the landslide with slope gradient until 45°–50°, reaching 0.044%, followed by a decay for slopes greater than 50°. This decay cannot be considered totally real for two main reasons. First, the coarse resolution of the DEM is not able to pick up large slope gradients (maximum values of the DEM is 65° whereas steeper slopes exist). Second, the quasi-nadir observations of the SPOT5 images (Table 1) do not allow detecting the landslides in steeper slopes (see also Section 4.1).
We also compute the distribution of slope orientation, shown in Figure 13(a,b) for the landslide number and area respectively. An unusual aspect of the landslides triggered by this subduction earthquake is their orientation along the western cordillera relief, which is mostly perpendicular to the fault rupture (Figure 10). We should note that this particular orientation of the slopes affected by landslide does not follow the general distribution of the regional slope orientation (Figure 13(a)). Slopes oriented E/NE are more affected than the W/SW slopes, with 3 times more landslides in the N70 than in the N250 directions. We suggest that this distribution reflects the interaction of the wave field with the topography (e.g., [4,54]). Indeed, in a first order, the landslides can be considered fairly far from the earthquake source, therefore being situated approximately at an azimuth E/NE (N70) from the source. In other words, the observed landslide orientation distribution reflects that flanks facing away from the wave source are more affected than the others. This has previously been interpreted as a result of the amplifications of the ground movement in these slopes [9]. Another explanation of this particular orientation could be the co-seismic tilt, which can destabilize specific slopes. In subduction zones, the tilt associated with earthquakes is first oriented away and then toward the subduction zone as we move away for this zone [55]. This change of tilt orientation appears at a distance close to the coastline [55]. This is notably the case for the Pisco earthquake [37]. Therefore, on the area where landslides occur, the tilt is oriented toward the trench. This mechanism can thus be excluded to explain the triggering of the Pisco triggered landslides.
Our inventory also shows that larger landslides are encountered in slopes facing N320 (mean size of 7,000 m2) with secondary peaks in slopes N160 (mean size of 5,000 m2) and N80 (mean size of 4,800 m2) (Figure 13(c)). We suggest that this effect comes from the geomorphology of the area. Indeed, the main (deepest) valleys near the seismic source, collecting and channeling the erosion down to the nearby coastal region, are oriented SW-NE (Figure 10). The walls of these valleys are thus trending NE-SW. Largest denivellation are therefore encountered on these particular slope orientations. We propose that the largest mean size of the landslides in the NE-SW orientation is directly due to this systematic pattern and topographic singularity, trending perpendicular to the topography and parallel to the subduction megathrust.

6. Discussion on the Method

The method used has been developed to work with simple datasets: panchromatic images coupled with a freely and globally available DEM. Thus, this method can also be used to realize landslide inventories in others parts of the globe subjected to earthquake activity, including for past events where only panchromatic images are available. It has been successfully applied to derive the first comprehensive landslide inventory triggered by the Pisco earthquake over a wide area (27,000 km2). The detection errors compare favorably to other studies (e.g., [21,23]) using other types of datasets (VHR, multi-spectral) in various areas. Commission errors are however larger with our method (50%) than for other studies using multispectral data (20%–25%). However, another study using panchromatic data [34] shows branching factors (ratio between the false alarms and the detected landslides) up to 3.6, corresponding to commission error of 78%. Therefore the 50% of commission rate obtained here can certainly be explained by various factors including the use of panchromatic versus multi-spectral data and the coarse resolution of the DEM [23]. Results of higher quality could probably be obtained if a higher resolution DEM was available.
A major lesson from this application is that pixel-based methods coupled with contextual analysis are working fine on panchromatic images in arid environment. Compared with other methods using panchromatic images (e.g., [34]), scores of detection obtained here are better (80% against 73% for the number of landslides, 93% against 88% of the landslide extent), and false alarm rates are similar (50%). It must be noted that the latter method has been applied to other climatic context and with DEM of better resolution (10 m). It is thus difficult to compare directly the scores obtained. A comparative study of the various methods with similar datasets would be very interesting for this purpose.
The region affected by the Pisco earthquake is arid with limited vegetation. We can therefore wonder what is the applicability of the developed method to other climatic and vegetation contexts. It is first important to note that the precise co-registration is a key step to obtain high score of detection. The co-registration is here realized with a precision better than half a pixel, thanks to the selection of successive images with similar incidence angles (less than 7° or B/H ratio below 0.15 (Table 1)). This prevents the histogram of the difference image from stretching out and thus allows a better definition of the detection thresholds (A). This enables both the higher detection rate of false alarms and the better definition of the limits of the landslides. This good detection of the landslide limits is also enhanced by the climatic context in the area studied. Indeed the arid environment keeps the areas around the landslides unchanged in the time interval between the two image acquisitions. Therefore the landslide limits are very clearly marked in the difference image. On areas covered by vegetation, the method would certainly be of lower efficiency. For these areas, the method must be modified for use with multispectral data. Indeed, multispectral data enable the computation of Normalized Difference Vegetation Index (NDVI) images, and changes of vegetation-cover can be detected between the two acquisitions (e.g., [22]). Therefore, only the change detection step (Section 3.3) would differ. The process of false alarm removal would be unchanged.
The method developed can be summarized in 5 steps: (1) orthorectification, co-registration; (2) estimation of the method parameters on a subset of the images where a manual inventory has been realized; (3) application to the whole images; (4) visual removal of the remaining false positives; (5) gathering of pixel groups belonging to the same landslide. Application of this method to other context would thus require the estimation of various parameters described in Table 2. The sensitivity analysis realized here (Section 4.2) shows that only 2 of them have large effects on the obtained scores (A and α). Therefore, the step (2) can be realized by running the sensitivity analysis over a representative sample of the total images, letting only A and α evolve. Due to its definition (Equation (5)), the range of evolution of A does not depend on the histogram of the difference image. Values of A between 3 and 5 are good first estimations. On the contrary, the range of evolution of α must be chosen depending on the DEM used and the regional slope characteristics. Indeed, as pointed out in Section 4.2, a main limitation of the method comes from the coarse resolution of the DEM, which precludes the precise computation of the local slopes. Therefore taking a large α will remove some landslides but decrease the false alarm rate. In this study we preferred using a smaller α to keep the score of detection the largest as possible. The drawback of our choice impacts on the false alarm rate (50%).
The two latest steps of the method (removal of the remaining false positives and gathering of pixel groups) have been realized manually. The false alarm removal is not a big task, as only 900 objects over 1800 must be removed, and these false alarms are easily recognized visually. This step can become tedious when more objects are kept (for instance the Wenchuan earthquake triggered about 60,000 landslides [5]). However, checking visually the quality of the obtained inventory is a necessary step. Concerning the merging of pixel groups, our experience on the Pisco earthquake teaches us that this step cannot be automated. For instance Figures 6 and 9(d) show two examples of instabilities where merging was necessary (groups of pixels of the same landslide are not connected). The three groups of pixels of the Figure 9(d) could be merged automatically because they are all aligned and situated along the same thalweg. However, the resolution of the DEM does not enable to detect precisely the small undulations of the topography. On Figure 6, the merging is less obvious for several reasons: (1) a criterion of distance must be defined to know which groups must be gathered; (2) we do not really know if all the groups of pixels were triggered exactly at the same time, and therefore if they can be gathered or not. We thus see that the step of merging is subjective and cannot be automated.

7. Conclusions

We developed a method for the detection of landslides in optical satellite images, with the purpose to cover extensive areas in arid environment. This method is adapted to the following dataset: high resolution (5 m) panchromatic SPOT5 images, coupled with coarse resolution (posting of 30 m) but freely available DEM (GDEM). The method enables to extract both the number and the area of the landslides. Score of detection reaches 80% of the landslides representing 93% of their area. The detection score, evaluated over two sub-regions where a detailed inventory has been manually obtained, is similar to that of previous studies using satellite images and DEM of better resolutions. The drawback of the method is the resulting high level of false alarms (50%), mainly due to the low resolution of the DEM that prevents picking the small undulations of the topography. However, these false alarms are easily identified and removed visually.
The application of this method to eight pairs of images enables us to realize the first comprehensive inventory of the landslides triggered by the Pisco earthquake. A total of 866 landslides larger than 100 m2 were detected covering an area of 1.29 km2. This inventory does not encompass a band of about 10 km along the coast of Peru, where no images were available or clouds masked the coastal relief. However, this number is far lower than predicted by generic laws [3]. The size/number relation is found to follow a power law with slope exponent of 1.63, lower than previously found for other case studies. The triggering by earthquakes in the coastal part of Peru seems therefore to be fairly particular, with low level of triggering but relatively bigger landslides than in other seismic/climatic contexts worldwide. We suggest that this low slope exponent of the power law is related to the type of the soil cover of the study area. Another interesting characteristic of our inventory is the greater average size of landslides oriented parallel to the earthquake rupture. We argue that this particularity highlights the main role of the geomorphology in the propagation of landslides.
A first analysis of the obtained database in relation with the topography shows an increase of the susceptibility of the slope to landslide with the slope gradient, and a site effect with flanks facing away from the wave source more susceptible than other flanks directions. Future work will include analysis of this database with the ground movements and with the geological settings, to decipher all these particularities. Because of the simple dataset required to use the proposed methodology, the method can be used for a large set of applications, mostly for covering wide areas, which is of particular interest for co-seismic landslides studies.

Acknowledgments

We particularly thank the ISIS program from the French Space Agency (CNES), who permits the acquisitions of all the SPOT5 images. We also acknowledge support from the CNES/TOSCA. The earthquake source model has been kindly provided by Anthony Sladen. Helpful discussions with Bernard Pirletta also benefits this study. We thank four anonymous reviewers and the scientific editor for very useful comments and suggestions.
  • Conflict of InterestThe authors declare on conflict of interest.

References

  1. Keefer, D.K. Landslides caused by earthquakes. Geol. Soc. Amer. Bull 1984, 95, 406–421. [Google Scholar]
  2. Petley, D. Global patterns of loss of life from landslides. Geology 2012, 40, 927–930. [Google Scholar]
  3. Keefer, D.K. Investigating landslides caused by earthquakes—A historical review. Surv. Geophys. 2002, 23, 473–510. [Google Scholar]
  4. Sato, H.P.; Hasegawa, H.; Fujiwara, S.; Tobita, M.; Koarai, M.; Une, H.; Iwahashi, J. Interpretation of landslide distribution triggered by the 2005 Northern Pakistan earthquake using SPOT 5 imagery. Landslides 2006, 4, 113–122. [Google Scholar]
  5. Gorum, T.; Fan, X.; van Westen, C.J.; Huang, R.Q.; Xu, Q.; Tang, C.; Wang, G. Distribution pattern of earthquake-induced landslides triggered by the 12 May 2008 Wenchuan earthquake. Geomorphology 2011, 133, 152–167. [Google Scholar]
  6. Sepulveda, S.A.; Murphy, W.; Jibson, R.W.; Petley, D.N. Seismically induced rock slope failures resulting from topographic amplification of strong ground motions: The case of Pacoima Canyon, California. Eng. Geol. 2005, 80, 336–348. [Google Scholar]
  7. Sidle, R.; Ochiai, H. Landslides: Processes, prediction, and land use. Water Resour. Monogr. 2006, 18, 1–312. [Google Scholar]
  8. Lin, G.W.; Chen, H.; Chen, Y.H.; Horng, M.J. Influence of typhoons and earthquakes on rainfall-induced landslides and suspended sediments discharge. Eng. Geol. 2008, 97, 32–41. [Google Scholar]
  9. Meunier, P.; Hovius, N.; Haines, J.A. Topographic site effects and the location of earthquake induced landslides. Earth Planet. Sci. Lett. 2008, 275, 221–232. [Google Scholar]
  10. Tang, C.; Zhu, J.; Qi, X.; Ding, J. Landslides induced by the Wenchuan earthquake and the subsequent strong rainfall event: A case study in the Beichuan area of China. Eng. Geol. 2011, 122, 22–33. [Google Scholar]
  11. Del Gaudio, V.; Coccia, S.; Wasowski, J.; Gallipoli, M.R.; Mucciarelli, M. Detection of directivity in seismic site response from microtremor spectral analysis. Nat. Hazards Earth Syst. Sci. 2008, 8, 751–762. [Google Scholar]
  12. Bozzano, F.; Lenti, L.; Martino, S.; Paciello, A.; Mugnozza, G.S. Self-excitation process due to local seismic amplification responsible for the reactivation of the Salcito landslide (Italy) on 31 October 2002. J. Geophys. Res. 2008. [Google Scholar] [CrossRef]
  13. Rodriguez, C.; Bommer, J.; Chandler, R. Earthquake-induced landslides: 1980–1997. Soil Dyn. Earthq. Eng. 1999, 18, 325–346. [Google Scholar]
  14. Galli, M.; Ardizzone, F.; Cardinali, M.; Guzzetti, F.; Reichenbach, P. Comparing landslide inventory maps. Geomorphology 2008, 94, 268–289. [Google Scholar]
  15. Hervas, J.; Barredo, J.I.; Rosin, P.L.; Pasuto, A.; Mantovani, F.; Silvano, S. Monitoring landslides from optical remotely sensed imagery: The case history of Tessina landslide, Italy. Geomorphology 2003, 54, 63–75. [Google Scholar]
  16. Barlow, J.; Franklin, S.; Martin, Y. High spatial resolution satellite imagery, DEM derivatives, and image segmentation for the detection of mass wasting processes. Photogramm. Eng. Remote Sensing 2006, 72, 687–692. [Google Scholar]
  17. Lee, S.; Lee, M.J. Detecting landslide location using KOMPSAT 1 and its application to landslide-susceptibility mapping at the Gangneung area, Korea. Adv. Space Res. 2006, 38, 2261–2271. [Google Scholar]
  18. Borghuis, A.M.; Chang, K.; Lee, H.Y. Comparison between automated and manual mapping of typhoon-triggered landslides from SPOT-5 imagery. Int. J. Remote Sens. 2007, 28, 1843–1856. [Google Scholar]
  19. Martha, T.R.; Kerle, N.; Jetten, V.; van Westen, C.J.; Kumar, K.V. Characterising spectral, spatial and morphometric properties of landslides for semi-automatic detection using object-oriented methods. Geomorphology 2010, 116, 24–36. [Google Scholar]
  20. Mondini, A.; Guzzetti, F.; Reichenbach, P.; Rossi, M.; Cardinali, M.; Ardizzone, F. Semi-automatic recognition and mapping of rainfall induced shallow landslides using optical satellite images. Remote Sens. Environ. 2011, 115, 1743–1757. [Google Scholar]
  21. Lu, P.; Stumpf, A.; Kerle, N.; Casagli, N. Object-oriented change detection for landslide rapid mapping. IEEE Geosci. Remote Sens. Lett. 2011, 8, 701–705. [Google Scholar]
  22. Martha, T.R.; Kerle, N.; van Westen, C.J.; Jetten, V.; Kumar, K.V. Segment optimization and data-driven thresholding for knowledge-based landslide detection by object-based image analysis. IEEE Trans. Geosci. Remote Sens. 2011, 49, 4928–4943. [Google Scholar]
  23. Stumpf, A.; Kerle, N. Object-oriented mapping of landslides using random forests. Remote Sens. Environ. 2011, 115, 2564–2577. [Google Scholar]
  24. Holbling, D.; Fureder, P.; Antolini, F.; Cigna, F.; Casagli, N.; Lang, S. A semi-automated object-based approach for landslide detection validated by persistent scatterer interferometry measures and landslide inventories. Remote Sens. 2012, 4, 1310–1336. [Google Scholar]
  25. Guzzetti, F.; Mondini, A.C.; Cardinali, M.; Fiorucci, F.; Santangelo, M.; Chang, K.T. Landslide inventory maps: New tools for an old problem. Earth-Sci. Rev. 2012, 112, 42–66. [Google Scholar]
  26. McKean, J.; Roering, J. Objective landslide detection and surface morphology mapping using high-resolution airborne laser altimetry. Geomorphology 2004, 57, 331–351. [Google Scholar]
  27. Czuchlewski, K.R.; Weissel, J.K.; Kim, Y. Polarimetric synthetic aperture radar study of the Tsaoling landslide generated by the 1999 Chi-Chi earthquake, Taiwan. J. Geophys. Res. 2003. doi: 200310.1029/2003JF000037. [Google Scholar]
  28. Yonezawa, C.; Watanabe, M.; Saito, G. Polarimetric decomposition analysis of ALOS PALSAR observation data before and after a landslide event. Remote Sens. 2012, 4, 2314–2328. [Google Scholar]
  29. Delacourt, C.; Allemand, P.; Berthier, E.; Raucoules, D.; Casson, B.; Grandjean, P.; Pambrun, C.; Varel, E. Remote-sensing techniques for analysing landslide kinematics: A review. Bull. Soc. Geol. France 2007, 178, 89–100. [Google Scholar]
  30. Riedel, B.; Walther, A. InSAR processing for the recognition of landslides. Adv. Geosci. 2008, 14, 189–194. [Google Scholar]
  31. Malamud, B.D.; Turcotte, D.L.; Guzzetti, F.; Reichenbach, P. Landslide inventories and their statistical properties. Earth Surf. Process. Landf. 2004, 29, 687–711. [Google Scholar]
  32. Tatard, L.; Grasso, J.R.; Helmstetter, A.; Garambois, S. Characterization and comparison of landslide triggering in different tectonic and climatic settings. J. Geophys. Res. 2010. doi: 201010.1029/2009JF001624. [Google Scholar]
  33. Rosin, P.L.; Hervas, J. Remote sensing image thresholding methods for determining landslide activity. Int. J. Remote Sens. 2005, 26, 1075–1092. [Google Scholar]
  34. Martha, T.R.; Kerle, N.; van Westen, C.J.; Jetten, V.; Vinod Kumar, K. Object-oriented analysis of multi-temporal panchromatic images for creation of historical landslide inventories. ISPRS J. Photogramm. 2012, 67, 105–119. [Google Scholar]
  35. Nichol, J.; Wong, M.S. Satellite remote sensing for detailed landslide inventories using change detection and image fusion. Int. J. Remote Sens. 2005, 26, 1913–1926. [Google Scholar]
  36. Zavala, B.; Hermanns, R.; Valderrama, P.; Costa, C.; Rosado, M. Procesos geologicos e intensidad macrosismica Inqua del sismo de Pisco del 15/08/2007, Peru. Revista de la Asociacion Geologica Argentina 2009, 65, 760–779. [Google Scholar]
  37. Sladen, A.; Tavera, H.; Simons, M.; Avouac, J.P.; Konca, A.O.; Perfettini, H.; Audin, L.; Fielding, E.J.; Ortega, F.; Cavagnoud, R. Source model of the 2007 Mw 8.0 Pisco, Peru earthquake: Implications for seismogenic behavior of subduction megathrusts. J. Geophys. Res. 2010. doi: 201010.1029/2009JB006429. [Google Scholar]
  38. Tavera; Bernal, I. The pisco (Peru) earthquake of 15 August 2007. Seismol. Res. Lett. 2008, 79, 510–515. [Google Scholar]
  39. Tavera, H.; Bernal, I.; Strasser, F.O.; Arango-Gaviria, M.C.; Alarcon, J.E.; Bommer, J.J. Ground motions observed during the 15 August 2007 Pisco, Peru, earthquake. Bull. Earthq. Eng. 2009, 7, 71–111. [Google Scholar]
  40. Berthier, E.; Vadon, H.; Baratoux, D.; Arnaud, Y.; Vincent, C.; Feigl, K.; Remy, F.; Legresy, B. Surface motion of mountain glaciers derived from satellite optical imagery. Remote Sens. Environ. 2005, 95, 14–28. [Google Scholar]
  41. Leprince, S.; Barbot, S.; Ayoub, F.; Avouac, J.P. Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Trans. Geosci. Remote Sens. 2007, 45, 1529–1557. [Google Scholar]
  42. ASTER GDEM Validation Team. ASTER Global DEM Validation, Summary Report; Technical Report; ASTER GDEM Validation Team (METI/ERSDAC, NASA/LPDAAC, USGS/EROS); 2009. [Google Scholar]
  43. Singh, A. Digital change detection techniques using remotely-sensed data. Int. J. Remote Sens. 1989, 10, 989–1003. [Google Scholar]
  44. Im, J.; Jensen, J.R.; Tullis, J.A. Object-based change detection using correlation image analysis and image segmentation. Int. J. Remote Sens. 2008, 29, 399–423. [Google Scholar]
  45. Larsen, I.J.; Montgomery, D.R.; Korup, O. Landslide erosion controlled by hillslope material. Nat. Geosci. 2010, 4, 247–251. [Google Scholar]
  46. Zhang, Q.; Wang, J.; Peng, X.; Gong, P.; Shi, P. Urban built-up land change detection with road density and spectral information from multi-temporal Landsat TM data. Int. J. Remote Sens. 2002, 23, 3057–3078. [Google Scholar]
  47. Krieger, G.; Moreira, A.; Fiedler, H.; Hajnsek, I.; Werner, M.; Younis, M.; Zink, M. TanDEM-X: A satellite formation for high-resolution SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2007, 45, 3317–3341. [Google Scholar]
  48. Dorbath, L.; Cisternas, A.; Dorbath, C. Assessment of the size of large and great historical earthquakes in Peru. Bull. Seismol. Soc. Amer. 1990, 80, 551–576. [Google Scholar]
  49. Stark, C.P.; Hovius, N. The characterization of landslide size distributions. Geophys. Res. Lett. 2001, 28, 1091–1094. [Google Scholar]
  50. Plafker, G.; Ericksen, G.E.; Concha, J.F. Geological aspects of the May 31, 1970, Peru earthquake. Bull. Seismol. Soc. Amer. 1971, 61, 543–578. [Google Scholar]
  51. Kampherm, T.S. Landslides Triggered by the 1946 Ancash Earthquake (Peru) and Geologic Controls on the Mechanisms of Initial Rock Slope Failure. M.Sc. Thesis, University Waterloo, Waterloo, ON, Canada. 2009. [Google Scholar]
  52. Parise, M.; Jibson, R.W. A seismic landslide susceptibility rating of geologic units based on analysis of characteristics of landslides triggered by the 17 January, 1994 Northridge, California earthquake. Eng. Geol. 2000, 58, 251–270. [Google Scholar]
  53. Dadson, S.J.; Hovius, N.; Chen, H.; Dade, W.B.; Hsieh, M.L.; Willett, S.D.; Hu, J.C.; Horng, M.J.; Chen, M.C.; Stark, C.P.; et al. Links between erosion, runoff variability and seismicity in the Taiwan orogen. Nature 2003, 426, 648–651. [Google Scholar]
  54. Qi, S.; Xu, Q.; Lan, H.; Zhang, B.; Liu, J. Spatial distribution analysis of landslides triggered by 2008.5.12 Wenchuan Earthquake, China. Eng. Geol. 2010, 116, 95–108. [Google Scholar]
  55. Kanda, R.V.; Simons, M. Practical implications of the geometrical sensitivity of elastic dislocation models for field geologic surveys. Tectonophysics 2012, 560–561, 94–104. [Google Scholar]
Figure 1. Summary of the SPOT5 image pairs used in this study: (a) temporal baseline and (b) corresponding images coverage. The 1 m contour of slip distribution on the fault during the Pisco earthquake (orange line) is extracted from Sladen et al.[37]. The white dots represent the landslides identified with the field inventory [36].
Figure 1. Summary of the SPOT5 image pairs used in this study: (a) temporal baseline and (b) corresponding images coverage. The 1 m contour of slip distribution on the fault during the Pisco earthquake (orange line) is extracted from Sladen et al.[37]. The white dots represent the landslides identified with the field inventory [36].
Remotesensing 05 02590f1
Figure 2. (a) Diagram showing the radiometry of each pixel for the pre- and post-event images (black points). The color lines represent the iso-level of pixel density. We notice that most pixels are situated along the line x = y. This shows that most pixels have similar radiometric values in the pre- and post-image. (b) Histogram of the difference image once the images have been normalized.
Figure 2. (a) Diagram showing the radiometry of each pixel for the pre- and post-event images (black points). The color lines represent the iso-level of pixel density. We notice that most pixels are situated along the line x = y. This shows that most pixels have similar radiometric values in the pre- and post-image. (b) Histogram of the difference image once the images have been normalized.
Remotesensing 05 02590f2
Figure 3. Different steps of the detection of large scale changes. First, the image of correlation C(x, y) is calculated (a). We note that low correlation corresponds to area with a high density of changes (black points). (b) Then, the objects situated on low slopes are detected (red points). (c) The area is finally expanded to connected points with similar or lower correlation than the objects detected in step (b). The remaining blue points are candidates for being landslides.
Figure 3. Different steps of the detection of large scale changes. First, the image of correlation C(x, y) is calculated (a). We note that low correlation corresponds to area with a high density of changes (black points). (b) Then, the objects situated on low slopes are detected (red points). (c) The area is finally expanded to connected points with similar or lower correlation than the objects detected in step (b). The remaining blue points are candidates for being landslides.
Remotesensing 05 02590f3
Figure 4. Valleys detected by the automatic procedure (red points). Objects identified as inside these valleys are removed (yellow points). Remaining blue points are candidates for being landslides.
Figure 4. Valleys detected by the automatic procedure (red points). Objects identified as inside these valleys are removed (yellow points). Remaining blue points are candidates for being landslides.
Remotesensing 05 02590f4
Figure 5. Results of the automatic classification for one subset image (area A2) where an expert based inventory has been realized.
Figure 5. Results of the automatic classification for one subset image (area A2) where an expert based inventory has been realized.
Remotesensing 05 02590f5
Figure 6. Image subsets over the same area (of 1.2 × 1.4 km2) at different dates surrounding the Pisco earthquake. The red contours correspond to the 2 new landslides detected by the automatic procedure described in this paper: one occurring between 26 May 2005 and 26 July 2007 (middle image), and the largest one occurring between 26 July 2007 and 30 May 2011 (right image). The yellow point corresponds to a 150 m3 landslide assumed to be triggered by the earthquake in the field inventory (13.17°S 75.64°W), which fits well with the first landslide detected by the automatic method but before the earthquake.
Figure 6. Image subsets over the same area (of 1.2 × 1.4 km2) at different dates surrounding the Pisco earthquake. The red contours correspond to the 2 new landslides detected by the automatic procedure described in this paper: one occurring between 26 May 2005 and 26 July 2007 (middle image), and the largest one occurring between 26 July 2007 and 30 May 2011 (right image). The yellow point corresponds to a 150 m3 landslide assumed to be triggered by the earthquake in the field inventory (13.17°S 75.64°W), which fits well with the first landslide detected by the automatic method but before the earthquake.
Remotesensing 05 02590f6
Figure 7. Effect of the minimum slope angle α (top) and the threshold parameter A (bottom) on the error of omission (black) and commission (red) for landslide number (left) and area (right). Error bars are calculated by the standard deviation of the omission and commission when parameters other than α or A vary in the range defined by the Table 2.
Figure 7. Effect of the minimum slope angle α (top) and the threshold parameter A (bottom) on the error of omission (black) and commission (red) for landslide number (left) and area (right). Error bars are calculated by the standard deviation of the omission and commission when parameters other than α or A vary in the range defined by the Table 2.
Remotesensing 05 02590f7
Figure 8. Multi-temporal inventories of the pairs #1. The 2005–2008 inventory is in blue and the 2008–2010 inventory is in yellow.
Figure 8. Multi-temporal inventories of the pairs #1. The 2005–2008 inventory is in blue and the 2008–2010 inventory is in yellow.
Remotesensing 05 02590f8
Figure 9. Example of landslides detected by the methodology. Left and right subplots correspond to pre- and post-event images respectively.
Figure 9. Example of landslides detected by the methodology. Left and right subplots correspond to pre- and post-event images respectively.
Remotesensing 05 02590f9
Figure 10. Field (white circle) and remotely-sensed (red square) landslide inventories. The 1 m contour of slip distribution on the fault during the Pisco earthquake is shown with dashed dot yellow line. The red contour is the coverage of the SPOT5 images. Iso-altitudes are lined in black at 1,000 m and 4,000 m levels.
Figure 10. Field (white circle) and remotely-sensed (red square) landslide inventories. The 1 m contour of slip distribution on the fault during the Pisco earthquake is shown with dashed dot yellow line. The red contour is the coverage of the SPOT5 images. Iso-altitudes are lined in black at 1,000 m and 4,000 m levels.
Remotesensing 05 02590f10
Figure 11. Probability density function (pdf) of landslide area.
Figure 11. Probability density function (pdf) of landslide area.
Remotesensing 05 02590f11
Figure 12. (a) Distribution of landslide numbers as a function of slopes; (b) Distribution of landslide area as a function of slopes; (c) Proportion of slopes affected by landslide, i.e., the percentage of area of each slopes covered by landslides.
Figure 12. (a) Distribution of landslide numbers as a function of slopes; (b) Distribution of landslide area as a function of slopes; (c) Proportion of slopes affected by landslide, i.e., the percentage of area of each slopes covered by landslides.
Remotesensing 05 02590f12
Figure 13. (a) Distribution (given in %) of landslide numbers as a function of slope orientation (solid line). The slope distribution of the study area is shown as reference with a dashed line. (b) Distribution of landslide area (given in %) as a function of slope orientation. (c) average landslide area (m2) as a function of the slope orientation.
Figure 13. (a) Distribution (given in %) of landslide numbers as a function of slope orientation (solid line). The slope distribution of the study area is shown as reference with a dashed line. (b) Distribution of landslide area (given in %) as a function of slope orientation. (c) average landslide area (m2) as a function of the slope orientation.
Remotesensing 05 02590f13
Table 1. SPOT5 images index table.
Table 1. SPOT5 images index table.
Pair #Pre-Event DatePost-Event DateB/HMean Incidence Angle
12005/04/042008/03/050.03072.6 °
2008/03/052010/04/240.00762.2°
22005/11/242009/07/120.00942.3°
2009/07/122010/04/240.01282.1°
32003/10/272008/08/290.00632.5°
42005/05/262011/05/190.0021
52003/06/192007/10/220.02486.5°
62004/04/202008/05/060.157824.3°
72005/05/262011/05/190.0021
82007/07/262011/05/300.001422°
Table 2. Definition and values of the method parameters.
Table 2. Definition and values of the method parameters.
NameDescriptionMin ValueMax ValueSectionOptimum
Aparameter related to the detection threshold34.53.33.5
wsize (in pixel) of the correlation window161283.4.2.32
αminimum slope of landslides10°25°3.4.2.15°
Bmaximum angle between road and slope90°3.4.4.20°

Share and Cite

MDPI and ACS Style

Lacroix, P.; Zavala, B.; Berthier, E.; Audin, L. Supervised Method of Landslide Inventory Using Panchromatic SPOT5 Images and Application to the Earthquake-Triggered Landslides of Pisco (Peru, 2007, Mw8.0). Remote Sens. 2013, 5, 2590-2616. https://doi.org/10.3390/rs5062590

AMA Style

Lacroix P, Zavala B, Berthier E, Audin L. Supervised Method of Landslide Inventory Using Panchromatic SPOT5 Images and Application to the Earthquake-Triggered Landslides of Pisco (Peru, 2007, Mw8.0). Remote Sensing. 2013; 5(6):2590-2616. https://doi.org/10.3390/rs5062590

Chicago/Turabian Style

Lacroix, Pascal, Bilberto Zavala, Etienne Berthier, and Laurence Audin. 2013. "Supervised Method of Landslide Inventory Using Panchromatic SPOT5 Images and Application to the Earthquake-Triggered Landslides of Pisco (Peru, 2007, Mw8.0)" Remote Sensing 5, no. 6: 2590-2616. https://doi.org/10.3390/rs5062590

APA Style

Lacroix, P., Zavala, B., Berthier, E., & Audin, L. (2013). Supervised Method of Landslide Inventory Using Panchromatic SPOT5 Images and Application to the Earthquake-Triggered Landslides of Pisco (Peru, 2007, Mw8.0). Remote Sensing, 5(6), 2590-2616. https://doi.org/10.3390/rs5062590

Article Metrics

Back to TopTop