SAR image change detection based on equal weight image fusion and adaptive threshold in the NSST domain

ABSTRACT

In order to improve the accuracy of change detection and reduce the running time, a change detection method based on equal weight image fusion and adaptive threshold in the NSST domain is proposed. First, the logarithmic transformation is used to transform images and the mean filter is applied to the transformed images. The log-ratio method and the mean ratio method are adopted to generate two kinds of difference images. The final difference image is achieved by equal weight image fusion method. Then, an adaptive threshold denoising method based on non-subsampled shearlet transform (NSST) is used to achieve noise reduction. Finally, the k-means clustering algorithm is utilized to get the change detection results. The experimental results show that the proposed algorithm has better change detection performance than the reference algorithms in visual effect and objective parameters.

KEYWORDS:

• Non-subsampled shearlet transform
• image fusion
• change detection
• difference map
• adaptive threshold
• k-mean algorithm

Introduction

Change detection is a process which uses multi-temporal remote sensing images to obtain the change information of the same area, that is to say, to qualitatively analyze and confirm the features and processes of the change of the earth’s surface from the remote sensing data taken at different times (Radke, Andra, Al-Kofahi, & Roysam, 2005; Singh, 1989). There are supervised methods and unsupervised methods for image change detection (Coppin, Jonckheere, Nackaerts, Muys, & Lambin, 2004; Wang, Du, & Dai, 2016). The former requires the known classifier to train the samples, and then use the trained samples to detect the change area of the images. The known classifier requires training set obtained from the ground reference data. The latter does not require a prior knowledge and can be applied directly to images (Celik & Ma, 2010; Hu & Ban, 2014). Obtaining ground reference data is expensive for supervised detection methods. Therefore, the lack of ground reference data makes efficient unsupervised detection methods more widely used in practical applications (Gopal & Woodcock, 1996).

At present, a lot of research work has focused on multiscale decomposition of images. Among them, the typical multiscale decomposition method is the wavelet transform (WT) (Gupta, Chauhan, & Sexana, 2004). As a typical decomposition tool, it can only describe the singularity of points. Nevertheless, representation the directionality of images is impossible for the WT. For the more sparsely representation of the image, the contourlet transform was proposed by Martin Vetterli in 2002(Do & Vetterli, 2005). However, the contourlet transform does not have translation invariance, and it is prone to pseudo Gibbs phenomenon, resulting in image distortion and information loss. Non-subsampled contourlet transform (NSCT) (Da, Zhou, & Do, 2006; Li, Fang, & Yin, 2012) solves this problem. Owing to the subdivision of the frequency space, the sparse representation of the image (Olshausen & Field, 1996) is weakened to a certain extent, and the computation process is relatively complex. Therefore, Guo et al (Guo & Labate, 2007) proposed the construction of the shearlet transform (ST) using the synthetic expansion affine system. But ST does not have translation invariance. Non-subsampled shearlet transform (NSST) (Kong, 2014) is an extension of the ST. NSST avoids the down sampling process of the ST and has translation invariance. And the problem of pseudo Gibbs phenomenon is solved effectively. NSST is a tool for multiscale decomposition. It can not only sparsely represent high-dimensional signals, but also have the best approximation order for edge curves in images. And NSST has been successfully applied in image denoising. Hou et al proposed SAR image denoising based on non-subsampled shearlet (Hou, Zhang, Bu, & Feng, 2012). The coefficients of the non-subsampled shearlet can be regarded as the projection of the image on an anisotropic basis function, while the noise is isotropic, and the coefficients are uniformly distributed in all directions. In this article, the adaptive threshold denoising algorithm based on non-subsampled shearlet domain is proposed to denoise the difference image.

The subtraction method and the ratio method are often used to produce the difference map (Jia et al., 2016). For example, difference operator and log-ratio operator are applied to obtain two kinds of difference images in the literature (Zheng, Zhang, Hou, & Liu, 2014). The probabilistic patch-based (PPB) filtering method is utilized to suppress noise in the images before obtaining the two difference maps. The final difference image is clustered by k-means (PPB-K). The algorithm improves the precision of change detection. However, the objective indicators need to be further improved. In the literature (Celik, 2009), the difference image is obtained by using the ratio operator for synthetic aperture radar (SAR) images. The subtraction operator is used to obtain the difference image for optical images. Owing to the single type of the difference image, the detection result is not ideal for some images. Gong et al proposed a Neighborhood-based ratio method to generate difference map (Gong, Cao, & Wu, 2012). Although the algorithm is superior to the traditional method of generating difference maps, the detection results are seriously affected by noise and the denoising effect is not very satisfactory. Compared with the subtraction method, the ratio method can effectively eliminate the speckle noise in SAR images. The subtraction method usually describes quantitatively whether the target area has changed or not, and it is difficult to determine the nature of the change of the target area, resulting in a higher error detection rate. The ratio operator includes the log-ratio operator and the mean ratio operator (Inglada & Mercier, 2007). Gong et al fuse the mean ratio difference image and the log-ratio difference image, and obtain the better detection results (Gong, Zhou, & Ma, 2012). Ma et al. proposed an unsupervised image change detection method based on data fusion and fuzzy clustering (Ma, Li, Wu, Jiao, & Xing, 2014). The context-independent variable behaviour method is used to fuse two difference images obtained the mean ratio operator and the log-ratio operator (CIVB-F), respectively. Although the algorithm reduces the error rate of the detection results, it requires a long running time to obtain the change detection results.

According to the above analysis, a single type and the quality of difference map affect the results of image change detection. To solve the above problems, a novel unsupervised SAR remote sensing image change detection algorithm is proposed in the article. The algorithm uses logarithmic transformation to convert images into logarithmic domain. And mean filtering is used to the transformed images to suppress background information. Log-ratio operator and mean ratio operator are applied to obtain two kinds of difference images in logarithmic domain. The final difference image is obtained by simple image fusion algorithm. Then the adaptive threshold algorithm in NSST domain is adopted to reduce image noise. Finally, the k-means clustering algorithm is employed to obtain the final detection results.

Production of the final difference image

The difference map based on the source image contains more background information, which has a negative impact on the subsequent detection work. Therefore, before obtaining the difference images, logarithm-2 transform is employed to convert the SAR images into the logarithmic domain, and the transformed image is processed by mean filtering to suppress background information. We will analyze the advantages of the difference map based on transformed image in the next section. Then we use the mean ratio operator and log-ratio (Bazi, Bruzzone, & Melgani, 2005; Carincotte, Derrode, & Bourennane, 2006; Hou, Wei, Zheng, & Wang, 2014) operator to generate two kinds of difference images. Finally, the final difference image is obtained by equal weight fusion algorithm.

Logarithmic transformation

Suppose that I₁ = {I₁(i, j), 1 < i < M, 1 < j < N} and I₂ = {I₂(i, j), 1 < i < M, 1 < j < N} are two temporal SAR images acquired over the same area. The size of the each image is M × N pixels. We can acquire U₁ and U₂ in the logarithmic domain transformed by logarithm-2 transform (Wang, Jia, Yang, & Kasabov, 2017). And the pixel values in the images I₁ and I₂ are compressed between 0 and 8. The transformation process is expressed as follows:

(1)

(2)

where I₁(i, j)+ 1 and I₂(i, j)+ 1 are used instead of X₁(i, j) and X₂(i, j), respectively, to avoid a case in which the pixel values in X₁(i, j) or is zero.

The two kinds of difference images

The window size for the mean filter is set as 3 × 3. The change image F₁ and the change image F₂ are obtained by the following formulas, respectively.

(3)

(4)

Where μ₁(i, j) and μ₂(i, j) represent the local mean values of U₁ and U₂ at the pixel point (i, j), respectively. To avoid a case in which the pixel values in U_i (i = 1,2) is zeros, U₁(i, j)+ eps and U₂(i, j)+ eps are applied instead of U₁(i, j) and U₂ (i, j), respectively. After using the mean ratio operator and the log-ratio operator, we normalize F₁ and F₂ to the range [0, 8].

The final difference image

Many of scholars have fused the difference maps produced by the various operators, and achieved good results. However, it is relatively complicated to implement and may increase the running time. Ma et al. proposed a change detector to fuse the log-ratio (LR) and the mean-ratio (MR) images by a context independent variable behavior operator (Ma et al., 2014.). Although this method improves detection precision, the running time is increased. Combined difference image method was proposed by Zheng et al (Zheng et al., 2014). However, the method requires many tests to determine the value of the weight, which results in a decrease in detection efficiency. We use equal weight image fusion method to get the final difference image to preserve more image change information. In order to make the gray value of the final difference image between 0 and 8, the sum of the weight of the two kinds of difference image should be 1. The final difference image R obtained by formula (5).

(5)

The implementation of denoising algorithm

We use the proposed adaptive threshold denoising algorithm in NSST domain to suppress noise for the final difference image. The high-frequency coefficients and the low-frequency coefficients can be obtained by conducting NSST to the difference image. The low frequency sub-band coefficient contains outline information of the image, and the low coefficient is usually larger and more evenly distributed. So we do not process the low frequency coefficient. The high frequency coefficients contain detail information and a lot of noise. The higher the absolute value of the high frequency coefficient, the more details of the image. Therefore, the traditional classical single threshold method is easy to weaken the detailed information and always ignores the neighborhood information’s relationship to the non-subsampled shearlet coefficients of the image. In order to keep the detail information of the difference image, an adaptive threshold denoising algorithm based on NSST domain is proposed in the article.

Adaptive threshold algorithm in NSST domain

The threshold denoising algorithm, also known as the Shrinkage, was first proposed by Donoho et al (Donoho & Johnstone, 1994). The hard threshold method is described:

(6)

where c(m, n, s, d) is a noisy sub-band coefficient dependent on scale index s and direction index d. c’(m, n, s, d) denotes the estimates for the denoised coefficients. T(s, d) is the threshold function and the coefficients K is used as the threshold coefficient in this paper T(s, d) = Kα (Ghosh, Mishra, & Ghosh, 2011). We make K equal to the decomposition scale by experiments on a large number of images. The standard deviation σ of noise is:

(7)

Consider that A₁(s, d, m, n) and A₂(s, d, m, n) are sub-band coefficients generated by mean filtering and median filtering, respectively.

(8)

(9)

where B is the neighborhood window of the high frequency sub-band coefficients c(s, d, m, n). m, n is spatial coordinate index. S is the number of coefficients in the neighborhood window. The neighborhood window size of the mean filter and median filter are set as 3 × 3.

A₁(s, d, m, n) and A₂(s, d, m, n) reflect the distribution of the high frequency sub-band coefficients. The bigger the A₁(s, d, m, n) is, the higher the coefficient energy in the region and the richer the edge and detail information contained within the region. Therefore, the threshold should be set with a smaller value, and the original coefficients are preserved as much as possible. Otherwise, the threshold should be set with a larger value, so as to eliminate noise as much as possible. A₂(s, d, m, n) determines noise standard deviation σ to a certain extent. The closer the two values are, the closer the coefficient values in the region are and the more uniform the distribution is. Otherwise, the coefficient values in this region differ greatly and the distribution is uneven. Considering the correlation between A₁(s, d, m, n) and A₂(s, d, m, n), an adaptive threshold denoising method is proposed in this article, which is defined as:

(10)

The formula (10) can adaptively adjust the threshold, and effectively overcome the defect that the hard threshold algorithm excessively loses the original information.

Description of denoising process

The final difference image is denoised by using the proposed adaptive threshold denoising algorithm based on NSST. The steps are as follows:

1. The low frequency coefficient and the high frequency coefficients{C, C^s,d}(s ≥ 1) can be obtained by conducting NSST to the difference image R. C is the low frequency coefficient. C^s,d are the high frequency coefficients.

2. The denoised high frequency sub-band coefficients are obtained by the proposed adaptive threshold denoising algorithm. The formula is:

   (11)

3. The reconstructed difference image is obtained by inverse NSST using the low frequency sub-band coefficient C and the high frequency sub-band coefficients

K-means clustering algorithm

Extracting the change region from a difference image can be regarded as dividing the difference image into two categories: changed area and unchanged area. The common segmentation algorithms include threshold algorithms (Patra, Ghosh, & Ghosh, 2011) and clustering algorithms (Ghosh et al., 2011; Mishra, Ghosh, & Ghosh, 2012) when extracting the change region from the difference image. Threshold algorithms need to establish a statistical model for the difference image. The complexity of statistical features makes it difficult to build an accurate model. Compared with the threshold algorithm, the clustering algorithms do not need to consider the statistical characteristics of the image and can effectively complete the classification for the image. The k-mean clustering algorithm (Celik, 2010; Yetgin, 2012) requires a short computation time, with simple operation. In this article, we use the k-mean algorithm to cluster the obtained difference image.

According to the Euclidean distance based on formula (12), we can acquire the change detection image C. The process is performed:

(12)

where V_c represents the mean feature vectors for the changed class and V_d represents the mean feature vectors for the changed class. ‘1’ and ‘0’ represent the changed pixel and the unchanged pixel, respectively.

Experimental study

First, we analyze the advantages of the difference image based on logarithmic domain to explain why we use the difference image based on logarithmic domain instead of using the difference image based on source images. Then, two different SAR image data sets are used for experiments to verify the effectiveness of the proposed algorithm. The proposed algorithm and reference algorithms are tested, and the detection results of different algorithms are obtained. We compare the results obtained by different algorithms from two aspects: visual effects and objective quantitative indices. The quantitative indices include the number of false negatives (FN), the number of false positives (FP), the overall error (OE), the Kappa index (K), and the running time (T) (Gong et al., 2012). The K is usually adopted to express the similarity. When K = 1, the detection result map is in complete agreement with the reference map (Ma et al., 2014).

The analysis of difference image

Taking the Yellow River Estuary data set as an example, we use the mean ratio method and the log-ratio method to produce two kinds of difference maps based on the source images. Then, the two methods are used to generate two kinds of difference images based on logarithmic domain. The various difference maps are shown in the following Figure 3.

Figure 1. The Ottawa data set (a) 1997,05; (b) 1997,08; (c) the reference image.

Figure 2. The Yellow River Estuary data set (a) 2008,06; (b) 2009,06; (c) the reference image.

Figure 3. The various difference maps (a) log-ratio difference map based on source image (b) mean ratio difference map based on source image (c) log-ratio difference map based on logarithmic domain (d) mean ratio difference map based on logarithmic domain.

As can be seen in Figure 3, the difference maps based on the source images contain a lot of background information, which cover some details and edges of the difference map. Obviously, this may reduce the accuracy of the detection result. The background information in (c) and (d) have been significantly reduced. Therefore, the method for generating difference maps can effectively improve the recognition of the change region, and suppress the background information of the non-change region in the difference maps.

The description of experimental data

The experimental images we used are all registered. A radiometric calibration is not applied in the paper. The first data set used in the experiment is a pair of two SAR images (290 × 350 pixels) from Ottawa, Canada. They are shown in Figure 1(a) and Figure 1(b). They are obtained by Radarsat in May 1997 and August 1997, respectively. In the two SAR images, the change area selected as a test site was caused by floods. Figure 1(c) shows change reference image. The white pixels represent the changed area in Figure 1(c).

The second data set consists of two SAR images, which were acquired by Radarsat-2 at the region of Yellow River Estuary in China. The two SAR images was taken in June 2008 and June 2009, respectively. Figure 1(c) shows the reference image. The original size (7666 × 7692 pixels) is too large to use them in the pages. The size of selected images are 257 × 289 pixels. They are shown in Figure 2(a) and Figure 2(b), respectively. Figure 2(c) shows the change reference image in which the white pixels show the changed area. The image in Figure2(a) is single look and the image in Figure2(b) were four looks, which increases the difficulty of detection due to the great difference in noise level.

Experimental results and comparisons

The proposed algorithm is compared with CIVB-F and PPB-K. The qualitative results are given in Figure 4 and Figure 5. As can be see from the results of visual effect in Figure 4 and Figure 5, Figure 4(c) shows that the detection result map effectively preserves the details of change. Figure 5 shows the superiority of the proposed algorithm in denoising performance. There is no isolated clutter in the Figure 5(c). However, a lot of isolated noise are distributed over non-change region of detection result map shown in Figure 5(a) and Figure 5(b) achieved by PPB-K algorithm and CIVB-F algorithm, respectively. By comparing the detection results of the reference algorithms and the proposed algorithm, the noise in the detection results achieved by the proposed algorithm is obviously reduced. Most of the change areas can be detected correctly and more similar to the change reference maps.

Figure 4. Change detection results of the Ottawa data set obtained by various algorithms: (a) PPB-K (b) CIVB-F (c) Proposed.

Figure 5. Change detection results of the Yellow River Estuary data set obtained by various algorithms: (a) PPB-K (b) CIVB-F (c) Proposed.

In order to make a comparison more intuitively, the PPB-K algorithm, CIVB-F and the proposed algorithm are quantitatively analyzed. As shown in Table 1, the proposed algorithm has the shortest change detection time and the highest detection accuracy compared with reference algorithms. The FN and FP in the detection results of the proposed algorithm method are greatly reduced. The K of Proposed is higher than that of the PPB-K algorithm and the CIVB-F algorithm. Although the number of false positives of PPB-K algorithm is lower, the algorithm has a higher OE. All objective quantitative indices of the proposed algorithm except time are higher. For example, in objective quantitative indices achieved by the various algorithm for Yellow River Estuary data, the OE pixels in this paper is reduced by 1108 as compared with the PPB-K algorithm, and the OE pixels is decreased by 2791 as compared with CIVB-F algorithm. As can be seen in Table 1, The objective indices of the Ottawa data indicate that the OE pixels is reduced by 135 and 864 as compared with the PPB-K algorithm and CIVB-F algorithm, respectively. As we can see from quantitative analysis for the Ottawa and Yellow River Estuary, our proposed algorithm has better detection performance. We use confusion matrix like Table 2 as percentages as well to show the detection results of the proposed approach for Ottawa image and Yellow River Estuary image.

Table 1. Change detection results of the two data sets obtained by various methods.

Data set	Algorithm	FP	FN	OE	K	T(s)
Ottawa	PPB-K	336	1368	1704	0.9353	72
	CIVB	241	2129	2433	0.9063	69.7
	proposed	720	849	1569	0.9418	8.5
Yellow	PPB-K	1530	2451	3981	0.8155	61.5
River	CIVB	880	4784	5664	0.7169	72.6
Estuary	proposed	818	2055	2873	0.8659	7.8

Table 2. Confusion matrix.

	Reference change result
	Type	Changed	Unchanged
Change detection result	Changed	n11	n12
	Unchanged	n21	n22

In the table, n₁₁ indicate the ratio of the number of changed pixels correctly detected to the total number of pixels. n₁₂ is the ratio of FP to the total number of pixels. n₂₁ is the ratio of FP to the total number of pixels. n₁₁ indicate the ratio of the unchanged number of pixels correctly detected to the total number of pixels. The confusion matrixes of Ottawa image and Yellow River Estuary image set are shown in (Tables 3 and 4), respectively.

Table 3. Confusion matrix of Ottawa.

	Reference change result
	Type	Changed	Unchanged
Change detection result	Changed	15%	0.8%
	Unchanged	0.8%	83.4%

Table 4. Confusion matrix of Yellow River Estuary.

	Reference change result
	Type	Changed	Unchanged
Change detection result	Changed	15.3%	1.1%
	Unchanged	2.7%	80.8%

In order to further illustrate the effectiveness of our proposed algorithm, The reference algorithms and the proposed algorithm are used to process the SAR image of 20 groups to obtain the average of quantitative indicators. The results are shown in Table 5. The average running time of the proposed algorithm is 11 s, much faster than those of PPB-K and CIVB-F of 71 s and 65 s, respectively. The kappa coefficient of the proposed algorithm are also the highest compared with the reference algorithms. Comprehensive analysis, the proposed algorithm has higher detection efficiency and higher detection accuracy.

Table 5. Change detection results obtained by various methods for detecting the average change of SAR images.

Algorithm	FP	FN	OE	K	T(s)
PPB-K	794	1607	2401	0.9097	71
CIVB	241	2092	2333	0.9163	65
proposed	681	864	1545	0.9426	11

We also compare the proposed algorithm with the another algorithm (AA) in our previous paper (Chen, Zhang, Jia, Jie, & Kasabov, 2017) that has been published. We use the AA algorithm to detect the 20 sets of images used in the article and get the average of the objective index. The average objective index K obtained by AA algorithm is 0.9497 and the average T is 24 s. The average objective index K obtained by AA algorithm is 0.9426 and the average T is 11 s. It shows that there is not much difference between the two algorithms in performance of keeping the consistent with the reference image. The proposed algorithm greatly improves the detection efficiency under the premise that the detection accuracy of the two algorithms is almost equal. In practical application, we need to process massive images. The higher detection efficiency can save us a lot of time. Therefore, the algorithm in the article has a better application value in detecting changes of SAR images.

The analysis of performance against noise of algorithms

Different levels of speckle noise are added to remote sensing image X₁ to test algorithm’s performance against noise. The peak signal-to-noise ratio (PSNR) is used to measure the level of noise (Celik, 2009). Let us consider an input image X₁ of size M × N and add speckle noise image X’₁. The PSNR is calculated as follows:

13

Firstly, the two phase remote sensing images X₁ and X₂ are used to generate the binary image CM₁. Then, the speckle noise is added to image X₁ to obtain the noisy image X’₁. The X’₁ and X₂ are used to generate the binary image CM₂. Finally, the robustness of the change detection algorithm against the speckle noise is measured by measuring the difference between CM₁ and CM₂ (Celik, 2009). The formula is:

14

where k reflects the robustness of the algorithm against noise, when the CM₁ and CM₂ are the same, and the value of k is 1. The closer the value of k is to 1, the better the robustness of algorithm’s performance against noise.

We add different intensity speckle noise to the phase 1 of Ottowa image and the phase 1 of Yellow River Estuary image to test robustness of the proposed algorithm and reference algorithms against different levels of noise. The performance results is shown in Figure 6. The PSNR value decreases as the noise increases. The k value obtained by PPB-K algorithm is smaller and away from 1 when PSNR≤ 40dB. The k value got by the CIVB-F algorithm is smaller than the k value achieved by the proposed algorithm when PSNR≤ 36dB. Through the analysis of Figure 6, the range of the k obtained by the proposed algorithm is small and the k value close to 1, which proves that the proposed algorithm has better robustness.

Figure 6. Change detection performance comparisons of various algorithms against different levels of noise. (a) Result on Ottowa image and (b) result on Yellow River Estuary image.

Conclusions

A new SAR image change detection algorithm is proposed in this article. In order to suppress the background information in the difference map, logarithmic transform is carried out to the SAR images. Then, the final difference image is obtained by equal weight image fusion method, so that the information in the two kinds of difference image is mutually complementary. NSST is applied to the final difference image. And the de-noised difference map can be acquired by using the proposed adaptive threshold method in the NSST domain. Finally, the simple k-means clustering is taken to divide the de-noised difference map into two classes. Experiments on real SAR image data sets show that our approach can improve the change detection results and has better performance in suppressing noise. Compared with the detection results obtained by the reference methods, the detection results obtained by the proposed algorithm are greatly improved in terms of the visual effect and objective quantitative indices.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the International Cooperative Research and Personnel Training Projects of the Ministry of Education of the People’s Republic of China [No. 2014-2029 and 2016-219] and the Natural Science Foundation of Xinjiang [No. 2015211C288].