Noncontact detection of earthquake-induced landslides by an enhanced image binarization method incorporating with Monte-Carlo simulation

Abstract

Detecting landslides using remote sensing images involves converting gray level images into binary images. Given the complex background and non-uniform illumination in the regional remote sensing image, the commonly-used global thresholding methods are limited to interpret some landslides those hidden in the shadow. In this paper, we report on an enhanced image thresholding method for the noncontact detection of earthquake-induced landslides. The proposed method incorporates Monte-Carlo simulation into the local thresholding method, and essential issues, regarding complex illumination condition and uncertainties of determining block size in local thresholding method, are addressed. To better separate landslide candidate objects from the background, we incorporate Digital Surface Model (DSM) into the binary image, such that the interferences by built-up areas, terraces, and rivers can be significantly reduced considering the indicator of slope gradient. The described method has been tested using two benchmark tests, showing that the proposed method performs well dealing with the complex background and illumination condition. As a case study, we use the proposed method to detect landslides in a remote sensing image near Beichuan area after the 2008 Ms8.0 Wenchuan earthquake. Results demonstrate that the presented method interprets landslides resembling those of manually visual delineation.

Keywords:

• Landslides
• Image binarization
• Local thresholding method
• Monte-Carlo simulation
• The Wenchuan Earthquake

1. Introduction

Landslides are often triggered by strong earthquakes (Zhang et al. 2015). These earthquake-induced landslides are a kind of dynamic process reshaping landscape, endangering human lives and infrastructure facilities in mountainous regions. The 1999 Ms7.6 Chi-Chi earthquake and the 2008 Ms8.0 Wenchuan earthquake for instance, respectively 9,272 and 60,104 landslides were reportedly induced (Chang et al. 2005; Yin et al. 2009; Gorum et al. 2011). These landslides could be remobilized by the heavy precipitations and transit into debris flows, posing a great threat to the local society in the following decades (Han et al. 2015a; Han et al. 2018). In this context, detecting and delineating these earthquake-induced landslides is therefore of great relevance for disaster prevention and mitigation (Rosin and Hervás 2005). Several lines of evidence have also supported the importance of landslides detection (Parker et al. 2011; Han et al. 2015b; Roback et al. 2018).

With the boom of high spatial resolution optical imagery from multi-platforms (e.g., orbital satellites, aerial crafts), remote sensing technique is increasingly applied not only in land management but also to assist in delineating landslide surface boundaries. Synthetic Aperture Radar (SAR) is one optional solution for this purpose. Differential SAR interferometry (DInSAR) has proved able to detect centimetric displacements for earthquake-induced landslides (Vietmeier et al. 1999; Refice et al. 2000). However, the application of DInSAR technique in landslide detection is often constrained because of the limitation with respect to inadequate spatial resolution, atmospheric effects, and especially signal phase decorrelation due to vegetation and soil moisture changes (Rosin and Hervás 2005).

As an indispensable component of remote sensing technique, digital aerophotogrammetry has proved feasible to detect the change of earth surface because of their high resolution and low cost. Usually, panchromatic (PAN) images are available for a catastrophic earthquake event. Landslides could be distinguished by their gray levels, which are used as an evidence for delineating these landslides. Traditional visual interpretation of these remote sensing photographs is the major solution for delineating landslides. However, artificially interpretation is rather experience-based, and the interpretation process is tedious and time-consumed (Li et al. 2014). In this context, it will be beneficial to pre-process these remote sensing images by highlighting the landslides before interpretation.

Binary image representation in textual image techniques is one of the preferred formats for image analysis (Stathis et al. 2008). It aims at distinguishing foreground from background. Kefali et al. (2010) illustrated that the binarization process decreases the presence of unwanted data and preserves only the desired data in the image. This technique could be applied in landslides detection in remote sensing images, by converting all pixels into two categories, i.e. landslides, and background, with separated colors (e.g. black and white). Pre-processing of remote sensing images removes the unnecessary information to increase the visibility of the useful information (Bataineh et al. 2011; Lech et al. 2015) and therefore accelerates the extraction and recognition process of the landslides.

A robust, automatic, and efficient binarization algorithm is required to attain a satisfactory binary result. Following Sezgin and Sankur (2004) and Chou et al. (2010), Wen et al. (2013) classified the up-to-date binarization methods into three major categories, clustering-based methods, threshold-based methods, and hybrid methods. In general, the threshold-based methods permit to rapidly binarize the images with one or many thresholds. Common wisdom holds that the threshold-based methods include two sub-categories: global thresholding methods and local thresholding methods (Stathis et al. 2008; Kefali et al. 2010; Bataineh et al. 2011; Wen et al. 2013). The global thresholding methods derive an overall threshold for the entire image. Each pixel of the image is binarized according to this threshold. Remarkable work is Otsu’s algorithm. As introduced in Pai et al. (2010), Otsu’s algorithm decides a global threshold based on the maximum of the between-class variances. The global thresholding method is advanced in its simple implementation and high computation efficiency (Poletti et al. 2012).

Logically, the simple thresholding methods will be better when using pre-processing and post-processing (Bataineh et al. 2011). In fact, in our previous studies (Li et al. 2014), we classified the histogram of the images into three categories and claimed that the global thresholding methods perform well when the histogram approximates a balanced bimodal type (document image for instance in Figure 1a). However, regional remote sensing image commonly covers a wide range of landscape, containing various background patterns, multiple noises and illumination levels, shadow areas also increases the interference. As shown in Figure 1b, the histogram of the remote sensing image approximates a mixed unimodal type. In this case, the global thresholding methods are often limited to interpret objects from the image (Gatos et al. 2006).

Figure 1. Histogram of different image sources. (a) document image with homogeneous background and uniform illumination (Image from Wen et al. 2013). (b) remote sensing image with inhomogeneous background and complex illumination (Image was taken after Ms 7.0 Lushan earthquake in April 20, 2013).

To solve the existing limitations of the global thresholding method, the local thresholding methods become a means of solution. The local thresholding methods compute many local thresholds based on the information contained in the neighborhood of each pixel, or a small region of the image independently. Hence, they are also referred as the local adaptive thresholding methods (Pai et al. 2010). Some remarkable works have been done by previous studies (e.g., Bernsen 1986; Niblack 1986; Sauvola and Pietikainen 2000; Yanowitz and Bruckstein 1989; Taxt et al. 1989; Eikvil et al. 1991; Gatos et al. 2011; Wen et al. 2013). Following studies also proposed some compound algorithms to improve the performance of the local threshold method (e.g., Chou et al. 2010; Tong et al. 2009; Zhang and Yang 2010). By adding some preparatory steps to an existing simple thresholding method, effect of varying levels of noise and complex illumination on the binarization can be minimized.

However, these sophisticated compound methods are often costly and complex to design. Moreover, the value of some required factors should be determined manually. Some studies (e.g., Trier and Jain 1995; Khurshid et al. 2010; Pai et al. 2010) have evaluated well-established local thresholding methods, indicating that Niblack’s method, Sauvola’s method, and NICK’s method may perform better. Unfortunately, as illustrated by Bataineh et al. (2011), the above-mentioned methods are based on adaptive thresholding algorithm in fixed square blocks, which means they are scale-dependent and require empirical tuning of the block size. Nevertheless, block size has long been substantiated having a significant influence on the binarization results (Wu and Amin 2003). Uncertainties of block size, therefore, challenges the application of the local threshold method. Firstly, the best-fitting block size varies case by case, lacking a widely accepted criterion. Secondly, local binarization in the fixed block is likely to generate interference left between two neighbored blocks (Wen et al. 2013), because the binarization result in the neighbored block may be opposite. We call this kind of interference ‘the boundary effect’ and find that it is difficult to eliminate in post-processing. As such, although the concept of the intelligent block has been proposed by Pai et al. (2010), challenges still remain in the field of image binarization, especially in the application of landslides detection using regional remote sensing image.

The object of this study is to present an enhanced local binarization method for landslides detection using remote sensing image. The proposed method incorporates Monte-Carlo simulation into the local thresholding method, and two essential issues, regarding complex illumination condition and uncertainties of determining block size in local thresholding method, are addressed. To better separate landslide candidate objects from the background, we incorporate Digital Surface Model (DSM) into the binary image, such that the interferences by built-up areas and rivers can be significantly reduced considering the indicator of slope gradient. Two benchmark tests are introduced to evaluate the performance of the presented method, and the application of landslides detection in Beichuan area after Ms 8.0 Wenchuan earthquake is conducted. Results show that it will be an alternative solution for dealing with the images of complex background and will be suitable for detecting landslides using remote sensing image.

2. Methodology

The presented method is based on Otsu’s thresholding algorithm; consequently, we begin our study with a review of Otsu’s algorithm.

2.1. Review of Otsu’s thresholding algorithm

Otsu’s thresholding algorithm (Otsu 1979) is a clustering-based automatic thresholding method. It has been regarded as one of the best global techniques and measures region homogeneity by variance.

Generally, pixels of a given image can be represented in $L$ gray levels $\left[1,\mathrm{ }2,\dots ,L\right]$. The purpose of Otsu’s method is to derive a threshold $T_o$ to separate the gray levels into two classes, $C_1=\left[1,\mathrm{ }2,\dots ,T_o\right]$ representing objects, and $C_2=\left[T_o+1,T_o+2,\dots ,L\right]$ representing background. The number of pixels at gray level $i$ is denoted by $z_i$. The total number of the pixels of the image can be expressed as (1) $$Z=\sum _{i=1}^L{z}_i$$(1)

The gray-level histogram is normalized and regarded as a probability distribution (2) $$p\left(i\right)=\frac{z_i}{Z},z_i\ge 0,\sum _{i=1}^{L}p\left(i\right)=1$$(2)

As demonstrated above, pixels are separated into two classes $C_1$ and $C_2$. The probabilities of each class occurrence can be expressed as (3) $${\omega }_1=\sum _{i=1}^{T_o}p\left(i\right)=\omega \left(T\right)$$(3)

(4) $${\omega }_2=\sum _{i=T_o+1}^{L}p\left(i\right)=1-\omega \left(T\right)$$(4)

and the class mean levels of each class can be written as (5) $${\mu }_1=\sum _{i=1}^{T_o}\frac{ip\left(i\right)}{{\omega }_1}$$(5)

(6) $${\mu }_2=\sum _{i=T_o+1}^L\frac{ip\left(i\right)}{{\omega }_2}$$(6)

(7) $${\mu }_T=\sum _{i=1}^{L}ip\left(i\right)={\omega }_1{\mu }_1+{\omega }_2{\mu }_2$$(7)

where ${\omega }_1$ and ${\omega }_2$ denote the probabilities of object class and background classes, respectively. ${\mu }_1$, ${\mu }_2$, and ${\mu }_T$ refer to the average gray levels of the object, the background, and the global image. Given these variables, the class variances are (8) $${\sigma }_1^2=\sum _{i=1}^{T_o}p\left(i\right)\frac{{\left(i-{\mu }_1\right)}^2}{{\omega }_1}$$(8)

(9) $${\sigma }_2^2=\sum _{i=T_o+1}^{L}p\left(i\right)\frac{{\left(i-{\mu }_2\right)}^2}{{\omega }_2}$$(9)

where ${\sigma }_1^2$ and ${\sigma }_2^2$ are the variances of the object class and the background respectively. Equation (8–9) and Equation (3–4) reduce to the within-class variance ${\sigma }_W^2$, the between-class variance ${\sigma }_B^2$, and the total variance of levels ${\sigma }_T^2$, which equate (10) $${\sigma }_W^2={\omega }_1{\sigma }_1+{\omega }_2{\sigma }_2$$(10)

(11) $${\sigma }_B^2={{\omega }_1\left({\mu }_1-{\mu }_T\right)}^2+{{\omega }_2\left({\mu }_2-{\mu }_T\right)}^2$$(11) (12) $${\sigma }_T^2=\sum _{i=1}^L{\left(i-{\mu }_T\right)}^{2}p\left(i\right)$$(12)

In the discrimination analysis, the separable degree $\eta$ of the class is (13) $$\eta =\underset{1\le T\le L}{\max }\frac{{\sigma }_B^2}{{\sigma }_W^2}$$(13)

Finally, a threshold $T_{ot}$, resulting in maximal $\eta$, is selected as the optimal threshold in the Otsu’s global thresholding method.

Comparing to Niblack’s and Sauvola’s local thresholding method as introduced in Section 1, although Otsu’s method is much easier in selecting a simple threshold for the entire image, avoiding uncertainties of block size, it is limited to yield acceptable binary result when images contain complex background and illumination.

2.2. The proposed thresholding method

The aim of the presented method is to address the binarization challenge due to complex background and illumination when dealing with regional remote sensing image. In line with the inspirations of previous local thresholding methods, we propose an enhanced local thresholding method integrating with Monte-Carlo simulation. Notably, the proposed method is beneficial to avoid the uncertainties of local block size and minimize the interference due to the so-called boundary effect. The proposed method is composed of four major steps as described below.

2.2.1. Configuration of Monte-Carlo simulation

The first step of the proposed method is to configure Monte-Carlo simulation. The Monte-Carlo method uses repeated random sampling to simulate data for a given mathematical model and evaluate the outcome. In this sense, the core of the Monte-Carlo simulation is to create a large, random data set such that uncertainties of inputs can be covered by this data set (Han et al. 2015). The total number of this data set is named of the Monte-Carlo steps $MCS$. For instance, $MCS=100$ refers to that the procedure will repeat 100 times. Ideally, the greater $MCS$ generates better data set for the following analysis. However, it linearly increases the amount of computation.

2.2.2. Image partitioning with random block size

In each step of Monte-Carlo simulation, the original gray-level image is imported. $f\left(x,y\right)$ refers to the gray-level image with size $M\times N$, where $\left(x,y\right)$ represents the coordinate of each pixel in the image. We use square blocks with uniform size to partition the entire image. As mentioned in Section 1, block size has a significant influence on the binarization results (Wu and Amin 2003), but it is scale-dependent and requires empirical tuning. In our study, the block size notated as $D$ is regarded as the only random variable in the Monte-Carlo simulation. It obeys the uniform distribution ranging between $D_{min}$ and $D_{max}$.

As shown in Figure 2, for most of images, the rows $M$ does not strictly equate the columns $N$, which means the square blocks could not completely cover the entire image. In this case, we use rectangular blocks of different sizes to fill the remaining fragmented part. The sizes of the remaining blocks are calculated by (14) $$D^{\mathrm{\prime }}=M-⌊M/D⌋D$$(14)

(15) $$D^{\mathrm{″}}=N-⌊N/D⌋D$$(15)

where $D^{\mathrm{\prime }}$ and $D^{\mathrm{″}}$ are the remainders of the term $\left(M/D\right)$ and $\left(N/D\right)$, respectively. In this sense, four types of blocks, with sizes of$\mathrm{ }\left(D,D\right)$, $\left(D^{\mathrm{\prime }},D\right),\mathrm{ }\left(D,D^{\text{″}}\right)$, and $\left(D^{\mathrm{\prime }},D^{\mathrm{″}}\right)$, are used to partition the entire image in a Monte-Carlo step.

Figure 2. Schematic illustration for image partitioning using blocks in a Monte-Carlo step. Blue region represents the area that partitioned by square blocks, while red region denotes the remaining fragmented area that partitioned by the rectangular blocks.

2.2.3. Image binarization

In each Monte-Carlo step $k$, the sub-images of the blocks are independently binarized using Otsu’s algorithm. The whole binary image is subsequently obtained by integrating the amount of the binary images of blocks. In order to reduce the computational complexity, the gray level histograms of the sub-images are preliminarily checked. We use the criterion by Pai et al. (2010) to determine whether the block is required to evaluate a local threshold. If the standard deviation $s$ and the mean $u$ of the sub-image satisfy the following condition (16) $$s>s_T\mathrm{ }\&\&\mathrm{ }u<u_T$$(16)

where $s_T$ and $u_T$ denote the standard deviation and the mean of the whole image, respectively, the block is assumed as background-dominated, such that the sub-image of the block contains only background pixels. In this case, the pixels in the block are substituted by background. In contrast, we used Otsu’s thresholding method as described in Section 2.1 to calculate optimal local thresholds for each block. In each block, the pixel that with a higher gray level than the local threshold is labeled as object, otherwise it is labeled as background. Assume that $f_k^{\mathrm{\prime }}\left(x,y\right)$ refers to the binary image at the step $k$, we assign the object pixel a value of $f_k^{\mathrm{\prime }}\left(x,y\right)=1$, while the background pixel $f_k^{\mathrm{\prime }}\left(x,y\right)=0$.

Notice that the representation of the object and background differs by purposes. For character recognition, the object (i.e., characters) is often highlighted in black, while for landslide detection black color often represents the background.

2.2.4. Monte-Carlo iteration

In each Monte-Carlo step $k$, the original image $f_k\left(x,y\right)$ is transited to the binary image $f_k^{\mathrm{\prime }}\left(x,y\right)$ using square blocks with size $D$. However, as illustrated by the previous studies, binary result is quite sensitive to the block size, as well as the line-approximating interferences left between blocks that observed by Wen et al. (2013). As such, owing to that the block size is randomly changed during Monte-Carlo iteration, a better solution is to comprehensively consider the binary images of all Monte-Carlo steps. An indicator $P\left(x,y\right)$ is introduced, which is represented by (17) $$P\left(x,y\right)=\frac{1}{MCS}\sum _{k=1}^{MCS}{f}_k^{\mathrm{\prime }}\left(x,y\right)$$(17)

where $P\left(x,y\right)$ denotes the probability of the pixel $\left(x,y\right)$ detecting as object. The value of $P\left(x,y\right)$ theoretically ranges from 0.0 to 1.0. The higher value demonstrates that the pixel is more likely to be a component of object. The probability image is comprehensively evaluated by the following criterion (18) $$P^{\mathrm{\prime }}\left(x,y\right)=\left\{\begin{array}{l}\mathit{\text{object}}, P\left(x,y\right)\ge T_P\\ \mathit{\text{background}}, P\left(x,y\right)\hspace{0.17em}<\hspace{0.17em}{T}_P\end{array}\right.$$(18)

where $P\mathrm{\prime }\left(x,y\right)$ refers to the final binary image after the Monte-Carlo iteration; $T_P$ denotes probability threshold to evaluate the result in Monte-Carlo iteration. In order to find an appropriate probability threshold, trial-and-error adjusting can be considered, because it does not increase computational complexity. In this paper, a better fitting probability threshold $T_P=0.8$ is suggested. In order to better illustrate the improved local thresholding method, the procedure described above is represented in the flowchart (Figure 3).

Figure 3. Procedure of the proposed Monte-Carlo iteration algorithm.

2.2.5. Post-processing

The above procedure aims to binarize images with complex illumination condition. However, for the purpose of interpreting landslides using remote sensing images, background is much more complex. Owing to the high gray level of build-up areas, farmlands, and rivers in the remote sensing image, these regions are often falsely recognized as landslide candidate objects in the binary image. Currently, the interferences due to this region is difficult to deal with. Except of the changing detecting requiring images before and after the earthquake, intelligent classification based on machine learning technique, e.g., Li et al. (2017), may be a mean of solution. However, it is an ongoing work that requires a lot calibration.

In this paper, we incorporate the DSM with the binary image to minimize the interferences by build-up areas, farmlands, and rivers. Commonly, these areas are low slope gradient, low relief regions on plain. As such, we use the indicator of slope gradient $\theta \left(x,y\right)$ to exclude these areas, which is (19) $$P^{\mathrm{″}}\left(x,y\right)=\left\{\begin{array}{l}{P}^{\mathrm{\prime }}\left(x,y\right), \theta \left(x,y\right)\ge T_{\theta }\\ \mathit{\text{background}}, \theta \left(x,y\right)<T_{\theta }\end{array}\right.$$(19)

In Equation (19)(19) $$P^{\mathrm{″}}\left(x,y\right)=\left\{\begin{array}{l}{P}^{\mathrm{\prime }}\left(x,y\right), \theta \left(x,y\right)\ge T_{\theta }\\ \mathit{\text{background}}, \theta \left(x,y\right)<T_{\theta }\end{array}\right.$$(19) , $\theta \left(x,y\right)$ refers to the slope gradient that analyzed using DSM data. $P^{\text{″}}\left(x,y\right)$ denotes the final output binary image after post-processing.$\mathrm{ }{T}_{\theta }$ represents the gradient threshold to classify the plains and slopes. As indicated by previous studies (e.g., Wysocki and Zanner 2006; Jahn et al. 2006; Chabala et al. 2013; Mokarram and Hojati 2016), a threshold $T_{\theta }=5.0\mathrm{\%}$ is used in this paper, which means the areas with the slope gradient less than 5.0% are classified as flat, gently sloping plains. These areas are then excluded in the binary image as background (as shown in Figure 4). Notably, the value 5.0% in Equation (19)(19) $$P^{\mathrm{″}}\left(x,y\right)=\left\{\begin{array}{l}{P}^{\mathrm{\prime }}\left(x,y\right), \theta \left(x,y\right)\ge T_{\theta }\\ \mathit{\text{background}}, \theta \left(x,y\right)<T_{\theta }\end{array}\right.$$(19) is a rather conservative criterion because common wisdom holds that landslides usually have slope gradients greater than 20° (Dai and Lee 2002; Shiels and Walker 2013; Lee 2014).

Figure 4. In the post-processing, Digital Surface Model (DSM) is incorporated to exclude the plain and rivers.

3. Results

The quality of the binary results must be of priority concern. In this study, two types of application are conducted. The first application is benchmark test using some typical images in previous studies. It provides a clear illustration of the advantages of the proposed method to deal with the complex background and illumination. The second application is a case study, where we use the proposed method to binarize the remote sensing image of Beichuan area, China, after the 2008 Ms8.0 Wenchuan earthquake. The triggered landslides are interpreted by the proposed method, and results are compared to those of manual delineation.

3.1. Benchmark tests

Figure 5a displays an original document image from Pai et al. (2010). The image was taken from a paper that involves gradations of shadows. Figure 5b shows the binary image using Otsu’s thresholding method. Owing to the shadows in the upper part, the characters in that part could not be well separated, and then merged with the background. Although Niblack’s method performs better as shown in Figure 5d, the binary image is rather sensitive to the block size. We test a smaller block size, and the binary result in Figure 5c shows that the characters are still broken. We configure a 50-steps Monte-Carlo simulation in our proposed method. Figure 5e represents the probability image. Using a probability threshold $T_P=0.8$ generates the final binary image (Figure 5f). It is obvious that the proposed method yields a better outcome resembling those of Niblack’s, without any broken characters.

Figure 5. Benchmark test 1. (a) The imported image; (b) Otsu’s method; (c) Niblack’s method with block size of 15; (d) Niblack’s method with block size of 50; (e) Probability image by the proposed method; (f) Final binary image by the proposed method.

Figure 6a shows a captured photo of car license plate at night (Martín-Rodríguez 2013). High brightness of the two car bulbs significantly increases the overall gray level of the image, such that the license plate in the middle of the car could not be well separated by Otsu’s method (Figure 6b). Figure 6c and d show the binary results by Niblack’s method. Although the segmentation quality is improved, the results depend upon the block size, that a large block generates a better result in spite of increasing the computational complexity. In the proposed method, a 50-steps Monte-Carlo simulation is used to produce the probability image (Figure 6e), and a probability threshold $T_P=0.8$ is selected to obtain the final binary image (Figure 6f). It is obvious that the license plate can be well recognized.

Figure 6. Benchmark test 2. (a) The source image; (b) Otsu’s method; (c) Niblack’s method with block size of 15; (d) Niblack’s method with block size of 50; (e) Probability image by the proposed method; (f) Final binary image by the proposed method.

The above benchmark tests demonstrate that Otsu’ global thresholding method is limited to deal with complex illumination condition in the image. Niblack’s local thresholding method performs better, however, the binary result is rather sensitive to the block size, which is currently empirically determined. In this sense, the proposed method yields the best results among of the methods.

3.2. Case study

The proposed method is applied to interpret the landslides triggered by the 2008 Ms 8.0 Wenchuan earthquake in Beichuan region. The study area is located at Qushan Town, 138 km northeast of Wenchuan earthquake epicenter. Tongkou River flows across the area as shown in Figure 7a. The study area is close to the main active fault rupture zone and experienced high seismic shaking levels during the Wenchuan earthquake on May 12, 2008. A large number of landslides were induced in the region (Tan et al. 2011).

Figure 7. (a) Location of the study area. Image from Google Earth, taken at May 2001; (b) PAN SPOT5 image used for landslides detection. Image was taken at September 2008.

A PAN SPOT 5 image was available on September 24, 2008 (as shown in Figure 7b). The image has a 2.5 m resolution, covering 4.75 km² landscape region. The image had been pre-processed, including georeferencing, orthorectification, and atmospheric correction. The proposed method is applied to interpret landslides mixed in the image. A 50-steps Monte-Carlo simulation is configured to reduce the influence of complex illumination condition owing to the shadow areas. Figure 8a presents the generated probability map after Monte-Carlo simulation. It is shown that most of the discrepancies in the binary image of each step are concentrated in Tongkou River and the slopes at the right bottom of the image. To reduce the discrepancies, a probability threshold $T_P=0.8$ is used to obtain the output binary image. $T_P=0.8$ denotes that a pixel can be classified as a component of the landslide if it is in more than 80% steps during the Monte-Carlo simulation. Figure 8b shows the final output binary image after Monte-Carlo simulation.

Figure 8. Binarizing results using the proposed method. (a) Probability map generated by 50-steps Monte Carlo iteration; (b) binary image using a probability threshold $T_P=0.8$.

Nevertheless, the binary image in Figure 8b still contains some interferences. These interferences are commonly produced due to the high gray level of the rivers, roads, and built-up regions on plains. The interferences have similar gray-level features that close to the landslides, and they are difficult to exclude only by a PAN image. We use DSM data to identify these regions. The DSM data was converted from 1:50,000 contour lines in the topographic vector map provided by Chinese Geological Survey (Figure 9a). Based on the DSM data, slope gradient map is calculated using the TopoToolbox in MATLAB environment (Schwanghart and Scherler 2014), as shown in Figure 9b. However, as limited by the low resolution of the DSM data, water surface at river bend is illogically inclined. It has been manually modified in the following analysis. Notably, because the slope gradient map in Figure 9b and the binary image in Figure 8b have different resolution, both two images have to be processed by resampling to the same resolution. In order to reduce the interferences in Figure 8b, gentle slopes, plains, and rivers are excluded using the criterion in Equation (19)(19) $$P^{\mathrm{″}}\left(x,y\right)=\left\{\begin{array}{l}{P}^{\mathrm{\prime }}\left(x,y\right), \theta \left(x,y\right)\ge T_{\theta }\\ \mathit{\text{background}}, \theta \left(x,y\right)<T_{\theta }\end{array}\right.$$(19) .

Figure 9. Ten meter resolution DSM data of the study area for excluding plains and rivers. Data was provided by Chinese Geological Survey in 2001. (a) hillshade and contour map; (b) slope gradient map generated from DSM data, using TopoToolbox (Schwanghart and Scherler 2014).

Figure 10 presents the final binary image after 50-steps of Monte-Carlo simulation and post-processing by DSM data. White components highlight the landslides induced by the earthquake, and black parts refer to the background. The binary image in Figure 10 interprets totally 33 landslides in the study area. To illustrate the performance of the proposed method, we compare the binary result with that of manually visual delineation. Indicators proposed is adopted for the empirically-dependent visual delineation for landslides, including shape, size, color, texture, and tone. Totally 28 landslides are interpreted, as delineated by the red polygons in Figure 11.

Figure 10. Final detection of landslides based on the binary image and the DSM data.

Figure 11. Comparison of the detection result between the proposed method and manually visual delineation.

As revealed in Figure 11, all the 28 landslides by manual delineation could be identified in the final binary image. However, there remain 5 slopes that falsely classified as landslides. The false delineations are majorly the alluvial flats of the bend of Tongkou River, where the gray levels of those parts are close to the neighboring landslides. As well as the newly built road at the right part of the image. It can be explained in part that the 10 m resolution DSM data is limited to distinguish the terrain relief at the junction area between slopes and alluvial flats. The delineation results in Figure 10 could be improved if DSM data with higher resolution is available. In this context, LiDAR-based high-resolution DSMs, and DSMs derived by Airborne laser scanning (ALS), terrestrial laser scanning (TLS), and photogrammetric topography reproduction may provide a powerful supplement to the proposed method.

4. Discussions

4.1. Computational complexity analysis

However, the proposed method in this paper requires repeatedly binarization during Monte-Carlo simulation, which significantly increases the computational cost comparing to Otsu’s method. Following Pai et al. (2010), we compare the arithmetic computational complexity of the proposed method with previous methods. Here, the commonly used global and local thresholding methods, i.e., Otsu’s method and Niblack’s method, are considered owing to their lower computation.

For an image of size $M\times N$, and gray level $L$, the computational complexity of Otsu’s method and Niblack’s method have been deduced by Pai et al. (2010) as shown in Table 1. In Otsu’s method, total computational complexity includes the term $M\times N$ for calculating the image histogram, and $7L^2+5L-12$ for estimating the probability, mean and variance in order to dichotomize the histogram into two classes. While in Niblack’s method, the threshold is computed by calculating the local mean and standard deviation inside a fixed block with size $b\times b$ for each pixel, hence the complexity is reduced to $\left({4b}^2+3\right)\left(MN\right)$. In the proposed method, the total amount of the segmented block is $\left(M/D\right)\times \left(N/D\right)$. Each individual block is binarized using Otsu’s procedure, such that the complexity of each Monte-Carlo is $MN\left(7L^2+5L-12\right)/D^2$. This procedure is repeatedly executed in n steps of Monte-Carlo iteration; hence the total complexity of the proposed method is $n\left(MN+MN\left(7L^2+5L-12\right)/D^2\right)$.

Table 1. Computational complexity of the proposed and the previous thresholding method.

Method	Procedure	Computational complexity
Otsu’s method	Histogram	$MN$
	Separated class mean	$\left(3L-2\right)\left(L-1\right)$
	Separated variance	$\left(4L+3\right)\left(L-1\right)$
	Maximal separable degree	$11\left(L-1\right)$
	Total	$MN+7L^2+5L-12$
Niblack’s method	Mean of neighborhood	$\left(b^2\right)\left(MN\right)$
	Variance of neighborhood	$\left({3b}^2+1\right)\left(MN\right)$
	Denote threshold	$2\left(MN\right)$
	Total	$\left({4b}^2+3\right)\left(MN\right)$
The proposed method	Histogram	$MN$
	Otsu’s operation	$7L^2+5L-12$
	Local thresholding	$\frac{MN}{D^2}\left(7L^2+5L-12\right)$
	Total	$n\left(MN+\frac{MN}{D^2}\left(7L^2+5L-12\right)\right)$

Table 2 compares the results of the above three methods with respect to arithmetical computation. We used images with four different resolutions. For general images, the gray level $L$ is 256. In Niblack’s method, a typical block size $b=15$ is used, while in the proposed method, we assume that the block size is $D=N/2$, which means the image is partitioned into four square blocks and fragmented regions. Results in Table 2 demonstrates that general complexity of the proposed method is between that of Otsu’s and Niblack’s method.

Table 2. Comparison of computational complexity of the methods.

Method	Image resolution
	$640\times 480$	$1280\times 960$	$2560\times 1920$	$3072\times 2304$
Otsu	767	1,689	5,375	7,538
NiBlack	277,402	1,109,606	4,438,426	6,391,333
Proposed^#1	27,606	36,822	73,686	95,313
Proposed^#2	55,213	73,645	147,373	190,627
Proposed^#3	138,032	184,112	368,432	476,566

Notice: Proposed^#1, Proposed^#2, and Proposed^#3 refer to that the Monte-Carlo steps $n=10$, $n=20$, and $n=50$, respectively. $D=N/2$. The unit of the data is 1000.

However, as shown in Figure 12, for the proposed method, the transition of the computational complexity from “large block size” to “small block size” is drastic. For a general case, where $D=N/4$ and $D=N/10$, the complexity significantly increases in a 50-steps Monte-Carlo iteration. Figure 12 also reveals that the proposed method performs better in a high-resolution image, because the computational complexity is not quite sensitive to the image resolution. For the worst case, when $D=N/10$, the computational cost of the proposed method is still only 53.5% that of Niblack’s method. It implies that the proposed method is a better solution for binarizing large-scale, high-resolution images, such as regional remote sensing images.

Figure 12. The tendency of computation complexity upon various image resolutions among the proposed method and Niblack’s method ($n=50$).

4.2. Limitations of the current method

4.2.1. Dealing with the clumping of individual landslides

In this paper, we present an alternative approach for rapid delineation of landslide perimeter. However, as indicated by Marc and Hovius (2015), amalgamation, i.e., the mapping of several adjacent landslides as a single object, often exists and leads to potentially severe distortion of the statistics of landslides inventories. It is one major concern of the automatic landslides mapping, because the overestimated area of the landslide perimeter results an illogical volume of the landslide according to the power-law scaling relationship between area and volume of a landslide (Larsen et al., 2010). In this context, the presented approach in this paper is limited, because it only identifies the pixels affected by landslides, rather than delineating individual landslide objects.

4.2.2. Empirical adjusting on the probability threshold $T_P$

Except of the computational complexity, the proposed method also has some limitations regarding some parameters in the algorithm. Although the proposed method addresses the essential issues with respect to the block size and the so-called boundary effect, it still involves some parameters that need a further discussion.

The first parameter that needs a further discussion is the total steps of Monte-Carlo runs $MCS$. Ideally, the greater $MCS$ generates better data set for the probability analysis. We have performed sensitivity analysis on the generated probability map $P(x,y)$ against different $MCS$ (as shown in Figure 13, benchmark test 1 for example). It demonstrated that $MCS$ has no obvious influence on the overall distribution of the generated probability map $P(x,y)$. However, as implied in Figure 13(a), insufficiency of the total steps of Monte-Carlo runs commonly results in a rough distribution of $P(x,y)$, limiting to overcome the boundary effect. Considering the fact that greater $MCS$ linearly increases the amount of computation, generally we suggest a total number of Monte-Carlo steps $MCS=25\sim 100$ for practical work.

Figure 13. The variation of binarization $P\left(x,y\right)$ against different values of $MCS$.

Empirical adjusting also exists while determining the thresholds. Equations (18)(18) $$P^{\mathrm{\prime }}\left(x,y\right)=\left\{\begin{array}{l}\mathit{\text{object}}, P\left(x,y\right)\ge T_P\\ \mathit{\text{background}}, P\left(x,y\right)\hspace{0.17em}<\hspace{0.17em}{T}_P\end{array}\right.$$(18) and (19)(19) $$P^{\mathrm{″}}\left(x,y\right)=\left\{\begin{array}{l}{P}^{\mathrm{\prime }}\left(x,y\right), \theta \left(x,y\right)\ge T_{\theta }\\ \mathit{\text{background}}, \theta \left(x,y\right)<T_{\theta }\end{array}\right.$$(19) involve two thresholds, the probability threshold $T_P$ and the gradient threshold $T_{\theta }$. As illustrated in Section 2.2, the value $T_{\theta }=5.0\mathrm{\%}$ is derived according to the previous geomorphological studies. With regard to the probability threshold $T_P$, it has a significant influence of the binarization results as shown in Figure 14. It is currently empirically manipulated and needs trial-and-error adjusting when generating the output binary image. $T_P$ only influences the binary image, having no relevance to compute the probability map during the Monte-Carlo simulation. A greater probability threshold $T_P$ is beneficial to decrease the presence of unwanted data, but on the contrary, lost some fine details. Unfortunately, the probability threshold $T_P$ currently needs empirically manipulation and adjusting for a better output binary image. We use a scrollbar in the MATLAB-based program to check in time the output binary image as a function of $T_P$, because adjusting on $T_P$ does not increase computational complexity. However, the uncertainty of $T_P$ somewhat limits the proposed method to a robust automatic one. As such, our ongoing work focuses on improving the here presented method by addressing this issue.

Figure 14. The variation of the probability map results against different values of $T_P$.

5. Conclusion

In this paper, we present an improved image binarization method for noncontact detection of earthquake-induced landslides. In general, the presented method first partitions the image into a number of square blocks with uniform size that reflects local characteristics in the entire image. Following image partitioning, the binary image of each block can be obtained by using Otsu’s thresholding method. To reduce the influence of the block size and interferences due to the so-called boundary effect, we incorporate the Monte-Carlo simulation into the binarization algorithm. Varying block sizes during the Monte-Carlo iteration generate a probability map. A final output binary image can be obtained with a pre-defined probability threshold. The presented method integrates the advantages of global and local methods. Benchmark tests show the effectiveness of the presented method in performing well to binarize the image with complex illumination condition.

To improve the binarization result for separating landslides from the background, we furthermore incorporate the digital surface model into the binary image, such that the interferences by build-up areas, farmlands, and rivers can be significantly reduced using the indicator of slope gradient. As a case study, we use the proposed method to detect landslides in the remote sensing image near Beichuan area, China, after the Ms8.0 Wenchuan earthquake. Results demonstrate that the presented method interprets landslides resembling those of manually visual delineation.

However, the proposed method is somewhat limited by introducing two extract thresholds during post-processing. Although criterions for determining the two thresholds are available, it still needs empirical adjusting. Efforts addressing such limitation will be the next step toward making the presented method a robust and automatic solution for rapid landslide detection.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was financially supported by the National Key R&D Program of China (Grant No. 2017YFB1201204, Z. Han); the National Natural Science Foundation of China (Grant No. 41702310, Z. Han); the Foundation of State Key Laboratory of Geo-Hazard Prevention and Geo-Environment Protection (Grant No. SKLGP2017K014, Z. Han); and the Natural Science Foundation of Hunan Province, China (Grant No. 2018JJ3644, Z. Han). These financial supports are gratefully acknowledged. The authors also extend their gratitude to editor-in-chief Ramesh Singh, and two nominated reviewers for their insightful comments.