An Approach for Unsupervised Change Detection in Multitemporal VHR Images Acquired by Different Multispectral Sensors

Abstract

This paper proposes an approach for the detection of changes in multitemporal Very High Resolution (VHR) optical images acquired by different multispectral sensors. The proposed approach, which is inspired by a recent framework developed to support the design of change-detection systems for single-sensor VHR remote sensing images, addresses and integrates in the general approach a strategy to effectively deal with multisensor information, i.e., to perform change detection between VHR images acquired by different multispectral sensors on two dates. This is achieved by the definition of procedures for the homogenization of radiometric, spectral and geometric image properties. These procedures map images into a common feature space where the information acquired by different multispectral sensors becomes comparable across time. Although the approach is general, here we optimize it for the detection of changes in vegetation and urban areas by employing features based on linear transformations (Tasseled Caps and Orthogonal Equations), which are shown to be effective for representing the multisensor information in a homogeneous physical way irrespectively of the considered sensor. Experiments on multitemporal images acquired by different VHR satellite systems (i.e., QuickBird, WorldView-2 and GeoEye-1) confirm the effectiveness of the proposed approach.

1. Introduction

The use of Remote Sensing (RS) in the analysis and evaluation of environmental processes evolution is a valuable tool whose relevance has increased in conjunction with the use of digital image processing techniques. Due to the improvement of both acquisition sensor technology and data processing algorithms, it is possible to get an accurate and automatic identification and extraction of features for understanding the environmental changes occurring on the ground due to natural and anthropogenic interactions. The technological evolution resulted in the availability of multitemporal and multispectral satellite images with Very High spatial Resolution (VHR) acquired by passive sensors (e.g., QuickBird, WorldView-2, GeoEye, Pleiades). In the context of this paper, we define VHR images as those acquired with a metric to sub-metric spatial resolution in the panchromatic channel (i.e., spatial resolution lower than one meter). These images allow a detailed geometrical analysis when compared to medium or high spatial resolution data [1,2,3,4]. When considering VHR satellite systems, it is difficult to define Time Series (TS) made of images from one single sensor that satisfy the application temporal resolution constraints and show homogeneous acquisition conditions characteristics (e.g., similar light conditions, similar acquisition angle). This is mainly due to the satellite revisit period, the possible competing orders of different users on the satellite pointing, the limited life of a satellite mission, and weather conditions. However, since a considerable number of satellites have been launched in the last decades, the above-mentioned limitations can be mitigated by building TS where images acquired by different multispectral VHR sensors are considered.

The use of multisensor TS increases the probability of having sequences of multitemporal images with a proper time sampling but poses some challenges. In addition to real changes occurring on the ground, multisensor multitemporal images are affected by differences induced by the acquisition conditions (e.g., atmospheric conditions and acquisition system). On the one hand, some of the differences in atmospheric conditions (e.g., cloud cover), and the differences in acquisition system (e.g., view angle and shadows) affect single-sensor multitemporal image processing as well [5]. On the other hand, multisensor TS poses the big challenge of having system induced differences due to the type of sensor and the sensor acquisition modes. They result in spectral and geometric differences. Such differences make state-of-the-art change detection methods (here based on image comparison) less effective since they usually assume that multitemporal images are acquired by the same sensor and under similar acquisition conditions. Similar observations hold for the use of quantities like Digital Numbers (DN) that become non-comparable, thus their use may increase Change Detection (CD) error. Moreover, standard methods for converting DN in reflectance can be insufficient when using images acquired by different sensors [6].

An appropriate multitemporal image homogenization is therefore required to reduce both spectral and geometrical differences, and to ensure that differences in the multitemporal images can be associated to real changes occurring on the ground. In the literature, multitemporal image homogenization methods are available for medium and low spatial resolution data [7,8,9,10,11,12]. However, their adaptation to VHR images is still an open issue because of the higher within-class spectral variability induced by the very high spatial resolution, when compared to lower spatial resolution sensors. In [13] and [14], the authors made a first attempt to adapt existing methods to radiometric homogenization of a pair of VHR images (IKONOS and QuickBird). Spatial resolution differences were mitigated by resampling to the images with the lowest spatial resolution [13,14]. Other works transform DNs to physical information [6,15]. Pacifici et al. [6] showed the importance of working with physical quantities that are homogenous across time when using VHR optical images acquired by the same sensor. They pointed out the improvement achieved while increasing the level of information abstraction in the feature space, and its consequence for the final image homogenization. In [15], a first attempt to create a fully automated system for CD in VHR and multisensor images is presented. Attention is paid to the pre-processing steps where the goal is to standardize images across the sensors. To do so, many steps are carried out, among them the conversion of DN to Top Of Atmosphere (TOA) reflectance values, which allows the direct comparison of images acquired by different sensors. In other words, some work on adaptation of homogenization methods is present in the literature. However, and to the best of the authors’ knowledge, no formal approach exists that properly accounts for the complexity of data acquired by similar, yet different VHR multispectral sensors [15,16,17,18], nor does an approach exist that guides the user from data pre-processing until the final CD process.

In the literature, methods are available for multitemporal VHR optical images information extraction in the context of CD [5,19,20,21]. Both supervised and unsupervised CD techniques have been widely used in several RS applications (e.g., flood detection, damage assessment, environmental monitoring). The main drawback of supervised methods lies in the need for collecting and constructing ground reference data for the system-training phase. However, they lead in general to more accurate results. On the other hand, unsupervised techniques have the advantage of not requiring any ground reference information. This mainly results in a faster and more operational processing, but lower accuracy in part due to the use of non-homogenous features. Accordingly, in the literature a large attention has been devoted to improving the accuracy of unsupervised approaches [5,6,20,21,22,23]. Among CD techniques, the most widely used are the ones based on Principal Component Analysis (PCA), Iteratively Reweighted Multivariate Alteration Detection (IR-MAD) and Change Vector Analysis (CVA) [20,24,25,26].

CVA has been widely applied to the original spectral bands of multitemporal images from low to very high spatial resolution [22,23,27,28,29]. Other authors used CVA on the vegetation index or the Tasseled Cap (TC) feature space, but with low and medium spatial resolution images only, and/or for the detection of specific/single vegetation changes [28,29]. In the case of VHR images, original spectral bands are traditionally used, and the CD mainly focuses on the separation of change and no change classes, or the identification of a single type of change, without considering the nature (type) of the changes [5,30,31,32,33]. In low, medium, high and VHR cases, features are mostly selected according to the possible changes occurring on the ground. However, less attention is devoted to explicitly mitigating the differences induced by the acquisition system, both from the homogenization and feature selection perspectives [30,31,32,33], which in turn would lead to an improvement on the CD accuracy. In the specific case of IR-MAD, little attention is devoted to the possibility of using these features for CD in multisensor images [33]. Mostly qualitative visual analysis about the possibility to detect changes is conducted, but no quantitative change information extraction approach is provided [24]. In this context, the need arises for defining proper techniques and operational strategies for homogenization of multitemporal images acquired by different VHR multispectral sensors that mitigates differences in both atmospheric conditions and acquisition system parameters, thus reducing their impact on the multitemporal information extraction, and on the CD accuracy. The technique should also be able to extract changes automatically and distinguish among different types.

This paper presents an approach for unsupervised CD in multitemporal VHR images acquired by different multispectral sensors. The proposed approach takes inspiration from two existing frameworks presented in [5] and [34], but improves them by introducing guidelines and solutions for how to deal with multisensor VHR images in the context of CD. It exploits some of the concepts in [5] and extends and integrates them with a strategy for mitigating the effects of the non-homogeneous properties of multitemporal images acquired by different VHR multispectral satellite systems at both pre-processing and feature extraction/change detection level. In conclusion, this paper aims at illustrating the entire processing chain from image acquisition to CD map using multi-sensor images. Attention is focused on the possible sources of differences, particularly on the ones becoming more critical when multisensor images are used and may result in change detection errors.

The main steps of the proposed approach are: (i) mitigation of differences induced by the use of VHR multitemporal images acquired by different sensors; and (ii) detection of multiple changes occurring on the ground by means of high level physical features. The first step is conducted by defining homogenization procedures that address radiometric, spectral and geometrical differences. Thus, multitemporal multisensor images become more comparable (i.e., more homogeneous) across time. Homogenization is further improved by extracting physical features from multisensor images that allow for an effective multitemporal comparison across sensors at a given level of abstraction. Features are designed to detect multiple changes relevant to the user. Although the approach is general, here we concentrate on the selection of high-level physical features (i.e., features that account for the physical quantities) suitable to detect changes in vegetation and urban areas. Linear/Orthogonal transformation features are selected for the detection step. They are shown to be effective for representing the multisensor information in a coherent physical way versus the considered sensor. Other physical features should be selected depending on the application. The second step is conducted by means of CVA in spherical coordinates, where the changes are represented by magnitude and direction variables. Separation among the changes is carried out in an automatic way. In the case of the magnitude variable, changed and non-changed pixels are separated by means of a Bayesian decision rule [35], whereas along direction variables, an adaptation of the Two-Stage Multithreshold Otsu (TSMO) method [36] is used. Experiments carried out on bi-temporal VHR image pairs acquired by different sensors confirm the effectiveness of the proposed approach.

The remainder of this paper is structured as follows. Section 2 presents an overview and an analysis of the properties of a CD system for VHR RS images, the proposed approach for the mitigation of differences induced by the use of VHR multisensor multitemporal images and the CD process. Section 3 introduces the multisensor datasets, describes the design of experiments, and illustrates and discusses the experimental results. Finally, Section 4 draws the conclusions and future developments.

2. Methods

In this section, an overview of the conceptual framework for analyzing the properties of a CD system in the context of VHR RS images is presented. The theoretical structure of the problem in the context of multisensor VHR images is drawn. This is the model for the design of solutions for different specific CD problems. After that, the proposed approach for the mitigation of differences induced using VHR multisensor images is introduced. Last, the description of the CD approach for detecting changes of interest in bi-temporal images is provided.

2.1. Conceptual Framework: CD Systems for VHR Remote Sensing Images

When considering VHR RS images, the CD problem becomes complex due to the heterogeneous spatial and spectral characteristics. Further, standard CD techniques often do not account for the semantic meaning of changes of interest. Standard CD methods often assume that unchanged pixels have similar signatures on the two dates, whereas changed ones do not. Unfortunately, this assumption is often not satisfied when we consider multitemporal VHR images, since additional differences may appear due to spectral and spatial heterogeneity. This becomes more critical when multisensor images are considered [5,37]. In accordance with [5], a proper understanding and modeling of changes is fundamental for the development of effective techniques that can mitigate the intrinsic differences in multitemporal VHR data, and accurately detect multiple changes. Further, most of the current methods for CD in VHR images focus on: (i) handling images acquired by the same sensor [38,39], or (ii) detecting specific type of changes by using multisensor images (i.e., deforestation, burned areas, buildings detection) [30,31,32,33]. Thus, their applicability to CD in multisensor image pairs is limited.

In this paper, we aim at developing a CD approach for multisensor VHR optical images. For this purpose, the framework in [5] is used as a baseline. The general flowchart proposed in [5] consists of two main steps: (i) definition of the tree of radiometric changes; and (ii) detection of changes.

2.1.1. Definition of the Tree of Radiometric Changes

In this step, possible classes of radiometric changes are analyzed and their taxonomy is defined. The resulting tree of radiometric changes is specific for the considered CD problem. To this end, a categorization of the different possible radiometric changes that may be present in a multitemporal VHR dataset is required.

Figure 1 shows the tree that models radiometric changes for a generic CD problem in multitemporal multisensor VHR images. Let $\mathsf{\Omega }$ be the symbol representing the concept of change class. Thus, each specific type of change will be represented by the symbol $\mathsf{\Omega }$ and a sub-script that refers to the meaning of the type of change. According to [5], two main types of radiometric changes $\left({\mathsf{\Omega }}_{\mathrm{Rad}}\right)$ originate because of the complexity of VHR images: (i) changes due to acquisition conditions (${\mathsf{\Omega }}_{\mathrm{Acq}}$); and (ii) changes occurring on the ground $({\mathsf{\Omega }}_{\mathrm{Grd}})$. The former corresponds usually to changes of no interest for the end user. The latter includes the changes relevant from the user’s viewpoint. ${\mathsf{\Omega }}_{\mathrm{Acq}}$ changes are the ones associated to differences in atmospheric conditions $\left({\mathsf{\Omega }}_{\mathrm{Atm}}\right)$ and in the acquisition system $\left({\mathsf{\Omega }}_{\mathrm{Sys}}\right)$. The latter results in the appearance of undesired change patterns, differences in the geometry and in shadows. ${\mathsf{\Omega }}_{\mathrm{Sys}}$ changes, such as the ones due to the type of sensor (${\mathsf{\Omega }}_{\mathrm{Sen}}$), can be mitigated by working with proper homogeneous features.

Among ${\mathsf{\Omega }}_{\mathrm{Sen}}$, we find changes due to differences in the spectral resolution (${\mathsf{\Omega }}_{\mathrm{Spe}}$) and spatial resolution (${\mathsf{\Omega }}_{\mathrm{Spa}}$). ${\mathsf{\Omega }}_{\mathrm{Sys}}$ includes changes due to the sensor view angle (${\mathsf{\Omega }}_{\mathrm{Ang}}$) or the solar rays incidence angle (${\mathsf{\Omega }}_{\mathrm{Sun}}$). ${\mathsf{\Omega }}_{\mathrm{Ang}}$ can be related to the topography (${\mathsf{\Omega }}_{\mathrm{Top}}$) and the relief (${\mathsf{\Omega }}_{\mathrm{Rel}}$). For example, they may generate changes in shadows, whereas ${\mathsf{\Omega }}_{\mathrm{Sun}}$ refers to effects such the ones induced by seasonal variations of the solar ray incidence angle, which generates shadows differences that are not associated to changes on the ground. On the other hand, ${\mathsf{\Omega }}_{\mathrm{Grd}}$ changes can be divided into four main categories: natural disasters $({\mathsf{\Omega }}_{\mathrm{Dis}}$—e.g., earthquake damages), vegetation phenology $({\mathsf{\Omega }}_{\mathrm{Veg}}$—e.g., leaves lost during winter), environmental conditions $({\mathsf{\Omega }}_{\mathrm{Env}}$—e.g., variation in soil moisture levels) and anthropogenic activities $({\mathsf{\Omega }}_{\mathrm{Ant}}$—e.g., harvested crop fields, new buildings). The tree structure illustrated in Figure 1 has a general validity and can be used to model most of the CD problems [5]. However, depending on the specific CD problem, the tree can be adapted and optimized. On the one hand, some leaves/nodes might be irrelevant and thus can be removed. On the other hand, some leaves may require to be further specified in accordance with the specific study case (see Section 3 for an example of how to define the tree of radiometric changes for a concrete CD problem). To extract the changes of interest, it is necessary to select effective features and to count on prior knowledge about the problem.

2.1.2. Detection of Changes

In this step, the type of changes identified in the previous step are detected by selecting the strategy for the design of the CD method. According to [5], two possible strategies can be adopted based on: (i) direct extraction of the radiometric changes of interest; or (ii) detection by cancellation of non-interesting radiometric changes. Most of the current available CD methods for VHR RS images make use of the direct extraction strategy since their goal is to extract a specific type of change. Nevertheless, sometimes it is easier to detect the radiometric changes that are of no interest, and therefore to detect relevant changes by cancellation.

2.2. Proposed Approach to Unsupervised CD in VHR Multispectral Images Acquired by Different Sensors

In this sub-section, details on the proposed approach to solve the CD problem in VHR multisensor optical images are given. For handling the problem, we focus on two issues: (i) mitigation of multisensory induced changes ${\mathsf{\Omega }}_{\mathrm{Sys}}$ by the homogenization of multispectral data acquired by different VHR sensors; and (ii) detection of ${\mathsf{\Omega }}_{\mathrm{Grd}}$ changes by mitigation of residual ${\mathsf{\Omega }}_{\mathrm{Sys}}$ at the level of feature extraction. Figure 2 depicts the block scheme of the proposed approach. In order to accomplish the ${\mathsf{\Omega }}_{\mathrm{Sys}}$ mitigation, two main steps are considered: (1) spectral; and (2) geometric differences mitigation. ${\mathsf{\Omega }}_{\mathrm{Grd}}$ detection is accomplished in two steps: (3) multisensor feature extraction; and (4) change detection.

Let us consider two VHR optical images, $X_{1,a}$ and $X_{2,b}$, acquired by different sensors $S_1$ and $S_2$_, where $a=1,2,\dots ,A$ and $b=1,2,\dots ,B$ represent the multispectral bands for $S_1$ and $S_2$. The number $A$ and $B$ of spectral channels in $S_1$ and $S_2$ can be equal or different depending on the sensor. Given the use of different VHR sensors, $X_{1,a}$ and $X_{2,b}$ are likely to show different number of acquisition bands with slightly different bandwidth and spatial resolution, and/or different view angle. In other words, different spectral and geometrical properties. Let us assume that $X_{1,a}$ and $X_{2,b}$ sizes are $I_1\times J_1$ and $I_2\times J_2$, respectively, and that the images are acquired over the same geographical area at different times $t_1$ and $t_2$.

2.2.1. ${\mathsf{\Omega }}_{\mathrm{Sys}}$ Differences Mitigation

When dealing with multitemporal images acquired by different sensors $S_1$ and $S_2$, one of the critical issues is to identify and remove acquisition system induced changes (${\mathsf{\Omega }}_{\mathrm{Sys}}$). Handling the differences due to ${\mathsf{\Omega }}_{\mathrm{Sys}}$, contributes to mitigating issues on the left side of the tree of radiometric changes (Figure 1). In single sensor VHR images, ${\mathsf{\Omega }}_{\mathrm{Sys}}$ are mainly due to differences given by the sensor view angle (${\mathsf{\Omega }}_{\mathrm{Ang}}$), and are accentuated by the topography (${\mathsf{\Omega }}_{\mathrm{Top}}$) and relief (${\mathsf{\Omega }}_{\mathrm{Rel}}$). All of them contribute to the geometrical differences and result in radiometric distortions. When multisensor VHR images are considered, additional problems arise due to the type of sensor (${\mathsf{\Omega }}_{\mathrm{Sen}}$) and thus the differences in the spectral (${\mathsf{\Omega }}_{\mathrm{Spe}}$) and spatial (${\mathsf{\Omega }}_{\mathrm{Spa}}$) resolution. ${\mathsf{\Omega }}_{\mathrm{Spe}}$ can be mitigated by performing a radiometric normalization of the images, whereas ${\mathsf{\Omega }}_{\mathrm{Spa}}$ should be managed by means of geometric corrections, since they contribute to geometric differences.

To mitigate ${\mathsf{\Omega }}_{\mathrm{Spe}}$, two macro-categories of normalization methods exist in the literature: absolute and relative methods. The former involves the conversion of the DN values into the corresponding ground reflectance ones [6], while the latter performs image-to-image adaptation in the histogram feature space [9,10,11,12,13,14]. When data from two different sensors $S_1$ and $S_2$ are considered, the spectral information in the multisensor images is not comparable. Thus, absolute normalization is preferred with respect to the relative one, though not limited.

Absolute normalization estimates surface reflectance values providing information at physical level and mitigating spectral differences. The steps for spectral differences mitigation are shown in Figure 3, where $X_{1,a}$ and $X_{2,b}$ are converted from DN to At-Surface Reflectance (${\rho }^{ASR}$ [unitless]) images, $X_{1,a}^s$ and $X_{2,b}^s$ (where $s$ stands for spectrally corrected). Although this step might seem obvious in the CD process, several works can be found in the literature that use DN in the comparison step. Further, the mitigation of ${\mathsf{\Omega }}_{\mathrm{Spe}}$ becomes more critical when multisensor data are considered.

To get ${\rho }^{ASR}$, the digitalization process performed at the sensor during image formation must be inverted [6]. Parameters such as the mean exoatmospheric solar irradiance, solar zenith angle, Earth–Sun distance, radiance value and others are required. They can be retrieved from the metadata or from user guides, and are specific for each satellite. The resulting $X_{1,a}^s$ and $X_{2,b}^s$ have the same physical meaning. However, some differences cannot be mitigated. Thus, in addition to physical driven strategies, some data driven ones (feature extraction) are required (see next sub-section).

Satellites carrying on-board VHR sensors have the capability to acquire images with different view angles; this increases the probability of having multitemporal images with the required time resolution and cloud free on a given area. However, when multitemporal images are considered, differences in the acquisition view angle can induce misalignments because of the impact of topography (${\mathsf{\Omega }}_{\mathrm{Top}}$), small changes in relief of the terrain or the presence of buildings (${\mathsf{\Omega }}_{\mathrm{Rel}}$) [40,41,42,43]. Further, when $X_{1,a}$ and $X_{2,b}$ are acquired by different VHR sensors, small differences in the spatial resolution (${\mathsf{\Omega }}_{\mathrm{Spa}}$) are also expected. To achieve the geometric differences mitigation (step 2, Figure 2), the block scheme shown in Figure 4 is followed. $X_{1,a}^s$ and $X_{2,b}^s$ are the input images and $X_{1,a}^{s,g}$ and $X_{2,b}^{s,g}$ are the spectrally and geometrically mitigated ones ($g$ stands for geometrically corrected).

Geometric distortions due to the joint effect of topography, relief and sensor view angle, are mitigated by applying orthorectification with a high resolution Digital Elevation Model (DEM). In this way, most of the misalignments between multitemporal images due to ${\mathsf{\Omega }}_{\mathrm{Top}}$ and ${\mathsf{\Omega }}_{\mathrm{Rel}}$ are corrected. However, additional issues remain because of differences in spatial resolutions, thus co-registration should be applied so that pixels with the same coordinates in the images may be associated with the same area on the ground. This step is very critical since a poor co-registration may result in an unreliable final CD map [44]. On the other hand, it is important to clarify that neither orthorectification, nor co-registration solve the problems derived by the presence of vertical structures (i.e., the parallax problem). These types of changes are usually considered as sources of noise and are not of interest, thus they can be removed/mitigated by some feature extraction strategy applied during the CD stage (e.g., shadow detection and removal).

Pansharpening (PS) could be applied between orthorectification and co-registration as an optional step. It is meant to improve the spatial information by integrating the high spectral and low spatial resolution bands with the high spatial and low spectral resolution panchromatic band. Several PS methods exist in the literature, e.g., Intensity–Hue–Saturation (IHS), PCA, wavelets and Gram–Schmidt [45]. While co-registering, a resampling of the image with highest spatial resolution is performed in order to get the same common spatial resolution of the one with lower spatial resolution. The spatial resolution of VHR images is metric to sub-metric. Looking at existing satellite missions, the spatial resolution of VHR optical images ranges from the 0.3 m of WV-3 and -4 to the 1 m of Kompsat-2 in the panchromatic channel. Thus, when considering VHR multisensor image pairs the maximum resolution difference may rise up to about 0.7 m. Therefore, VHR images spatial resolutions are different yet similar and comparable. The outputs are multisensor VHR images, $X_{1,a}^{s,g}$ and $X_{2,b}^{s,g}$, showing the same spatial resolution and the same physical information, where ${\mathsf{\Omega }}_{\mathrm{Sys}}$ have been mitigated.

2.2.2. ${\mathsf{\Omega }}_{\mathrm{Grd}}$ Detection

Once the ${\mathsf{\Omega }}_{\mathrm{Sys}}$ have been mitigated, the proposed approach performs ${\mathsf{\Omega }}_{\mathrm{Grd}}$ detection. It extracts the changes of interest by selecting and extracting significant features for specific changes present in the study area. Standard CD approaches like Univariate Image Difference (UID) [46] and CVA [47] perform multitemporal comparison by means of the difference operator. The multispectral (or single spectral, in the case of UID) difference image $X_D$ is composed by Spectral Change Vectors (SCV).

The rationale behind the use of the difference operator is that unchanged samples show similar spectral signatures and thus result in SCVs with almost all zero components, whereas changed samples show SCVs with components far from zero. However, when multisensor images are considered, such an assumption is seldom satisfied, even after ${\mathsf{\Omega }}_{\mathrm{Sys}}$ mitigation. Thus, further homogenization is required in order to satisfy the a priori assumption for a successful employment of simple methods like UID and CVA. When multisensor images are considered, a proper feature space should be explicitly identified where pixel based comparison is meaningful. Here we propose the use of higher-level physical features, derived from the ${\rho }^{ASR}$ ones. While reducing residual sensor-induced differences in the unchanged areas, thus improving the homogenization level, a better highlighting of the changes of interest is achieved. In other words, working with higher-level physical quantities improves the level of abstraction while increasing the probability to detect the changes of interest and reducing the number of false alarms [6].

Since the proposed approach is general, any feature with high-level physical meaning can be used (e.g., radiometric indices). Further, since the approach is designed for VHR images, it may benefit from the use of spatial context information [5,48]. However here we are interested in understanding the performance of the approach for CD and thus of the multisensor homogenization procedure in mitigating the effects of sensor differences on the CD map. Accordingly, pixel-based features are considered such as radiometric indices. As an example, if changes due to vegetation phenology (${\mathsf{\Omega }}_{\mathrm{Veg}}$) are present, a radiometric index to detect vegetation can be used. In the case of natural disasters (${\mathsf{\Omega }}_{\mathrm{Dis}}$) in urban areas, a building index plus vegetation or soil index could be considered. Radiometric indices suitable to detect most of the relevant type of changes can be found in the literature.

Selection of a proper index becomes more and more complex when the ${\mathsf{\Omega }}_{\mathrm{Grd}}$ are coming from different sources. Since here we focus on vegetation and urban changes, we select features based on linear transformations such as TC or Orthogonal Equations (OrE) among the others [49,50]. TC features were designed originally as a linear transformation for the agricultural analysis on single date images [49], but they were further analyzed and used for CD analysis in multitemporal images. The literature works mainly used TCs in medium resolution and single sensor images (e.g., Landsat) [28,29]. Three main TCs have been studied (i.e., Brightness, Greenness and Wetness) because of their ability to detect and monitor soil content or transitions, vegetation and canopy moisture. Since TC is an invariant transformation in the physical feature space, its features are consistent between different scenes in a multitemporal time series [49] and therefore could be invariant between multisensor multitemporal images. Similarly, the OrE were derived following the TC philosophy and highlight information on crops, vegetation and soil. Additionally, OrE were designed as a substitution for the sensors which TC coefficients have not been derived in the literature, yet.

Given the above discussion, we use TC or OrE for the ${\mathsf{\Omega }}_{\mathrm{Grd}}$ detection process. These transformations compute a linear combination of the spectral bands. The number of features derived from TC is the same as the number of input features (in our case $A$ and $B$, respectively), but only 3 TCs are generally used for CD. In the case of OrE, only 3 features are derived. Equation (1) shows the general equation to calculate TC or OrE features ($F$), where $j$ represents each of the $F$ features and $C_{j,a}$ are the coefficients calculated for each $X_j^F$. The same equation applies for TC and OrE, though only the red, blue, green and NIR bands of the sensors are used for the latter. The following analysis applies for both TC and OrE:

$$X_j^F={{\displaystyle \sum }}_{a=1}^A{C}_{j,a}{X}_{1,a}^{s,g}.$$

(1)

Once physical level features have been extracted from different sensors, step 4 (Figure 2) applies CD. As mentioned above, CVA is employed. Since three features are considered, a 3-D representation is obtained [34]. Each SCV of $X_D^F$ is defined as in Equation (2):

$$X_D^F=X_2^F-X_1^F.$$

(2)

Each SCV component captures the multitemporal behavior of the corresponding feature (either TC or OrE). SCV components tend to assume small values when no change occurred, whereas if changes occurred, components assume large values (either negative or positive) depending on the type of change. This is true, even if single date features may have different ranges across each other since the difference operator (2) is applied to corresponding features computed at t₁ and t₂ showing similar range. To effectively perform CD in the multidimensional space defined by $X_D^F$, the information in $X_D^F$ vector is represented in spherical coordinates by computing its magnitude ($\rho$), azimuth angle ($\theta$), and elevation angle ($\phi$). The relationship between $X_D^F$ in Cartesian coordinates and Spherical coordinates is described by Equations (3)–(5):

$$\rho =\sqrt{X_{D,1}^F{}^2+X_{D,2}^F{}^2+X_{D,3}^F{}^2},$$

(3)

$$\theta =arctan\left(\frac{X_{D,2}^F}{X_{D,1}^F}\right),$$

(4)

$$\phi =arcos\left(\frac{X_{D,3}^F}{\rho }\right).$$

(5)

In the spherical representation, unchanged samples having small values in all the $X_D^F$ components assume small magnitude ($\rho$) values, whereas changed samples assume large values in one or more $X_D^F$ components thus showing a large magnitude ($\rho$) and a direction along $\theta$ and $\phi$ variables that depends on the ratios among values of the $X_D^F$ components (i.e., on the type of change). Therefore, the magnitude variable carries information about the presence/absence of changes, whereas the direction variables carry information about the possible type of changes [22,34,47]. According to these observations, a magnitude-direction domain (D) (Figure 5) can be defined as in Equation (6) that includes all SCVs:

$$D=\left\{\rho \in \left[0,{\rho }_{max}\right], 0\le \theta <2\pi and 0\le \phi <\pi \right\},$$

(6)

where ${\rho }_{max}$ is the maximum magnitude of $X_D^{TC}$.

Three subsets of D are of interest: (i) the sphere ($S_n$) that includes unchanged pixels, i.e., the ones with small magnitude values; and (ii) the spherical shell ($S_c$) that includes changed pixels, i.e., the ones with large magnitude values. $S_n$ and $S_c$ are complementary, the radius T of $S_n$ is the inner radius of $S_c$, and their union provides D. T separates changed from unchanged samples along the magnitude. Since the magnitude is a compressed 1-dimensional representation of the change problem, T is obtained as a trade-off among the effects of the types of change. Yet, the definition of T along the magnitude only is a simple and effective solution [22,34]; and (iii) truncated cone sectors ($S_k$) of changed pixels associated to preferred directions (${\theta }_k,{\phi }_k$) in $S_c$ (gray shaded truncated cone in Figure 5). Each preferred direction is associated to a specific type of change ${\mathsf{\Omega }}_{\mathrm{Grd}}$ (see Figure 1). The volume $S_k$ associated with the k-th change is defined as:

$$S_k=\left\{\rho ,\theta ,\phi :T\le \rho <{\rho }_{max},{\theta }_{k_1}\le \theta <{\theta }_{k_2},{\phi }_{k_1}\le \phi <{\phi }_{k_2}\right\}.$$

(7)

The upper and lower bounds ${\theta }_{k_1}$, ${\theta }_{k_2}$, ${\phi }_{k_1}$ and ${\phi }_{k_2}$, as well as T, can be calculated manually or automatically [51]. Once the angular thresholds have been estimated, the magnitude threshold T can be refined for $S_k$ to account for the behavior of each specific type of change [52]. Finally, the change detection map ($C{D}_{map}$) is built by including the following labels $\mathsf{\Omega }=\left\{{\omega }_n,{\mathsf{\Omega }}_{\mathrm{Grd}}\right\}$, with ${\mathsf{\Omega }}_{\mathrm{Grd}}=\left\{{\mathsf{\Omega }}_{\mathrm{Dis}},{\mathsf{\Omega }}_{\mathrm{Veg}},{\mathsf{\Omega }}_{\mathrm{Env}},\dots ,{\mathsf{\Omega }}_{\mathrm{Ant}}\right\}$, where ${\omega }_n$ refers to unchanged areas. As last step, non-relevant changes (e.g., misregistration, shadows), usually associated with the left side of the tree of radiometric changes shown in Figure 1, are removed from the $C{D}_{map}$.

3. Experimental Results and Discussion

3.1. Dataset Description and Design of Experiments

The proposed approach was validated over different areas located in the Trentino region in the north of Italy (Figure 6). These areas show interesting properties from the point of view of orographic conformation and environmental variety. Over a relatively small region it is possible to find: (i) flat regions including precious apple and vineyard fields, and urban, sub-urban and industrial areas with different density and structure; and (ii) hill and mountain environments with a variety of tree species.

Three bi-temporal data sets made up of different couple of images among QuickBird (QB), two WorldView-2 (WV-2) and one GeoEye-1 (GE-1) were constructed over the sample areas (yellow squares in Figure 6). The three datasets were selected such that different types of change are represented. Therefore, dataset 1 shows the transition from forest area to several types of vegetation, dataset 2 shows transitions among different phenological states of crop areas; and dataset 3 shows transitions from vegetation and bare soil (and vice versa) and changes in roofs and roads around an urban area. These datasets allow us to evaluate the complexity of working with multisensor VHR optical images. Details about the three datasets are provided in

Table 1. Datasets description. QB: QuickBird; WV: WorldView; GE: GeoEye.

	Dataset 1 and 2		Dataset 3
	${\mathit{t}}_1$	${\mathit{t}}_2$	${\mathit{t}}_1$	${\mathit{t}}_2$
Sensor	QB	WV-2	WV-2	GE-1
Acquisition date	July 2006	August 2010	May 2011	September 2011
Off-nadir angle	14.1°	19.3°	7.8°	14.4°

. The three satellites show some remarkable differences due to differences in the view angle, the spectral resolution, the number of bands and the spatial resolution.

The QB and GE-1 images have four multispectral bands, whereas WV-2 has eight. The spatial resolution of the QB image is 0.6 m for the panchromatic band and 2.4 m for multispectral bands, whereas WV-2 and GE-1 offer a higher spatial resolution in both panchromatic and multispectral bands with 0.5 m and 2 m, respectively.

Table 2. Main characteristics of QuickBird, WorldView-2 and GeoEye-1 optical sensors [53].

Satellite	QuickBird	WorldView-2	GeoEye-1
Bands (nm)	445–900 (pan)	450–800 (pan)	450–800 (pan)
		400–450 (coastal)
	450–520 (blue)	450–510 (blue)	450–510 (blue)
	520–600 (green)	510–580 (green)	510–580 (green)
		585–626 (yellow)
	630–690 (red)	630–690 (red)	655–690 (red)
		705–745 (red edge)
	760–900 (NIR)	770–895 (NIR 1)	780–920 (NIR)
		860–1040 (NIR 2)
Spatial Resolution (m)	0.61	0.46	0.41
	2.44	1.84	1.65

summarizes the characteristics of QB, WV-2 and GE-1 images from the spectral and spatial point of view. The spatial resolution differences imply that the sizes $I_1\times J_1$ and $I_2\times J_2$ of $X_1$ and $X_2$ images, respectively, are different despite they cover the same surface. The size of QB image in Figure 6 is 8674 × 6300 pixels, whereas the size of WV-2 and GE-1 images is 10297 × 7139 pixels. Thus, pixel-by-pixel comparison cannot be directly applied since the same pixel coordinates in the two images do not correspond to the same position on the ground. Concerning spectral resolution, we can observe that the four primary multi-spectral bands of QB, GE-1 and WV-2 are acquired over similar spectral ranges (e.g., red), but not fully identical (e.g., blue). Similar considerations hold for green and NIR bands.

In order to apply the proposed approach, we define the tree of radiometric changes (${\mathsf{\Omega }}_{\mathrm{Rad}}$) specific for the considered study areas (one single joint tree is provided for the three datasets), apply mitigation, extract suitable features and perform ${\mathsf{\Omega }}_{\mathrm{Grd}}$ detection. As a first step, we define the specific tree of radiometric changes $\left({\mathsf{\Omega }}_{\mathrm{Rad}}\right)$ for the considered problem by starting from the general tree given in Figure 1 [5]. The three datasets show changes due to acquisition conditions (${\mathsf{\Omega }}_{\mathrm{Acq}}$) and changes occurring on the ground (${\mathsf{\Omega }}_{\mathrm{Grd}}$). With regard to ${\mathsf{\Omega }}_{\mathrm{Grd}}$, there are changes in the phenological state of the vegetation $({\mathsf{\Omega }}_{\mathrm{Veg}})$ (e.g., radiometry of some crops yards, trees, small roads between crop yards (${\mathsf{\Omega }}_{\mathrm{Cro}}$), re-vegetation) and changes due to anthropogenic activities $({\mathsf{\Omega }}_{\mathrm{Ant}})$ (e.g., changes in road ${\mathsf{\Omega }}_{\mathrm{Roa}}$, deforestation ${\mathsf{\Omega }}_{\mathrm{Def}}$, roofs ${\mathsf{\Omega }}_{\mathrm{Bui}}$ and crop planting). It is important to clarify that even though the types of changes can be visually separated by photointerpretation, we do not have enough information to give a precise “from-to” label to them.

Concerning ${\mathsf{\Omega }}_{\mathrm{Acq}}$, both atmospheric conditions $\left({\mathsf{\Omega }}_{\mathrm{Atm}}\right)$ and acquisition system $\left({\mathsf{\Omega }}_{\mathrm{Sys}}\right)$ effects are present. ${\mathsf{\Omega }}_{\mathrm{Atm}}$ is mitigated by means of atmospheric corrections [6]. ${\mathsf{\Omega }}_{\mathrm{Sys}}$ is related to the type of sensor (${\mathsf{\Omega }}_{\mathrm{Sen}}$) and to the sensor view angle (${\mathsf{\Omega }}_{\mathrm{Ang}}$). The ${\mathsf{\Omega }}_{\mathrm{Ang}}$ changes generate small differences in the appearance of objects, leading to geometric distortions, thus to residual misregistration even after proper alignment of images, and to spectral differences when tall buildings are present. We can also see some differences in shadows, which become more critical when high buildings, structures or reliefs are present. ${\mathsf{\Omega }}_{\mathrm{Ang}}$ and ${\mathsf{\Omega }}_{\mathrm{Sen}}$ changes are non-relevant from the application viewpoint. Therefore, they are explicitly handled before proceeding to detect ${\mathsf{\Omega }}_{\mathrm{Grd}}$ and tuning the final CD map. ${\mathsf{\Omega }}_{\mathrm{Sys}}$ such as the ones due to sensor acquisition mode (${\mathsf{\Omega }}_{\mathrm{Mod}}$) are not considered since we are working with passive sensors. ${\mathsf{\Omega }}_{\mathrm{Grd}}$ like the ones due to natural disasters (${\mathsf{\Omega }}_{\mathrm{Dis}}$) or environmental conditions (${\mathsf{\Omega }}_{\mathrm{Env}}$) are ignored, since such events did not occur in the considered study area. According to this analysis and to the general taxonomy in Section 2, the tree of radiometric changes for the considered problem becomes the one in Figure 7.

Once the radiometric tree was defined, we performed spectral and geometric differences mitigation. All images were provided by DigitalGlobe Foundation in the context of the “MS-TS—Analysis of Multisensor VHR image Time Series” project [54]. Conversion from DNs to At-Surface Reflectance (ASR) was conducted before delivery by means of the Atmospheric Compensation (AComp) algorithm [6,55,56]. Given the orography of the study area and the possible distortions, we applied orthorectification by using a DEM obtained from LiDAR data [57]. Further distortions appear in dataset 1, since it is located in mountain area, and dataset 3 because of the presence of buildings. Additional pixel-to-pixel problems are also observed due to ${\mathsf{\Omega }}_{\mathrm{Ang}}$, and co-registration should be applied.

In order to achieve a better co-registration, PS was applied by means of the Gram–Schmidt method. Here ENVI software package was employed [58]. After PS, the spatial resolution for QB is 0.6 m, and 0.5 m for GE-1 and WV-2 multispectral bands. Co-registration was carried out over the QB–WV-2 and WV-2–GE-1 pairs, covering the whole study area in Figure 6, by using a polynomial function of second order. For the QB 2006 and WV-2 2010 couple, 79 uniformly distributed Ground Control Points (GCP) were selected, whereas 68 uniformly distributed GCP were selected for the WV-2 2011 and GE-1 2011 couple. The WV-2 2010 image was resampled during co-registration. Resampling was performed by means of the nearest neighbor interpolation. Figure 8 shows the pansharpened multisensor VHR QB, GE-1 and WV-2 images after applying spectral and geometric mitigation in the first and second column, respectively.

Datasets 1 and 2 show a common spatial resolution of 0.6 m and a size of 640 × 640 pixels, whereas dataset 3 shows a common spatial resolution of 0.5 m and a size of 1800 × 1800 pixels. Datasets 1, 2 and 3 appear in first, second and third row of Figure 8, respectively. In order to perform qualitative and quantitative analysis, a reference map for datasets 1 and 2 was defined by photointerpretation and exploiting prior knowledge on the scene as no ground truth was available (see Figure 8c,f, showing 332,414 and 280,928 unchanged pixels (white color) and 77,186 and 128,672 changed pixels (different colors), respectively). For dataset 3, considering the extent of the area and the fact that we have no complete knowledge of the changes occurring on the ground, it was not possible to derive a complete reference map. Thus, quantitative analysis was based on 62808 pixels marked as changed, and 6263 as unchanged, selected by photointerpretation. For comparison purposes, a false color composition of the two acquisitions is provided, green and fuchsia areas represent changes (Figure 8i). Changed pixels in the reference map include ${\mathsf{\Omega }}_{\mathrm{Grd}}$, only. For dataset 1, changes from (i) forest to bare soil; (ii) forest to grass; (iii) base soil to grass and (iv) bare soil to some road are identified. For dataset 2, changes from (i) dense vegetation to sparse or light vegetation and vice versa; and (ii) bare soil to vegetation are present. In addition, for dataset 3, changes are from (i) bare soil to vegetation, both dense and sparse; (ii) one to another color of the roofs; and (iii) old to new roads.

Once spectral and geometric mitigation was achieved, mitigation of residual ${\mathsf{\Omega }}_{\mathrm{Sys}}$ was performed at the level of feature extraction. The selection of the features is therefore designed to detect the residual ${\mathsf{\Omega }}_{\mathrm{Sys}}$ and the ${\mathsf{\Omega }}_{\mathrm{Grd}}$. According to the tree of radiometric changes (Figure 7), residual ${\mathsf{\Omega }}_{\mathrm{Sys}}$ might be related to ${\mathsf{\Omega }}_{\mathrm{Ang}}$, resulting in possible shadows and/or registration noise, whereas ${\mathsf{\Omega }}_{\mathrm{Grd}}$ includes three type of changes: ${\mathsf{\Omega }}_{\mathrm{Cro}},{\mathsf{\Omega }}_{\mathrm{Roa}}$ and ${\mathsf{\Omega }}_{\mathrm{Def}}$. Residual ${\mathsf{\Omega }}_{\mathrm{Sys}}$ due to shadows, due to vegetation or buildings, were detected by means of the method in [59], whereas ${\mathsf{\Omega }}_{\mathrm{Sys}}$ due to registration noise are negligible. ${\mathsf{\Omega }}_{\mathrm{Grd}}$ were detected by employing TC and OrE features.

Three main TCs (i.e., Brightness, Greenness and Wetness) and three OrE (Crop mark, Vegetation and Soil) have been studied because of their sensibility to soil content or transitions from-to soil, vegetation, canopy moisture and anthropogenic activities have been studied. Thus, we expect them to properly detect the different changes in the study area, with the exception of some transitions between green areas that do not show up in TC features. These are the cases of datasets 1 and 3, where transitions from forest to grass and crop to grass are misdetected. This is due to the fact that the difference is more in texture rather than in TC features. In order to evaluate and compare the proposed approach, experiments carried on features such as ASR, are also conducted. The detection of changes is obtained by applying CVA to a 3D TC and OrE feature space. However, other features such as IR-MAD [24] could be also used, despite no approaches for automatic detection with these features in multisensor VHR images exist yet. The drawback with IR-MAD is that it requires an end-user interaction to select the most specific components that represent the specific change of interest and to separate among changed and un-changed samples. This is time consuming and makes the approach not fully automatic.

Two experiments were designed: (i) experiment 1 (exp. 1) applies CVA to the transformed ASR bands features; and (ii) experiment 2 (exp. 2) applies CVA to TC (for datasets 1 and 2) or OrE (for dataset 3 which has a large presence or urban areas) features. In exp. 1, the first two selected bands are the Near-IR (NIR) and Red (R) given their high spectral sensibility in the analysis of vegetation and anthropogenic activities. The R bands of QB, GE-1 and WV-2 have quite similar spectral range, whereas the NIR ones do not (see Table 2). Therefore, between NIR1 and NIR2 WV-2 bands, NIR1 (770–995 nm) was selected given that its spectral range better matches to the QB NIR (760–900 nm) and GE-1 NIR (757–853 nm) band spectral range. Another ASR feature to be selected could be the Green or Blue band. Empirical experiments showed a slightly improvement in the final CD accuracy while using Green band instead of Blue one. Therefore, Green band was selected as third feature. In exp. 2, for datasets 1 and 2, TC features were selected based on: (i) the maximum number of TC features that can be derived for each specific sensor, (ii) the radiometric tree of changes; and (iii) the possibility to compare between multisensor TC features. Thus, 4 TC features were derived, bounded by QB properties. In accordance with the radiometric tree, features should be selected that are able to highlight ${\mathsf{\Omega }}_{\mathrm{Veg}}$ and ${\mathsf{\Omega }}_{\mathrm{Ant}}$. Based on the level of comparison between the multisensor TC coefficients, and according to the state-of-the-art, only the first three TC features of each sensor show similar physical meaning. For dataset 3, only three OrE features exist in the literature and are derived by means of the 4 main spectral bands of the sensors.

TC and OrE features were derived directly from the spectrally mitigated data and by using the coefficients in

Table 3. Tasseled Cap (TC) coefficients for QuickBird Digital Numbers (DN) values [60].

Bands	Brightness (TC1)	Greenness (TC2)	Wetness (TC3)
B1	0.319	−0.121	0.652
B2	0.542	−0.331	0.375
B3	0.490	−0.517	−0.639
B4	0.604	0.780	−0.163

,

Table 4. TC coefficients for WorldView-2 TOA values [61].

Bands	Brightness (TC1)	Greenness (TC2)	Wetness and shadows (TC3)
B1	−0.060	−0.140	−0.271
B2	0.012	−0.206	−0.316
B3	0.126	−0.216	−0.317
B4	0.313	−0.314	−0.243
B5	0.412	−0.411	−0.256
B6	0.483	0.096	−0.097
B7	−0.161	0.601	−0.743
B8	0.673	0.504	0.202

and

Table 5. Orthogonal coefficients for WorldView-2 and GeoEye-1 [50].

Orthogonal Component		Blue	Green	Red	NIR
WV-2	Crop Mark	−0.38	−0.71	0.20	−0.56
	Vegetation	−0.37	−0.39	−0.67	0.52
	Soil	0.09	0.27	−0.71	−0.65
GE-1	Crop Mark	−0.39	−0.73	0.17	−0.54
	Vegetation	−0.35	−0.37	−0.68	0.54
	Soil	0.08	0.27	−0.71	−0.65

. For TC, only coefficients corresponding to the first three TC feature are present. Here the set of coefficients in [60] was applied to the QB image. The coefficients are derived from the DN feature space (Table 3). There are no TC coefficients derived from TOA values for QB images in the literature. Thus, we applied the QB TC coefficients over the QB DN features, and compared the derived TC features as a higher level primitive. For the WV-2 image, coefficients are applied as given in [61], but to the ASR features instead of TOA ones, as originally derived (Table 4). The OrE features were derived by means of the coefficients shown in Table 5 and as per Equation (1).

3.2. Results and Discussion

In order to assess the effectiveness of the ${\mathsf{\Omega }}_{\mathrm{Sys}}$ mitigation and the ${\mathsf{\Omega }}_{\mathrm{Grd}}$ detection approach based on higher-level physical features, CVA was applied by considering the 3D feature space defined above and by means of Equations (2)–(5). We first extracted all the areas that correspond to radiometric changes in the image, by thresholding the magnitude variable. The selection of a threshold T over $\rho$ showed to be a simple and fast way to separate among changed and non-changed pixels. A good separation, among the three dimensions of CVA was guaranteed in average. T was automatically selected by applying a Bayesian decision rule [35] and by following the implementation presented in [62]. In [62], the statistical distribution of the magnitude as a mixture model representing the classes of unchanged and changed pixels is approximated by a Gaussian mixture. Parameters for the statistical distribution are derived by means of the EM algorithm and decision of the threshold is made using a Bayesian minimum cost criterion [46]. The T values for each of the datasets in the two experiments are shown in

Table 6. Magnitude ($T$) threshold values for the three datasets in Exp. 1 and 2.

Dataset	Exp.	$\mathit{T}$
1	1	0.060
	2	0.080
2	1	0.025
	2	0.030
3	1	0.090
	2	0.060

.

Figure 9 shows the multispectral difference image 3D histogram for the dataset 1, where bigger circles represent higher density (corresponding to the un-changed samples) and small circles represent lower densities (related to the different types of changes). In Figure 9a,c, it is possible to see how $X_D^{ASR}$ and $X_D^{TC}$ are distributed in a coplanar and sparse way, respectively. This kind of distribution leads to an easier visual interpretation of the $X_D^{TC}$, if compared to $X_D^{ASR}$. In fact, from the 3D histograms in Figure 9, we can see that the coplanar distribution of ASR features does not allow good separation between changes of interest and changes of no interest. A similar behavior is observed between the ASR and OrE features. Right column of Figure 9, presents the changed samples after removing the unchanged ones. Different clusters fit to different $S_k$ sectors and are associated with different changes. As we move from $X_D^{ASR}$ to $X_D^{TC}$, or $X_D^{ASR}$ to $X_D^{OrE}$, it becomes evident how different clusters locate around preferred directions and how the number of clusters increases, making their detection and separation more effective. Separation among different clusters can be performed by automatic or manual methods.

The last step builds the multiple $C{D}_{map}$ by means of the extraction by cancellation strategy. To this end, selection of each of the clusters in the changed region was conducted automatically by means of an adaptation of the TSMO proposed in [36]. The TSMO method yields the same set of thresholds as the conventional Otsu method, but it greatly decreases the required computation time. This is achieved by introducing an intermediate step based on histogram valley estimations, which is used to automatically determine the appropriate number of thresholds for an image. Final multiple CD maps were built by cancelling the remaining ${\mathsf{\Omega }}_{\mathrm{Sys}}$ and the ${\mathsf{\Omega }}_{\mathrm{Grd}}$ that are out of interest (e.g., cars in road and parking areas, Figure 10). As expected, some vegetation changes that affect more texture rather than TC features are misdetected in datasets 1 and 3. In other words, the selected higher-level physical features are not optimized for those changes, but different higher-level physical features could be selected to properly model most of the types of changes.

In order to perform a qualitative analysis, a comparison of the CD maps with their reference maps (for datasets 1 and 2) and the set of points collected by photointerpretation (dataset 3) was carried out. The comparison pointed out the improvement achieved when working with higher-level physical quantities, specifically for transitions from and to bare soil and different types of crops. This is confirmed when we analyze dataset 1 where changes are mainly from-to vegetation and bare soil. TC features outperform the results obtained when using ASR features (

Table 7. Confusion matrices for dataset 1 in Exp. 1 and Exp. 2.

Exp.	Changes Found	Actual Changes
		C1	C2	C3	C4	C5	Reliability
1	C1	39775	596	91	0	2129	93.38%
	C2	2206	9500	1501	44	28428	22.79%
	C3	9	5	2427	128	2077	52.24%
	C4	0	0	0	0	0	0.00%
	C5	3765	13090	2718	1327	299780	93.48%
	Accuracy	86.93%	40.96%	36.02%	0.00%	90.18%
Overall Accuracy		85.81%
2	C1	36558	250	92	0	910	96.69%
	C2	2412	8153	1355	94	6571	43.87%
	C3	0	2	1464	42	1152	55.03%
	C4	0	0	0	0	0	0.00%
	C5	6785	14786	3829	1363	323781	92.36%
	Accuracy	79.9%	35.15%	21.72%	0.00%	97.40%
Overall Accuracy		90.32%

). As expected, both TC and ASR features have similar problems to identify the change from forest to vegetation. An analogous situation occurs on dataset 2, with changes from-to different types of crops. In this case, TC outperforms ASR being able to separate among changes C2 and C3 (which cannot be discriminated when using ASR features—see Figure 10d–f).

For datasets 1 and 2, the major improvement is related to the decrease of False Alarms (FA) (see Table 7 and

Table 8. Confusion matrices for dataset 2 in Exp. 1 and Exp. 2.

Exp.	Changes Found	Actual Changes
		C1	C2	C3	C4	C5	C6	Reliability
1	C1	5381	6	1	10387	12286	16408	12.10%
	C2	142	40242	19330	3560	4	47406	36.35%
	C3	0	0	0	0	0	0	0.00%
	C4	18	0	50	9472	10	2688	77.40%
	C5	0	0	0	0	0	0	0.00%
	C6	1198	7973	3643	13800	1169	214426	88.53%
	Accuracy	79.85%	83.45%	0.00%	25.48%	0.00%	76.33%
Overall Accuracy		65.80%
2	C1	3269	8	3	10447	12423	12197	8.52%
	C2	66	37971	291	0	22	29301	56.13%
	C3	116	773	19534	3713	2	15080	49.81%
	C4	59	3	196	15476	51	4156	77.61%
	C5	0	0	0	0	0	0	0.00%
	C6	3229	9466	3000	7583	971	220194	90.08%
	Accuracy	48.51%	78.74%	84.84%	41.58%	0.00%	78.38%
Overall Accuracy		72.37%

). In dataset 1, the FA correspond to the main road passing through the area and to some of the remaining shadows generated by the tree lines. Even though an index was used to remove the shadows, a small percentage of them remained. In dataset 2, the FA correspond mainly to the linear structures like roads in between the different crops. Moreover, it is possible to observe improvements in terms of detection and separation of the different types of changes. The number of Missed Alarms (MA) decreased as well when working with TC, leading to a better detection of changes, especially when a higher number of changes exist, such as the case of dataset 2. From the quantitative viewpoint, the reduction in both FA and MA reflects in the Overall Accuracy (OA) that increases of about five percent and seven percent for dataset 1 and 2, respectively (see Table 7 and Table 8).

In the case of dataset 3, ASR features poorly separates among C3 and C4 classes that correspond to transitions from bare soil to sparse and dense vegetation, respectively. The number of MA for the class C3 by ASR features is clearly larger than for OrE ones (see

Table 9. Confusion matrices for dataset 3 in Exp. 1 and Exp. 2.

Exp.	Changes Found	Actual Changes
		C1	C2	C3	C4	C5	C6	C7	C8	C9	Reliability
1	C1	905	0	0	0	0	0	0	0	0	100%
	C2	0	206	0	0	0	0	0	0	0	100%
	C3	0	0	1642	334	0	0	0	0	11	82.64%
	C4	0	0	27	450	0	0	16	0	0	91.28%
	C5	0	0	0	0	240	0	0	0	0	100%
	C6	0	0	0	0	0	400	0	0	12	97.09%
	C7	3	0	0	68	0	0	195	0	0	73.31%
	C8	0	0	0	0	0	0	0	0	0	0.0%
	C9	725	0	203	392	0	198	89	213	62785	97.18%
	Accuracy	55.42%	100%	87.71%	36.17%	100%	66.89%	65.00%	0.00%	99.96%
Overall Accuracy		96.75%
2	C1	1294	0	0	0	0	0	0	0	91	93.43%
	C2	0	206	0	0	0	0	0	0	0	100%
	C3	0	0	1792	134	0	0	0	0	22	91.99%
	C4	0	0	27	1052	0	0	16	0	0	96.07%
	C5	0	0	0	0	240	0	0	0	0	100%
	C6	0	0	0	0	0	582	0	0	143	80.27%
	C7	0	0	0	0	0	0	158	0	205	43.53%
	C8	0	0	0	0	0	0	0	66	0	100%
	C9	339	0	26	58	0	16	110	147	62347	98.90%
	Accuracy	79.24%	100%	97.13%	84.57%	100%	97.32%	55.63%	30.98%	99.27%
Overall Accuracy		98.07%

). Moreover, OrE features are able to better detect the classes C2 and C7 than ASR features. C2 and C7 correspond to changes occurring on building roofs and roads renewed in the studied period. Finally, ASR features do not detect class C8 (change in roof color of a building), whereas OrE features do. Other transitions from-to vegetation and bare soil can be seen in the change classes with less MA and FA when the OrE features are used. It is worth noting that changes due to moving cars on roads or in parking areas were considered as non-relevant for this study and thus neglected. From the quantitative perspective, the OA obtained by OrE features increased the overall accuracy of about 2% over that of ASR features. This improvement can be considered relevant given the complexity of the scene with the presence of more types of changes. Proper higher-level physical features, such as some radiometric indexes or texture features, may provide better results. Table 7, Table 8 and Table 9 show the confusion matrices obtained for each dataset in the two experiments.

4. Conclusions and Future Developments

In this paper, an approach for CD in VHR multispectral multisensor optical images has been presented. The proposed approach aims at defining and illustrating a data flow for effectively handling differences due to acquisition sensors. It is based on a general framework for the design of CD systems for VHR multitemporal images presented in [5]. In order to deal with multispectral and multitemporal images acquired by different sensors, it integrates in the general approach the following two concepts: (i) spectral, radiometric and geometric homogenization between images acquired by different sensors; and (ii) detection of multiple changes by means of features that guarantees homogeneity over time and across sensors. Experimental results on real datasets, made-up of VHR bi-temporal and multisensor optical images, confirmed the effectiveness of the proposed block scheme and the improvement achieved by the use of higher-level physical features (i.e., TC and OrE) over the traditional features (i.e., TOA or ASR). In the specific cases of datasets 1 and 2, a major improvement is observed when changes from-to vegetation and bare soil, and different types of crops are considered. This given that the higher-level physical feature (i.e., TC) were selected to highlight such type of changes. In the case of dataset 3, the use of OrE for the detection of changes described above, as well as for changes from-to vegetation and bare soil (i.e., forest to grass, crop to grass, bare soil to grass/forest) and small changes on roads and roofs, resulted in better CD accuracy than that obtained by using ASR. In general, both TC and OrE features allow a better separation and interpretation of ${\mathsf{\Omega }}_{\mathrm{Grd}}$ by guaranteeing that these changes are distributed in compact and well separated clusters (when in the 3D feature space).

As future developments, further analysis should be carried out to determine which cluster represents a specific type of change, and to define appropriate features for other types of changes. For the mitigation of remaining ${\mathsf{\Omega }}_{\mathrm{Sys}}$ and the better detection of the ${\mathsf{\Omega }}_{\mathrm{Grd}}$, the use of additional features, either in the physical feature space or in the spatial feature space, could help to make the separation and distinction better, thus improving the final OA. Additional improvements from the ${\mathsf{\Omega }}_{\mathrm{Sys}}$ mitigation process point of view in both spectral and geometric perspective should be considered. For the spectral differences, and given that some of the VHR multisensor optical images have different number of bands and different spectral ranges, the use of regression methods for predicting bands that match from the spectral viewpoint could be considered. For the geometric differences, improvements on the co-registration process by the use of co-registration methods designed for multisensor images could be also explored. Improvements on the detection of building changes could be also integrated.