Urban Change Detection Based on Dempster – Shafer Theory for Multitemporal Very High-Resolution Imagery

Fusing multiple change detection results has great potentials in dealing with the spectral variability in multitemporal very high-resolution (VHR) remote sensing images. However, it is difficult to solve the problem of uncertainty, which mainly includes the inaccuracy of each candidate change map and the conflicts between different results. Dempster–Shafer theory (D–S) is an effective method to model uncertainties and combine multiple evidences. Therefore, in this paper, we proposed an urban change detection method for VHR images by fusing multiple change detection methods with D–S evidence theory. Change vector analysis (CVA), iteratively reweighted multivariate alteration detection (IRMAD), and iterative slow feature analysis (ISFA) were utilized to obtain the candidate change maps. The final change detection result is generated by fusing the three evidences with D–S evidence theory and a segmentation object map. The experiment indicates that the proposed method can obtain the best performance in detection rate, false alarm rate, and comprehensive indicators.


Introduction
Remote sensing imagery has been one of the major data sources for studying human development and monitoring ecosystems [1][2][3][4][5][6].When multitemporal remote sensing images covering the same areas are able to be obtained, it is a matter of course to analyze where the landscape conditions have changed and what types the changes are.Therefore, change detection is one of the earliest and most important interpretation approaches for remote sensing technology [7][8][9].Dynamic analysis of landscape variation by change detection is extremely useful in many applications, such land-use/land-cover change, urban expansion, ecosystem monitoring, and resource management [10][11][12][13][14][15][16][17][18].
In the literatures, numerous studies have been proposed for change detection, which can be categorized into three classes [19]: (1) direct comparison, which calculates the difference between the original spectral features, such as image difference, image ratio, and change vector analysis [9,[20][21][22]; (2) image transformation, which obtains the change information by transforming the original data into new discriminant features, such as principle component analysis, multivariate alteration detection (MAD), and slow feature analysis (SFA) [23][24][25]; (3) classification based method, which analyzes the "from-to" change types by comparing the classification maps or classifying the multitemporal combinations, such as post-classification method and direct classification [26,27].
With the development of very high-resolution (VHR) remote sensing technology, multitemporal VHR images are much easier to obtain and are widely used in building change detection, disaster assessment, and detailed land-cover change monitoring [28][29][30][31][32][33].In order to take advantage of the abundant textual information and detailed spatial relationship [34][35][36], various studies have been proposed for multitemporal VHR data [32].We can group them into two categories: (1) object-oriented change detection, which utilizes the object as the process unit to improve the completeness and accuracy of the final result [19,37,38]; (2) change detection fusing spatial features, which takes into consideration the texture, shape, and topology features in the process of change analysis [28,39,40].
In most previous works for multitemporal VHR data, one specific change detection method was utilized to obtain the change intensity, measuring the probability to be changed.However, with the increase of spatial resolution, the high reflectance variability within individual objects in urban areas also increases, so that no single method can obtain a satisfactory performance.When fusing multiple change detection results, the uncertainty will be the main problem, which contains the inaccuracy of each result and the conflict between different decisions.
Dempster-Shafer Theory (D-S) is a decision theory effectively fusing multiple evidences from different sources [41].One important advantage of D-S theory is that it can provide explicit estimations of uncertainty between different evidences from different data sources [42,43].Thus, the fusion of change detection results by D-S theory aims to improve the reliability of decisions by taking advantage of complementary information while decreasing their uncertainty [44].
In this paper, we proposed an urban change detection method for VHR remote sensing images by fusing multiple change detection methods with D-S theory and segmentation objects.Firstly, the multitemporal VHR images were stacked for object segmentation.Since change vector analysis (CVA) [20], iteratively reweighted multivariate alteration detection (IRMAD) [24], and iterative slow feature analysis (ISFA) [23] are all classical unsupervised methods and perform well in urban change detection, they were then implemented to obtain three candidate change maps, respectively.Finally, these three change maps were fused based on the object map and D-S theory to obtain the final result.
The contribution of this paper is in proposing a framework of fusing multiple change detection methods by D-S theory considering the uncertainty of each binary candidate result.Although there have been numerous studies about change detection, it is hard to obtain a satisfactory result with a single method due to the complexity of the urban environment in VHR images.When combining multiple change detection results, few studies take the uncertainty, which mainly includes the inaccuracy of each candidate change map and the conflicts between different results, into consideration.In our paper, the proposed fusion method by D-S theory has the ability to explore the uncertainty among different results and to generate a more accurate change map.
The rest of this paper is organized as follows.Section 2 details the proposed change detection method.The experiments are addressed in Section 3. Section 4 shows the detailed discussion.Finally, Section 5 draws the conclusion.

Methodology
In this section, we elaborate how to fuse the results of different change detection methods with D-S theory to improve the performance for multitemporal VHR images.The whole procedure of the proposed method is shown in Figure 1.The main steps are as follows: (1) Obtain the object map by the segmentation of stacked multitemporal VHR images; (2) Implement three change detection methods (CVA, ISFA, and IRMAD) to get change intensity, and utilize OTSU thresholding method to obtain the candidate change maps; (3) Fuse the three candidate change maps by the object map and D-S theory to obtain the final object-oriented change map.

Segmentation
VHR remote sensing imagery can provide detailed and abundant spatial information.However, the high spectral variability within individual objects have challenged the traditional pixelwise image interpretation [19].Object-oriented image analysis (OOIA) is a good idea to deal with the problem of spectral variability and has the ability to keep the object completeness.
For change detection with multitemporal VHR images, it is important to maintain the temporal consistency of objects, which means the pixels within one object are all changed or unchanged.In this paper, we firstly stacked the multitemporal images and implemented multiresolution segmentation to obtain the object map [45].Multiresolution segmentation (MRS) is an effective segmentation algorithm, which forms homogenous objects by a bottom-up region merging.The software of eCognition Developer was utilized to employ MRS in this study.

Change Detection
Due to the complexity of multitemporal VHR images, no single change detection can obtain a very accurate result.Therefore, in this paper, we utilized three effective change detection methods and fused their results by D-S theory.We call them "component methods".The three component methods are CVA, IRMAD, and ISFA, since they are all effective unsupervised change detection methods.
CVA is a classical change detection method and has been the foundation of numerous researches [20].The change vector was obtained by differencing the spectral values of corresponding bands, and the change intensity was calculated with the Euclidean distance of the change vector as follows: where and indicate the spectral vectors and means the change intensity.IRMAD is an improved version of the MAD algorithm [24].It assigns high weights to the unchanged pixels during the iterations so as to reduce the negative effect of changed pixels in the feature space learning [24,46,47].After the convergence of IRMAD is achieved, the chi-squared distance of the transformed features is used as the change intensity [24].For IRMAD, the transformation vectors and should satisfy the following optimization objective: The change intensity of IRMAD is calculated by the chi-squared distance as

Segmentation
VHR remote sensing imagery can provide detailed and abundant spatial information.However, the high spectral variability within individual objects have challenged the traditional pixelwise image interpretation [19].Object-oriented image analysis (OOIA) is a good idea to deal with the problem of spectral variability and has the ability to keep the object completeness.
For change detection with multitemporal VHR images, it is important to maintain the temporal consistency of objects, which means the pixels within one object are all changed or unchanged.In this paper, we firstly stacked the multitemporal images and implemented multiresolution segmentation to obtain the object map [45].Multiresolution segmentation (MRS) is an effective segmentation algorithm, which forms homogenous objects by a bottom-up region merging.The software of eCognition Developer was utilized to employ MRS in this study.

Change Detection
Due to the complexity of multitemporal VHR images, no single change detection can obtain a very accurate result.Therefore, in this paper, we utilized three effective change detection methods and fused their results by D-S theory.We call them "component methods".The three component methods are CVA, IRMAD, and ISFA, since they are all effective unsupervised change detection methods.
CVA is a classical change detection method and has been the foundation of numerous researches [20].The change vector was obtained by differencing the spectral values of corresponding bands, and the change intensity was calculated with the Euclidean distance of the change vector as follows: where x and y indicate the spectral vectors and d means the change intensity.
IRMAD is an improved version of the MAD algorithm [24].It assigns high weights to the unchanged pixels during the iterations so as to reduce the negative effect of changed pixels in the feature space learning [24,46,47].After the convergence of IRMAD is achieved, the chi-squared distance of the transformed features is used as the change intensity [24].For IRMAD, the transformation vectors a and b should satisfy the following optimization objective: The change intensity of IRMAD is calculated by the chi-squared distance as where σ k indicates the variance of k th feature band.
ISFA is an iterative improvement of the SFA change detection method [23].SFA can extract the temporal invariant features from multitemporal images and reduce the spectral variance between unchanged landscapes [48].In this way, the real changes are anomalous from the invariant features and can be well separated in the feature differences [49].ISFA is proposed to consider more unchanged pixels in feature learning during the iterations by assigning high weights to unchanged pixels and almost zero weights to the changed pixels [23].For ISFA, the transformation vectors w should satisfy the following optimization objective: where the transformed features should be under the constraints of zero mean, unit variance, and decorrelation.The chi-squared distance are also used to calculate the change intensity [23]: where σ k also indicates the variance of k th SFA feature band.
After the change intensity were obtained, the automatic thresholding method was utilized to generate the binary change maps.OTSU was selected as the thresholding method due to its effectiveness and simplicity [23,48].In order to test the effectiveness of OTSU automatic thresholding algorithm, we also evaluated the results with different thresholding methods, including K-means two-class clustering [25] and expectation maximization (EM) thresholding [50].The maximum accuracies by traversing all possible thresholds are also shown as comparisons.

Fusion with D-S Theory
Due to the complexity of VHR images, it is evident that no single method can provide a consistent means of detecting landscape changes.In order to improve the performance and decrease the uncertainties, it is feasible to set a criterion combining complementarities of the three candidate change maps.To this end, D-S evidence theory was utilized in the proposed method.D-S theory incorporates the basic probability assignment function (BPAF) by fusing the probability of each evidence to measure the event probability [51,52].Assuming that there is a space of hypotheses, denoted as Θ, in change detection applications, Θ is the set of hypotheses about change/nonchange, and its power set is 2 Θ .Supposing A is a nonempty subset of 2 Θ , thus m(A) indicates the BPAF of subset A, representing the degree of belief.The BPAF m: 2 Θ → [0, 1], and it must follow the following constraints: Supposing that we have q independent evidences, m i (B i ) indicates the BPAF computed from the evidence i (1 ≤ i ≤ n, n ≥ 3) and B i ∈ 2 Θ , B i = ∅.Therefore, the computation of BPAF m(A) is shown as follows: which denotes the probability of A by fusing the probabilities of difference evidences.
In the change detection problem, the space of hypotheses Θ equals to {h 0 , h 1 }, where h 0 indicates nonchange and h 1 indicates change.Therefore, the three nonempty subsets of 2 Θ are {h 0 }, {h 1 }, and {h 0 , h 1 }, which means nonchange, change, and uncertainty.The three evidences are the binary change detection results from CVA, IRMAD, and ISFA.The BPAF for each evidence is computed considering the binary change maps and object maps.For each object j in VHR images, the BPAF of {h 0 }, {h 1 }, and {h 0 , h 1 } for the evidence i is defined as where N j U indicates the number of unchanged pixels in object j in evidence i, N j C indicates the number of changed pixels in object j in evidence i, N j T indicates the total number of pixels in object j, and p i (0 ≤ p i ≤ 1) measures the certainty weight of this evidence.If p i is very large, the BPAF of uncertainty {h 0 , h 1 } will be very small.In this way, the three candidate change maps are fused by considering the segmentation object information and uncertainty of each evidence.
The final decision of each object to be changed is assigned when the following rule is satisfied: where the m(A) (A is {h 0 }, {h 1 } or {h 0 , h 1 }) is calculated by ( 8) and ( 9).

First Dataset
The first experiment of change detection was carried out based on a pair of GaoFen-2 (GF-2) multispectral images, as shown in Figure 2. The multispectral images contained three visible spectral bands and one near-infrared band.The image size was 3000 × 3000, with the resolution of 4 m.Accurate coregistration was employed with ground control points (GCPs), and the residual mis-registration was less than 1 pixel.The study scene covered the Hanyang urban area of Wuhan city, Hubei province, China.Due to the rapid development of Wuhan, the study area showed obvious land-cover changes.
where indicates the number of unchanged pixels in object in evidence , indicates the number of changed pixels in object in evidence , indicates the total number of pixels in object , and 0 1 measures the certainty weight of this evidence.If is very large, the BPAF of uncertainty , will be very small.In this way, the three candidate change maps are fused by considering the segmentation object information and uncertainty of each evidence.
The final decision of each object to be changed is assigned when the following rule is satisfied: where the ( is , or , ) is calculated by ( 8) and ( 9).

First Dataset
The first experiment of change detection was carried out based on a pair of GaoFen-2 (GF-2) multispectral images, as shown in Figure 2. The multispectral images contained three visible spectral bands and one near-infrared band.The image size was 3000 3000, with the resolution of 4 m.Accurate coregistration was employed with ground control points (GCPs), and the residual misregistration was less than 1 pixel.The study scene covered the Hanyang urban area of Wuhan city, Hubei province, China.Due to the rapid development of Wuhan, the study area showed obvious land-cover changes.The reference map was selected by manual interpretation, which is shown in Figure 3.The changed pixels are labeled as red with the number of 106,851, and the unchanged pixels are shown The reference map was selected by manual interpretation, which is shown in Figure 3.The changed pixels are labeled as red with the number of 106,851, and the unchanged pixels are shown as labeled with the number of 1,269,986.The change detection was regarded as a two-class classification problem, where change and nonchange were the two classes.Then, with the selected change and nonchange samples, the confusion matrix for a two-class classification could be generated, and the evaluation values, including Kappa, overall accuracy (OA), detection rate (DR), false alarm rate (FAR), and F1-Score, were calculated.In the GF-2 dataset, the value range of the pixels was from 0-1000.If the pixels were overexposed, the spectral features would not follow the tendency of spectral differences caused by different observation environments at different times.These overexposed pixels would lead to pseudo changes that decreased change detection accuracy.Therefore, the pixels, the values of which were higher than 980 in both multitemporal images, were determined as overexposed and excluded from the calculation and evaluation.Besides, since the unchanged water areas in the image showed quite different spectral features, the water areas were also excluded from the calculation and evaluation by thresholding the NDWI features.The NDWI features of the multitemporal images were calculated separately, and 0.35 was used as a threshold to extract water areas in both images.The areas belonging to water in both images were masked.
The change intensity of CVA, IRMAD, and ISFA are shown in Figure 4a-c.After the OTSU automatic thresholding, the binary change maps are shown in Figure 4d-f, where white indicates change.It can be observed that CVA, IRMAD, and ISFA detected most obvious changes.The Kappa coefficients for the binary results of CVA, IRMAD, and ISFA were 0.8444, 0.8680, and 0.8656, respectively.In the GF-2 dataset, the value range of the pixels was from 0-1000.If the pixels were overexposed, the spectral features would not follow the tendency of spectral differences caused by different observation environments at different times.These overexposed pixels would lead to pseudo changes that decreased change detection accuracy.Therefore, the pixels, the values of which were higher than 980 in both multitemporal images, were determined as overexposed and excluded from the calculation and evaluation.Besides, since the unchanged water areas in the image showed quite different spectral features, the water areas were also excluded from the calculation and evaluation by thresholding the NDWI features.The NDWI features of the multitemporal images were calculated separately, and 0.35 was used as a threshold to extract water areas in both images.The areas belonging to water in both images were masked.
The change intensity of CVA, IRMAD, and ISFA are shown in Figure 4a-c  In the GF-2 dataset, the value range of the pixels was from 0-1000.If the pixels were overexposed, the spectral features would not follow the tendency of spectral differences caused by different observation environments at different times.These overexposed pixels would lead to pseudo changes that decreased change detection accuracy.Therefore, the pixels, the values of which were higher than 980 in both multitemporal images, were determined as overexposed and excluded from the calculation and evaluation.Besides, since the unchanged water areas in the image showed quite different spectral features, the water areas were also excluded from the calculation and evaluation by thresholding the NDWI features.The NDWI features of the multitemporal images were calculated separately, and 0.35 was used as a threshold to extract water areas in both images.The areas belonging to water in both images were masked.
The change intensity of CVA, IRMAD, and ISFA are shown in Figure 4a-c In order to evaluate the performances of different thresholding methods, the quantitative assessments are shown in Table 1.For quantitative assessment, Kappa coefficient and OA were used as comprehensive indicators, where the Kappa coefficient was more reliable due to the imbalance between changed reference samples and unchanged reference samples."DR" and "FAR" were the mean detection rate and false alarm rate, respectively, while "F1-Score" was a comprehensive evaluation used in the detection problem of computer vision [48].The best accuracies for indicators are highlighted in bold.It can be seen that, among the three automatic thresholding methods, OTSU obtained the highest accuracies, which were very close to the maximum accuracies.For the purpose of making the proposed adaptation, the binary results by OTSU automatic thresholding method were selected for the fusion of D-S theory.Combining Figure 4 and Table 1, it can be observed that CVA with OTSU can obtain a higher detection rate, while IRMAD and ISFA with OTSU can obtain a lower false alarm rate.This indicates that these three binary candidate maps can offer different information for change detection, while CVA provides more information for detecting changes, and IRMAD and ISFA provide more information for avoiding false alarms.The segmentation object map is shown in Figure 5a, with the scale of 50, the shape weight of 0.5, and the compactness weight of 0.3.With this object map, the final binary change map by D-S fusion is shown in Figure 5b, where the weight for CVA, IRMAD, and ISFA are 0.7, 0.1, and 0.1, respectively.In order to evaluate the performances of different thresholding methods, the quantitative assessments are shown in Table 1.For quantitative assessment, Kappa coefficient and OA were used as comprehensive indicators, where the Kappa coefficient was more reliable due to the imbalance between changed reference samples and unchanged reference samples."DR" and "FAR" were the mean detection rate and false alarm rate, respectively, while "F1-Score" was a comprehensive evaluation used in the detection problem of computer vision [48].The best accuracies for indicators are highlighted in bold.It can be seen that, among the three automatic thresholding methods, OTSU obtained the highest accuracies, which were very close to the maximum accuracies.For the purpose of making the proposed adaptation, the binary results by OTSU automatic thresholding method were selected for the fusion of D-S theory.Combining Figure 4 and Table 1, it can be observed that CVA with OTSU can obtain a higher detection rate, while IRMAD and ISFA with OTSU can obtain a lower false alarm rate.This indicates that these three binary candidate maps can offer different information for change detection, while CVA provides more information for detecting changes, and IRMAD and ISFA provide more information for avoiding false alarms.The segmentation object map is shown in Figure 5a, with the scale of 50, the shape weight of 0.5, and the compactness weight of 0.3.With this object map, the final binary change map by D-S fusion is shown in Figure 5b, where the weight p for CVA, IRMAD, and ISFA are 0.7, 0.1, and 0.1, respectively.Since the study region was very large, for a better visual comparison, the zoomed multispectral images and binary change maps are shown in Figure 6.Compared with the pixelwise results, the changed objects were more complete, and most false alarms were avoided in results of the proposed method in Figure 6c.This result took advantage of the three change detection methods and decreased the uncertainty by D-S theory.The object map was considered in the fusion to keep the homogeneity of landscape objects.In order to evaluate the performance of the proposed method, we chose the pixelwise methods and object-oriented methods for comparison.The comparison is shown in Table 2. "CVA", "IRMAD", and "ISFA" indicate the pixelwise change detection results shown in Figure 4.The automatic thresholding method utilized the OTSU method, since it showed the best performance among the Since the study region was very large, for a better visual comparison, the zoomed multispectral images and binary change maps are shown in Figure 6.Compared with the pixelwise results, the changed objects were more complete, and most false alarms were avoided in results of the proposed method in Figure 6c.This result took advantage of the three change detection methods and decreased the uncertainty by D-S theory.The object map was considered in the fusion to keep the homogeneity of landscape objects.In order to evaluate the performance of the proposed method, we chose the pixelwise methods and object-oriented methods for comparison.The comparison is shown in Table 2. "CVA", "IRMAD", and "ISFA" indicate the pixelwise change detection results shown in Figure 4.The automatic thresholding method utilized the OTSU method, since it showed the best performance among the In order to evaluate the performance of the proposed method, we chose the pixelwise methods and object-oriented methods for comparison.The comparison is shown in Table 2. "CVA", "IRMAD", and "ISFA" indicate the pixelwise change detection results shown in Figure 4.The automatic thresholding method utilized the OTSU method, since it showed the best performance among the automatic Remote Sens. 2018, 10, 980 9 of 18 thresholding methods."CVA_MajorVote", "IRMAD_MajorVote", and "ISFA_MajorVote" indicate the object-oriented results, where the pixelwise results were fused according to the segmentation object map in Figure 5a by major vote."MajorVote_Fusion" indicates the result that are obtained by major vote between the three object-oriented results.Finally, "D-S_Fusion" denotes the result by the proposed method.In Table 2, compared with the pixel-wise methods, the object-oriented methods show improvements in most cases.Especially for FAR, the objected-oriented methods can get obvious fewer false alarms than their pixelwise versions.It can prove that the object-oriented process has the ability to avoid the "salt and pepper" noise in change detection results.The "MajorVote_Fusion" by fusing the three object-oriented results did not show a better performance than its components, which is because the uncertainty between different results will negatively affect the fusion.Among all the results, the proposed method obtained the largest Kappa coefficient and overall accuracy.It showed the third largest detection rate with the fourth lowest false alarm rate, thus obtaining the largest F1-score.Therefore, Table 2 illustrates the effectiveness of the proposed fusion method.
The results with different segmentation parameter settings are evaluated in Figure 7.During the evaluation, the weight p in D-S fusion were fixed to be 0.9, 0.5, and 0.5.The results with different scales are evaluated in Figure 7a, when the compactness weight and shape weight are fixed to be 0.5.Figure 7a illustrates that the scale of 30 can obtain the best performance.Figure 7b shows the accuracies for different shape weights with different scales, when the compactness weight is fixed to be 0.5.When the scale is 60, the change of shape weight will lead to great differences.With other scales, the shape weight of 0.3 and 0.5 provide better performances.Figure 7c shows the different results with various compactness weights, when the shape weight is fixed to be 0.5.It illustrates that when the scale will affect the tendency of the performances.When the scale is 30 or 40, the curves are all above the others.It is worth noting that in almost all the cases, the accuracies of the proposed method are higher than 0.91, which outperforms the other methods.
Sens. 2018, 10, x FOR PEER REVIEW 9 of 18 automatic thresholding methods."CVA_MajorVote", "IRMAD_MajorVote", and "ISFA_MajorVote" indicate the object-oriented results, where the pixelwise results were fused according to the segmentation object map in Figure 5a by major vote."MajorVote_Fusion" indicates the result that are obtained by major vote between the three object-oriented results.Finally, "D-S_Fusion" denotes the result by the proposed method.In Table 2, compared with the pixel-wise methods, the object-oriented methods show improvements in most cases.Especially for FAR, the objected-oriented methods can get obvious fewer false alarms than their pixelwise versions.It can prove that the object-oriented process has the ability to avoid the "salt and pepper" noise in change detection results.The "MajorVote_Fusion" by fusing the three object-oriented results did not show a better performance than its components, which is because the uncertainty between different results will negatively affect the fusion.Among all the results, the proposed method obtained the largest Kappa coefficient and overall accuracy.It showed the third largest detection rate with the fourth lowest false alarm rate, thus obtaining the largest F1score.Therefore, Table 2 illustrates the effectiveness of the proposed fusion method.
The results with different segmentation parameter settings are evaluated in Figure 7.During the evaluation, the weight in D-S fusion were fixed to be 0.9, 0.5, and 0.5.The results with different scales are evaluated in Figure 7a, when the compactness weight and shape weight are fixed to be 0.5.Figure 7a illustrates that the scale of 30 can obtain the best performance.Figure 7b shows the accuracies for different shape weights with different scales, when the compactness weight is fixed to be 0.5.When the scale is 60, the change of shape weight will lead to great differences.With other scales, the shape weight of 0.3 and 0.5 provide better performances.Figure 7c shows the different results with various compactness weights, when the shape weight is fixed to be 0.5.It illustrates that when the scale will affect the tendency of the performances.When the scale is 30 or 40, the curves are all above the others.It is worth noting that in almost all the cases, the accuracies of the proposed method are higher than 0.91, which outperforms the other methods.Figure 8 shows the performances with different certainty weights p, when the segmentation is employed with the scale of 50, the shape weight of 0.5, and the compactness weight of 0.3.The three certainty weights involved in the D-S fusion were tested with traversal combinations from 0.1 to 0.9 by step of 0.2.The Kappa coefficient surface in Figure 8 shows the Kappa values generated by a certain weight p of ISFA with all the possible combinations of the other two weights.It can be observed that with the increase of p for CVA, the accuracy increases obviously.Therefore, the results of CVA plays an important role in the fusion.High weight of IRMAD can also lead to good performance when the weight of CVA is very high.The increase of p for ISFA leads to more robust performance.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 18 Figure 8 shows the performances with different certainty weights , when the segmentation is employed with the scale of 50, the shape weight of 0.5, and the compactness weight of 0.3.The three certainty weights involved in the D-S fusion were tested with traversal combinations from 0.1 to 0.9 by step of 0.2.The Kappa coefficient surface in Figure 8 shows the Kappa values generated by a certain weight of ISFA with all the possible combinations of the other two weights.It can be observed that with the increase of for CVA, the accuracy increases obviously.Therefore, the results of CVA plays an important role in the fusion.High weight of IRMAD can also lead to good performance when the weight of CVA is very high.The increase of for ISFA leads to more robust performance.

Second Dataset
The second experiment dataset was also acquired by GF-2 with the resolution of 4 m.The study scene covered the Hanyang area of Wuhan city.The multitemporal images with the size of 1000 1000 are shown in Figure 9.

Second Dataset
The second experiment dataset was also acquired by GF-2 with the resolution of 4 m.The study scene covered the Hanyang area of Wuhan city.The multitemporal images with the size of 1000 × 1000 are shown in Figure 9.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 18 Figure 8 shows the performances with different certainty weights , when the segmentation is employed with the scale of 50, the shape weight of 0.5, and the compactness weight of 0.3.The three certainty weights involved in the D-S fusion were tested with traversal combinations from 0.1 to 0.9 by step of 0.2.The Kappa coefficient surface in Figure 8 shows the Kappa values generated by a certain weight of ISFA with all the possible combinations of the other two weights.It can be observed that with the increase of for CVA, the accuracy increases obviously.Therefore, the results of CVA plays an important role in the fusion.High weight of IRMAD can also lead to good performance when the weight of CVA is very high.The increase of for ISFA leads to more robust performance.

Second Dataset
The second experiment dataset was also acquired by GF-2 with the resolution of 4 m.The study scene covered the Hanyang area of Wuhan city.The multitemporal images with the size of 1000 1000 are shown in Figure 9.The reference map of this study area is shown in Figure 10.The changed samples contained 15,975 pixels and the unchanged samples contained 484,546 pixels.The multitemporal images were preprocessed by the same approach as the first dataset.The water areas and overexposure areas were also masked from the calculation and evaluation.The reference map of this study area is shown in Figure 10.The changed samples contained 15,975 pixels and the unchanged samples contained 484,546 pixels.The multitemporal images were preprocessed by the same approach as the first dataset.The water areas and overexposure areas were also masked from the calculation and evaluation.The reference map of this study area is shown in Figure 10.The changed samples contained 15,975 pixels and the unchanged samples contained 484,546 pixels.The multitemporal images were preprocessed by the same approach as the first dataset.The water areas and overexposure areas were also masked from the calculation and evaluation.Table 3 shows the quantitative evaluation of different thresholding methods.The "Max" indicates the threshold corresponding to the maximum Kappa coefficient.It can be observed that OTSU can obtain the highest Kappa coefficient for CVA, and the second highest Kappa coefficients for IRMAD and ISFA.Especially, for IRMAD and ISFA, OTSU shows the highest OA.It worth noting that OTSU can obtain good performances for all methods, while K-means show extremely low accuracy for CVA.Therefore, in order to make the proposed method automatic and robust, OTSU was selected as the automatic thresholding method to produce the binary candidate change map. Figure 12a shows the segmentation object map, with the scale of 50, the shape weight of 0.5, and the compactness weight of 0.3.With this object map, Figure 12b shows the final binary change map by D-S fusion, where the weight for CVA, IRMAD and ISFA are 0.9, 0.1, and 0.3, respectively.Compared with the binary candidate change maps in Figure 11d-f, the fused change map shows more complete detected areas and most false alarms were removed.Table 3 shows the quantitative evaluation of different thresholding methods.The "Max" indicates the threshold corresponding to the maximum Kappa coefficient.It can be observed that OTSU can obtain the highest Kappa coefficient for CVA, and the second highest Kappa coefficients for IRMAD and ISFA.Especially, for IRMAD and ISFA, OTSU shows the highest OA.It worth noting that OTSU can obtain good performances for all methods, while K-means show extremely low accuracy for CVA.Therefore, in order to make the proposed method automatic and robust, OTSU was selected as the automatic thresholding method to produce the binary candidate change map. Figure 12a shows the segmentation object map, with the scale of 50, the shape weight of 0.5, and the compactness weight of 0.3.With this object map, Figure 12b shows the final binary change map by D-S fusion, where the weight p for CVA, IRMAD and ISFA are 0.9, 0.1, and 0.3, respectively.Compared with the binary candidate change maps in Figure 11d-f, the fused change map shows more complete detected areas and most false alarms were removed.Table 4 shows the quantitative evaluation of the proposed method.It illustrates that the proposed method can obtain the best accuracy with Kappa, OA, and F1-Score.The proposed method can obtain the second highest detection rate and the second lowest false alarm rate at the same time.Therefore, it shows the best performance in comprehensive evaluators.All the results with major vote got better performances than their original results, which indicates the effectiveness of the objectoriented process.However, it is worth noting that the major vote fusion with three methods cannot obtain an accurate result, since the uncertainty of each result is not taken into consideration.Figure 13 illustrates the evaluation of segmentation parameters.The weights for CVA, IRMAD, and ISFA are fixed to be 0.9, 0.1, and 0.3.It can be found that the scale of 30 is more suitable for this dataset in Figure 13a.Figure 13b indicates that when the scale is large and the shape weight increases, the accuracy will be decreased.In Figure 13c, when the scale is 30, the accuracy with various compactness weights will be higher and more stable.Table 4 shows the quantitative evaluation of the proposed method.It illustrates that the proposed method can obtain the best accuracy with Kappa, OA, and F1-Score.The proposed method can obtain the second highest detection rate and the second lowest false alarm rate at the same time.Therefore, it shows the best performance in comprehensive evaluators.All the results with major vote got better performances than their original results, which indicates the effectiveness of the object-oriented process.However, it is worth noting that the major vote fusion with three methods cannot obtain an accurate result, since the uncertainty of each result is not taken into consideration.Figure 13 illustrates the evaluation of segmentation parameters.The weights for CVA, IRMAD, and ISFA are fixed to be 0.9, 0.1, and 0.3.It can be found that the scale of 30 is more suitable for this dataset in Figure 13a.Figure 13b indicates that when the scale is large and the shape weight increases, the accuracy will be decreased.In Figure 13c, when the scale is 30, the accuracy with various compactness weights will be higher and more stable.Table 4 shows the quantitative evaluation of the proposed method.It illustrates that the proposed method can obtain the best accuracy with Kappa, OA, and F1-Score.The proposed method can obtain the second highest detection rate and the second lowest false alarm rate at the same time.Therefore, it shows the best performance in comprehensive evaluators.All the results with major vote got better performances than their original results, which indicates the effectiveness of the objectoriented process.However, it is worth noting that the major vote fusion with three methods cannot obtain an accurate result, since the uncertainty of each result is not taken into consideration.Figure 13 illustrates the evaluation of segmentation parameters.The weights for CVA, IRMAD, and ISFA are fixed to be 0.9, 0.1, and 0.3.It can be found that the scale of 30 is more suitable for this dataset in Figure 13a.Figure 13b indicates that when the scale is large and the shape weight increases, the accuracy will be decreased.In Figure 13c, when the scale is 30, the accuracy with various compactness weights will be higher and more stable.Figure 14 shows the accuracies with various weights for CVA, IRMAD, and ISFA.The parameter settings for segmentation are fixed as: the scale is 50, the shape weight is 0.5, and the compactness weight is 0.3.Figure 14 illustrates that the weight of CVA plays a key role in the fusion, since the high accuracies are mostly obtained with high weights of CVA. Figure 14 shows the accuracies with various weights for CVA, IRMAD, and ISFA.The parameter settings for segmentation are fixed as: the scale is 50, the shape weight is 0.5, and the compactness weight is 0.3.Figure 14 illustrates that the weight of CVA plays a key role in the fusion, since the high accuracies are mostly obtained with high weights of CVA.

Discussion
Since VHR multitemporal images show obvious variability within objects, it is hard to obtain a good performance by a single change detection method.When fusing different change detection results, the uncertainty is the main problem.Often, two types of uncertainty are associated with the change detection problem: one is due to the inaccuracy in each change detection result, and the other is due to the conflict between different results [41].In the proposed method, the inaccuracy is measured by the weight of in (9).The conflict is solved by calculating the BPAF in (8).D-S evidence theory utilizes the basic probability assignment function to combine the probability of each evidence and reduce the uncertainty.
The binary change maps for different methods can be obtained by thresholding the change intensities.In order to make the proposed method adaptive, automatic thresholding methods were employed.In Tables 1 and 3, K-means method, EM method, and OTSU method are evaluated for three change intensities and compared with the maximum accuracies.It can be found that the OTSU method outperformed the other two automatic thresholding methods in most cases and was comparatively robust.Therefore, the binary change maps by OTSU method were used for the candidate maps in the proposed D-S fusion.
The performance of the proposed method is evaluated in Tables 2 and 4. For each component change detection result, the corresponding object-oriented results show more accurate and complete change maps.It illustrates that the object-oriented process is an effective approach to deal with multitemporal VHR imagery analysis.Tables 2 and 4 also illustrates that the proposed method outperforms its component results (CVA, IRMAD, and ISFA).Even compared with the objectoriented fusion, the proposed D-S fusion shows higher accuracies in most cases.
According to the proposed method, the component methods with higher accuracy should be assigned as a larger , since is measured as the opposite of uncertainty.As shown in Figure 7,

Discussion
Since VHR multitemporal images show obvious variability within objects, it is hard to obtain a good performance by a single change detection method.When fusing different change detection results, the uncertainty is the main problem.Often, two types of uncertainty are associated with the change detection problem: one is due to the inaccuracy in each change detection result, and the other is due to the conflict between different results [41].In the proposed method, the inaccuracy is measured by the weight of p in (9).The conflict is solved by calculating the BPAF in (8).D-S evidence theory utilizes the basic probability assignment function to combine the probability of each evidence and reduce the uncertainty.
The binary change maps for different methods can be obtained by thresholding the change intensities.In order to make the proposed method adaptive, automatic thresholding methods were employed.In Tables 1 and 3, K-means method, EM method, and OTSU method are evaluated for three change intensities and compared with the maximum accuracies.It can be found that the OTSU method outperformed the other two automatic thresholding methods in most cases and was comparatively robust.Therefore, the binary change maps by OTSU method were used for the candidate maps in the proposed D-S fusion.
The performance of the proposed method is evaluated in Tables 2 and 4. For each component change detection result, the corresponding object-oriented results show more accurate and complete change maps.It illustrates that the object-oriented process is an effective approach to deal with multitemporal VHR imagery analysis.Tables 2 and 4 also illustrates that the proposed method outperforms its component results (CVA, IRMAD, and ISFA).Even compared with the object-oriented fusion, the proposed D-S fusion shows higher accuracies in most cases.
According to the proposed method, the component methods with higher accuracy should be assigned as a larger p, since p is measured as the opposite of uncertainty.As shown in Figure 7, most performances of the proposed method with the weight for CVA, IRMAD, and ISFA are 0.9, 0.5, and 0.5 are higher than the comparative methods.It can be seen in Figure 8 that with the high weight of ISFA, the proposed method can obtain good performance, which is mostly higher than 0.92. Figure 8 indicates that the high weight of ISFA can lead to more robust performance, and the result of CVA plays an important role in the D-S fusion.The accuracies of most results in Figure 8 are higher than 0.9, which are much higher than those of the candidate change maps.It indicates that the proposed method is also robust to various parameter settings.In practical applications, CVA, IRMAD, and ISFA may be more effective in different datasets [23,24,48].The determination of weights is one of the important works in the proposed method.Figures 13 and 14 also show the similar phenomena.
The main limitation of the proposed method is the parameter setting.There are two groups of parameters in the experiment.The first one is the group of segmentation parameters.The optimum setting of these parameters is affected by the resolution, object scale, and landscape distribution.Figures 7 and 13 indicate that when the weights of D-S fusion are suitable, most parameter settings can lead to a better performance than the simple major vote fusion.The second group includes the weight p for different evidences.Since the weight p measures the uncertainty of each evidence, in most cases, the more accurate component methods should be associated with a larger p.However, according to the experiments, the determinations of weights are not independent, and there are interactions among different weights.It can be observed that CVA shows higher detection rate, while IRMAD and ISFA show lower false alarms.Therefore, CVA provides more information about changed areas, and IRMAD and ISFA play a part in avoiding false alarms.Accordingly, the result of CVA is assigned as high weights to better detect changes, and one of the results of IRMAD and ISFA is assigned with comparatively low weights to provide complementary information.
It is hard to determine the weight p automatically.For this problem, we can introduce several effective ways.Firstly, the weight p can be determined according to interpretation and expert knowledge.With the interpretation of the binary candidate change map, the weight p can be estimated manually.Secondly, the weight p can be evaluated according to experiment in a smaller classical study area.Some statistical theories, such as orthogonal experiment design (OED), can be utilized to evaluate the sensitivity of each parameter and the interaction effects among all the parameters [52,53].Also, the best parameter setting can be estimated by OED and used for the whole experimental dataset.

Conclusions
Multitemporal VHR remote sensing images have shown great potentials in the applications of change detection.However, due to the spectral variability within landscape objects, it is hard to obtain an accurate result with a single method.When fusing multiple change detection results, the uncertainty between different evidences will be a big problem.Therefore, in this paper, we have proposed a novel change detection method by fusing multiple change detection results with D-S evidence theory.CVA, IRMAD, and ISFA with OTSU thresholding are the component methods providing candidate binary change maps that are regarded as the evidences.The BPAF for each evidence is calculated with the candidate change map, segmentation objects, and the weight for uncertainty.Finally, the decision for each object to be changed is made by fusing these three evidences with D-S theory.
The experiments with two VHR datasets were undertaken for evaluation.The results indicate that the proposed method can provide a more accurate change detection result.It shows the best performances in comprehensive indicators.The parameter setting for segmentation is determined according to the dataset.The weights of p are related to the uncertainties of their corresponding evidence and are also affected by the correlations among different evidences.
In our future work, we will mainly focus on the parameter setting, especially on the automatic determination of weights p. Their optimum values and interaction effects between different evidences should be considered and analyzed.The statistical method, such as orthogonal experimental design, will be utilized for further analysis [53].

Figure 1 .
Figure 1.Flowchart of the proposed method.

Figure 1 .
Figure 1.Flowchart of the proposed method.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 18 as labeled with the number of 1,269,986.The change detection was regarded as a two-class classification problem, where change and nonchange were the two classes.Then, with the selected change and nonchange samples, the confusion matrix for a two-class classification could be generated, and the evaluation values, including Kappa, overall accuracy (OA), detection rate (DR), false alarm rate (FAR), and F1-Score, were calculated.

Figure 3 .
Figure 3.The reference map where red indicates change and green indicates nonchange.

Figure 3 .
Figure 3.The reference map where red indicates change and green indicates nonchange.
. After the OTSU automatic thresholding, the binary change maps are shown in Figure 4d-f, where white indicates change.It can be observed that CVA, IRMAD, and ISFA detected most obvious changes.The Kappa coefficients for the binary results of CVA, IRMAD, and ISFA were 0.8444, 0.8680, and 0.8656, respectively.Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 18 as labeled with the number of 1,269,986.The change detection was regarded as a two-class classification problem, where change and nonchange were the two classes.Then, with the selected change and nonchange samples, the confusion matrix for a two-class classification could be generated, and the evaluation values, including Kappa, overall accuracy (OA), detection rate (DR), false alarm rate (FAR), and F1-Score, were calculated.

Figure 3 .
Figure 3.The reference map where red indicates change and green indicates nonchange.

Figure 5 .
Figure 5. (a) The segmentation object map; and (b) binary change map after D-S fusion.

Figure 5 .Figure 5 .
Figure 5. (a) The segmentation object map; and (b) binary change map after D-S fusion.

Figure 10 .
Figure 10.The reference map where red indicates change and green indicates nonchange.

Figure 10 .
Figure 10.The reference map where red indicates change and green indicates nonchange.

Figure 10 .
Figure 10.The reference map where red indicates change and green indicates nonchange.

Figure 12 .
Figure 12.(a) The segmentation object map; and (b) binary change map after D-S fusion.

Figure 12 .
Figure 12.(a) The segmentation object map; and (b) binary change map after D-S fusion.

Figure 12 .
Figure 12.(a) The segmentation object map; and (b) binary change map after D-S fusion.

Table 1 .
Quantitative evaluation of different thresholding methods in test dataset 1.

Table 1 .
Quantitative evaluation of different thresholding methods in test dataset 1.

Table 2 .
Quantitative evaluation of the proposed method and comparative methods in test dataset 1.

Table 2 .
Quantitative evaluation of the proposed method and comparative methods in test dataset 1.

Table 3 .
Quantitative evaluation of different thresholding methods in test dataset 2.

Table 3 .
Quantitative evaluation of different thresholding methods in test dataset 2.

Table 4 .
Quantitative evaluation of the proposed method and comparative methods in test dataset 2.

Table 4 .
Quantitative evaluation of the proposed method and comparative methods in test dataset 2.

Table 4 .
Quantitative evaluation of the proposed method and comparative methods in test dataset 2.