Laplacian pyramid-based change detection in multitemporal SAR images

ABSTRACT

The Synthetic Aperture Radar (SAR) has many attractive and desired features to cope with uncertain weather conditions and also provides superlative resolution in images. In this paper, a unique multiscale approach of change detection (CD) that integrates the preprocessing and CD technique of SAR imagery is discussed. Condensed anisotropic diffusion method is employed as a preprocessing tool for despeckling, and the resultant image is processed by multiscale approach in the Laplacian pyramid (LP) domain. In each layer of LP, the difference between the pyramid coefficient of two images of the same area obtained at different periods is estimated through log ratio operation. Each layer of pyramid coefficient has explicit features at particular scale. Fine scale related to a band pass layer localizes edges of the change area and the coarse scale related to a low pass layer preserves more detail information of the change region. Reconstruction of LP provides exact difference image with perfect boundary. Otsu's algorithm is used to generate binary change map to distinguish changed and unchanged pixels. The efficiency of the proposed method is tested by using simulated and real SAR image datasets. The confusion matrix parameters are used to prove competence of the proposed CD method.

KEYWORDS:

• SAR image
• change detection
• condensed anisotropic diffusion
• Laplacian pyramid
• Otsu’s algorithm
• performance analysis

Introduction

In remote sensing, satellite images are acquired by both active and passive sensors. Synthetic Aperture Radar (SAR) comes under active microwave sensors which uses a broad range of radiometric spectrum. It can work in all weather conditions, even at day and night, which enables the SAR to extract correct information about the targets. Since the microwave sensors are self-reliant to atmosphere and a cloud cover condition, it encourages the researcher to use SAR images in various fields. The remote sensing satellites such as ENVI-ASAR, TerraSAR-X, ALOS and Radarsat-2 effectively replicate the active changes about spatial and temporal information of earth’s surface. Especially ENVI-ASAR provides constant data sets for geological monitoring (soil-moisture, agriculture, coastal (Ghosh, Kumar, & Roy, 2015), forest (John & Woodcock, 1996; Kennedy, Cohen, & Schroeder, 2007) and ship tracking) and mapping applications (topographic (Shalaby & Tateishi, 2007), ice and snow applications). The above ENVI-ASAR has broad angle coverage with excellent resolutions during every 35-days revisit cycle using ScanSAR technique. The ASAR high-resolution images have 30 m of precision and works from 15° to 45° of incident angle. The ASAR wave mode creates 5 km x 5 km sketch spaced with 100 km of azimuth resolution (VV-polarization), whereas the ASAR swath mode provides 400 km of swath coverage with 150 m resolution. The Radarsat satellite is most leading surveillance system developed by Canada centre for mapping, and earth observation is used to track atmospheric changes of flood (Hazarika, Apurba, & Borah, 2015), oil spill, earthquake, land use (Hadjimitsis et al., 2010) and change map (Ban & Yousif, 2012) applications. The Radarsat also provides high-quality multitemporal images with a steerable beam and single frequency to offer fine image scale and resolution. This satellite operates at 5.3 GHz/C-band, 5.6 cm wavelength with HH-polarization. The Radarsat image mode has a swath width of 100 km, incident angle span of 20– 49° with 30 m resolution.

Radar sensors play essential role in SAR imagery, particularly the realization of change map applications. Change detection (CD) is the measure of identifying the differences of same geological area observed at different times (Radke, Andra, Al-kofahi, & Roysam, 2005; Singh, 1989; Baselice et al., 2014). Generally, CD methods are categorized into two: unsupervised and supervised. Unsupervised method is an analysis of change and no-change classes without prior knowledge about the target image, whereas supervised method required former information. CD process consists of three stages (Dogan & Perissin, 2014; Mura, Benediktsson, Bovolo, & Bruzzone, 2008): (1) preprocessing of multitemporal SAR images; (2) image comparison to identify the changes; (3) analysis of CD by appropriate thresholding or clustering methods (Celik, T. 2009b).

In preprocessing stage, SAR images undergo speckle reduction, geometric correction, digital elevation model correction, co-registration, etc., and these operations are preferred according to the application area. Out of all the preprocessing steps, the speckle reduction is essential for CD application. SAR images are very sensitive to the coherent nature of multiplicative noise called speckle (Gong, Cao & Wu 2012). Speckle noise is a granular noise which degrades the quality of SAR image consequently. Various adaptive filtering methods are available to reduce speckle noise like Lee (Lee, 1980), Frost (Frost, Stiles, Shanmugan, & Holtzman, 1982) and Kuan (Kuan et al., 1985, Kuan, Sawchuk, Strand, & Chavel, 1987) methods; they do speckle reduction by simply converting the multiplicative noise model into additive model using a logarithmic operator. In addition, the sensitivity of the filter window size causes smear edges and blur spatial information of an image. Recently, diffusion method shows improved merit compared with adaptive filtering techniques. Diffusion filter removes speckle from an image by modifying the image by means of a partial differential equation (PDE) (you & Kaveh, 2000; Feddern et al., 2006; Krissian et al., 2007). Condensed Anisotropic Diffusion (CAD) (Kalaivani et al., 2012) method is used, which incorporates energy condenser into nonlinear coherent diffusion. It can reduce the speckle content with preservation of edges and fine spatial details (Masoomi, Hamzehyan, & Shirazi, 2012).

Then, speckle-free multitemporal SAR images are to be interpreted to identify the difference image (DI) (Bruzzone & Prieto, 2000) in the subsequence stage. The DI is achieved by comparing two co-registered images by various mathematical functions. In the literature, many comparison methods (Dekker, 1998; Lu, Mausel, Brondizio, & Moran, 2004; Singh, 1989; yousif & Ban, 2014) are available to find the DI such as algebra-based approach (Masroor, Dongmei, Angela, Hui, & David, 2013), transformation-based approach (Li & Yeh, 1998), classification-based approach (Jesus, Arie, & Joost, 2012) and similarity measures (Chesnokova & Erten, 2013). Out of these methods, algebra-based approach is a simple method to find DI which utilizes image differencing (Ma et al., 2012), image rationing, image regression, image correlation and change vector analysis (Sartajvir & Rajneesh, 2014). Existing methods are based on the analogy of the local neighborhood window to estimate the changes. Generally irrelevant window size would cause wrong and unwanted detection. In case of using small window size for correlation analysis, no detection will be carried on the homogeneous regions; on the other hand, larger window size will detect changes only in larger area or high-intensity regions, such cases produce a coarse resolution DI (Inglada & Mercier, 2007). The image differencing method and image ratio methods augment the foreground pixels by confining the background pixel information (Gong, Zhou & Jingjing 2012). The ratio operations have the following limitations: the mean ratio operator may fail to detect changes because the DI is emerged from the conversion of the local mean value of the input images. Image ratio operator and its distribution depend on relative change of average intensity between a pair of multitemporal SAR images. The logarithmic characteristic of log ratio (LR) method strengthens the low-intensity pixels, while deteriorating the pixels in the area of high intensity. Hence, it may not give the real changed trends to its maximum level. To overcome the above issues, multiscale approach is introduced in the CD analysis.

In image processing, the multiscale approach is a significant concept for retaining spatial information. Multiscale images are computed by orthogonal and non-orthogonal basis. Laplacian Pyramid (LP) is well-known non-orthogonal method. The representation of multiscale image is in decreased order of the image resolution and the sample density, obtained through regular steps. The Gaussian pyramid (GP) consists of a low pass version of the image. On the wavelet pyramid, decomposition is based on the outcome of the orthogonal method. Also, wavelet pyramid transforms are more computationally complex, and the decomposition is to be done recursively on the low pass band to obtain other levels.

Techniques in the literature,Bovolo and Bruzzone (2005) developed CD technique by LR operation and multiresolution analysis done by stationary wavelet transform (SWT) to decompose multiple resolutions of the difference coefficient. Adaptive scale-driven approach is used to fuse the final change map. This algorithm provides considerable CD result with accuracy, but the decision threshold was elected based on manual selection of trial and error procedure. Celik (2009a) proposed undecimated discrete wavelet transform to decompose the DI. Each pixel of feature vector is sampled with multiscale difference data. The binary change map is derived by k-means clustering algorithm. Though variety of sensor images acquired in multiple geometric conditions, guaranteed results are not obtained for multiple selection of sensor data. Celik (2010b) adapted unsupervised CD by Bayesian interference with expected maximization algorithm. Also, dual tree-complex wavelet transform exploits multiscale approach for DI. However, this multiscale approach has not provided efficient number of scales to prefer, and also speckle filter failed to preserve the boundaries of the changed region exactly. Wang et al. (2013) used multiclass segmentation method by triplet markov field model. This method obtained a better solution for CD with low miss rate and high precision at the cost of long computational time. yousif and Ban (2013) developed non-local means (NLM (Coupe et al., 2009)) algorithm for speckle noise reduction. Principal component analysis (PCA) is used to reduce the neighborhood dimensionality of the feature vector. This NLM-PCA model takes high computational complexity. Hu and Ban (2014) proposed that the satellite images are compared by ratio method and then classified by thresholding algorithm. If the DI is bimodal histogram, the Kittler Illingworth algorithm is used to classify the change map. Otherwise, multiscale region selection is executed for unimodal pixel classification. But this method has not provided the dependency information associated with parallel neighborhood pixels to increase the accuracy of change map. Ajadi et al. (2016) carried CD in satellite images by multiscale-driven approach. Fast non-local filter is applied for speckle noise reduction, and SWT is employed to decompose the DI. Bayesian threshold is estimated in each decomposition layer. The localized features are fused by measurement-level fusion to obtain the change map, though this algorithm shows blur change result with false detection.

In order to improve upon previous work, pyramid-based multiscale approach is proposed. The pyramid-based multiscale approach has a great deal of significance in the representation that contains spatial localization in the spatial-frequency domain. The pyramid transformation may separate critical modules of the image pattern, hence they are more feasible to access and also more compact for data storage and transfer. Image pyramid is a data architecture designed to aid dynamic scaled convolution over the reduced sequence of the original image. It consists of recursive copies of the original image which have self-similarity form decreased in both sample density and resolution. In GP, the original image is convoluted with low pass filter ensueing with sub sample by a factor of two to earn the later pyramid levels. Most of the image processing applications require medium band resolution (band pass and low pass images), which demands LP (Adelson, Anderson, Bergen, Burt, & Ogden, 1984) structure. This may be achieved by differencing each GP level from immediate downstream level. Due to this, interpolation has to take place before balancing the sample density. Band pass copies of the pyramid layers provide valid details and related components residing in each decomposing level. LP have been used to analyze images at various image processing applications such as image compression (Adelson et al., 1984; Perry, 1997), texture analysis and synthesis, image mosaics and image filtering techniques. In image compression technique by LP, less significance of compression ratio with loss of information is provided (Michael, 1992). In texture analysis technique, the pyramid structure could not recognize complete perceptual information of natural textures (Heeger & Bergen, 1995). Pyramid domain on image mosaics has not supported an image with movable object and different exposures (Burt & Adelson, 1983a). LP-based noise filtering technique enables only with fine scale layers and does not comfort to use with course scale (Zhang, Yoo, Koh, & Kim, 2007). However, pyramid structure regards poor choice of selection of above-existing applications. In our paper, we exploit the strength of the LP domain for CD application to overcome the existing limitation of LR operation. Moreover, the band pass copies obtained in LP enhance the image structures, linear features and edges which are essential for image interpretation. The essential property of LP is that it gives an entire image portrayal: the degenerative stage recoups the original image absolutely, without any lose of information.

Finally, the decision theory is utilized to make perfect automatic classification of the DI. The realistic CD approach is enhanced by thresholding algorithm to provide an automatic binary change map. Threshold is a simple and effectual technique to segregate objects (changed) and the background (unchanged) pixels (Moser & Serpico, 2006; Rosin & Ioannidis, 2003; Gong et al., 2014) of an image. Most of the thresholding algorithms are categorized according to changed and unchanged classes in the assumptions of the specific mode (unimodal or bimodal) of the histogram of an image. In change map application, unidentified tiny portion may lead to enormous detriment of environment changes, particularly in terrain control changes such as land use, land cover, agricultural monitoring, crop area observation, deforestation, urban expansion in the agricultural area and flood risk assessment. On the above consideration, several threshold algorithms are referred to classify the DI using Kittler-Illingworth (KI) algorithm (Kittler & Illingworth, 1986), Otsu algorithm (Otsu, 1979; Yang, Shen, Long, & Chen, 2012), Kapur algorithm (Kapur, Sahoo, & Wong, 1985), Bayesian thresholding (Bovolo, Camps-Valls, & Bruzzone, 2010), Tsai method (Tsai, 1995), Rosin algorithm (Rosin, 2001) and Expectation-Maximization algorithm (Hao, Shi, Zhang, & Li., 2014). In this paper, Otsu’s thresholding algorithm is utilized for classification. The optimal threshold is chosen by either minimum value of within-class variance or maximum value of between-class variance of the histogram of an image. Otsu algorithm is preferred for the following merits. It has an unsupervised quality of threshold selection, quite simple and applies zeroth- and first-order cumulative distribution to elect gray-level thresholding. An optimal threshold is analyzed through the probability of class mean levels, which will offer automatic and constant threshold.

In order to consider the above standard processing stages, we proposed speckle-free, unsupervised CD for SAR images by following aspects: (1) preprocessing by CAD; (2) Image difference by LP-based CD; (3) Otsu thresholding algorithm for effective classification of change map. The rest of the paper is structured as follows. The Background section presents basics of the preprocessing technique. Methodology section elaborates the details of the proposed method. Result and discussion section describes experimental significance and related discussion of the proposed technique. Outcome of proposed method with future work is given in the conclusion part.

Background

Diffusion in speckle reduction

The fundamental speckle reduction technique is developed from diffusion process. Diffusion is termed as a physical process (Rajan & Kaimal, 2006) that equilibrates two unequal concentrations without creating or destroying mass. The physical process of diffusion is expressed by Flick’s law as,

(1) $j=-D.\left(\mathrm{\nabla }I\right)$(1)

The homogeneity between the flux $j$ and gradient $\mathrm{\nabla }I$ is represented by the diffusion tensor $D$, which is a positive scalar symmetric matrix. Diffusion related to the image processing is the concentration of gray-level value in a particular location. Perona & Malik (1990) proposed anisotropic diffusion using nonlinear PDE. Nonlinear diffusion uses scalar valued diffusivity instead of the diffusion tensor.

(2) $\left\{\begin{array}{c}\frac{\mathrm{\partial }I}{\mathrm{\partial }t}=div\left[D\left(∥\mathrm{\nabla }I∥.\mathrm{\nabla }I\right)\right]\\ I_{\left(t=0\right)}=I_0\end{array}\right.$(2)

where $div$ is the divergence operator, $∥\mathrm{\nabla }I∥$ is the gradient, $D∥\mathrm{\nabla }I∥$ is diffusion function, $D\left(x\right)$ is scalar function and $I_0$ is the initial image. This diffusivity provides maximum likelihood value. Generally, the second-order PDE-based approach is subjected to feature blocky effects (Kalaivani et al., 2012b). To overcome the above issue, fourth-order PDE-based speckle reduction method (you & Kaveh, 2000), which gives improved result for speckle noise removal and reduces edge sensitivity formation, is used. However, fourth-order PDE method takes higher number of iterations to converge.

Methodology

This paper aims a unique multiscale approach for SAR imagery in an unsupervised CD strategy, which integrates the preprocessing and CD. The proposed method is demonstrated in sections below.

Preprocessing

A CAD method is suggested as a preprocessing tool to eliminate speckle noise in SAR images. CAD method (Kalaivani et al., 2012a) consists of two terms: a diffusion term and a condenser. The term of diffusion coefficient is weakening the speckle by diffusion tensor. The second term of energy condenser enhances the performance of diffusion tensor and enables to correct the feature broadening effect. Blending energy condenser into diffusion tensor leads to preserve the edges of change area, enhance the local details of lines and bar features, and also rectifies the feature broadening defect. In this method, the scalar function $D\left(x\right)$ acts as a tensor and controls diffusion process in both gradient and contour direction. The CAD is stated by,

(3) $\frac{\mathrm{\partial }I\left(x,y,t\right)}{\mathrm{\partial }t}=\mathrm{\nabla }\left(D.\mathrm{\nabla }I\right)+\beta \gamma (I_c{)}^{\gamma }{I}^{\gamma -1}\left(\vec{x},t\right)$(3)

The anisotropic diffusion term is expressed in discrete form as,

(4) $\frac{\mathrm{\partial }I\left(x,y,t\right)}{\mathrm{\partial }t}=\left[\begin{array}{c}\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}{I}_x\\ I_y\end{array}\right)\\ diffusivitymatrixGradientvector\end{array}\right]$(4)

With reference to semi-explicit form,

(5) $\frac{I_{i,j}^{n+1}-I_{i,j}^n}{\tau }=D_N\left({\mathrm{\nabla }}_N{I}_{i,j}^n\right)+D_S\left({\mathrm{\nabla }}_S{I}_{i,j}^n\right)+D_E\left({\mathrm{\nabla }}_E{I}_{i,j}^n\right)+D_W\left({\mathrm{\nabla }}_W{I}_{i,j}^n\right)+\beta \gamma {\left(I_{ci,j}^n\right)}^{\gamma }{\left(I_{i,j}^n\right)}^{\gamma -1}$(5)

where $D_{N,}{D}_S,D_E$ and $D_W$ are diffusion coefficients summarized as given,

(6) $\begin{array}{c}\begin{array}{c}{\mathrm{\nabla }}_N{I}_{ij}^n=I_{i-1,j}^n-I_{i,j,}^n\\ {\mathrm{\nabla }}_S{I}_{ij}^n=I_{i+1,j}^n-I_{i,j,}^n\end{array}\\ {\mathrm{\nabla }}_W{I}_{ij}^n=I_{i,j-1}^n-I_{i,j,}^n\\ {\mathrm{\nabla }}_E{I}_{ij}^n=I_{i,j+1}^n-I_{i,j.}^n\end{array}$(6)

where $I_{ij}^n$ is a value of ${\left(i,j\right)}^{th}$ pixel at time “ $n$ ” and $I_{ij}^{n+1}$ is pixel value at time $n+1.I_{ci,j}^n$ is an average of $I_{i,j}^n.$ The $\tau$ is the time interval between two consecutive iterations. The time interval is set at 0.25 to maintain $0\le \gamma \le \frac{h^2}{4},$ and these values are chosen by $\gamma \ll 1$ and neighborhood pixel distance is $h=1$ for the image. In case of weight factor $\beta =0,$ the image region is subjected to nonlinear coherent diffusion. When the value $\beta \ne 0$, it prevents the coarsening of high-frequency structure in predominant image regions.

(7) $I_{ij}^{n+1}=I_{ij}^n+D_N\left({\mathrm{\nabla }}_N{I}_{i,j}^n\right)+D_S\left({\mathrm{\nabla }}_S{I}_{i,j}^n\right)+D_E\left({\mathrm{\nabla }}_E{I}_{i,j}^n\right)+D_W\left({\mathrm{\nabla }}_W{I}_{i,j}^n\right)+\beta \gamma {\left(I_{ci,j}^n\right)}^{\gamma }{\left(I_{i,j}^n\right)}^{\gamma -1}$(7)

The diffusion is an iterative process. The stopping criterion is decided by mean absolute error (MAE) between two successive steps (Yu & Acton, 2002).

(8) $MAE\left(I\right)=\frac{1}{mn}{\sum }_{i,j=1}^{m,n}\sqrt{{\left(I\left(i,j,t\right)-I\left(i,j,t-1\right)\right)}^2}$(8)

where $\left(I\left(i,j,t\right),I\left(i,j,t-1\right)\right)$ are filtered images at times $t$ and $t-1$ and $m,n$ are columns and rows of the image, respectively. MAE value decides finest diffusion process to earn desirable despeckling. For the implementation of CAD as a preprocessing tool, the parameters are chosen as $\alpha =1,s=70,\beta =0.05,\gamma =0.75$ and $I_c=⟨I⟩$. In each iteration, the energy of $I$ is updated by a factor of $I_0/⟨I⟩$, where $⟨\cdot ⟩$ is the mean value. Hence, the proposed image has similar energy like input $I_0$. This method provides impressive speckle noise removal with condensed feature broadening effect.

Laplacian pyramid

A basic structure of image pyramid is a hierarchical representation of multiscale images. The pyramid structure consists of decomposition and a reconstruction stage, which are designed by approximation and interpolation filter (Burt & Adelson, 1983b; Yu, 2009). In the decomposition stage, the signal of the input image is consecutively decomposed into reduced copies of approximation and residual signal information. The residual signal information is estimated by the difference between the signal on the fine scale and the interpolated signal from the coarser scale. A fine scale with respect to lower pyramid layer has its size equal to the input image. Pyramidal structure is designed by its specific decimation factor, approximation and interpolation filter. The pyramid structure decomposition is done by two techniques: GP and LP. In GP, the input image $G_0$ is continually filtered and sub-sampled to achieve the sequence of low-resolution images $G_1,G_2\dots ,G_{d-1}.,$ where $d$ is number of decomposition layers. This consists of a set of low pass filter layers of the input image whose bandwidth reduces by one octave steps.

(9) $\begin{array}{rl}& G_0=I\\ & G_l=Reduce\left[G_{l-1}\right]\\ & L_l=G_l-Expand\left[G_{l+1}\right]\end{array}$(9)

LP constructed by two standard operators: Reduce and Expand The Reduce operator in LP (Michael, 1992) is realised by 2D low pass filter followed by sub sampling by a factor of 2 in both directions. The Expand operator is performed by upsample with low pass filter followed by a factor of 4; the multiplication factor is essential to retain the average intensity being decreased by the insertion of zeros. This operator increases the image to double in size on both directions. For an original image I, LP at layer l is $L_l$, where $l=0,1,2,3\dots d-1$ and $d$ is number of decomposition layers. The decomposition of GP and LP structures are shown in Figure 1. The GP consists of different size of low pass layers of the original image and the LP decomposes the original image into the layers of band pass images and a unique low pass image. Reconstruction stage is done by simply regenerating the image from LP; it can be realized by reversing the decomposition steps.

Figure 1. Flowchart of Gaussian (G_l) and Laplacian pyramid (L_l)

Log ratio (LR) operation

Consider two despeckled SAR images acquired in the same geographical area at different time sequences, $I_1,I_2$, respectively. These input images are decomposed by LP transform and are represented as $L_1^{i,j}$ and $L_2^{i,j}$, respectively.

(10) $\begin{array}{c}{Y}_1=\left\{L_1^{i,j}\right\}wherei=1\dots m.\\ Y_2=\left\{L_2^{i,j}\right\}j=1\dots n.\end{array}$(10)

where $m,n$ represents the size of the each decomposition layer. Ratio operators are available to predict the differences between two images (Ashok, Varshney, & Arora, 2007). In that, most frequently used method is a logarithmic ratio method (Bazi, Bruzzone, & Melgani, 2005). Change differences can be identified by pixelwise comparison of ratio operator. Log operator is used to constrict the range of differences of the ratio image and to stabilize the values of the pixel either lesser or greater than 1. The LR operator is defined below and is expressed in a natural-logarithmic scale.

(11) $LR=log\left(\frac{Y_2}{Y_1}\right)=log\left(Y_2\right)-log\left(Y_1\right)$(11)

Significance of the proposed method

The aim of the LP is to decompose the multitemporal SAR images into different resolution layers in order to improve essential signal components like edges, homogeneous region and geometrical detail information of the change map image. LP decomposition of SAR images is embed with LR operator and its reconstructed image is the DI. The flow graph of the Laplacian Pyramid based change detection (LPCD) process is depicted in Figure 2. The $\mathfrak{D}{\stackrel{\sim }{\mathrm{I}}}_{decomp}$ is LP decomposion of the DIs that is expressed below,

(12) $\begin{array}{rl}& \mathfrak{D}\tilde{\mathrm{I}}{ }_{decomp}=\left[log\left(L_{I_{2}l}\right)-log\left(L_{I_{1}l}\right)\right];\\ & \begin{array}{c}{I}_1,I_2-\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{s}\\ l=0,1,2,3-\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}\mathrm{s}\end{array}\end{array}$(12)

Figure 2. Block diagram of LP-based change detection

The prime objective (i.e. CD in this study) and the essential signal components (i.e. edges and change pixel information) of the DI are typically lying in different band layers as a result of different spatial components. Mainly, spatial resolution of the DI provides change information from coarse to fine scale. Coarse to fine search is a process of consecutively applying Reduce and Expand operators to the given image. They repetitively zoom in on the accurate position of an area depicted at the finest scale. The edges of the objects are enhanced by fine scale layer and detailed information is predicted in coarse scale layer. Thus, the estimation of spatial difference in medium band resolution (band pass and low pass) can efficiently detect changes without losing primary information. We have preferred four decomposition $L_l:\left(l=0,1,2,3\right)$ layers of the pyramid structure, each make explicit features at specific scale. The significance of different layers is given below;

First layer $\begin{array}{c}{L}_l\\ l=0\end{array}$ : This is fine scale related to a low pyramid layer and has a size equal to the input image. The Reduce and Expand operators are used in original image to obtain band pass layer of the image. LR operator is utilized in high-frequency component layer to yield exact edges, lines and bar features of the image.

Second layer $\begin{array}{c}{L}_l\\ l=1\end{array}$ : This is second band pass layer of the DI. The operators are applied in $\left(l=0\right)$ layer and transition takes place from finest scale to next layer. Generating DI in this medium band resolution restrains unchanged pixels and enhances changed pixels exactly.

Third layer $\begin{array}{c}{L}_l\\ l=2\end{array}$ : This is third band pass layer of the DI. The Reduce and Expand operators are applied in $\left(l=1\right)$ layer to maintain symmetry among the distribution of two (changed and unchanged) classes of the DI.

Fourth layer $\begin{array}{c}{L}_l\\ l=3\end{array}$ : This is coarse scale related to a high pyramid layer. Only Reduce operator is used in $\left(l=2\right)$ layer to get low pass image. It can verify the hypothesized features tentatively identified at coarse scale. More detail information presents in this layer, they computationally efficient by means for applying valuable matching operator only in area of interest. .

Hence, only four layers $l=0,1,2,3$ of pyramidal structures are considered in the proposed method. The amount of information obtained from the individual layers is significantly increased when LR operator is combined with each layer of pyramid decomposition. Finally, the DI is generated by reconstructing the LP through degeneration steps. The flattening of background information is a common problem in LR method, which is significantly reduced by combining with LP.

The proposed LP-based DI is $\mathfrak{D}\tilde{\mathrm{I}}$, which is expressed below,

(13) $\mathfrak{D}\stackrel{\sim }{\mathrm{I}}=reconst(\mathfrak{D}{\stackrel{\sim }{\mathrm{I}}}_{decomp})$(13)

Optimal thresholding algorithm

The CD is to detect the meaningful changes, thereby increasing the CD accuracy and reducing error rate (ER). Thresholding is one among the finest ways to deal with. This threshold algorithm is proposed by Otsu (Deng-Yuan, Ta-Wei, & Wu-Chih, 2011; Otsu, 1979) and is a proven method for selecting threshold for gray-level distribution of the image. The normalized gray-level histogram of the DI from $\mathfrak{D}\tilde{\mathrm{I}}$ ($L$ is gray levels, $n_i$ is number of pixels at level $i$ and total number of pixels $N$) in the form of probability distribution is given by,

(14) $P_i=\frac{n_i}{N}{P}_i\ge 0;\sum _{i=1}^L{P}_i=1$(14)

The total class mean of a gray-level image is,

(15) $M_T={\sum }_{i=1}^{L}i{P}_i$(15)

The effectiveness of optimal threshold is followed by discriminant criterion measures either from the minimum value of the weighted within-class variance or maximum value of the between-class variance of the image. The optimized threshold by maximum between-class variance ${\sigma }_B^2$ can be represented as,

(16) ${\sigma }_B^2\left(k\right)=\frac{{\left[M_{T}S\left(k\right)-M\left(k\right)\right]}^2}{S\left(k\right)\left[1-S\left(k\right)\right]}$(16)

$S\left(k\right)$ is zeroth-order cumulative distribution of the gray-level histogram assigned to object (background), $M\left(k\right)$ is first-order cumulative distribution of the gray-level histogram up to level $k\left(k\in \left[0,L\right]\right)$.

(17) ${\sigma }_B^2\left(t^{\ast }\right)={}_{0\le k\le L}^{max}{\sigma }_B^2\left(k\right)$(17)

Table

Algorithm 1 Laplacian pyramid-based change detection Algorithm.
Step 1	:	Use multitemporal SAR images $I_1$ and $I_2$ as input images.
Step 2	:	Preprocess the images by condensed anisotropic diffusion (CAD) method.
		Estimate the absolute gradient in $x,y$ direction $\left|\mathrm{\nabla }I\right|=\sqrt{I_x^2+I_y^2}$. Evaluate D, ${\lambda }_1$ and ${\lambda }_{2.}$ $D=\frac{1}{{\left|\mathrm{\nabla }I\right|}^2}{\left(\begin{array}{cc}{I}_x& I_y\\ -I_y& I_x\end{array}\right)}^T\left(\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right)\left(\begin{array}{cc}{I}_x& I_y\\ -I_y& I_x\end{array}\right)$ ${\lambda }_1=\left\{\begin{array}{ccc}\alpha \left(1-{\left(\frac{A}{S}\right)}^2\right)& ,& ifA=\left({\mu }_1-{\mu }_2\right)\le S^2\\ 0& ,& else\end{array}\right.;$ ${\lambda }_2=\alpha$ Calculate $I_{ij}^{n+1}$ from $I_{ij}^n$ using semi implicit scheme Equation (5)(5) $\frac{I_{i,j}^{n+1}-I_{i,j}^n}{\tau }=D_N\left({\mathrm{\nabla }}_N{I}_{i,j}^n\right)+D_S\left({\mathrm{\nabla }}_S{I}_{i,j}^n\right)+D_E\left({\mathrm{\nabla }}_E{I}_{i,j}^n\right)+D_W\left({\mathrm{\nabla }}_W{I}_{i,j}^n\right)+\beta \gamma {\left(I_{ci,j}^n\right)}^{\gamma }{\left(I_{i,j}^n\right)}^{\gamma -1}$(5) . Continue iterations until converge for time step $\tau .$
Step 3	:	Generate Gaussian and Laplacian pyramid images $G_l$, $L_l$ from despeckled images $I_1$ and $I_2$ using Equation (9)(9) $\begin{array}{rl}& G_0=I\\ & G_l=Reduce\left[G_{l-1}\right]\\ & L_l=G_l-Expand\left[G_{l+1}\right]\end{array}$(9) .
Step 4	:	Involving LR operator in each Laplacian pyramid decomposition layers using Equation (12)(12) $\begin{array}{rl}& \mathfrak{D}\tilde{\mathrm{I}}{ }_{decomp}=\left[log\left(L_{I_{2}l}\right)-log\left(L_{I_{1}l}\right)\right];\\ & \begin{array}{c}{I}_1,I_2-\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{s}\\ l=0,1,2,3-\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}\mathrm{s}\end{array}\end{array}$(12) .
Step 5	:	Reconstruct the degenerated images by reversing the construction steps to obtain the difference image using Equation (13)(13) $\mathfrak{D}\stackrel{\sim }{\mathrm{I}}=reconst(\mathfrak{D}{\stackrel{\sim }{\mathrm{I}}}_{decomp})$(13) .
Step 6	:	Apply automatic Otsu thresholding algorithm to get change map using Equations (14)(14) $P_i=\frac{n_i}{N}{P}_i\ge 0;\sum _{i=1}^L{P}_i=1$(14) –(17(17) ${\sigma }_B^2\left(t^{\ast }\right)={}_{0\le k\le L}^{max}{\sigma }_B^2\left(k\right)$(17) ).
Step 7	:	Analyze the change map image.
Step 8	:	Evaluate the performance metrics.

Data set

In this work, simulated images and three real SAR datasets are utilized for experiments;

1. Simulated Images: This is a Bolivia deforestation image, monitored by the Sentinal-1B spacecraft system under the operation of ESA during 2016. This observation has considered for large-scale deforestation for sugarcane cultivation. The simulated image is generated with consideration of agricultural field area patches placed to observe changes.

2. Athabasca Dataset: The Athabasca is the largest oil sand valley that contains the severe crude oil deposition in a mixed form. This is the largest basin of tar sand in the world and also applicable to large-scale open mining. These images are acquired by ENVI-ASAR satellite during 2000 and 2005 surveillance.

3. Alberta Dataset: This is High-river region of Alberta caused by flood in June 2013. This dataset has been tracked by RADARSAT-2 satellite in space and is major disaster to analyze the extent level of the flood.

4. Aare River Dataset: The Aare River flooded utterly and entered into the parts of cities and airports. This dataset was captured by ERS-2 satellite before and after flood in April and May 1999.

Quality measures

This section evaluates the performance metrics in various stages. In the first stage, the efficiency of the preprocessing method is evaluated through the typical measures like Structural Similarity Index Measure (SSIM) (Pratt, 1978; Michailovich & Tannenbaum, 2006; kalaivani et al., 2012b), Equivalent Number of Looks (ENL) (Patel, Easley, Chellappa, & Nasrabadi, 2014) and Speckle Suppression Index (SSI) (Sheng & Xia, 1996) that are used to examine the effectiveness of speckle and are tabulated in Table 1.

Table 1. Quality measures for speckle reduction assessment

Quality Measure	Definition	Description
SSIM	$\mathrm{S}\mathrm{S}\mathrm{I}\mathrm{M}\left(\mathrm{x},\mathrm{y}\right)=\frac{\left(2{\mu }_{\mathrm{x}}{\mu }_{\mathrm{y}}+{\mathrm{c}}_1\right)\left(2{\sigma }_{\mathrm{x}\mathrm{y}}+{\mathrm{c}}_2\right)}{\left({\mu }_{\mathrm{x}}^2+{\mu }_{\mathrm{y}}^2+{\mathrm{c}}_1\right)\left({\sigma }_{\mathrm{x}}^2+{\sigma }_{\mathrm{y}}^2+{\mathrm{c}}_2\right)}$, where ${\mathrm{c}}_1,{\mathrm{c}}_2$ are constant and chosen lesser than one.	Structural similarity index is used to characterize luminosity, disparity and structural variations, where mean intensity ${\mu }_x=\frac{1}{N}{\sum }_{i=1}^N{x}_i$  Standard deviation ${\sigma }_x={\left(\frac{1}{N-1}{\sum }_{i=1}^N{\left(x_i-{\mu }_x\right)}^2\right)}^{1/2}$  Covariance ${\sigma }_{xy}=\frac{1}{N}{\sum }_{i=1}^N\left(x_i-{\mu }_x\right)\left(y_i-{\mu }_y\right)$
ENL	$\mathrm{E}\mathrm{N}\mathrm{L}=\frac{\mathrm{E}{\left[{\mathrm{f}}_{\mathrm{I}}\right]}^2}{\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathrm{f}}_{\mathrm{I}}\right)}$	Equivalent number of looks. The amplitude of an image is calculated by the ratio of mean to variance. Maximum value corresponds to the better speckle reduction. Here, ${\mathrm{f}}_{\mathrm{I}}$ is filtered image.
SSI	$\mathrm{S}\mathrm{S}\mathrm{I}=\frac{\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathrm{f}}_{\mathrm{I}}\right)}{\mathrm{E}{\left[{\mathrm{f}}_{\mathrm{I}}\right]}^2}\cdot \frac{\mathrm{E}{\left[{\mathrm{s}}_{\mathrm{I}}\right]}^2}{\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathrm{s}}_{\mathrm{I}}\right)}$	Speckle Suppression Index. It is defined by the coefficient of variance of the filtered image and normalized by that of the speckle image. SSI provides better results for lesser than one. ${\mathrm{s}}_{\mathrm{I}}$ is a speckled image.

In the second stage, the performance of CD method is qualified by various quality measures (Lal & Margret Anouncia, 2015; Shang, Qi., Jiao, Stolkin, & Li, 2014) to assess the CD image with respect to ground truth image. False positive (FP) defines the number of unchanged pixels that are wrongly classified as changed pixels. False negative (FN) defines the number of changed pixels that are wrongly classified as unchanged pixels. True negative (TN) is the number of pixels that are correctly classified as unchanged pixels. True positive (TP) is the number of pixels that are correctly classified as changed pixels. False alarm rate (FAR) is the probability of unchanged pixels that are wrongly classified as changed pixels and it is also termed as FP rate, which is close to a minimal value for ideal case. $F1$ score measures the harmonic mean of detection rate and precision. $F_1$ closer to one shows better score. The overall error (OE) is the sum of FP and FN pixels, and less number of OE represents better result. The ER measures the ratio between OE and total number of pixels $\left(N\right)$. It provides best result for minimum rate. Kappa coefficient (KC) measures the effectiveness of classification and is to be closed to 1 for ideal case. Percentage correct classification (PCC) measures accuracy of prediction. On the third stage, the critical analysis is carried on a threshold algorithm using the Region nonuniformity (NU) measure (Mehmet & Sankur, 2004) and Area overlap measure (AOM) (Rosenfield et al., 1986; Xu, Xu, Jin, & Song, 2011) to justify the meaningful change map.

Results and discussion

We have tested the proposed algorithm with one set of simulated image and three sets of the real SAR dataset images. The experimental results and its discussions are given below,

Preprocessing by CAD algorithm

In preprocessing stage, CAD algorithm is used as a despeckling tool to perform speckle noise reduction in SAR images. Speckle-free images can able to give improved response in post-processing stages. The CAD provides excellent tool for speckle noise reduction, edge preservation and feature enhancement without broadening effect. The quality of the filtering technique is verified by effective performance metrics such as SSIM, ENL and SSI. These metric values are evaluated and tabulated in Table 2 and it shows precise values for all three real SAR datasets. An additional, image profile is used to give pictorial representation to validate the intensity profile of an image, which gives visual performance to compare the original speckle image and the filtered image. In Figure 3, the performance of image profile is depicted in distance along profile with respect to intensity on 330th row of Athabasca dataset, 33rd row of Alberta dataset and 67th row of the Aare river dataset. Thus, validation by above representations’ CAD algorithm exhibits excellent speckle noise reduction without disturbing structural originality of the given image. Importantly, CAD is applied as a preprocessing tool for all dataset images to maintain uniformity among the various comparison techniques and proposed method to define exact change map analysis.

Table 2. Performance of various datasets using SSIM, ENL and SSI

Dataset	SSIM	ENL	SSI
Athabasca	0.8767	56.3572	0.2406
Alberta	0.8807	24.2520	0.1391
Aare River	0.9200	57.7557	0.2622

Figure 3. Experimental results of various dataset (a1) Athabasca image with marginal on 330th row, (a2) Image profile beside on 330th row of original image, (a3) Image profile beside on 330th row of despeckled image, (b1) Alberta image with marginal on 33rd row, (b2) Image profile beside on 33rd row of original image, (b3) Image profile beside on 33rd row of despeckled image, (c1) Aare image with marginal on 67th row, (c2) Image profile beside on 67th row of original image, (c3) Image profile beside on 67th row of despeckled image

The performance of proposed LPCD method

Flattening effect is a common problem of LR operation, which is characterized as strengthening the low-intensity pixels at the same time as reducing the strength of pixels in the area of high intensity. In order to prevent this flattening effect of high-intensity pixels and to improve the linear features and edges, LP-based CD is introduced. The proposed LP-based DI gives excellent pixel distribution between changed and unchanged classes. The LP is decomposed into $l=4$ layers and each decomposition layer is compared by LR operation to work under medium band resolution. The LP-based multiscale decomposition is constructed and verified by LR operation, and the performance of band pass and low pass difference coefficients is depicted in Figure 4 along with reconstructed image of LP. Reconstruction of LP is the DI. This multiple scale analysis enhances boundaries, lines and bars features in the first three band pass layers and verified hypothesized feature tentatively identified by a low pass layer. The proposed LP-based DI is a gray-level image and it is converted into binary change map by optimized Otsu thresholding algorithm.

Figure 4. Experimental results on LP-based multiscale decomposition combined with LR operator. Athabasca dataset (a₁) L₁₀−L₂₀, (a₂) L₁₁−L₂₁, (a₃) L₁₂−L₂₂, (a₄) L₁₃−L₂₃ and (a₅) Reconstructed image; Aare River dataset (b₁) L₁₀−L₂₀, (b₂) L₁₁−L₂₁, (b₃) L₁₂−L₂₂, (b₄) L₁₃−L₂₃ and (b₅) Reconstructed image; Alberta dataset (c₁) L₁₀−L₂₀, (c₂)) L₁₁−L₂₁, (c₃) L₁₂−L₂₂, (c₄) L₁₃−L₂₃ and (c₅) Reconstructed image

Two time series SAR images are considered as the input images to experiment for CD and the available ground truth image is manipulated with reference to the input images. The performance of the proposed LPCD method is compared with various change difference methods like Mean Ratio (MR), LR, and MR with Kittler-Illingworth algorithm (MRKI), LR with Kittler-Illingworth algorithm (LRKI), MR with Otsu algorithm (MRO), LR with Otsu algorithm (LRO), Fusing MR and LR methods with kittler-Illingworth algorithm (FMLKI) and Fusing MR and LR method with Otsu algorithm (FMLO). Comparisons with above-existing methods are depicted in Figures 5–8 for (i) simulated dataset (ii) Athabasca dataset (iii) Alberta dataset and in (iv) Aare River dataset, respectively. It is clearly understandable that the figures mentioned above have the following observations, and the MR method finds differences by local mean value of the input images. Change map obtained by MR method represents unchanged pixels are wrongly identified as changed pixels, which carry more hot spots on the CD image. The LR method fails to detect change class pixels since the characteristics of logarithmic operator influenced to fade away the high-intensity pixels. MRKI and LRKI methods used KI algorithm for thresholding, and due to that, the error pixels are reduced when compared with direct analysis. MRO and LRO utilizes Otsu threshold algorithm, and it gives better pixel distributions of changed and unchanged classes and shows lesser hot spots. FMLKI and FMLO methods fuse MR and LR operator using principal component analysis (PCA) fusion rule. Fusing MR and LR method used thresholding techniques of KI and Otsu algorithm providing asymmetrical detection between two classes of pixels.

Figure 5. Experimental results of simulated images (a) I₁ (b) I₂ (c) Ground truth (d) MR (e) LR (f) MRKI (g) LRKI (h) MRO (i) LRO (j) FMLKI (k) FMLO (l) LPCD

Figure 6. Experimental results of Athabasca Dataset (a) I₁ (b) I₂ (c) Ground truth (d) MR (e) LR (f) MRKI (g) LRKI (h) MRO (i) LRO (j) FMLKI (k) FMLO (l) LPCD

Figure 7. Experimental results of Alberta Dataset (a) I₁ (b) I₂ (c) Ground truth (d) MR (e) LR (f) MRKI (g) LRKI (h) MRO (i) LRO (j) FMLKI (k) FMLO (l) LPCD

Figure 8. Experimental results of Aare River Dataset (a) I₁; (b) I₂; (c) Ground truth; (d) MR; (e) LR; (f) MRKI; (g) LRKI; (h) MRO; (i) LRO; (j) FMLKI; (k) FMLO; (l) LPCD

The proposed method is compared with above CD methods and evaluated by performance measures using confusion matrix parameters and are tabulated in Tables 3–6 with respective (i), (ii), (iii) and (iv) datasets. The comparison results are proved that the proposed method is competent and afford excellent KC value of 0.9574 for (i), 0.8099 for (ii), 0.9018 for (iii) and 0.9064 for (iv) datasets and provide highest PCC of 99.21% for (i), 99.02% for (ii), 99.05% for (iii) and 99.78% (iv) datasets, respectively. ER 0.0079 for (i), 0.0098 for (ii), 0.0095 for (iii) and 0.0022 for (iv) are acquired, and it shows that the proposed method illustrates most of the changed pixels are identified correctly with minimum misclassification error. The rate of FAR and Miss rate (MsR) are lesser, which implicit the better class distribution of LPCD method. F₁ scores obtained by LPCD method are 0.9618 for (i), 0.8148 for (ii), 0.9068 for (iii) and 0.9075 for (iv) dataset, which are comparatively higher than above-existing methods. The Detection Rate (DR) is more for LPCD, which shows the perfect classification of the proposed CD method. The pixel classifications by FN, FP and OE are comparatively less which proves the effective performance of the proposed method. Comparison plot for KC and PCC (%) are portrayed in Figure 9.

Table 3. Evaluation of change detection for simulated images with different methods

Methods	FP	FN	FAR	F1 Score	OE	ER	Kappa	MsR	DR	PCC (%)
MR	5889	457	0.1479	0.5448	6346	0.144	0.4743	0.1074	0.8926	85.60
LR	3372	315	0.0869	0.7284	3687	0.0836	0.6819	0.0599	0.9401	91.64
MRKI	5091	522	0.1271	0.5559	5613	0.1273	0.4927	0.1294	0.8706	87.27
LRKI	1897	948	0.0491	0.7595	2845	0.0645	0.7225	0.1743	0.8257	93.55
MRO	4217	605	0.1049	0.5758	4822	0.1094	0.5202	0.156	0.844	89.06
LRO	1450	961	0.0371	0.7703	2411	0.0547	0.7393	0.1921	0.8079	94.53
FMLKI	1188	289	0.0305	0.8678	1477	0.0335	0.8488	0.0562	0.9438	96.65
FMLO	1030	280	0.0264	0.8782	1310	0.0297	0.8614	0.056	0.944	97.03
LPCD	237	111	0.0028	0.9618	348	0.0079	0.9574	0.0513	0.9487	99.21

Table 4. Evaluation of change detection for Athabasca oil sand dataset with different methods

Methods	FP	FN	FAR	F1 Score	OE	ER	Kappa	MsR	DR	PCC (%)
MR	23,383	522	0.1992	0.2273	23,905	0.1969	0.1799	0.1293	0.8707	80.31
LR	16,142	192	0.1354	0.1991	16,334	0.1345	0.1721	0.0864	0.9136	86.55
MRKI	17,304	644	0.1474	0.2744	17,948	0.1478	0.2317	0.1595	0.8405	85.22
LRKI	11,283	192	0.0947	0.2613	11,475	0.0945	0.2374	0.0864	0.9136	90.55
MRO	12,753	778	0.1087	0.3252	13,531	0.1114	0.2873	0.1927	0.8073	88.86
LRO	6910	192	0.058	0.3637	7102	0.0585	0.3445	0.0864	0.9136	94.15
FMLKI	4033	106	0.0342	0.6209	4139	0.0341	0.6055	0.0303	0.9697	96.59
FMLO	2080	378	0.0178	0.7763	2458	0.0202	0.7659	0.0814	0.9186	97.98
LPCD	1123	64	0.0095	0.8148	1187	0.0098	0.8099	0.0239	0.9761	99.02

Table 5. Evaluation of change detection for Alberta flood dataset with different methods

Methods	FP	FN	FAR	F1 Score	OE	ER	Kappa	MsR	DR	PCC (%)
MR	7783	162	0.1453	0.305	7945	0.1432	0.2628	0.085	0.915	85.68
LR	2914	304	0.0554	0.6175	3218	0.058	0.5893	0.1048	0.8952	94.20
MRKI	4313	351	0.0803	0.3791	4664	0.0841	0.3472	0.1977	0.8023	91.59
LRKI	1964	273	0.038	0.7603	2237	0.0403	0.739	0.0715	0.9285	95.97
MRO	2559	410	0.0475	0.4391	2969	0.0535	0.4158	0.2608	0.7392	94.65
LRO	1008	124	0.0191	0.8227	1132	0.0204	0.8121	0.0451	0.9549	97.96
FMLKI	1016	165	0.0189	0.7203	1181	0.0213	0.7097	0.0979	0.9021	97.87
FMLO	867	199	0.0161	0.7361	1066	0.0192	0.7265	0.118	0.882	98.08
LPCD	456	73	0.0086	0.9068	529	0.0095	0.9018	0.0276	0.9724	99.05

Table 6. Evaluation of change detection for Aare River dataset with different methods

Methods	FP	FN	FAR	F1 Score	OE	ER	Kappa	MsR	DR	PCC (%)
MR	13,481	269	0.1602	0.1096	13,750	0.1613	0.0874	0.2413	0.7587	83.87
LR	8224	132	0.0985	0.2767	8356	0.098	0.2508	0.0763	0.9237	90.2
MRKI	10,469	281	0.1244	0.1343	10,750	0.1261	0.1132	0.252	0.748	87.39
LRKI	3889	155	0.048	0.6681	4044	0.0474	0.6451	0.0367	0.9633	95.26
MRO	9170	298	0.1089	0.1319	9468	0.1111	0.1127	0.293	0.707	88.89
LRO	2560	275	0.0314	0.701	2835	0.0333	0.6845	0.0764	0.9236	96.67
FMLKI	1922	161	0.0229	0.5372	2083	0.0244	0.5266	0.1175	0.8825	97.56
FMLO	1433	207	0.0171	0.5865	1640	0.0192	0.5776	0.1511	0.8489	98.08
LPCD	148	40	0.0018	0.9075	188	0.0022	0.9064	0.0416	0.9584	99.78

Figure 9. Comparison Result of PCC for: (a) Athabasca dataset, (b) Alberta dataset, (c) Aare River dataset; Comparison Result of Kappa Coefficient for: (d) Athabasca dataset, (e) Alberta dataset, (f) Aare River dataset

The effectiveness of the proposed method is verified by comparing with other multiscale approaches like wavelet-based multiscale decomposition method with scale-driven approach (SDA) (Bovolo & Bruzzone, 2005), multiscale change analysis of similarity measure (KL-MCP) (Inglada & Mercier, 2007), multiscale CD by bayesian approach (BA-MCD) (Celik, 2010a) and multiscale SAR image preparation by the kennaugh element framework (K matrix-MS) (Schmitt, Wendleder, & Hinz, 2015), and the comparison result is shown in Figure 10. The analysis of above various multiscale representations give the following facts that, in Figure 10(a), wavelet-based scale-driven approach showed better accuracy based on trial-and-error threshold selection, also this technique is unable to provide correct classification of change map upto the maximum extent. In Figure 10(b), KL-MCP is statistical similarity measure and its extension is considered for multiscale representation, which is wrongly detecting unchanged class pixels as a change class pixels and also shows false prediction on class distribution. Figure 10(c) shows BA-MCD approach, which has not identified number of scale selections for effective CD, and the computational complexity is more. In Figure 10(d), K matrix-MS is Kennaugh element framework for multiscale SAR image preparation to detect change map, which gives CD result in concern with number of looks of the SAR image. This framework gives blurring effect on the edges of each scale and takes more computation complexity. By examination of other multiscale approaches, it is obvious that the proposed LPCD method yields better performance, easy execution, exact class distribution among changed and unchanged pixels, edge enhancement and accurate area detection. The comparison of various performance measures is evaluated with above-existing approaches in terms of accuracy, FAR, ER, DR and KC, and are tabulated in Table 7. This performance measure shows that the LPCD method has better performance of high accuracy, less FAR and ER, high DR, high KC and less computation time.

Table 7. Comparison of different multiscale methods

Multiscale methods	FAR	ER	DR	Kappa	PCC (%)	CT (sec)
SDA	0.4328	0.0188	0.9454	0.6999	98.12	2.3323
KL-MCP	0.6446	0.5806	0.6071	0.0251	41.94	1.8937
BA-MCD	0.0180	0.0203	0.9222	0.7671	97.97	3.1409
K matrix-MS	0.0244	0.0237	0.8321	0.6983	97.63	4.8321
Proposed method	0.0095	0.0098	0.9761	0.8099	99.02	0.9639

Figure 10. Experiment results of comparison with multiscale methods (a) SDA (b) KL-MCP (c) BA-MCD (d) K matrix-MS (e) Proposed LPCD method

Class distribution by Otsu thresholding algorithm

The perfect classification of the CD is carried out by Otsu thresholding algorithm. The optimal threshold value $\left(t\right)$ depends on the maximum class probability and class mean, where the class mean is the first-order statistical criterion (cumulative probability). In our proposed technique, Otsu method selects the threshold by exploiting the between-class variance $({\sigma }_B^2)$. The Figure 11 shows the pictorial representation of between-class variance ${\sigma }_B^2$ for various dataset images (Fan & Lei, 2012). Figure 11(a,c) represents the plot of ${\sigma }_B^2$ value with corresponding gray-level value. Figure 11(d,f) represents the histogram of an image with corresponding threshold value adapted by various CD methods (Muthukumaran, Gopalakrishnan, Purna, Soumitra, & Saravanan, 2017) and shows effective class distribution and automatic thresholding of the proposed method.

Figure 11. Otsu thresholding results: ${\sigma }_{\mathrm{B}}^2$ Vs gray level (a) Athabasca image, (b) Alberta image, (c) Aare river image; Histogram with threshold values: (a) Athabasca image, (b) Alberta image, (c) Aare river image

Thresholding method offers fine distribution of changed and unchanged pixels of the image. Here, Otsu algorithm is capable to distinguish the changed/unchanged criterion automatically, even the number of changed pixels is lesser than 5% of the whole image pixels. In Figure 12, small change area (tiny portion with <5% of changed pixels) of the image is considered for testing.

Figure 12. Experiment results of Alberta dataset for Otsu thresholding (a) ${\sigma }_{\mathrm{B}}^2$ (b) ${\sigma }_{\mathrm{W}}^2$

Otsu’s approach is experimented with ${\sigma }_w^2$ (within-class variance) and ${\sigma }_B^2$ (between-class variance). The maximum value of between-class variance ${\sigma }_{B\left(max\right)}^2$ can detect number of changed pixels closely and provides less OE. The simulation result of ${\sigma }_w^2$ and ${\sigma }_B^2$ are experimented for small percentage of changed pixels. This experiment shows ${\sigma }_B^2$ detected number of 3020 changed pixels with respect to ground truth of 3033 pixels, but ${\sigma }_w^2$ detected only a number of 2150 pixels. This shows minimum value of within-class variance ${\sigma }_{w\left(min\right)}^2$ is fails to detect foreground pixels exactly and spatial information are not considered into an account. Hence, the within-class variance ${\sigma }_{w\left(min\right)}^2$ is not suitable to classify less than 5% of change area of the whole image (Alexander et al., 2011) which is verified by NU and AOM, and is tabulated in Table 8.

Table 8. Experimental results for small percentage of changed pixels

Dataset	Image size	Ground Truth (changed pixels)	% of change	Class Separability	Change area pixels	OE	NU	AOM
Alberta Images	$175$ x $317$	$3033$	≤5%	${\sigma }_W^2$ (Within-class)	$2150$	$1385$	$0.0392$	$0.5783$
				${\sigma }_B^2$ (Between-class)	$3020$	$433$	$0.0551$	$0.8665$

The above experiment results show ${\sigma }_B^2$ provides valid information of changed and unchanged class distribution compared to ${\sigma }_w^2$. The proposed LPCD method is classified by the maximum value of between-class variance at ${\sigma }_B^2=270.5766$ and is obtained at threshold value of $t=58$. The consideration of ${\sigma }_B^2$ dataset (ii) is used to experiment for various values of between-class variance, and the test image is classified by minimum to maximum threshold values in Figure 13. The Figure 13(e) represents that the foreground and background pixels are effectively classified compared to other ${\sigma }_B^2$ values. Table 9 carries the performance of AOM evaluated for different values of between-class variance ${\sigma }_B^2$.

Table 9. Area overlap measure (AOM) for Otsu thresholding for different between class variance value in values. Fig.13

	Otsu segmentation for different between-class variance values.
	${\sigma }_B^2=124.76$	${\sigma }_B^2=164.33$	${\sigma }_B^2=221.26$	${\sigma }_B^2=256.84$	${\sigma }_B^2=270.58$	${\sigma }_B^2=251.27$	${\sigma }_B^2=225.65$	${\sigma }_B^2=201.78$	${\sigma }_B^2=175.15$	${\sigma }_B^2=150.70$	${\sigma }_B^2=124.19$	${\sigma }_B^2=99.68$
AOM	0.0479	0.0750	0.2223	0.5829	0.8542	0.7386	0.5682	0.4579	0.3569	0.2832	0.2138	0.1307

Maximum value of between class variance${\sigma }_B^2=270.58$ yields accurate binary classification experimented by AOM.

Figure 13. Otsu segmentation result for Athabasca image (a) ${\sigma }_{\mathrm{B}}^2$ = 124.76; Threshold = 7, (b) ${\sigma }_{\mathrm{B}}^2=164.33$; Threshold = 14, (c) ${\sigma }_{\mathrm{B}}^2=221.26$; Threshold = 27, (d) ${\sigma }_{\mathrm{B}}^2=256.84$; Threshold = 39, (e) ${\sigma }_{\mathrm{B}}^2=270.57$; Threshold = 58, (f) ${\sigma }_{\mathrm{B}}^2=251.27$; Threshold = 88, (g) ${\sigma }_{\mathrm{B}}^2=225.65$; Threshold = 104, (h) ${\sigma }_{\mathrm{B}}^2=201.78$; Threshold = 115, (i) ${\sigma }_{\mathrm{B}}^2=175.15$; Threshold = 126, (j) ${\sigma }_{\mathrm{B}}^2=150.70$; Threshold = 135, (k) ${\sigma }_{\mathrm{B}}^2=124.19$; Threshold = 145, (l) ${\sigma }_{\mathrm{B}}^2=99.68$; Threshold = 154

The region NU and AOM are performance measurements for evaluating the effective classification of the thresholding algorithm. The assessment of effective classification of LPCD and various change map techniques are measured by AOM and NU, shown in Figure 14(a,b) respectively. The performance measures of NU and AOM are itemized in Tables 10 and 11 respectively.

Table 10. Performance measures (region nonuniformity)

Dataset	MRKI	MRO	LRKI	LRO	FMRKI	FMLO	LPCD
Athabasca	0.1705	0.1097	0.1319	0.0736	0.0611	0.0523	0.0296
Alberta	0.1326	0.0934	0.116	0.069	0.0367	0.0304	0.0126
Aare River	0.1034	0.0993	0.0671	0.0655	0.0557	0.0524	0.0446

Table 11. Performance measures (area overlap measure)

Dataset	MRKI	MRO	LRKI	LRO	FMRKI	FMLO	LPCD
Athabasca	0.159	0.1503	0.1942	0.2223	0.4503	0.6343	0.8542
Alberta	0.2307	0.5016	0.3316	0.5397	0.3673	0.4149	0.8306
Aare River	0.2339	0.6132	0.2813	0.6988	0.5629	0.5825	0.8296

Figure 14. Comparison results of (a) area overlap measure (AOM), (b) region nonuniformity (NU)

From an analysis of the proposed LPCD method it can be seen that the overall performance perceives best change map detection. The LP-based multiscale DI incorporated with speckle reduction and thresholding algorithm provide appreciable CD results. It enables speckle-free correct classification result and less computational time of 0.9639s compared to other existing multiscale methods. These significance makes LP multiscale analysis is apt for CD application and this method is suitable to work under crisis situation in disaster management. The entire experiments are executed on MATLAB in windows platform on the system configuration of Intel(R) Core(TM) i5-4210U [email protected] 2.40GHz, 4GB of RAM for SAR images.

Conclusions

In this paper, a novel CD method based on LP multiscale representation and the effective classification by thresholding algorithm is proposed. The significance of the proposed method is to restrain the unchanged class pixels ideally and to enhance the changed class pixels effectively with perfect boundaries. As specified in the literature, the characteristic of LR operation is flattening the high-intensity pixels and enhancing the pixels in the area of low intensity, which provides asymmetrical class distribution. Hence, the multiscale representation by LPCD is introduced to detect changes in medium band resolution and it provides changed and unchanged class distributions exactly. The proposed LP-based CD method is experimented with a simulated image and three real SAR image datasets.

The SAR images are highly degraded by speckle noise, hence preprocessing is an essential task for CD. In the first step, CAD is employed as a preprocessing tool for removing the speckle noise. CAD effectively suppresses the speckle content and preserve edges, enhance sharp details and reduce feature broadening effects. Profile for despeckled image shows the effective speckle reduction ability of CAD, which is proved by the computation of quality measures like SSIM, ENL and SSI. In second stage, speckle-free images are processed by LP domain, where the input images are decomposed into three $l=0,1,2$ band pass layers and a unique $l=3$ low pass layer of the image. The proposed method preserve edges from the fine scale layers and the detailed information tentatively presents in coarse scale layer. The decomposed pyramid coefficients are manipulated by a LR operation and the reconstructed image is a DI. In the third stage, a threshold algorithm of Otsu’s method is applied on the DI to realize the automatic classification of the change map image. With the aid of thresholding experiment by two different variance approaches, it has proved that, our proposed method ${\sigma }_{B\left(max\right)}^2$ can detect even less than 5% of changed class pixels (tiny area) of the whole image.

The quality of the LPCD method is analyzed by various standard measures like PCC, KC, F₁- score, FAR, DR and so on. The proposed method affords perfect CD map with the preservation of required details. The correct classification of threshold algorithm is verified by the region NU measure and AOM. The efficacy of the proposed method is compared with various multiscale approaches. LP-based multiscale CD method provides exact change information from coarse to fine scale layer and can deal with a large selection of changed pixels with less computing time of 0.9639 s compared to other existing methods, which makes the proposed method is suitable to work under emergency conditions in disaster management. The future development of this work could be carried out by considering the similarity measures in pyramid layers which may provide higher number of dependency pixels between two successive layers of the image. Similarity measures are insensitive to sensor parameters, hence would make fast evaluation of changed and unchanged pixel separation framework. Also, it would be a concern to extend the proposed multiscale approach for multiband and polarimetric SAR images.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Disclosure statement

No potential conflict of interest was reported by the authors.