Evaluation of Unsupervised Change Detection Methods Applied to Landslide Inventory Mapping Using ASTER Imagery

Natural hazards include a wide range of high-impact phenomena that affect socioeconomic and natural systems. Landslides are a natural hazard whose destructive power has caused a significant number of victims and substantial damage around the world. Remote sensing provides many data types and techniques that can be applied to monitor their effects through landslides inventory maps. Three unsupervised change detection methods were applied to the Advanced Spaceborne Thermal Emission and Reflection Radiometer (Aster)-derived images from an area prone to landslides in the south of Mexico. Linear Regression (LR), Chi-Square Transformation, and Change Vector Analysis were applied to the principal component and the Normalized Difference Vegetation Index (NDVI) data to obtain the difference image of change. The thresholding was performed on the change histogram using two approaches: the statistical parameters and the secant method. According to previous works, a slope mask was used to classify the pixels as landslide/No-landslide; a cloud mask was used to eliminate false positives; and finally, those landslides less than 450 m2 (two Aster pixels) were discriminated. To assess the landslide detection accuracy, 617 polygons (35,017 pixels) were sampled, classified as real landslide/No-landslide, and defined as ground-truth according to the interpretation of color aerial photo slides to obtain omission/commission errors and Kappa coefficient of agreement. The results showed that the LR using NDVI data performs the best results in landslide detection. Change detection is a suitable technique that can be applied for the landslides mapping and we think that it can be replicated in other parts of the world with results similar to those obtained in the present work.


Introduction
The instability of land slopes is a natural risk with a high destruction ability; landslides and other related complex movements occur every day around the world and the victims and damage caused are substantial [1][2][3][4], and the risks from natural disasters have become a growing concern [5][6][7].Significant efforts have been made to collect, record, and analyze information on the occurrence and impacts of natural disasters around the world [8].The development of natural disasters databases is essential for risk management, since it facilitates the assessment of their natural and social impacts and the vulnerability of regions on different scales [9].In particular, geodatabases that manage information on landslides, including inventories and thematic data, represent a powerful tool at local, regional, national, or continental management and organization levels [10,11].
A fundamental component of landslide geodatabases is the landslide inventory maps, since they provide historical information of past events, that, when analyzed in combination with information about geoenvironmental conditioning or landslides triggers factors, can produce essential and useful information expressed in spatial models, such as landslide susceptibility maps or landslide risk maps.
Landslide inventory maps have multiple purposes: showing the location and type of landslides in a region, demonstrating the effects of trigger events (earthquake, heavy rain, or rapid thaw), integrating statistics about fault frequency into slopes, or providing relevant information to build models of susceptibility or slip risk [5,[12][13][14][15][16][17][18][19].Thus, this kind of information must be available to decision-makers to be considered in territorial management and the planning of the development of civil protection plans and policies to face the wide range and accelerated change of the landscapes that are currently underway in several regions around the world.This can be achieved through the application of appropriate strategies aimed to assess the transformation of the space over the time, to prevent and mitigate losses due to disasters [9], eventually to support the management of emergencies related to this phenomenon [1], and develop suitable landscape and urban planning, environmental resources management, risk assessment, or disaster management [20][21][22].
The geographical location of the Guerrero State in México makes it susceptible to hydrometeorological events that, in combination with topographic and geological conditions, increase the threat and vulnerability of populations to landslides [4].
The availability of new remote sensing technologies can facilitate the production of landslide inventory maps, where the application of change detection techniques using multispectral data can be a useful approach for mapping construction and for defining criteria to evaluate its quality [23].There are some new approaches for landslides identification, where visual exploration and field surveys are the traditional methods used for the generation of inventory maps [24].However, the visual exploration of images for assessing the quality of results in remote sensing studies is a subjective evaluation that can be inconvenient because significant work and monetary investment are required.Also detailed field studies are time-consuming and challenging to perform on inaccessible terrain [25].
Some previous studies state that the use of remote sensing has not been adequately exploited in landslide studies, and only a limited number of researchers have used multispectral images to assess landslide activity [26,27].Stereoscopic air-photo interpretation continues to be the most frequent remote sensing application technique in the mapping of landslide and monitoring its characteristics [27].Singhroy and Molch [28] and Singhroy [29] mentioned two typical approaches to determine the characteristics of landslides from remotely sensed data.The initial approach can be achieved with satellite or aerial images of the visible and infrared regions to determine qualitative characteristics, such as the spatial distribution of the landslides.The second approach complements the initial characterization in estimating the dimensions such as the length, width, depth, and local slope along the detected mass movement, usually using stereo Synthetic Aperture Radar (SAR), interferometric SAR (InSAR), and topographic information.
In other studies, satellite imagery was used to conduct landslide detection and inventory mapping [3,[30][31][32][33][34][35][36][37].Singhroy et al. [30] used RadarSat and Landsat Tematic Mapper (TM) to apply interferometry, georeference, and data fusion to obtain a good representation of the changes in elevation and slope to identify landslides.Gupta and Saha [31] integrated maps of landslide inventories with small spatial extent areas through high-resolution Indian Remote Sensing (IRS) images, by applying georeferencing and visual interpretation.Cheng et al. [32] used multi-temporal Satellite Pour l'Observation de la Terre (SPOT) images to perform spectral ratioing (IR/R), image differencing, and image thresholding, aiming to achieve change detection and assess land use changes pre-and post-landslide occurrence.
Aksoy et al. [33] developed a segmentation process using Landsat Enhanced Thematic Mapper Plus (ETM+) images to identify landslides.In their analyses, they considered NDVI, slope angle, curvature, brightness, mean band blue, asymmetry, shape index, length/width ratio, gray level co-occurrence matrix, and mean difference of infrared band.Rau et al. [34] used high spatial pan-sharpened FORMOSAT-2 imagery and a Digital Elevation Model (DEM) with an object-oriented analysis to develop a scheme for landslide identification using multilevel segmentation and a hierarchical semantic network.The threshold values were determined semi-automatically by statistical estimation from a few training samples.Pour and Hashim [35] used Landsat Operational Land Imager (OLI) and Phased Array type L-band Synthetic Aperture Radar/Advanced Land Observation Satellite (PALSAR-2/ALOS-2) data to map geologic structural and topographical features for identification of high potential risk and susceptible zones for landslides and flooding areas by analyzing major geological structures and characterizations of lineaments, drainage patterns and lithology.The analytical hierarchy process approach was used for landslide susceptibility mapping.Hashim et al. [36] identified landslides determined by tracking changes in vegetation pixel data using Landsat OLI images that were acquired before and after a flooding event.Si et al. [37] proposed a semi-automatic detection method for landslides at a regional scale.They applied the Change Vector Analysis (CVA) technique to Landsat OLI images.Then, the image difference was thresholded by statistical parameters to obtain areas of change and possible landslides.Based on the candidate regions, an analysis of susceptibility to landslides was carried out where the topographic, geomorphological, geological, and land cover factors were considered.This analysis showed the possibility of a landslide occurring in each candidate region, thus providing a criterion for the identification of landslides.
Change detection is a technique applied in the landslide identification, by comparing its spatial representation at two points in time and detecting differences between the state of specific characteristics of a phenomenon or elements [22,[38][39][40].Typical approaches for change detection based on multi-spectral image data analysis have been developed: image subtraction, principal component analysis, ratioing, change vector analysis, linear image regression, and spectral features variance [41][42][43][44][45].These approaches show that change detection techniques are suitable for understanding, evaluating and addressing the effect of changes in land use and land coverage, and so we think it would be a useful and practical approach to detect past landslide events.
In change detection using unsupervised techniques, the spectral or derivative information for each date is compared without relying on any additional information as ground-truth [38,46].This is an advantage given that the ground-truth integration is a difficult, slow, and expensive task and the information is often not available.
The classical method for thresholding and categorizing change and no-change pixels is based on the assumption that only a few changes occur between the analyzed dates.Therefore, the unchanged pixels are distributed around the mean of the frequency distribution, whereas the change pixels are distributed in the tails [47] and separated from the mean by a certain number of standard deviations of the distribution.Usually, the definition of the threshold values is based on the statistical parameters (mean and standard deviation) of the difference image, interactively adjusting the number of standard deviations and evaluating the results until they meet the established accuracy and error values [43,47,48].In previous works focused on detecting landslides, Hervás and Rosin [49] proposed an automatic threshold definition through the intensity value of the pixel in the difference image, which corresponds to the point in the frequency histogram (usually unimodal) of the maximum distance to the secant line between the highest and lowest points of the frequency distribution [40].
This work aims to perform three different of change detection techniques: Chi-Square Transformation (CST), Linear Regression (LR), and Change Vector Analysis (CVA), on Principal Components (PC) and NDVI Aster-derived images and using slope, clouds, and isolated pixels masks to automatically generate a cartographic inventory of landslides.
To fulfill the proposed objective, the following stages are proposed: the first stage is the application of the CST, LR, and CVA change detection methods to two kinds of input data: PC and NDVI, to obtain a single difference image for each data type between the dates included in the study.The second stage is the analysis of the change images using a fully automatic detection of changes, carried out using the statistical parameters of change images and by applying the secant technique previously mentioned, to obtain a binary change/no-change map.In this stage, masks of slopes, clouds, and isolated pixels were used to mask the resulting binary maps and thereby obtain the final cartographic inventory of landslides.The final stage is assessing the obtained results and validating them with the collected ground-truth information through confusion matrices.

Study Area
The study area was located in the central part of the Guerrero State in México, and covered an area of about 3300 km 2 in a mountainous zone with elevations ranging from 280 m to 3540 m above the mean sea level, with slopes over 40 • (Figure 1).To fulfill the proposed objective, the following stages are proposed: the first stage is the application of the CST, LR, and CVA change detection methods to two kinds of input data: PC and NDVI, to obtain a single difference image for each data type between the dates included in the study.The second stage is the analysis of the change images using a fully automatic detection of changes, carried out using the statistical parameters of change images and by applying the secant technique previously mentioned, to obtain a binary change/no-change map.In this stage, masks of slopes, clouds, and isolated pixels were used to mask the resulting binary maps and thereby obtain the final cartographic inventory of landslides.The final stage is assessing the obtained results and validating them with the collected ground-truth information through confusion matrices.

Study Area
The study area was located in the central part of the Guerrero State in México, and covered an area of about 3300 km 2 in a mountainous zone with elevations ranging from 280 m to 3540 m above the mean sea level, with slopes over 40° (Figure 1).The temperature ranged between 14.3 °C and 28.3 °C.The precipitation ranged from 2100 mm as the maximum average and 800 mm as the minimum, which was recorded from June to September.The area was 74.8% covered by forest (coniferous, mesophilic and mixed), 14.1% by deciduous forest, 7.8% by agricultural use, 3.2% by induced vegetation, and 0.1% by human settlements and urban areas [50].According to the 2010 Population and Housing Census, there are 187 cities and towns, in which there are 15,230 homes inhabited by 59,098 people [51].
Geologically, the area is located physiographically in the Sierra Madre del Sur [52] and has a variety of metamorphic rock composition consisting of schists and gneisses of biotite and quartzite, outcrops of deposited limestone, metavolcanic rocks with sedimentary influence, siltstones, sandstones, conglomerates, and carbonate rocks.Rhyolitic rocks are also found as a result of Oligocene-Miocene volcanism.The youngest rocks correspond to alluvial deposits present in the margins and river beds [53,54] (Figure 2).The temperature ranged between 14.3 • C and 28.3 • C. The precipitation ranged from 2100 mm as the maximum average and 800 mm as the minimum, which was recorded from June to September.The area was 74.8% covered by forest (coniferous, mesophilic and mixed), 14.1% by deciduous forest, 7.8% by agricultural use, 3.2% by induced vegetation, and 0.1% by human settlements and urban areas [50].According to the 2010 Population and Housing Census, there are 187 cities and towns, in which there are 15,230 homes inhabited by 59,098 people [51].
Geologically, the area is located physiographically in the Sierra Madre del Sur [52] and has a variety of metamorphic rock composition consisting of schists and gneisses of biotite and quartzite, outcrops of deposited limestone, metavolcanic rocks with sedimentary influence, siltstones, sandstones, conglomerates, and carbonate rocks.Rhyolitic rocks are also found as a result of Oligocene-Miocene volcanism.The youngest rocks correspond to alluvial deposits present in the margins and river beds [53,54] (Figure 2).
The area is of interest due to the presence of extraordinary hydrometeorological phenomena in recent years, which have triggered massive landslides that have severely affected the population and infrastructures.The events of September 2013 were particularly notable, when tropical depression No. 13 in the Pacific ocean, and subsequent simultaneous hurricanes: Manuel in the Pacific and Ingrid in the Gulf of México caused significant floods and landslides on the coast of the Guererro state [55].
Particularly noteworthy is the landslide that occurred in the community La Pintada in the municipality of Atoyac, where there were 70 deaths, 379 victims, and 20 damaged buildings [4] (Figure 3).The area is of interest due to the presence of extraordinary hydrometeorological phenomena in recent years, which have triggered massive landslides that have severely affected the population and infrastructures.The events of September 2013 were particularly notable, when tropical depression No. 13 in the Pacific ocean, and subsequent simultaneous hurricanes: Manuel in the Pacific and Ingrid in the Gulf of México caused significant floods and landslides on the coast of the Guererro state [55].Particularly noteworthy is the landslide that occurred in the community La Pintada in the municipality of Atoyac, where there were 70 deaths, 379 victims, and 20 damaged buildings [4] (Figure 3).

Dataset
Two Aster images level AST_L1T from Path 26, Row 48, of the World Reference System (WRS-2), courtesy of the U.S. Geological Survey, were used, corresponding to 10 December 2012 and 13 December 2013, respectively.The dates are before and after and extraordinary hydrometeorological event registered in September 2013 in the study area.The bands used were green, red and near infrared and whose spatial resolution was 15 m.
The image available on Google Earth dated 12 August 2014 was also used, which is the closest after the described meteorological events of September 2013.Also, a 15-m grid cell size DEM, based on a topographic map scale of 1:50,000, courtesy of the National Institute of Statistics and Geography of México (INEGI), was used.
For this work, Dinamica EGO and ArcGIS software were used to develop the models and processes described.The area is of interest due to the presence of extraordinary hydrometeorological phenomena in recent years, which have triggered massive landslides that have severely affected the population and infrastructures.The events of September 2013 were particularly notable, when tropical depression No. 13 in the Pacific ocean, and subsequent simultaneous hurricanes: Manuel in the Pacific and Ingrid in the Gulf of México caused significant floods and landslides on the coast of the Guererro state [55].Particularly noteworthy is the landslide that occurred in the community La Pintada in the municipality of Atoyac, where there were 70 deaths, 379 victims, and 20 damaged buildings [4] (Figure 3).

Dataset
Two Aster images level AST_L1T from Path 26, Row 48, of the World Reference System (WRS-2), courtesy of the U.S. Geological Survey, were used, corresponding to 10 December 2012 and 13 December 2013, respectively.The dates are before and after and extraordinary hydrometeorological event registered in September 2013 in the study area.The bands used were green, red and near infrared and whose spatial resolution was 15 m.
The image available on Google Earth dated 12 August 2014 was also used, which is the closest after the described meteorological events of September 2013.Also, a 15-m grid cell size DEM, based on a topographic map scale of 1:50,000, courtesy of the National Institute of Statistics and Geography of México (INEGI), was used.
For this work, Dinamica EGO and ArcGIS software were used to develop the models and processes described.

Dataset
Two Aster images level AST_L1T from Path 26, Row 48, of the World Reference System (WRS-2), courtesy of the U.S. Geological Survey, were used, corresponding to 10 December 2012 and 13 December 2013, respectively.The dates are before and after and extraordinary hydrometeorological event registered in September 2013 in the study area.The bands used were green, red and near infrared and whose spatial resolution was 15 m.
The image available on Google Earth dated 12 August 2014 was also used, which is the closest after the described meteorological events of September 2013.Also, a 15-m grid cell size DEM, based on a topographic map scale of 1:50,000, courtesy of the National Institute of Statistics and Geography of México (INEGI), was used.
For this work, Dinamica EGO and ArcGIS software were used to develop the models and processes described.

Topographic Normalization
An overview of the process for evaluating the unsupervised detection of landslides proposed to build inventory maps is shown in Figure 4.

Topographic Normalization
An overview of the process for evaluating the unsupervised detection of landslides proposed to build inventory maps is shown in Figure 4. To prevent false results in the change detection analysis, it is recommended to perform a precise co-registration, radiometric calibration, and atmospheric and topographic corrections between the multitemporal images [56,57].Thus, to reduce the effect caused in the reflectance values by the slope changes, terrain orientation and solar geometry at the time of image acquisition, the topographic correction was performed using the Canopy Sensor method + Correction (SCS + C) [58], which is recommended for forested mountain areas over other land-based methods, because it preserves the geotropic nature of the trees (normal to geoid growth) [59].
The C parameter used in the topographic correction to improve the correction by moderating the overcorrection of dimly illuminated pixels [60] was determined by linear regressions between the values of illumination and reflectance, according to a classification of the topographic slopes of the studied area [61].

Pre-Processing
Principal Components (PC) analysis [62] aims to summarize a large group of variables in a new smaller set without significantly losing the original information [43].This increases the possibility of detecting differences in land cover, since the reduction in the dimensionality of data increases the efficiency in the classification process.PC analysis facilitates the interpretation of the variability axes of the image, which allows identifying features that are present in most of the bands and others that are specific only to some of them.PC analysis generates new variables (the components) through a linear combination of the original  bands.Although  PC are ultimately required to reproduce the total variability, most of this variability is contained in a smaller number of  components.As such, by replacing the  bands with  components, the dimensionality is reduced, conserving almost all the information [63].In the present study, six PC were generated for each date study image from the algorithm included in ArcGIS, of which only the first three PC were used in the change detection stage since they concentrated most of the information from the bands (Figure 5).To prevent false results in the change detection analysis, it is recommended to perform a precise co-registration, radiometric calibration, and atmospheric and topographic corrections between the multitemporal images [56,57].Thus, to reduce the effect caused in the reflectance values by the slope changes, terrain orientation and solar geometry at the time of image acquisition, the topographic correction was performed using the Canopy Sensor method + Correction (SCS + C) [58], which is recommended for forested mountain areas over other land-based methods, because it preserves the geotropic nature of the trees (normal to geoid growth) [59].
The C parameter used in the topographic correction to improve the correction by moderating the overcorrection of dimly illuminated pixels [60] was determined by linear regressions between the values of illumination and reflectance, according to a classification of the topographic slopes of the studied area [61].

Pre-Processing
Principal Components (PC) analysis [62] aims to summarize a large group of variables in a new smaller set without significantly losing the original information [43].This increases the possibility of detecting differences in land cover, since the reduction in the dimensionality of data increases the efficiency in the classification process.PC analysis facilitates the interpretation of the variability axes of the image, which allows identifying features that are present in most of the bands and others that are specific only to some of them.PC analysis generates new variables (the components) through a linear combination of the original m bands.Although n PC are ultimately required to reproduce the total variability, most of this variability is contained in a smaller number of m components.As such, by replacing the m bands with n components, the dimensionality is reduced, conserving almost all the information [63].In the present study, six PC were generated for each date study image from the algorithm included in ArcGIS, of which only the first three PC were used in the change detection stage since they concentrated most of the information from the bands (Figure 5).Vegetation shows a sharp absorption peak caused by the photosynthetic pigments in the red wavelengths, in contrast to the strong reflection that they exhibit in the near infrared lengths.In contrast, bare soil is characterized by a smooth monotonic increase as the wavelength increases, whereas water bodies behave inversely, with significant absorption at all wavelengths (except for lengths corresponding to blue) [64].
Thus, it is possible to define an index that estimates green biomass density or chlorophyll density, based on radiometric data, as does the Normalized Difference Vegetation Index (NDVI), which is expressed as the difference between near infrared (NIR) and red channels divided by their sum of them:    Vegetation shows a sharp absorption peak caused by the photosynthetic pigments in the red wavelengths, in contrast to the strong reflection that they exhibit in the near infrared lengths.In contrast, bare soil is characterized by a smooth monotonic increase as the wavelength increases, whereas water bodies behave inversely, with significant absorption at all wavelengths (except for lengths corresponding to blue) [64].
Thus, it is possible to define an index that estimates green biomass density or chlorophyll density, based on radiometric data, as does the Normalized Difference Vegetation Index (NDVI), which is expressed as the difference between near infrared (NIR) and red channels divided by their sum of them: where NIR and Red are the normalized reflectance values of NIR and red bands, respectively.
The NDVI resulting values are between −1 and +1 in direct relation to the actual vegetation cover of each pixel image (Figure 6), [65,66].Vegetation shows a sharp absorption peak caused by the photosynthetic pigments in the red wavelengths, in contrast to the strong reflection that they exhibit in the near infrared lengths.In contrast, bare soil is characterized by a smooth monotonic increase as the wavelength increases, whereas water bodies behave inversely, with significant absorption at all wavelengths (except for lengths corresponding to blue) [64].
Thus, it is possible to define an index that estimates green biomass density or chlorophyll density, based on radiometric data, as does the Normalized Difference Vegetation Index (NDVI), which is expressed as the difference between near infrared (NIR) and red channels divided by their sum of them:

Chi-Square Transformation
The Chi-Square Transformation (CST) is a statistical technique applied to obtain a measure of divergence or distance between groups regarding multiple characteristics.The most-used measure is the Mahalanobis distance (Md), which plays a fundamental role in data analysis with multiple measurements, finding applications in statistical patterns recognition in areas such as medical diagnosis, archeology, or remote sensing [67,68].
The CST is applied in remote sensing to recognize patterns in a single global image, starting from the multivariate information stored in numerous multispectral bands [40].Following Ridd and Liu [69], the CST was applied to each group of PC images, in order to create group images for each date, a single global image that integrates the change represented by Md squared [70,71], determined by the equation: where Y is the square root of Md for each pixel in the change image, X is the difference vector of the n values between the two dates, M is the vector of the mean residuals of each band, and T is the transpose of the matrix (X − M), and ∑ −1 is the inverse covariance matrix of the PC images between the two dates.
Md considers the variance in each variable and the covariance between variables.Geometrically, this is performed by transforming the data into standardized uncorrelated data and computing the ordinary Euclidean distance for the transformed data, thus providing a method to measure distances that considers the scale of the data (standard deviation) [70].
Md is Chi-square distributed with degrees of freedom equal to the number of variables.In this case, this was the number of PC images used in the transformation for each date.Thus, it was possible to compare the distances just like comparing data with the normal distribution when the number of bands has a joint multivariate normal distribution with the covariance matrix ∑ and mean residual vector M [69,71].
The square root of Md (Y) can be displayed as a single image and the pixels represent the value of change through the dates, so theoretically a zero-value pixel means an absolute no-change.

Linear Regression
The linear regression (LR) method for change detection, assumes that the pixel values (Y) of a final date image f 2 results from a linear function of the pixel values (X) from the initial date image f 1 .Thus, it was possible to perform a regression from Y K I,J ( f 2 ) to X K I,J ( f 1 ) by least squares [38,69,72]  to obtain the slope m and ordinate b regression line parameters to find an equation in the form Y = mX + b to model it.
In this case, we applied the LR to the NDVI values and to the first PC (PC1) of the f 1 image.Thus, we obtained a new image Y K I,J , which corresponds to the expected values generated by the prediction model.With the expected values and the actual NDVI and PC1 values of f 2 , it was possible to obtain an image of the residuals of NDVI and PC1 values calculated by the equation: where R k i,j is the residual pixel value of i line and j column for the k band.The values obtained by the prediction model in the expected image will be the same as the actual values of the pixels in the subsequent date, provided no changes have been registered during the period analyzed (Residual = 0).On the other hand, if there are changes, they will be recorded in the corresponding residual value, whose magnitude will indicate the intensity of the change.

Change Vector Analysis
The Change Vector Analysis (CVA) model considers that it is possible to define a change vector for each co-registered pixel between the images of two dates, by using the magnitude and the direction of change [73].If the values of the pixels in two images corresponding to two different dates f 1 and f 2 , are given by G = (g 1 , g 2,..., g n ) F and H = (h 1 , h 2,..., h n ) F , respectively, where n is the number of bands, then the magnitude of the change ∆G can be obtained from the equation: where ∆G represents the absolute magnitude of the total difference between the characteristics of the same pixel between the two dates [74].Therefore, a pixel with a value of ∆G = 0 represents an absolute non-change, whereas high ∆G values represent a high intensity of the registered change.
The geometric concept of the CVA change detection method was applied to the PC images of both study dates.According to Equation ( 4), H represents the three PC of f 1 and G represents the three PC of f 2 .Finally, ∆G represents the changes map by CVA.

Threshold Definition
Any of the explained methods produces as a result, a difference image, where the pixel values correspond to the intensity of the change estimated according to the method applied.However, to complete the process, it was necessary to generate a categorical binary map that isolates the pixels into classes: change and no-change.Thus, we segmented the resulting difference images according to a value assigned as a threshold, thereby the dynamic areas were distinguished from those that remained stable during the period analyzed.
For the study two thresholding methods were tested, by statistical parameters and by the secant method.

Statistical Thresholding Method
Considering that in general there are few pixels in which relative changes occur in a period, the classical method to define the threshold value necessary to differentiate change from no-change is based on the fact that the density function of the change or difference image can be considered almost equal to the density function of unmodified pixels [75].As such, the unmodified pixels are distributed around the mean (µ), whereas the changes are distributed in the tails of the frequency distribution [38,47], separated from the value of the mean by a certain number of times of the value of the standard deviation (σ) of the distribution.
The thresholding by statistical parameters is defined by an analysis of the frequency histogram, which involves the calculation of the mean and standard deviation of the change image, and is represented by T = µ ± nσ, where T is the threshold value, µ is the average of the values of changes, σ represents the standard deviation and n is a parameter empirically set by the user, which can be adjusted with values ranging from 0.1 to 2. For this study, considering that this value would determine a radius of change of ≈4.6%, a value of n = 2 was chosen, which is consistent with previous works [43,49] (Figure 7a).where ‖∆‖ represents the absolute magnitude of the total difference between the characteristics of the same pixel between the two dates [74].Therefore, a pixel with a value of ‖∆‖ = 0 represents an absolute non-change, whereas high ‖∆‖ values represent a high intensity of the registered change.
The geometric concept of the CVA change detection method was applied to the PC images of both study dates.According to Equation (4),  represents the three PC of  and  represents the three PC of  .Finally, ‖∆‖ represents the changes map by CVA.

Threshold Definition
Any of the explained methods produces as a result, a difference image, where the pixel values correspond to the intensity of the change estimated according to the method applied.However, to complete the process, it was necessary to generate a categorical binary map that isolates the pixels into classes: change and no-change.Thus, we segmented the resulting difference images according to a value assigned as a threshold, thereby the dynamic areas were distinguished from those that remained stable during the period analyzed.
For the study two thresholding methods were tested, by statistical parameters and by the secant method.

Statistical Thresholding Method
Considering that in general there are few pixels in which relative changes occur in a period, the classical method to define the threshold value necessary to differentiate change from no-change is based on the fact that the density function of the change or difference image can be considered almost equal to the density function of unmodified pixels [75].As such, the unmodified pixels are distributed around the mean (), whereas the changes are distributed in the tails of the frequency distribution [38,47], separated from the value of the mean by a certain number of times of the value of the standard deviation () of the distribution.
The thresholding by statistical parameters is defined by an analysis of the frequency histogram, which involves the calculation of the mean and standard deviation of the change image, and is represented by  =  ±  , where  is the threshold value,  is the average of the values of changes,  represents the standard deviation and  is a parameter empirically set by the user, which can be adjusted with values ranging from 0.1 to 2. For this study, considering that this value would determine a radius of change of ≈4.6%, a value of  = 2 was chosen, which is consistent with previous works [43,49] (Figure 7a).Due to the nature of the equations models, for the CST and CVA methods the resulting change image started from 0 (pixels with no-change), and the rest of the pixels records absolute values depended on the change magnitude.Thus, the resulting changes from the CST and CVA methods were concentrated in the right tail of the histograms (Figure 4b).The LR method allowed the determination of the magnitude of the change in the values of the residual image and the direction Due to the nature of the equations models, for the CST and CVA methods the resulting change image started from 0 (pixels with no-change), and the rest of the pixels records absolute values depended on the change magnitude.Thus, the resulting changes from the CST and CVA methods were concentrated in the right tail of the histograms (Figure 4b).The LR method allowed the determination of the magnitude of the change in the values of the residual image and the direction of the change through using pixel values (positive/negative) in the resulting normal distribution histogram.Thus, it was possible to have pixels with changes in both tails of the histogram; however, according to the type of information (NDV or PC) analyzed with the LR method, some pixel changes were found in the right tail and others in the left tail.(Figure 7a).

Secant Thresholding Method
In the secant method the threshold is automatically defined by selecting the value of the pixel that corresponds to the point in the histogram distribution where the maximum perpendicular line intersects the secant line between the maximum and minimum points of the histogram [49].According to the distribution of the change image of the implemented methods, this procedure was performed only in the right tail for the CST and CVA detection methods (Figure 7b), and in both tails for the LR detection method (Figure 7a).
Any of the described thresholding methods automatically produces a binary (change/no-change) map that contain all the changes detected by the combination of the applied detection method (CST, RL, or CVA) and the thresholding method (statistical or secant).
That is, the resulting map records all changes detected by thresholding, including those identified changes caused by landslide events.To improve the landslide detection, the following additional criteria were implemented: for the LR detection method, only the tail of the histogram was considered first, where the changes caused by the landslides were identified (right for PC and left for NDVI).Then, following Alcántara-Ayala [76], the slope mask only considered the changed pixels with a slope over 5 • as a landslide, as were the changed pixels located outside the mask of clouds and cloud-shadow on any of the study images.Finally, the single isolated pixels detected as changed were eliminated in order to only consider landslides equal to or greater than 450 m 2 in the final maps, equivalent to two Aster joined pixels.

Accuracy Assessment
To assess the accuracy of the unsupervised change detection process, 617 polygons were sampled as ground-truth through the interpretation of color in aerial photo slides near to the Aster imagery dates.A total of 602 polygons (17,385 pixels) corresponded to verified landslides recorded in 2013 due to hydrometeorological events, and 15 polygons (17,632 pixels) corresponded to zones identified as verified no-landslide (Figure 8).These ground-truth samples were compared with each final thematic map produced by the unsupervised change detection method through confusion matrices to obtain the omission-commission errors for each case.
For each final thematic map, the Kappa concordance coefficient of agreement [77] was obtained to quantify the difference between the observed map-reality agreement and the map that would be randomly expected.The Kappa index attempts to define the degree of adjustment only due to the accuracy of the categorization, regardless of random causes [45,78].The Kappa coefficient was calculated by: where k is the Kappa coefficient of agreement, n is the sample size, X ii is the observed agreement, and X i+ X +i is the expected agreement in each category i.
The Kappa coefficient shows if the marked degree of agreement draws away or is not significantly different from the expected random agreement.The observed agreement highlights the diagonal of the confusion matrix, and the expected agreement was used to calculate the fit between the map and the reality due to randomness [45].

Accuracy Assessment
To assess the accuracy of the unsupervised change detection process, 617 polygons were sampled as ground-truth through the interpretation of color in aerial photo slides near to the Aster imagery dates.A total of 602 polygons (17,385 pixels) corresponded to verified landslides recorded in 2013 due to hydrometeorological events, and 15 polygons (17,632 pixels) corresponded to zones identified as verified no-landslide (Figure 8).These ground-truth samples were compared with each final thematic map produced by the unsupervised change detection method through confusion matrices to obtain the omission-commission errors for each case.

Integration of the Final Landslide Map Inventory
According to the assessment of the individual maps resulting from the applied change detection methods, different kinds of data were contained in the input images and different thresholding methods.To further confirm the spatial distribution of the landslides detected in the individual maps, we selected the three maps that reported the best results according to the information of the ground-truth to produce a single final map of landslides.For this task, the applied criterion was that if a pixel was considered as a landslide for at least two of the three best-valued maps, then it was also considered as a landslide in the final map.However, if only one of the three best-valued maps considered the same pixel as a landslide, then the pixel was considered as non-landslide.
We think that, by using this method, a pixel considered as a landslide in the final map will have a higher weight and, consequently, a higher confidence level than any of the pixels detected as a landslide in the resulting individual maps.

Change Detection Maps
Figure 9 shows the change images obtained from the application of CST, LR, and CVA change detection methods.These maps display histograms in which the non-change values are those near to zero.At a glance, the pixels in which a change occurred between the two study dates are those that were the highest (in dark blue) and the lowest (in red).Independent of the applied method, the maps show spatial coherency showing changes in the same areas, so the reliability of the maps appear to be high.

Change Detection Maps
Figure 9 shows the change images obtained from the application of CST, LR, and CVA change detection methods.These maps display histograms in which the non-change values are those near to zero.At a glance, the pixels in which a change occurred between the two study dates are those that were the highest (in dark blue) and the lowest (in red).Independent of the applied method, the maps show spatial coherency showing changes in the same areas, so the reliability of the maps appear to be high.

Threshold Definition
Figure 10 shows the frequency histograms corresponding to the change maps from the LR to PC1 (Figure 10a) and LR to NDVI (Figure 10b) detection method, between the images of dates 10 December 2012 and 13 December 2013.

Threshold Definition
Figure 10 shows the frequency histograms corresponding to the change maps from the LR to PC1 (Figure 10a) and LR to NDVI (Figure 10b) detection method, between the images of dates 10 December 2012 and 13 December 2013.The difference in the calculation of the threshold values calculated by both methods were observed, whereas considering statistical parameters for the LR-PC1 image, the resulting thresholds were -30.755 and +30.755; the secant method resulted in -16,058 and +18,615.As such, in the case of the LR-NDVI image, threshold values obtained using statistical methods resulted in -0.176 and +0.176, whereas the secant method resulted in -0.146 and +0.148.These values indicate that the range of pixels considered as possible landslides in the statistical threshold method was symmetric, i.e., same threshold values with opposite signs; whereas the secant threshold method determines the threshold independently for each tail of the histogram.
The model categorizes all the pixels lower and higher than the threshold values as changes (Figure 11).
Table 1 shows the threshold values and pixel change ratio obtained for each implemented method.
As mentioned earlier, any of the obtained binary (change/no-change) maps contain all the changes detected (including those identified changes caused by landslide events) by combining The difference in the calculation of the threshold values calculated by both methods were observed, whereas considering statistical parameters for the LR-PC1 image, the resulting thresholds were -30.755 and +30.755; the secant method resulted in -16,058 and +18,615.As such, in the case of the LR-NDVI image, threshold values obtained using statistical methods resulted in -0.176 and +0.176, whereas the secant method resulted in -0.146 and +0.148.These values indicate that the range of pixels considered as possible landslides in the statistical threshold method was symmetric, i.e., same threshold values with opposite signs; whereas the secant threshold method determines the threshold independently for each tail of the histogram.
The model categorizes all the pixels lower and higher than the threshold values as changes (Figure 11).Table 1 shows the threshold values and pixel change ratio obtained for each implemented method.As mentioned earlier, any of the obtained binary (change/no-change) maps contain all the changes detected (including those identified changes caused by landslide events) by combining applied detection and thresholding methods.
Figure 12 shows the final individual landslide maps produced for each detection method, type of analyzed data and threshold method, after the application of the additional criteria described in the methodology to improve only the detection of landslides.
Table 2 shows the pixel change ratio after the implementation of the additional criteria to detect only landslides.

Accuracy Assessment
Final maps with more or less pixels identified as landslides do not provide information beyond the agreement in the spatial distribution of the detected landslides in comparison with the spatial distribution of landslides in our ground-truth.Thus, to assess the accuracy achieved in the landslide detection, a comparison of each resulting landslide map with the ground-truth and confusion matrices was performed to obtain the omission and commission errors and the Kappa coefficient of agreement (Table 3).

Final Landslide Inventory Map
To produce a final map integrated from the generated information, the individual landslide maps produced by the LR-PC1-Secant, LR-NDVI-Statistic, and LR-NDVI-Secant methods were selected due to their reported the best results of Kappa concordance index at 84.6%, 81.89%, and 84.81% respectively.
The three selected methods produced acceptable results, with Kappa indexes over 80%.We think it would be possible to use some of them as a final map for subsequent processes.However, a debugging process was applied to generate a unique landslide map, in which, according to the methodology, a pixel is considered to be landslide if at least two of the three best individual maps analyzed, considered it as a landslide.The final landslide inventory map is shown in Figure 13.

Accuracy Assessment
Final maps with more or less pixels identified as landslides do not provide information beyond the agreement in the spatial distribution of the detected landslides in comparison with the spatial distribution of landslides in our ground-truth.Thus, to assess the accuracy achieved in the landslide detection, a comparison of each resulting landslide map with the ground-truth and confusion matrices was performed to obtain the omission and commission errors and the Kappa coefficient of agreement (Table 3).

Final Landslide Inventory Map
To produce a final map integrated from the generated information, the individual landslide maps produced by the LR-PC1-Secant, LR-NDVI-Statistic, and LR-NDVI-Secant methods were selected due to their reported the best results of Kappa concordance index at 84.6%, 81.89%, and 84.81% respectively.
The three selected methods produced acceptable results, with Kappa indexes over 80%.We think it would be possible to use some of them as a final map for subsequent processes.However, a debugging process was applied to generate a unique landslide map, in which, according to the methodology, a pixel is considered to be landslide if at least two of the three best individual maps analyzed, considered it as a landslide.The final landslide inventory map is shown in Figure 13.

Change Detection Maps
The information used as input for the LR detection method and their range of values, i.e., from -1 to +1 for the NDVI data and from -370.4 to +239.1 for the PC1 data, produced the extreme residual values set in the opposite histogram tails from the resulting change map (Figure 9a,b).In the CST and CVA detection methods, the three PC images were used as input.Due to the calculation method of the model, the changes will always be found in absolute values and the right histogram tail (Figure 9c,d).
The same can be seen in a zoom in to a particular detailed zone in the change images (Figure 14), where we can see the aspect of changing zones caused by clouds (Figure 14a

Change Detection Maps
The information used as input for the LR detection method and their range of values, i.e., from -1 to +1 for the NDVI data and from -370.4 to +239.1 for the PC1 data, produced the extreme residual values set in the opposite histogram tails from the resulting change map (Figure 9a,b).In the CST and CVA detection methods, the three PC images were used as input.Due to the calculation method of the model, the changes will always be found in absolute values and the right histogram tail (Figure 9c,d).
The same can be seen in a zoom in to a particular detailed zone in the change images (Figure 14), where we can see the aspect of changing zones caused by clouds (Figure 14a

Threshold Definition
The difference in the calculation of the threshold values calculated by both methods can be observed in Figure 10.Considering the statistical parameters for the LR-PC1 image, the resulting thresholds were -30.755 and +30.755; the secant method resulted in -16.058 and +18.615.As such, in the case of the LR-NDVI image, threshold values by statistical methods resulted in -0.176 and +0.176 whereas the secant method resulted in -0.146 and +0.148.These values indicate that the range in pixels considered as possible landslides in the statistical threshold method was symmetric, i.e., same threshold values with opposite sign; whereas the secant threshold method determines the threshold independently for each tail of the histogram.
The categorical maps of changes resulting from each detection method in combination with each thresholding method, according to the determined threshold values, include more or less pixels as changed, resulting in more or less visually intense maps (Figure 12).
Considering this argument, the categorical map of changes resulting from the analysis of the histogram of the residuals NDVI, detected a greater radius of changes, (including landslides) because the discretization of changes included a more significant portion of the tails of the frequency distribution.
According to the type of data analyzed by the detection and thresholding methods as the input image, the thresholds values varied (Table 1).Consequently, the radius of changed pixels also varied.Thus, in general, the secant thresholding method produces thematic maps with a higher number of pixels identified as changed than the maps produced using statistical-parameters thresholding.The maps obtained by the LR detection method applied to PC1 and NDVI images and thresholded by the secant method (Figure 8b,d), due to the higher number of identified changed pixels (6.220% and 6.257%) were more visually intense than the maps produced by applying the same method to the same images, but thresholded by statistical parameters (Figure 8a,c) with 2.395% and 4.239% as the change ratio values respectively.
The highest change ratio was reported in the map produced by the CVA detection method applied to PC images and thresholded by secants (7.279%).The lowest change ratio corresponds to the map generated by the same method but thresholded by statistical parameters (1.823%).
The consistency in the resulting maps can be observed in Figure 12.All maps identified landslides in the same zones of the study area, regardless of the detection method implemented, the type of information analyzed in images or the applied thresholding method.However, this combination of method-data-threshold, enabled identifying some differences in the resulting thematic maps.Through visual exploration, a larger or smaller number of pixels can be observed as categorized as landslides.Thus, we produced more or less visually intense maps.
Different land covers react differently to energy [79,80] according to the kind of data analyzed in the input images in the landslide detection process.Thus, as seen in Figure 12, the use of PC data as input images for the CVA detection method and thresholded by secants (Figure 12h) produced a map with a higher number of pixels identified as landslides, followed by the map produced by the LR detection method using NDVI data thresholding by secants (Figure 12d).The use of PC data in input images using the CVA detection model and thresholded by statistical parameters (Figure 12g), produces the map that detects fewer pixels as landslides and the least visually intense map.
As shown in Table 2, the change ratio significantly decreased when applying the additional criteria considered for the detection of only landslides, resulting in a landslide map produced by the CVA-PC method and thresholded by secants that maintains the highest number of pixels detected as landslides.The map produced by the same method but thresholded by statistical parameters detected fewer pixels as landslides.

Accuracy Assessment
According to the assessment of the resulting maps shown in Table 3, except for the CVA-PC detection method thresholded by statistical parameters, the other analyzed detection methods had errors under 20% with respect to to the ground-truth validation.
It can be argued that the values determined for the thresholds are related to the type of data contained in the analyzed images during the detection process.For this reason, the maps are not entirely comparable independently; however, they provide an idea about the similarity in the results obtained.
As mentioned above, the input images used in the application of the landslide detection process reacted differently, according to the information contained and the combination detection-thresholding methods.Thus, some processes were more or less sensitive in detecting landslides.
The higher mean omission/commission errors (20.54% and 14.41%, respectively) correspond to the landslide map resulting from the CVA-PC method, thresholded by statistical parameters, which was perceived in the Kappa concordance index of 59.09%, corresponding to the lowest value of the analyzed processes.
The lower mean omission/commission errors (7.64% and 6.56%, respectively) correspond to the LR-NDVI detection method, applying secants to define the threshold value, which was favorably noted in the Kappa concordance index of 84.81%.The next most accurate method was the LR-PC1 detection method, also applying the threshold secant method, with 7.75% mean omission error, 6.62% mean commission error, and a Kappa concordance index of 84.60%.Although the difference between these two Kappa indices is meager (0.21%), this is the same detection method; however, the type of data contained in the images analyzed was different (NDVI vs. PC1).
According to the results, the use of NDVI data by the LR detection method and thresholding by statistic and secant methods produced Kappa indexes with values higher than 80%, being two of the three best results obtained.Also, the implementation of the same detection method, but using PC1 as the input data and thresholding by the secant method, produced a Kappa index higher than 80% (Table 3).
For results obtained by the CST detection method using PC data as input images were very similar, producing Kappa indexes of 63.32% for thresholding by statistic parameters and 65.74% for thresholding by the secant method (Table 3).
Finally, the CVA detection method using PC data as input images and thresholding by statistical parameters produced the lowest Kappa index (59.09%) of all the analyzed detection processes.However, using this same method but applying the secant method for thresholding produced a Kappa index of 77.14%.

Final Landslide Inventory Map
As can be seen in the final map (Figure 13), the landslides were concentrated mainly in the north-east and central zone of the study area, and some others occurred in smaller quantities in the south-west zone.The landslide zones occurred in September 2013 (indicated in blue).Over these areas, we identified and digitized landslide polygons on aerial photos during the integration of the ground-truth used to assess the accuracy of maps.The presence of landslides in these areas were confirmed in the individual maps produced by the best detection methods and also in the final inventory map.
As for the individual maps resulting from the detection methods, the final landslide inventory map was also accuracy-evaluated with ground-truth data.The results show that the final map detected a landslide ratio of 1.48%.The mean omission error was 7.97%, and 6.79% was recorded as the mean commission error.The Kappa concordance index was 84.15%.
Notably, this Kappa index is less than the index of the two best individual maps analyzed: LR-PC1 and LR-NDVI thresholding by secants, which reached Kappa indexes of 84.60% and 84.81%, respectively (Table 3).However, this represents a minimal variation (0.45% and 0.66%, respectively) and places the final landslide inventory map in an acceptable rank according to the Kappa index, which was over 80%.Moreover, the final landslide inventory map is actually the final integration resulting from a combination of the different kinds of data analyzed and two automatic threshold definition methods applied.
Figure 15 shows zoomed in areas on the Aster images before (10 December 2012) and after (13 December 2013) the extraordinary hydrometeorological events that occurred in September 2013 in the study area.The absence of landslides can be seen before September 2013 (Figure 15a,d); this is the same area indicated in blue color in Figure 13.In this same area, the post-landslide image (Figure 15b,e) shows the scars caused by a significant number of landslides.In a visual exploration, the detected landslides through the method proposed, included in the final inventory map in yellow (Figure 15c,f), matches the scars caused by landslides.
This situation was also observed in the area indicated in green (Figure 13), where no identification or digitization of landslide polygons was carried out in the ground-truth integration.However, the resulting map from the detection identified pixels as landslides in this area, which was confirmed as scars caused by real landslides on the post-Aster image (Figure 15h).This same analysis was carried out in the south-west zone of the study area.In the Aster preimage (Figure 15j) the absence of landslides can be observed.Figure 15j shows La Pintada (red arrow).As previously described, this town was affected by a massive landslide in September 2013 (Figure 3).In the Aster post-image (Figure 15k), the scars caused by landslides can be seen, and particularly the massive landslide that affected La Pintada.The figure also confirms that this and other landslides were detected by the proposed method and are represented in the final landslide inventory map (Figure 15l).The absence of landslides can be seen before September 2013 (Figure 15a,d); this is the same area indicated in blue color in Figure 13.In this same area, the post-landslide image (Figure 15b,e) shows the scars caused by a significant number of landslides.In a visual exploration, the detected landslides through the method proposed, included in the final inventory map in yellow (Figure 15c,f), matches the scars caused by landslides.
This situation was also observed in the area indicated in green (Figure 13), where no identification or digitization of landslide polygons was carried out in the ground-truth integration.However, the resulting map from the detection identified pixels as landslides in this area, which was confirmed as scars caused by real landslides on the post-Aster image (Figure 15h).This same analysis was carried out in the south-west zone of the study area.In the Aster pre-image (Figure 15j) the absence of landslides can be observed.Figure 15j shows La Pintada (red arrow).As previously described, this town was affected by a massive landslide in September 2013 (Figure 3).In the Aster post-image (Figure 15k), the scars caused by landslides can be seen, and particularly the massive landslide that affected La Pintada.The figure also confirms that this and other landslides were detected by the proposed method and are represented in the final landslide inventory map (Figure 15l).
Other previous works use digital elevation models, semi-automatic analysis of satellite images (panchromatic, multispectral and radar) and field mapping tools [23].One of the most widely used methods involves calculating the differences between pre-and post-landslide images.The methods developed in this work are more complex than image differences; however, in concordance with previous works, the results of the modified zones can be directly observed [81].

Conclusions
In the creation of unsupervised change detection methods to produce automatic landslide cartography, the type of data analyzed, the detection methods and the automatic definition of thresholds implemented provide a theoretical and methodological framework suitable for the development of automated detection methods for landslides and automatic generation of map landslides inventory.
The proposed method for the integration of automatic landslide mapping automatically selects the threshold through a process that limits human participation, adding objectivity, and speed.The definition of the threshold is an important step that defines the overall capacity of the method, since an appropriate threshold in combination with additional criteria, maximizes the capacity to discriminate the areas of changes caused by landslide events of stable areas.
The suppression of isolated pixels from the individual landslide maps, which correspond to doubtful landslides, confirmed the partial previous results.We also verified that the landslide inventory map produced from the LR detection method using NDVI data and thresholding by secants was the best method reporting the best Kappa concordance index at 84.81%, in accordance with similar previous studies [49].The debugging process provided certainty to the spatial dimension of the detected landslides.
The main advantages of the proposed method of automatic generation of the landslide inventory map are: (1) the use of no-cost satellite images; (2) the implemented unsupervised methods for change detection do not require large amounts of antecedent information, which considerably reduces the time used in the integration of digitized landslide inventory cartography and the costs by campaigns to gather information in the field; (3) all the individual processes that include the automatic definition of thresholds and the production and debugging of thematic maps can be carried out by automated algorithms, independent of the possible subjectivity of the users integrating a fully unsupervised method to produce landslide inventory maps; and (4) we think it is possible to replicate the processes described in the proposed methodology in any area of the world where landslides have occurred, obtaining results similar to those shown here.

Figure 1 .
Figure 1.Study area: the central part of Guerrero State in México.

Figure 1 .
Figure 1.Study area: the central part of Guerrero State in México.

Figure 2 .
Figure 2. Geological map of study area.

Figure 2 .
Figure 2. Geological map of study area.
NDVI = (NIR − Red)/(NIR + Red) (1) where NIR and Red are the normalized reflectance values of NIR and red bands, respectively.The NDVI resulting values are between −1 and +1 in direct relation to the actual vegetation cover of each pixel image (Figure 6), [65,66].

3. 3 .
Change Detection3.3.1.Chi-Square TransformationThe Chi-Square Transformation (CST) is a statistical technique applied to obtain a measure of divergence or distance between groups regarding multiple characteristics.The most-used measure is the Mahalanobis distance (Md), which plays a fundamental role in data analysis with multiple measurements, finding applications in statistical patterns recognition in areas such as medical diagnosis, archeology, or remote sensing[67,68].
NDVI = (NIR − Red)/(NIR + Red) (1) where NIR and Red are the normalized reflectance values of NIR and red bands, respectively.The NDVI resulting values are between −1 and +1 in direct relation to the actual vegetation cover of each pixel image (Figure 6), [65,66].

3. 3 .
Change Detection3.3.1.Chi-Square TransformationThe Chi-Square Transformation (CST) is a statistical technique applied to obtain a measure of divergence or distance between groups regarding multiple characteristics.The most-used measure is the Mahalanobis distance (Md), which plays a fundamental role in data analysis with multiple measurements, finding applications in statistical patterns recognition in areas such as medical diagnosis, archeology, or remote sensing[67,68].

Figure 7 .
Figure 7. Definition of thresholds values by (a) Statistical parameters, and (b) Secant method.

Figure 7 .
Figure 7. Definition of thresholds values by (a) Statistical parameters, and (b) Secant method.

Figure 8 .
Figure 8. Location and distribution of ground-truth.Figure 8. Location and distribution of ground-truth.

Figure 8 .
Figure 8. Location and distribution of ground-truth.Figure 8. Location and distribution of ground-truth.

Figure 9 .
Figure 9. Change images resulting from the application of the detection models.(a) Linear Regression to NDVI images, (b) Linear Regression to first Principal Components PC1 images, (c) Chi-Square Transform to PC images, and (d) Change Vector Analysis (CVA) to PC images.

24 Figure 9 .
Figure 9. Change images resulting from the application of the detection models.(a) Linear Regression to NDVI images, (b) Linear Regression to first Principal Components PC1 images, (c) Chi-Square Transform to PC images, and (d) Change Vector Analysis (CVA) to PC images.

24 Figure 11 .
Figure 11.Categorical change maps resulting from the (a,c) statistical and (b,d) secant thresholding methods applied to change images (a,b) LR-PC1 and (c,d) LR-NDVI.

Figure 11 .
Figure 11.Categorical change maps resulting from the (a,c) statistical and (b,d) secant thresholding methods applied to change images (a,b) LR-PC1 and (c,d) LR-NDVI.

Figure 13 .
Figure 13.Final landslide inventory map.The blue highlighted areas are where photointerpretation and digitization was completed.The green highlighted are zones where no-photointerpretation or digitization work was carried out.

24 Figure 13 .
Figure 13.Final landslide inventory map.The blue highlighted areas are where photointerpretation and digitization was completed.The green highlighted are zones where no-photointerpretation or digitization work was carried out.
,b) and by landslides (Figure 14c,d) identified in different tails of the histogram.

Figure 14 .
Figure 14.Zoom to change zone caused by clouds in change images (a) LR-PC1 and (b) LR-NDVI; and caused by landslides in (c) CST-PC and (d) CVA-PC images.

Figure 14 .
Figure 14.Zoom to change zone caused by clouds in change images (a) LR-PC1 and (b) LR-NDVI; and caused by landslides in (c) CST-PC and (d) CVA-PC images.

Figure 15 .
Figure 15.Landslide zoomed areas.(a,d,g,j) Aster Images from 10 December 2012 (before the events of September 2013).(b,e,h,k) Aster images from 13 December 2013 (after the events of September 2013).(c,f,i,l) Final landslide inventory map (in yellow color) superimposed on the post-Aster image.

Figure 15 .
Figure 15.Landslide zoomed areas.(a,d,g,j) Aster Images from 10 December 2012 (before the events of September 2013).(b,e,h,k) Aster images from 13 December 2013 (after the events of September 2013).(c,f,i,l) Final landslide inventory map (in yellow color) superimposed on the post-Aster image.

Table 1 .
Threshold values and change ratio resulting from the threshold definition.

Table 2 .
Change ratio resulting from the landslide detection.

Table 1 .
Threshold values and change ratio resulting from the threshold definition.

Table 2 .
Change ratio resulting from the landslide detection.

Table 3 .
Omission/commission error and the kappa coefficient of agreement, obtained for each landslide map.

Table 3 .
Omission/commission error and the kappa coefficient of agreement, obtained for each landslide map.