Landslide mapping from multi-sensor data through improved change detection-based Markov random field

Accurate landslide inventory mapping is essential for quantitative hazard and risk assessment. Although multi- temporal change detection techniques have contributed greatly to landslide inventory preparation, it is still challenging to generate quality change detection images (CDIs) for accurate landslide mapping. The recently proposed change detection-based Markov random field (CDMRF) provides an effective approach for rapid mapping of landslides with minimum user interventions. However, when CDI is generated by change vector analysis (CVA) alone, the CDMRF method may suffer from noise especially when the pre- and post-event remote sensing images are acquired under different atmospheric, illumination, and phenological conditions. This paper improved such CDMRF approach by integrating normalized difference vegetation index (NDVI), principal component analysis (PCA), and independent component analysis (ICA) generated CDIs with MRF for landslide inventory mapping from multi-sensor data. To justify the effectiveness and applicability, the improved methods were applied to map rainfall-, typhoon-, and earthquake-triggered landslides from the pre- and post-event sa- tellite images acquired by very high resolution QuickBird, high resolution FORMOSAT-2, and moderate resolution Sentinel-2. Moreover, they were tested on pre-event Landsat-8 and post-event Sentinel-2 datasets, in- dicating that they are operational for landslide inventory mapping from combined multi-temporal and multisensor data. The results demonstrate that the improved δ NDVI-, PCA-, and ICA-based approaches perform much better than CVA-based CDMRF in terms of completeness, correctness, Kappa coefficient, and F -measures. To the best of our knowledge, it is the first time that NDVI, PCA, and ICA are integrated with MRF for landslide inventory mapping from multi-sensor data. It is anticipated that this research can be a starting point for developing new change detection techniques that can readily generate quality CDI and for applying advanced machine learning algorithms (e.g., deep learning) to automatic detection of natural hazards from multi-sensor time series data.

An accurate landslide inventory that typically documents the date of the event, spatial extent, and morphological features is essential for quantitative hazard and risk assessment, particularly the investigation of landslide characteristics such as magnitude/frequency relationships corresponding to different predisposing and triggering factors (Guzzetti et al., 2005;van Westen et al., 2008;Corominas et al., 2014). Traditional landslide inventory mapping approaches rely on 3D visual interpretation of aerial photos or satellite images in addition to intensive field surveys (Guzzetti et al., 2000Brardinoni et al., 2003;Malamud et al., 2004;Sato et al., 2007;Galli et al., 2008;Saba et al., introducing new remote sensing technologies (Metternicht et al., 2005;Guzzetti et al., 2012;Scaioni et al., 2014).
The previous landslide mapping methods using optical remote sensing images can be classified into two groups: pixel-and objectbased methods Li et al., 2016b). The commonly used pixel-based landslide mapping methods mainly involve change detection technique (Hervás et al., 2003;Park and Chi, 2008;Yang and Chen, 2010;Mondini et al., 2011;Li et al., 2016aLi et al., , 2016b, image correlation technique (Travelletti et al., 2012;Lucieer et al., 2014;Stumpf et al., 2017), and image classification algorithms such as maximum likelihood classification (Nichol and Wong, 2005;Mondini et al., 2017). However, pixel-based landslide mapping methods are comparatively sensitive to noise and may fail when the spectral information is limited. Thus, for further investigation of spatial and contextual features of landslides, object-oriented analysis (OOA) has been developed for landslide inventory mapping (Martha et al., 2010(Martha et al., , 2012(Martha et al., , 2013Lu et al., 2011;Stumpf and Kerle, 2011;Behling et al., 2014Behling et al., , 2016Kurtz et al., 2014;Rau et al., 2014). Although landslide mapping through OOA can achieve better accuracy than pixel-based approach (Keyport et al., 2018), the effectiveness of OOA is highly dependent on the quality of image segmentation. Over-and under-segmentation of landslides often occur especially with complex land cover types such as large urban areas inside images. This problem may be difficult to solve even with multi-scale optimization segmentation, thus to some extent reducing the efficiency and accuracy of landslide mapping (Drǎguţ et al., 2010). Mondini et al. (2011) used normalized difference vegetation index (NDVI)-, principal component analysis (PCA)-, and independent component analysis (ICA)-based change detection techniques to map rainfall-induced landslides caused by the 2009 Messina event in southern Italy. The NDVI approach is based on the assumption that landslides generally result in vegetation disturbances. PCA and ICA are nonparametric feature extraction methods and they transform the original dataset into new feature spaces, in which the desired information can be obtained more easily. Although such methods have been successfully used for landslide inventory mapping, they only took into account the spectral information of landslides, without considering the spatial contextual features between landslides and their surroundings. Li et al. (2016a) developed a semi-automatic change detectionbased Markov random field (CDMRF) method to compile a landslide inventory on western Lantau Island, Hong Kong. Although this method takes into account both the spectral and spatial contextual information of landslides, it is not perfect and suffers from limitations in some cases. For example, the change vector analysis (CVA) used in CDMRF generally requires careful image preprocessing such as relative radiometric correction to preserve the similar spectral and radiometric characteristics between the pre-and post-event images. Also, CVA generated CDIs may also potentially suffer from noise when the land cover types inside images are spectrally heterogeneous. In addition, only red, green, and blue bands were used to generate the CDI for landslide mapping in CDMRF. It is thus interesting to further examine the suitability of other spectral bands (e.g., near infrared) for landslide mapping. Finally, Li et al. (2016a) applied CDMRF only on single sensor data, without attempting to demonstrate its effectiveness in multi-sensor data applications. Thus research to test CDMRF on multi-sensor data and even on different pre-and post-event sensor data, establishing CDMRF as a generic method for landslide mapping, is appealing. This paper further develops previous research conducted by Mondini et al. (2011) andLi et al. (2016a). We improved the CDMRF approach by integrating NDVI, PCA, and ICA generated CDIs with MRF for landslide mapping from remote sensing images. To test the effectiveness and applicability, the improved methods were applied to map landslides triggered by rainfall, a typhoon, and an earthquake, respectively. The approaches were tested on multi-sensor data, including very high resolution QuickBird, high resolution FORMOSAT-2, and moderate resolution Sentinel-2 images, and also on pre-event Landsat-8 and post-event Sentinel-2 data. The results show that the improved CDMRF methods significantly outperform the CVA-based approach proposed in Li et al. (2016a).

Study sites and datasets
This study follows the landslide classification system proposed by Hungr et al. (2014). Thus, "slides", "flows" and "falls", the three typical types of mass movement in the study sites, are classified as landslides.

The rainfall-triggered Messina event
A catastrophic landslide event occurred on the night of 1 October 2009 in Messina due to an intense rainstorm (about 225 mm in 8 h with a peak of 115 mm in 3 h) (Ciampalini et al., 2015). The main types of triggered slope failures were shallow debris slides and flows that formed from landslide depletion zones up to several meters thickness. The debris slides chiefly developed rectilinearly along the slopes while the debris flows mainly initiated at upper slopes and propagated along the topographic plan concavities downhill into the drainage channels, consisting of materials from the weathered eluvial or colluvial layer (Lombardo et al., 2014). The landslides affected several villages in the south of Messina, with 31 people killed and 6 people missing. One hundred and twenty-two people were reported injured and 2019 people

The typhoon-triggered Taiwan event
On 7 August 2009, Typhoon Morakot swept Taiwan with a maximum wind speed of 40 m/s and nearly 3000 mm rainfall in 100 h, resulting in 673 deaths and at least 17,417 new landslides (Tsou et al., 2011;Lin et al., 2011;Mondini et al., 2013). The test site is located in the southern Huaguoshan watershed, an area often affected by rainfallinduced landslides. The main types of the triggered landslides in this area are shallow debris slides and flows. Major geological units in this site are composed of the Changchikeng Formation, the Chaujou Formation, and terrace deposits, with deeply fractured bedrocks caused by active tectonics (Weng et al., 2011;Mondini et al., 2013). The study site was abundantly covered by vegetation such as bamboo and forests that can be readily destroyed by slope failures. Pre-and post-event 8 bit FORMOSAT-2 images acquired on 3 July 2007 and 23 October 2009, respectively were used for landslide mapping. The 8 m FORMOSAT-2 multispectral bands were used in this study.

The earthquake-triggered Jiuzhaigou event
At 21:19:46 on 8 August 2017, a M s 7.0 earthquake with a depth of 20 km occurred in Jiuzhaigou. The epicenter is situated in Zhangzha Town (33.20°N 103.82°E). More than 1700 aftershocks occurred in the subsequent 36 h (Fan et al., 2018). The earthquake caused at least 25 deaths, with more than 70,000 buildings damaged, and triggered a large number of co-seismic landslides. Almost 1883 landslides were identified in the earthquake-affected areas, with the main types of mass movement as rockfalls and debris slides (Fan et al., 2018), which removed vegetation and exposed a large number of slopes. The test site is close to the extension line of the Zhongcha village secondary fault, where a tectonic surface rupture zone trending NW was identified (Yi et al., 2018). The average altitude is about 4000 m above sea level. Pre-and post-event 12 bit Sentinel-2 level-1C images were acquired on 29 July 2017 and 7 September 2017, respectively. The 10 m multispectral bands were utilized in this study. In order to test the effectiveness of CDMRF for different sensor data, the 30 m 12 bit Landsat-8   Lu, et al. Remote Sensing of Environment 231 (2019) 111235 OLI image acquired on 11 July 2017 was used as another pre-event image and the 30 m Sentinel-2 image aggregated from 10 m was used as another post-event image.

Methodology
The improved CDMRF methods for landslide mapping in this study primarily include four steps. The first step is image preprocessing. Then, three change detection techniques are used to generate landslide change detection images (CDIs). Next, a multi-threshold method is applied to CDIs to automatically generate training samples of landslides and non-landslides. Finally, the MRF approach is used to automatically map landslides from the post-event images.

Image preprocessing
For the Messina event, both pre-and post-event QuickBird images were first pan-sharpened to a spatial resolution of 0.6 m using the builtin Gram-Schmidt algorithm in ENVI. Then, ortho-rectification was implemented to the pan-sharpened images by using satellite rational polynomial coefficients (RPC), a 1 m LiDAR-derived digital elevation model (DEM), and two sets of ground control points (GCPs) that included 20 and 15 GCPs on pre-and post-event images, respectively. Next, a dark object subtraction proposed by Mondini et al. (2011) was used to reduce spectral radiance difference between bi-temporal images resulting from different atmospheric and illumination conditions. Such relative radiometric correction method adjusts the brightness of postevent image to pre-event image. Finally, in order to reduce change detection commission errors potentially existing in built-up areas caused by infrastructure construction, all the built-up areas in both bitemporal images were removed using the pre-mapped vectorized topographic data. The preprocessed pre-and post-event QuickBird images are shown in Figs. 1(a) and (b), respectively. The manually digitized reference landslides (red) superimposed on the post-event image are shown in Fig. 1(c). Fig. 1(d) shows the landslide binary reference map.
For the Taiwan event, the pre-and post-event FORMOSAT-2 images are the Level 4 products that had been radiometrically corrected and orthorectified using a DEM, GCPs, and a rigorous model supplied by Chen et al. (2005). The preprocessed pre-and post-event FORMOSAT-2 images are shown in Figs. 2(a) and (b), respectively. The manually digitized reference landslides (red) overlapped on the post-event image are displayed in Fig. 2(c). Fig. 2(d) illustrates the landslide binary reference map.
For the Jiuzhaigou event, the pre-and post-event Sentinel-2 Level-2A Bottom of Atmosphere (BOA) reflectance images were derived from the Level-1C Top of Atmosphere (TOA) reflectance images through Sen2cor, the official atmospheric correction algorithm for Sentinel-2 data (Roy et al., 2017). The preprocessed pre-and post-event 10 m Sentinel-2 images are shown in Figs. 3(a) and (b), respectively. The manually prepared landslide reference (red) is overlaid on the post- In addition, to test the applicability of the improved CDMRF on different pre-and post-event sensor data, the 30 m Landsat-8 L1TP image acquired on 11 July 2017 and atmospherically corrected using LaSRC algorithm was used. The post-event Sentinel-2 BOA image was then aggregated from 10 m to 30 m. Next, registration was completed by selecting two sets of 10 GCPs on pre-event Landsat-8 and post-event Sentinel-2 images, respectively. The preprocessed pre-event Landsat-8 and post-event Sentinel-2 images are shown in Figs. 4(a) and (b), respectively. The landslide reference (red) superimposed on the postevent image is presented in Fig. 4(c). Fig. 4(d) shows the landslide binary reference map.

Change detection images
Generating high quality CDI for multi-temporal change analysis of land surfaces has always been an important research topic (Lu et al., 2004). The change detection method used in Li et al. (2016a) for landslide mapping is based on CVA of the pre-and post-event images: where I t1 and I t2 are pixel values at the pre-event time t 1 and post-event time t 2 , b is the band number of remote sensing images, ρ(I) is the pixel value in CDI. However, such a method potentially suffers from noise existing in CDIs especially when the land cover types inside images are spectrally heterogeneous and pre-and post-event images are acquired under different atmospheric, illumination, and phenological conditions. This study instead attempts to use three other change detection methods (i.e., NDVI, PCA, and ICA) for landslide mapping. The bitemporal NDVI difference (δNDVI) is used to detect potential landslides based on the assumption that landslides generally alter vegetated land surfaces: where ρ NIR and ρ red are near-infrared (NIR) and red bands, respectively. Both PCA and ICA are widely used feature extraction approaches (Abdi and Williams, 2010;Jolliffe, 2011;Oja, 2000a, 2000b), which can reduce data redundancy and enhance discrepancies between different spectral information (Lu et al., 2004) and have great noise separation capabilities (Thomas et al., 2002). In this study, the red and NIR bands in both pre-and post-event images were used as where cov(•) means covariance and R pre , R post , NIR pre , and NIR post are the red and NIR bands in the pre-and post-event images with their means subtracted. The blue and green bands were not used given that they have a high correlation with red band, which may potentially result in data redundancy. Then, 4 eigenvectors and the corresponding 4 eigenvalues of the covariance matrix C were obtained. The eigenvector with the greatest eigenvalue normally represents the principal spectral information related to the main land cover types (e.g., vegetation) in the study areas, whereas the eigenvector with relatively smaller eigenvalue represents the minor spectral information (e.g., landslide). Next, the 4 eigenvectors were ordered by eigenvalues from greatest to lowest and an eigenvector matrix was constructed. Finally, the new images were obtained by multiplying the transposed eigenvector matrix with the transposed original column vectors. In this way, the original multi-temporal spectral bands were transformed by PCA into 4 new images.
ICA separates the multi-temporal spectral bands into 4 statistically independent components. Specifically, the multi-temporal spectral bands were formulated as a linear combination of the 4 independent components (Hyvärinen and Oja, 2000b): where X is the multi-temporal spectral bands mentioned before; A is the mixing matrix; and s is the independent component matrix. In this research, FastICA proposed by Hyvärinen and Oja (2000a) was used to estimate the 4 independent components.
The CDIs generated by CVA, δNDVI, PCA, and ICA from the pre-and post-event QuickBird images for the Messina event are shown in Fig. 5(a) to (d), respectively. Fig. S1 shows the CDIs generated from bitemporal FORMOSAT-2 images for the Taiwan event. The CDIs derived from the pre-and post-event Sentinel-2 images for the Jiuzhaigou event are shown in Fig. S2, while CDIs generated from the pre-event Landsat-8 and the post-event Sentinel-2 images are shown in Fig. S3. Visually, the landslides in CVA generated CDI are not even discernible but they are much brighter than the surroundings in δNDVI, PCA, and ICA generated CDIs. However, it is apparent that there are landslide-like errors in CDIs caused by phenology and illumination differences in preand post-event images. For accurate landslide mapping, a multithreshold method and MRF are further used.

Multi-threshold for training samples generation
In CDIs, brighter pixels mainly indicate larger land cover changes, which can be potentially considered as landslides. The darker pixels mainly represent less changes and thus could be considered as nonlandslides. In order to automatically obtain training samples of landslides and non-landslides from CDIs, as in Li et al. (2016a), a multithreshold method adapted from Chuvieco et al. (2002) where I CDI is the pixel value in CDI, μ CDI and σ CDI are the mean and standard deviation of CDI, respectively. The thresholds T and ΔT (T ∈ R + ) are used to determine the training samples of landslides and non-landslides, respectively. The values of T and ΔT determined by trial and error are summarized in Table 1. Figs. 6 to 9 show the training sample masks generated from CDIs using Eq. (5) for Messina, Taiwan, Jiuzhaigou (single sensor), and  Jiuzhaigou (two sensors), respectively. The red, green, and black pixels represent landslides, non-landslides, and uncertain areas, respectively. By superimposing these training sample masks onto the post-event images, the corresponding spectral information of landslides and nonlandslides in the post-event images was acquired and used as training samples. Next, MRF was applied to classify these uncertain areas in the post-event images into landslides or non-landslides over the study areas.

Markov random field
MRF is an established approach that has been used in computer vision and pattern recognition (Boykov and Kolmogorov, 2004;Rother et al., 2004;Szeliski et al., 2008;Li, 2009). However, its application in geoscience is relatively limited. Recently, CDMRF was used for landslide inventory mapping from bi-temporal aerial photos (Li et al., 2016a), in which landslide mapping was regarded as a pixel-labeling problem, assigning label 1 for landslide and 0 for non-landslide. This pixel-labeling problem can be realized by minimizing the following energy function in the MRF framework: where E r (L) and E b (L) are regional (mainly captures the spectral information) and boundary (mainly captures the spatial contextual information) energy in the MRF framework, respectively, and their weightings are balanced by the non-negative coefficient λ. L = {l 1 , l 2 , … , l i } is the label set in which l i ∈ {0, 1} represents the label of the ith pixel in the image. MRF achieves the globally optimized segmentation by finding the label set L that minimizes the energy function (6): Further details on the definitions of energy terms E r (L) and E b (L) in (6) and the implementation of (7) can be found in (Boykov and Funka-Lea, 2006) and (Boykov and Kolmogorov, 2004). The program in this paper was run under MATLAB R2014a 64 b in Windows 10 OS with a Dell workstation of Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10 GHz, 16 GB RAM. The source code is available upon request.

Results
All the results in Figs. 10 and S4 to S6 demonstrate that the improved methods especially PCA-and ICA-based CDMRF can be operational for mapping landslides triggered by rainfall, typhoon, and earthquake over different study sites from different satellite sensor data.
The landslide mapping results of the 2009 Messina event are shown in Fig. 10. By visual inspection, most landslides are shallow mass movements such as debris slides and flows that were induced by the extreme rainfall event. Such types of mass movements alter land surfaces, thus making it easier to identify them using multi-temporal change detection techniques. Compared with CVA-and δNDVI-based competitors (Figs. 10a to d), PCA-and ICA-based CDMRF can detect landslides more completely and accurately (Figs. 10e to h). Despite the phenological and illumination differences between the pre-and postevent images, PCA-and ICA-based CDMRF can still obtain decent landslide mapping results.
The landslide mapping results for Taiwan, Jiuzhaigou (single sensor), and Jiuzhaigou (two sensors) events are shown in Figs. S4, S5, and S6, respectively. The results in Figs. S4 and S5 further corroborate the applicability and robustness of the improved CDMRF methods to different satellite sensor data (i.e., FORMOSAT-2 and Sentinel-2) over different study sites (i.e., Taiwan and Jiuzhaigou). The results in Fig. S6 show that the improved CDMRF methods can map landslides from combined two different satellite sensor data (i.e., Landsat-8 and Sentinel-2), indicating that they are operational for landslide inventory mapping from multi-temporal and multi-sensor data.
All landslide mapping results were quantitatively compared with the manually digitized reference landslide maps. Four types of indices were used for accuracy assessment, i.e., Completeness = A m /A r , Correctness = A m /A lm , Kappa coefficient, and F-measure, in which A m is the total pixel number of mapped landslides that are matched with the reference map, A r is the total pixel number of the landslide reference map, and A lm is the total pixel number of mapped landslides. F-measure is defined as (Sokolova et al., 2006): where F 0.5 (β = 0.5), F 1 (β = 1), and F 2 (β = 2) were derived to represent lower, equal, and greater weight of correctness, respectively. The quantitative evaluation results for CVA-, δNDVI-, PCA-, and ICAbased CDMRF are summarized in Table 2. As shown in Table 2, δNDVI, PCA-, and ICA-based CDMRF methods have overwhelming advantages over CVA-based CDMRF in all the test sites, as indicated by the completeness, correctness, Kappa coefficient, and F-measures. The landslide mapping completeness and correctness of CVA-based CDMRF are less than 60% and the Kappa coefficients and F-measures are less than 0.6 in all the study sites. By contrast, the improved methods especially PCA-and ICA-based CDMRF in this study can achieve much more accurate landslide mapping results with correctness and Kappa coefficients greater than 0.75 in single sensor data tests and greater than 0.70 in two sensor data test.

Discussion
This study improved CDMRF approach by integrating δNDVI, PCA, and ICA generated CDIs with MRF for landslide inventory mapping from multi-sensor data in different areas. Despite the good performance in landslide mapping, the improved CDMRF methods are not perfect and there is room for further research.
One main error source may come from the selected change detection methods used for CDI generation. In particular, for PCA and ICA, cares should be taken in the choice of components that primarily represent landslides. Since PCA and ICA are the methods that extract global information from the entire image, namely they are driven by image statistics instead of landslides, components that contain the strongest landslide information may vary in different study areas. Landslides are heterogeneous from a spectroscopic point of view may be expressed in more than one component. Thus, the optimal PCA and ICA components that represent landslides were determined by manual inspection. In this respect, PCA-and ICA-based CDMRF are not fully automated. Nevertheless, they are promising to be used for rapid landslide mapping due to the following advantages. First, although the optimal components of PCA and ICA for landslide mapping are unknown in advance, they are not hard to determine in real applications. Land cover changes due to landslides are rarely found in the first component, which often represents the most unchanged land covers (Mondini et al., 2011). Also, minor components (e.g., beyond the fourth component) are mostly represented by noise (Lu et al., 2011). Second, PCA and ICA are more robust to noise than CVA because they transform noise into minor components, thus reducing data redundancy and Fig. 11. An example of commission errors (indicated by the yellow arrows) in Messina landslide mapping. (a) and (b) The landslide reference data and its binary map. Landslide mapping results of CVA-based CDMRF (c) and its binary map (d), δNDVI-based CDMRF (e) and its binary map (f), PCA-based CDMRF (g) and its binary map (h), and ICA-based CDMRF (i) and its binary map (j). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) enhancing discrepancies between different spectral information (Lu et al., 2004) and bringing great noise separation capabilities (Thomas et al., 2002). However, the image preprocessing requirement such as the radiometric correction to keep the similar spectral and radiometric characteristics between the pre-and post-event images for CVA is generally higher than PCA and ICA. Finally, from the perspective of landslide mapping accuracy (Table 2), PCA-and ICA-based CDMRF are much better than CVA-based competitor in terms of completeness, correctness, Kappa coefficient, and F-measures.
Training samples of landslides in this study were automatically generated from CDIs of δNDVI, PCA, and ICA to avoid manual selection. These training samples can be reliably produced when the assumption that higher CDI values indicate the higher probability of landslide occurrence is met. However, this assumption may not always be achieved. For instance, the relation of CDIs to landslide may be disturbed by other land cover changes such as harvesting, logging, and construction. Also such relations may be confused by changes resulting from clouds, haze, image misregistration, and so on. Moreover, this relation may be obscured by the natural variations of landslide process. For example, even for the same type of mass movement, landslides may occur over both dense and sparse vegetation covers, thereby rendering different CDI values. Under these circumstances, additional processing steps may be necessary, such as the removal of built-up areas, as conducted in Messina landslide mapping.
In terms of parameter tuning, T and ∆T are the only two parameters that were determined through trial and error in this study. First, T is the parameter that determines the lower threshold to differentiate training samples of landslides and non-landslides. However, δNDVI, PCA, and ICA are different methods in principle, thus the values of T for δNDVI, PCA, and ICA cannot be identical and were manually tuned in this study to ensure the proper thresholds used for training sample generation. Second, ∆T is the parameter that provides incremental value to determine uncertain pixels for further MRF analysis. The parameter ∆T is ideally to be specified through trial and error. The uncertainty caused by the parameter tuning for ∆T is considered tolerable given that those uncertain pixels can be further labeled by MRF. Third, the parameter value of λ, which is a weighting coefficient that balances the unary and pairwise potentials in the MRF framework, was set at a fixed value of 50 as optimized in (Rother et al., 2004;Szeliski et al., 2008;Li et al., 2016a). This study attempts to avoid manual tuning of λ because it primarily determines the level of smoothness and accordingly has less influential than T that generates the landslide training samples.
The implementation of CDMRF in Messina appears to have more challenges. Due to the long time interval between the pre-and postevent images, the buildings were masked out using the available topographic data in order to reduce the potential changes caused by building construction or demolition. In places with clustered habitations, such as Messina, these vectorized topographic data for built-up areas may often be available. However, in areas where habitations are scattered, such built-up topographic data may be difficult to acquire. Thus, the preference of CDMRF is to use pre-and post-event images with a short time interval, preferably a couple of days before and after the event, reducing the possible impacts of construction. If such images with a short time interval cannot be acquired, it is suggested to mask out buildings in advance using image processing techniques such as anomaly detection (Lu et al., 2011), GrabCut partitioning (Ok et al., 2013), level set methods (Li et al., 2014(Li et al., , 2015, or morphological building index (Huang et al., 2014). In this study, buildings were masked out in advance but the roads remained because no topographic road data are available in the study area. Thus, both PCA-and ICAbased CDMRF detected the spectrally similar roads as landslides, leading to over-detection in some sub-areas, as can be seen in Figs. 12(g)-(j).
The lower accuracy of CVA-based CDMRF is mainly due to the following reasons. CVA is essentially image differencing technique that subtracts pre-event image from post-event image to obtain land cover changes. However, image differencing results are typically contaminated by errors due to the discrepancies in phenological and illumination conditions between pre-and post-event images, thus generating commission errors. This problem may also exist in δNDVI-based CDMRF which is also based on image differencing technique. For example, in Messina landslide mapping, as shown in Fig. 11, shadows in post-event image (as indicated by the yellow arrows) cause spectral differences in bi-temporal NDVI images and consequently they were misclassified as landslides by CVA-and δNDVI-based methods. By contrast, the landslides occurred over the spectrally similar bare rocks (as indicated by the blue arrows) are hard to detect by image differencing techniques (Fig. 12), thereby resulting in omission errors in results of CVA-and δNDVI-based methods. By contrast, PCA-and ICAbased methods can achieve much better performance in Messina. In addition, both CVA and δNDVI approaches are essentially pixel-based image differencing methods, thus producing many salt and pepper noise, as shown in Figs. 11 and 12. These errors are also reflected by the lower Kappa coefficients of CVA-and δNDVI-based methods.
Other errors may arise from image preprocessing, particularly in image co-registration, ortho-rectification, pan-sharpening, and relative radiometric correction. These have already been elaborately discussed by Mondini et al. (2011) and Li et al. (2016a).
In this study, the improved methods combine change detection techniques and MRF for landslide mapping from multi-sensor data. They take advantage of the spectral and spatial contextual information of bi-temporal remote sensing images. Since landslide occurrences are often related to terrain, it is interesting to further incorporate topographic factors (e.g., main scarp, toes, and accumulation zones) into MRF or other advanced machine learning algorithms for landslide mapping. Li et al. (2016b) have proposed a CVA-based level set evolution (LSE) method for landslide mapping. The change detection techniques (i.e., δNDVI, PCA, and ICA) used in this study can also be integrated with LSE for landslide inventory preparation. Future work should focus on developing new change detection techniques that can readily generate quality CDIs and applying new and advanced machine learning algorithms (e.g., deep learning) to automatic natural hazards inventory mapping from remote sensing images. Future research examining the improved methods for mapping historical landslide inventories using multi-sensor time series data such as the studies conducted by Martha et al. (2013) and Behling et al. (2016) is also encouraged.
Future improvements can be directed at preparing landslide inventory that can be readily used for subsequent landslide hazard and risk assessment. Compiling a multi-temporal landslide inventory as discussed above is definitely helpful for decision makers because it can provide fundamental information on landslide spatial and temporal frequencies. Besides, landslide characteristics such as landslide activity, type of mass movement, linkage to specified predisposing and triggering factors, and landslide magnitude (or even intensity) contribute greatly to estimating magnitude-frequency relations (Corominas et al., 2014). If CDMRF can be used for extracting these landslide characteristics and producing a more detailed inventory, it may greatly facilitate quantitative landslide hazard and risk assessment.

Conclusion
This study improved the change detection-based Markov random field (CDMRF) approach by integrating normalized difference vegetation index (NDVI), principal component analysis (PCA), and independent component analysis (ICA) generated change detection images (CDIs) with MRF for landslide mapping from multi-sensor data. To corroborate the effectiveness and applicability, the improved methods were applied to the pre-and post-event very high resolution QuickBird, high resolution FORMOSAT-2, and moderate resolution Sentinel-2 images, for mapping landslides that occurred after three catastrophic events triggered by rainfall, typhoon and earthquake, respectively. In addition, the improved methods were also tested on the pre-event Landsat-8 and the post-event Sentinel-2 datasets, two encouraging and free moderate resolution satellite data that can be combined for more frequent multi-temporal inventory mapping. Quantitative results show that δNDVI-, PCA-, and ICA-based CDMRF perform much better than CVA-based CDMRF (Li et al., 2016a) in terms of completeness, correctness, Kappa coefficient, and F-Measures. Future research developing new change detection techniques that can readily generate quality CDIs and using advanced machine learning algorithms (e.g., deep learning) for natural hazard mapping are encouraged. Also, CDMRF deserves further attention for preparation of historical landslide inventories from multi-sensor time series data.