LANDSLIDE INVENTORY MAPPING FROM BITEMPORAL 10 m SENTINEL-2 IMAGES USING CHANGE DETECTION BASED MARKOV RANDOM FIELD

Landslide inventory mapping is essential for hazard assessment and mitigation. In most previous studies, landslide mapping was achieved by visual interpretation of aerial photos and remote sensing images. However, such method is labor-intensive and timeconsuming, especially over large areas. Although a number of semi-automatic landslide mapping methods have been proposed over the past few years, limitations remain in terms of their applicability over different study areas and data, and there is large room for improvement in terms of the accuracy and automation degree. For these reasons, we developed a change detection-based Markov Random Field (CDMRF) method for landslide inventory mapping. The proposed method mainly includes two steps: 1) change detection-based multi-threshold for training samples generation and 2) MRF for landslide inventory mapping. Compared with the previous methods, the proposed method in this study has three advantages: 1) it combines multiple image difference techniques with multi-threshold method to generate reliable training samples; 2) it takes the spectral characteristics of landslides into account; and 3) it is highly automatic with little parameter tuning. The proposed method was applied for regional landslides mapping from 10 m Sentinel-2 images in Western China. Results corroborated the effectiveness and applicability of the proposed method especially the capability of rapid landslide mapping. Some directions for future research are offered. This study to our knowledge is the first attempt to map landslides from free and medium resolution satellite (i.e., Sentinel-2) images in China. * Corresponding author


INTRODUCTION
In recent years, earthquake or rainfall-induced landslide hazard frequently occurred in China, leading to heavy casualties and property losses (Xu, 2012a).For landslide hazard assessment and mitigation, landslide inventory mapping recording geographic location, occurrence date, magnitude and type of landslide is required.In most previous studies, landslide inventory mapping was achieved by visual interpretation of aerial photos or remote sensing images, which however is laborintensive and time-consuming.Over the past few years, enormous methods have been proposed for landslide mapping from high resolution (HR) and even very high resolution (VHR) remote sensing images.However, limitations remain in terms of their applicability and accuracy in different study areas and data.

Previous Work
Previous work can be roughly divided into two parts, pixelbased and object-based.Pixel-based methods generally take advantage of the spectral characteristic of every pixel.Lu et al. (2004a) provided a comprehensive review on change detection approaches in remote sensing community.In Van Westen and Getahun (2003a), multi-temporal landslide maps were obtained by visual interpretation of sequential aerial photos analyse evolution of the Tessina landslide in Italy.In Mondini et al. (2011a), four change variables related to landslides occurrence were selected to calibrate landslide classification models.In Hervas et al. (2003a), automatic change detection from bitemporal aerial photographs was used to map landslides in Italy.In Li et al. (2016a), change detection and level set method using bitemporal aerial photographs were used to map landslides in Hong Kong.In addition to change detection methods, some other methods are proposed for landslide inventory mapping.In Yang et al. (2014a), multi-temporal satellite images were used to obtain the landslide distribution map in Sichuan earthquake area on May 12, 2008. In Chen et al. (2014a), random forest was applied to detect forested landslides.However, pixel-based methods only consider spectral information and ignore the similarity among adjacent pixels, which causes a lot of salt and pepper noise in the final mapping result (Hölbling, 2017a).
Object-based methods regard adjacent pixels with similar spectral, spatial, texture information as an object.A large number of object-based landslide inventory mapping methods were developed using eCognition software.In Martha et al. (2010a), an object-oriented method using multispectral data and a digital terrain model were used to detect landslides.In Lu et al. (2011a), a semi-automatic object-oriented approach using bitemporal VHR images was applied to map the Messina landslide in Italy.In Stump et al. (2011a), object-oriented image analysis was combined with the random forest algorithm for landslide mapping.In Rau et al. (2014a), a semi-automatic object-oriented technology with multi-sensor and multispectral imagery and DEM was proposed for landslide recognition.In Pradhan et al. (2015a), an object-oriented classification method with fused data between LiDAR and VHR imagery was used to detect landslides.
Although object-based method is widely used for landslide inventory mapping, limitations still remain, such as: 1) it is difficult to find general principles to delineate different landslides because of variant features existing; 2) the labourintensive parameter tuning causes the low automation degree of landslide mapping; and 3) the accuracy of landslide mapping can be further improved.Li et al. (2016a) proposed a change detection-based Markov random field (CDMRF) for landslide inventory mapping from 0.5 m bitemporal aerial orthophotos over Lantau island, Hong Kong.Experimental results demonstrated that this method can map landslides accurately and efficiently.However, this method is far from perfect.We found that it has difficulty in mapping complex and low-contrast landslides from some remote sensing images.To improve the applicability and accuracy of CDMRF method over large areas, we made a further development of the method proposed in (Li et al. 2016a) and proposed a new CDMRF method for landslide inventory mapping from bitemporal 10 m Sentinel-2 multispectral images.

Study Area
The study area is located on Jiuzhaigou Country, northeast of Tibetan Qiang Autonomous Prefecture of Ngawa, Sichuan Province, China, with an area about 16.5 km 2 .This area features a subtropical monsoon climate with an annual average temperature of 12.7 ℃ and an annual precipitation of 550 mm.The terrain is dominated by mountains covered with dense vegetation with an annual average altitude of 4000 m (Liu, 2017a).An earthquake of 7.0 magnitude occurred in Jiuzhaigou Country at 21:19:46 Beijing Time on August 8, 2017.The epicenter is located at 33.20°N and 103.82°E, and the depth of focus is 20 km.The major earthquake induced a large number of shallow landslides which caused vegetation destruction and mountains exposed over seismic hazard areas.The landslides in the study area were mainly distributed along National Highway 544 from Huanglong Airport to Zhangzha Town.

Data
Sentinel-2 satellite is a multi-spectral Earth observation system implemented by Global Monitoring for Environment and Security (GMES) and jointly organized by the European Commission (EC) and the European Space Agency (ESA) (Drusch, 2012a).This system includes two polar-orbiting satellites (i.e., Sentinel-2A and Sentinel-2B) in the same orbit but phased at 180°to each other.Together, the two satellites have a high revisit time of 5 days.The satellites carry the Multi Spectral Instrument (MSI) that samples 13 spectral bands (10 m, 20 m, 60 m bands) from visible and near-infrared (VNIR) bands to short wave infrared bands (SWIR).Level-1C and Level-2A product levels are freely available to users for land monitoring, emergency management and risk mapping (ESA, 2015).The cloud-free pre-and post-event Sentinel-2 Level-1C multispectral images were acquired at 3:45:41 UTC on July 29, 2017 and 3:45:31 UTC on September 7, 2017, respectively and they were downloaded from the ESA Sentinel Online website (The Copernicus Open Access Hub, 2018).Four bands, i.e., the near-infrared, red, green and blue bands with 10 m spatial resolution, were used in this study.
A landslide inventory map manually digitized from the pre-and post-event Sentinel-2 images was used as a reference for accuracy assessment.

METHODOLOGY
The proposed method includes four sub-steps.First, atmospheric correction was performed on Level-1C product to obtain Level-2A product, i.e., Bottom of Atmosphere (BOA) reflectance images.Then the near-infrared (NIR) and red bands of the pre-and post-event images were used to generate change detection images (CDI).Next, a multi-threshold method-based on CDI was used to generate training samples of landslides and non-landslides.Finally, the result of landslide inventory mapping is accomplished by using MRF.

Preprocessing of Remote Sensing Images
Radiometric and geometric corrections have been carried out on acquired Level-1C product, including ortho-rectification and spatial registration.For the subsequent bitemporal images process and analysis, the Level-2A product was derived from the Level-1C products through the Sentinel-2 Toolbox.

Change Detection Images
Two image difference techniques were selected to identify shallow landslides between the pre-and post-event images.
The first image difference technique is achieved using Normalized Difference Vegetation Index (NDVI) (Mondini, 2011a).NDVI is an effective indicator of vegetation growth and cover (Lu, 2004a).The occurrence of shallow landslides leads to vegetation destruction, thus NDVI is widely used to detect landslides.In this study, changes in the pre-and post-event NDVI (δNDVI) (Fig. 3(a)) were obtained by the following equation (Mondini, 2011a): The second image difference technique is based on principal component analysis (PCA).PCA converts multispectral bands from multi-temporal images into linear independence components, which can reduce noise and data redundancy and increase separability between bands, and the first principal component always concentrates the most information in images and other components account for lesser information (Lu, 2004a).Similar to NDVI, four bands, i.e., the pre-and postevent near-infrared and red bands, were used to generate four linear independence components (Mondini, 2011a).Through visual interpretation of PCA components, we found that the fourth component (PCA-4) (Fig. 3(c)) has a good match with landslides.

Multi-Thresholding for Training Samples Generation
δNDVI and PCA-4 had a high correlation with the occurrence of landslides, thus they are selected to generate landslides and non-landslides training samples.
In CDI, the brighter pixels indicate larger changes which can be considered as landslides, and the darker pixels often change a little which can be considered as non-landslides.Thus, multithresholding method were applied to generate training samples of landslides and non-landslides as the following equation (Xian, 2009a;Li, 2016a):

Markov Random Field for Landslide Mapping
First, landslide mapping is transformed into pixel labelling problem.MRF assigns a label 1 or 0 for every pixel in uncertain areas (1 for landslide and 0 for non-landslide), which forms a label set.And then an energy function for the label set is established using the joint probability distribution of MRF as the equation (3).Next the energy function is constructed as an undirected graph.Finally, the max-flow algorithm is used to segment the undirected graph to get the globally optimal minimum cut, which can achieve the optimal landslide mapping (Boykov, 2004a).
where Eu(L) = the unary potential Ep(L) = the pairwise potential L = (l1, l2, …, ln), a label set and li represents the label of the ith pixel λ = a coefficient reflecting the relative importance between the unary and pairwise potential where C1 = the single-site clique log(p(O|Ii)) = the posterior probability of pixels in uncertain areas belonging to the object O (landslide) log(p(B|Ii)) = the posterior probability of pixels in uncertain areas belonging to the background B (nonlandslide) Vi(li) = negative logarithmic of the posterior probability

3.4.2
The pairwise potential: Ep(L) represents the level of similarity between adjacent pixels (4-neighborhood).It is defined as the following equation ( 6), ( 7) and ( 8): where C2 = the pair-site clique (Ii-Ij) 2 = the spectral differences of red, green and blue of visible bands between adjacent pixels

II
, and represents the expectation operator on the whole image

Energy minimization:
The energy function is constructed as an undirected graph G = <V, E>.V and E represent the set of vertex and edge, respectively.In addition, the graph has two kinds of special vertex set, S and T, which represent landslide and non-landslide training samples, respectively.In this graph, the weight of edges connecting pixels in V with S and T is decided by the unary potential, and the weight of edges connecting neighbouring pixels in V is decided by the pairwise potential.A cut in the graph G refers to a set of edges, which cut them can cause the set V divided into two disjoint sets, VO for landslide and VB for non-landslide.If the sum of weights of edges reaches a minimum, it's called the minimum cut.To achieve the optimal landslide mapping, the max-flow algorithm (Boykov, 2004a) is used to segment the undirected graph to get the globally optimal minimum cut.

RESULTS AND ACCURACY ASSESSMENT
By combining change detection based multi-threshold method with MRF, two landslide inventory mappings were obtained as shown in Figure 4. To evaluate the accuracy of the proposed method, these two mappings were compared with the manually digitized reference mapping (Fig. 1(c)).And three accuracy assessment indices: Completeness, Correctness and Quality are used as follows (Li, 2016a): It has been calculated that in the study area, using bitemporal NDVI (δNDVI) for CDI and MRF for mapping landslides, the completeness, correctness, and quality of LIM achieved 79.67%, 97.12% and 76.07%, respectively.Using PCA-4 instead of δNDVI, the completeness of 95.93%, the correctness of 88.03%, and the quality of 76.08% were achieved.It shows clearly that in the study area, using δNDVI for CDI can achieve more correct mapping result but the omission error is more than PCA-4.And the PCA-4 would cause excessive landslide mapping.The quality of both methods is almost the same, which indicates the proposed new CDMRF can achieve mapping landslide inventory rapidly and accurately over large areas.

CONCLUSION
In this study, we integrated two image difference technologies into the CDMRF for automatic and accurate landslide inventory mapping from bitemporal 10 m multispectral images.First, a difference image based on bitemporal NDVI was generated and the fourth band image of principal component analysis (PCA-4) was generated using the pre-and post-event near-infrared (NIR) and red bands.Second, multi-threshold method based on generated change detection images (CDI) was used to generate training sample masks, which include landslides, non-landslides and uncertain pixels.Finally, using Markov random field (MRF) with training sample masks onto the post-event Sentinel-2 image achieved the optimal landslide inventory mapping (LIM).
The proposed method was used to map landslides in an area about 16.5 km 2 in Western China.δNDVI for CDI achieved completeness and correctness of 79.67% and 97.12% and PCA-4 for CDI achieved completeness and correctness of 95.93% and 88.03%.The quality of these two methods achieved 76.07% and 76.08%, respectively.
Experimental results verified the effectiveness of the proposed method: 1) it takes advantage of spectral characteristics of landslides; 2) it takes into account the spatial contextual information in post-event remote sensing image, which makes up for the deficiency of ignoring spatial characteristic of pixels in traditional pixel-based methods; 3) it has a high degree of automation due to little parameter tuning; and 4) it is a generic landslide inventory mapping method.This method has great potential to provide a technical support for future landslide inventory mapping over large areas in China.

Figure 2 .
Figure 2. Flowchart of the proposed landslide mapping method near-infrared band red = red band pre = the pre-event image post = the post-event image Figure 3. (a) Change detection image of bitemporal NDVI.(b) Training sample mask of δNDVI.(c) Fourth component of PCA.(d) Training sample mask of PCA-4.
unary potential: Eu(L) represents the level of similarity between the label set L and training samples.It is defined as the following equation (4) and (5):

Figure 4 .
Figure 4. Landslide mapping results of MRF.(a) Binary result of using δNDVI.(b) Landslides in (a) overlaid on the post-event image.(c) Binary result of using PCA-4.(d) Landslides in (c) overlaid on the post-event image.