Robust fusion algorithms for unsupervised change detection between multi-band optical images - A comprehensive case study

Rob ust fusion algorithms for unsupervised change detection between multi-band optical images – A comprehensi v e case study V inicius Ferraris a, , Nicolas Dobigeon a , Marie Chabert a a University of T oulouse, IRIT/INP-ENSEEIHT , 31071 T oulouse, F rance Abstract Unsupervised change detection techniques are generally constrained to two multi- band optical images acquired at dif ferent times through sensors sharing the same spa- tial and spectral resolution. While the optical imagery modality is largely studied in the remote sensing community , this scenario is suitable for a straight comparison of ho- mologous pix els such as pix el-wise dif ferencing. Howe ver , in some speciﬁc cases such as emergenc y situations, punctual missions, defense and security , the only available images may be those acquired through different kinds of sensors with dif ferent resolu- tions. Recently some change detection techniques dealing with images with dif ferent spatial and spectral resolutions, ha ve been proposed. Nev ertheless, the y are focused on a speciﬁc scenario where one image has a high spatial and low spectral resolution while the other has a low spatial and high spectral resolution. This paper addresses the prob- lem of detecting changes between an y two multi-band optical images disregarding their spatial and spectral resolution disparities. T o o vercome those limitations, state-of-the art methods consist in applying conv entional change detection methods after prepro- cessing steps applied independently on the tw o images, e.g. resampling operations intended to reach the same spatial and spectral resolutions. Nev ertheless, these prepro- cessing steps may waste relev ant information since they do not take into account the strong interplay existing between the two images. Con versely , in this paper , we propose Email addresses: vinicius.ferraris@enseeiht.fr (V inicius Ferraris), nicolas.dobigeon@enseeiht.fr (Nicolas Dobigeon), marie.chabert@enseeiht.fr (Marie Chabert) Pr eprint submitted to one journal April 10, 2018 a method that more ef fecti vely uses the a v ailable information by modeling the two ob- served images as spatially and spectrally de graded versions of tw o (unobserved) latent images characterized by the same high spatial and high spectral resolutions. Covering the same scene, the latent images are expected to be globally similar except for possi- ble changes in spatially sparse locations. Thus, the change detection task is en visioned through a rob ust fusion task which enforces the differences between the estimated la- tent images to be spatially sparse. W e sho w that this rob ust fusion can be formulated as an in verse problem which is iterativ ely solved using an alternate minimization strategy . The proposed framew ork is implemented for an exhaustiv e list of applicative scenar - ios and applied to real multi-band optical images. A comparison with state-of-the-art change detection methods e vidences the accuracy of the proposed rob ust fusion-based strategy . K eywor ds: Image fusion, change detection, different resolutions, h yperspectral imagery, multispectral imagery . 1. Introduction Remote sensing consists in collecting measurements, without an y ph ysical contact, about an object or phenomenon. This paper focusses on applications to Earth obser- vation and surface monitoring [ 1 , 10 , 29 , 5 ]. The type of acquired measurements, also referred to as modality , is intimately related to the sensor . Each modality provides a predeﬁned amount and type of information about the scene. The technological growth and new data policies directly impact the availability of multi-temporal data (i.e., ac- quired at dif ferent time instants) [ 3 ], which overcomes this limitation, yet, simulta- neously introducing new challenges. Notably , multi-temporal data acquired over the same geographical location can be used to detect changes or v ariations. Thus, analyz- ing multi-temporal data has culminated in the development of an extremely important area for the remote sensing community , namely , change detection (CD). CD refers to the techniques used to detect areas where potential changes hav e oc- curred between multiple multi-temporal and possibly multi-source (i.e., acquired by 2 different sensors) images acquired over the same scene (geographical location) [ 31 , 8 ]. CD is generally conducted within a supervised or unsupervised context [ 3 ]. The for - mer requires prior ground-truth knowledge in order to train algorithms maximizing the detection rate while minimizing the false alarm rate. Con versely the latter tries to infer changes after carefully designing a blind model-based distance operator . As ground- truth information is not easily av ailable, signiﬁcant ef forts ha ve been made so that un- supervised CD techniques reach the supervised CD performance. Ne vertheless, almost all unsupervised CD methods only focus on a particular scenario, actually the most fa vorable one which considers two multi-band optical images with same spatial and spectral resolutions [ 3 ]. There are two main reasons for considering such a scenario: i) multi-band optical images represent the most commonly used and largely studied re- mote sensing modality , according to [ 33 ], and ii) images with same spatial and spectral resolutions are pix el-wisely comparable, which makes easy the use of simple distance operators. Multi-band optical sensors provide a particular representation of the observed scene according to some intrinsic characteristics, particularly , the ability of sensing the re- ﬂected electromagnetic spectrum of the incoming light. W ell suited to map horizontal structures like land-cover type at large scales [ 7 ], its facility to be directly interpreted has contributed to popularize its use. Another important aspect of optical imaging is the widely admitted Gaussian modeling of the sensor noises, which has lead to a mas- siv e de velopment of least-square lik e methods, specially for CD. Indeed, the properties of the noise model, for instance the symmetry of the Gaussian probability distrib ution function, allow CD techniques to be implemented through a simple operation of im- age dif ferencing, as noticed in [ 31 ] and [ 3 ]. Although CD dif ferencing methods ha ve been adapted to handle multi-band images by considering spectral change vectors [ 2 , 4 ] and transform analysis [ 28 , 27 ], all of them rely on the crucial premise of a fav orable scenario which assumes that the observed images share the same spatial and spectral resolutions. Howe ver , the need for ﬂexible and reliable CD techniques that are able to handle more scenarios is real. In some situations, for instance consecutive to natural disas- ters or within punctual imagery missions, the limited av ailability of the sensors and the 3 time constraints may preclude the use of the same sensor at two distinct time instants. Thus, in these cases, observed images are possibly from different modalities and do not share the same spatial and/or spectral resolutions. T o make existing con ventional CD methods usable in these cases, one strategy , hereafter referred to as the worst-case (WC) method, consists in individually and independently , spatially and/or spectrally , resampling the images to reach the same spatial and spectral resolutions. Although this WC technique allows off-the-shelf CD techniques to be used directly , it may re- mains suboptimal since i) resampling operations independently applied to each image do not take into account their joint characteristics and thus crucial information may be missed and ii) these spatial and spectral operations are generally from a higher to a lower resolution, which results in a signiﬁcant loss of information. T o o vercome these limitations, [ 14 ] and [ 15 ] recently proposed two CD approaches speciﬁcally designed to deal with multi-band images with different spatial and spectral resolution. Both ap- proaches rely on the inference of a latent (i.e., unobserved) image which results from the fusion of the two observed images. Fusing information contained in remote sens- ing images has moti v ated a lot of research works in the literature [ 24 , 23 , 32 , 17 , 25 ]. W ithin a CD context, the underlying assumption is most of pixels of the fused image, which are supposed not to ha ve been changed during the time interv al, produce consis- tent information while the fe w remaining ones, locating in the change regions, produce aberrant information. More precisely , the method proposed in [ 15 ] is based on a 3 -step procedure (namely fusion, prediction and detection) which, instead of independently preprocess each observed image, recov ers a latent high spatial and spectral resolution image containing changed and unchanged regions by fusing observ ed images. Then, it predicts pseudo-observed images by artiﬁcially degrading this estimated latent image using forward models underlying the actually observed images. The result is two pairs, each composed of a predicted image and an observ ed image with the same spatial and spectral resolutions. Then, any classical multi-band CD method can be ﬁnally applied to estimate two change images, that can be thresholded to build the change maps. Con- versely , [ 14 ] propose a robust fusion-based CD technique which aims at recovering two high spatial and spectral resolution latent images related to the observed images via a double physically-inspired forward model. The rob ust fusion of multi-band images is 4 formulated as an inv erse problem. The difference between the two latent images is as- sumed to be spatially sparse, implicitly locating the changes at a high resolution scale. The resulting objective function is solv ed through the use of an alternating minimiza- tion algorithm, which iterati vely optimizes with respect to (w .r .t.) one latent image and the change image. The CD map can be ﬁnally generated from the recovered change image. Both methods of fer a way to conduct unsupervised CD between two multi-band optical images. Even if they hav e sho wn signiﬁcant impro vements in the detection rate when compared to WC method, they are also still limited to a single scenario: one high spatial lo w spectral resolution image and one lo w spatial high spectral resolution image. In this paper , capitalizing on the method proposed in [ 14 ], we show that the unsu- pervised CD task can be formulated in a general robust-fusion form for all multi-band optical image scenarios in volving two observed images. Contrary to the pre vious tech- nique, the proposed approach is not limited to one high spatial and low spectral res- olution observed image and one low spatial high spectral resolution observed one. It generalizes the robust fusion model to handle all possible conﬁgurations of two multi- band optical images. Note that the scenario and the solution proposed in [ 14 ] is a speciﬁc instance of the framework dev eloped in this paper . Therefore, the same as- sumptions re garding the two observed images is adopted. Namely the y can be jointly approximated by a standard linear decomposition model complemented with an out- lier term corresponding to the change image. The outlier term is still characterized by a spatial sparsity-inducing regularization. The resulting objecti ve function, re gardless the scenario proposed, is solv ed through the use of an alternate minimization algorithm. Remarkably , optimizing w .r .t. the latent image al ways relies on a closed-form solution, which ensures the conv ergence of the alternate minimization procedure. The CD map can be ﬁnally generated from the recov ered change image. The paper is organized as follows. Section 2 formulates the change detection prob- lem for multi-band optical image. Section 3 presents the solution for the formulated problem based on robust fusion for each possible scenarios. Experimental CD ex- amples are considered in Section 5 for each possible scenario described in Section 3 . Section 6 concludes the paper . 5 2. Problem formulation 2.1. Generic forwar d model for multi-band optical images In digital image processing, the image formation process inherent to multi-band op- tical sensors can be generally modeled as a sequence of successi ve transformations and degradations. These transformations are applied over the original scene and resulting in an output image, commonly referred as the observed image and denoted Y ∈ R m λ × m where m and m λ are the numbers of pix els and spectral bands in the observ ed image. It corresponds to the particular limited representation of the original scene according to the characteristics imposed by the image signal processor describing the sensor . The original scene cannot be exactly represented because of its continuous nature, but it can be conv eniently approximated by an (unkno wn) latent image of higher spatial and spectral resolutions, X ∈ R n λ × n , where n ≥ m and n λ ≥ m λ are the numbers of pixels and spectral bands, respecti vely . In what follows, as a well-admitted approxima- tion, the observ ed and latent images are assumed to be related according to the generic forward model [ 37 , 40 , 30 ] Y = LXR + N (1) where • L ∈ R m λ × n λ is a spectral degradation matrix, • R ∈ R n × m is a spatial degradation matrix, • N is an additiv e term comprising sensor noise and modeling errors. In ( 1 ), the left-multiplying matrix L and right-multiplying matrix R spectrally and spatially degrade the latent image, respecti vely , by combining some spectral bands of each pixel or by combining neighboring pixel measurements in each spectral band. More precisely , the spectral degradation L represents a spectral resolution reduction with respect to the latent image X , as already considered in [ 40 ], [ 30 ] and [ 38 ]. In practice, this matrix is fully deﬁned by spectral ﬁlters characterizing the optical sensors. When the speciﬁcations of the sensor are a vailable, these ﬁlters are known. Otherwise, they can be learned by cross-calibration, e.g., follo wing the strategies proposed in [ 30 ] 6 or [ 39 ]. On the other hand, the spatial degradation matrix R models the combination of different spatial transformations applied to the pixel measurements within each spectral band. These transformations are speciﬁc of the sensor architecture and include warp, blur , translation and decimation [ 39 , 38 ]. In this work, geometrical transformations such as warp and translation are assumed to have been pre viously corrected, e.g., using image spatial alignment techniques. Thus, similarly to the model considered in [ 38 ], the spatial degradation matrix R only stands for a spatially in v ariant blurring, followed by a decimation (i.e., do wnsampling) operation. Thus, in what follows, the spatial degradation matrix R will be assumed of the form R = BS . (2) The sparse symmetric T oeplitz matrix B ∈ R n × n in ( 2 ) operates a cyclic conv olution on each individual band to model a space-in v ariant blur associated with a symmetric con volution k ernel. The decimation operation, denoted by the n × m matrix S in ( 2 ), corresponds to a uniform do wnsampling operator 1 of factor d = d r × d c with m = n/d ones on the block diagonal and zeros elsewhere, such that S T S = I m [ 38 ]. The noise corrupting multi-band optical images is generally modeled as additive and Gaussian [ 3 , 10 , 26 , 38 ]. Thus the noise matrix N in ( 1 ) is assumed to be distrib uted according to the following matrix normal distrib ution (see Appendix Appendix A ) N ∼ MN m λ ,m ( 0 m λ × m , Λ , Π ) . The row co variance matrix Λ carries information regarding the between-band spectral correlation. In what follows, similarly to the approach in [ 38 ], this cov ariance matrix Λ will be assumed to be diagonal, which implies that the noise is spectrally indepen- dent and characterized by a speciﬁc variance in each band. Con versely , the column cov ariance matrix Π models the noise correlation w .r .t. to the pixel locations. Follo w- ing a hypothesis widely admitted in the literature, this matrix is assumed to be identity , Π = I m , to reﬂect the fact the noise is spatially independent. In real applications, both matrices Λ and Π can be estimated by calibration [ 39 ]. 1 The corresponding operator S T represents an upsampling transformation by zero-interpolation from m to n . 7 2.2. Pr oblem statement Let us consider two co-registered multi-band optical images Y 1 ∈ R m λ 1 × m 1 and Y 2 ∈ R m λ 2 × m 2 acquired by two sensors S 1 and S 2 at times t 1 and t 2 , respectively . It is not assumed any speciﬁc information about time ordering of acquisitions, either t 2 < t 1 or t 2 > t 1 are possible cases. The problem addressed in this paper consists in detecting signiﬁcant changes between these two multi-band optical images. This is a challenging task mainly due to the possible spatial and/or spectral resolution dis- similarity (i.e., m λ 1 6 = m λ 2 and/or m 1 6 = m 2 ), which prev ents any use of simple yet efﬁcient dif ferencing operation [ 31 , 3 ]. T o alleviate this issue, this work proposes to generalize the CD framew ork introduced in [ 14 ] to handle all possible combinations (scenarios) of tw o multi-band optical images. More precisely , following the widely ad- mitted forward model described in Section 2.1 and adopting consistent notations, the observed images Y 1 and Y 2 can be related to two latent images X 1 ∈ R n λ × n and X 2 ∈ R n λ × n with the same spatial and spectral resolutions Y 1 = L 1 X 1 R 1 + N 1 (3a) Y 2 = L 2 X 2 R 2 + N 2 . (3b) where L 1 and L 2 denote two spectral degradation operators and R 1 and R 2 are two spatial de gradation operators that can be decomposed according to ( 2 ). Note that ( 3a ) and ( 3b ) are a speciﬁc double instance of the model ( 1 ). In particular, the two multi- band latent images X 1 and X 2 share the same spectral and spatial resolutions, generally higher than those of the observed images: n λ ≥ max { m λ 1 , m λ 2 } and/or n ≥ max { m 1 , m 2 } . (4) Thereby , after inferring the latent images, any classical differencing technique can be subsequently implemented to compute a change image ∆ X = [∆ x 1 , . . . , ∆ x n ] ∈ R n λ ,n deﬁned by ∆ X = X 2 − X 1 (5) where ∆ x p ∈ R n λ denotes the spectral change vector in the p th pixel ( p = 1 , . . . , n ). It is worth noting that, under the assumptions ( 4 ), these changes can be identiﬁed at 8 a high spatial and spectral resolutions. Finally this change image can be further ex- ploited by conducting a pixel-wise change vector analysis (CV A) which exhibits the polar coordinates (i.e., magnitude and direction) of the spectral change vectors [ 22 ]. Then, to spatially locate the changes, a natural approach consists in monitoring the in- formation contained in the magnitude part of this representation, summarized by the change energy image [ 31 , 2 , 4 ] e = [ e 1 , . . . , e n ] ∈ R n with e p = k ∆ x p k 2 , p = 1 , . . . , n. When the CD problem in the p th pixel is formulated as the binary hypothesis testing    H 0 ,p : no change occurs in the p th pixel H 1 ,p : a change occurs in the p th pixel a pixel-wise statistical test can be written by thresholding the change energy image pixels e p H 1 ,p ≷ H 0 ,p τ . The ﬁnal binary CD map denoted d = [ d 1 , . . . , d n ] ∈ { 0 , 1 } n can be deriv ed as d p =    1 if e p ≥ τ ( H 1 ,p ) 0 otherwise ( H 0 ,p ) . As a consequence, to solv e the multi-band image CD problem, the k ey issue lies in the joint estimation of the pair of HR latent images { X 1 , X 2 } from the joint forward model ( 3 ) or , equi v alently , the joint estimation of one latent image and the difference image, i.e., { X 1 , ∆ X } . Finally , the next paragraph introduces the CD-driv en optimization problem to be solved. 2.3. Optimization pr oblem Follo wing a Bayesian approach, the joint maximum a posteriori (MAP) estimator n ˆ X 1 , MAP , ∆ ˆ X MAP o of the latent and change images can be deriv ed by maximizing 9 the posterior distribution p ( X 1 , ∆ X | Y 2 , Y 1 ) ∝ p ( Y 2 , Y 1 | X 1 , ∆ X ) p ( X 1 ) p (∆ X ) where p ( Y 2 , Y 1 | X 1 , ∆ X ) is the joint likelihood function and p ( X 1 ) and p (∆ X ) cor- respond to the prior distrib utions associated with the latent and change images, respec- tiv ely , assumed to be a priori independent. Because of the additive nature and statistical properties of the noise discussed in Section 2.1 , this boils do wn to solve the following minimization problem n ˆ X 1 , MAP , ∆ ˆ X MAP o ∈ argmin X 1 , ∆ X J ( X 1 , ∆ X ) (6) with J ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − L 2 ( X 1 + ∆ X ) R 2 )    2 F + 1 2    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F + λφ 1 ( X 1 ) + γ φ 2 (∆ X ) . (7) where k·k F denotes the Frobenius norm. The regularizing functions φ 1 ( · ) and φ 2 ( · ) can be related to the negati ve log-prior distributions of the latent and change images, respectiv ely , and the parameters λ and γ tune the amount of corresponding penaliza- tions in the ov erall objectiv e function J ( X 1 , ∆ X ) . These functions should be carefully designed to e xploit any prior knowledge regarding the parameters of interest. As dis- cussed in Section 3.1 , numerous regularizations can be advocated for the latent image X 1 . In this w ork, to maintain computational efﬁcienc y while providing accurate results [ 26 ], a T ikhonov regularization proposed in [ 37 ] has been adopted φ 1 ( X 1 ) =   X 1 − ¯ X 1   2 F where ¯ X 1 refers to a crude estimate of X 1 . Regarding the regularizing function φ 2 ( · ) , as already mentioned in the previous section, it should reﬂect the fact that most of the pixels are expected to remain un- changed i.e., most of the columns of the change image ∆ X are expected to be null vectors. Thus, the regularizing function φ 2 ( · ) is chosen as in [ 16 ] as the ` 2 , 1 -norm of the change image φ 2 (∆ X ) = k ∆ X k 2 , 1 = n X p =1 k ∆ x p k 2 . (8) 10 The ne xt section describes the general iterativ e algorithm scheme which solv es the minimization problem ( 6 ). 3. Robust multi-band image fusion algorithm: generic formulation Computing the joint MAP estimator of the latent image X 1 at time t 1 and of the change image ∆ X can be achie ved by solving the minimization problem in ( 6 ). How- ev er, no closed-form solution can be deriv ed for this problem for all the scenarios of interest. Thus this section introduces a minimization algorithm which iterati vely con- ver ges to this solution. This alternating minimization (AM) algorithm, summarized in Algo. 1 , consists in iterativ ely minimizing the objecti ve function ( 7 ) w .r .t. X 1 and ∆ X , within so-called fusion and corr ection discussed below . Algorithm 1 Algorithm for robust multi-band image fusion Input: Y 1 , Y 2 , L 1 , L 2 , R 1 , R 2 , Λ 1 , Λ 2 . 1: Set ∆ X 1 . 2: for k = 1 , . . . , K do 3: % Fusion step 4: X ( k +1) 1 = arg min X 1 J ( X 1 , ∆ X ( k ) ) 5: % Corr ection step 6: ∆ X ( k +1) = arg min ∆ X J ( X ( k +1) 1 , ∆ X ) 7: end for Output: ˆ X 1 , MAP , X ( K +1) 1 and ∆ ˆ X MAP , ∆ ˆ X ( K +1) 3.1. Fusion step As mentioned abov e, the forward model ( 3 ) relying on the pair { X 1 , X 2 } of latent images can be rewritten as a function of { X 1 , ∆ X } , i.e., Y 1 = L 1 X 1 R 1 + N 1 (9a) Y 2 = L 2 ( X 1 + ∆ X ) R 2 + N 2 . (9b) 11 Generalizing the strate gy proposed in [ 14 ], gi ven the change image ∆ X and the image Y 1 observed at time t 1 , a corrected image denoted ˜ Y 2 that would be acquired by the sensor S 2 at time t 1 can be deﬁned as ˜ Y 2 = Y 2 − L 2 ∆ XR 2 . (10) W ith this notation, the forward model ( 9 ) can be easily rewritten, leading to Y 1 = L 1 X 1 R 1 + N 1 (11a) ˜ Y 2 = L 2 X 1 R 2 + N 2 . (11b) Thus, the fusion step, at iteration k , consists in minimizing ( 7 ) w .r .t. X 1 , i.e., ˆ X ( k +1) 1 = argmin X 1 J 1 ( X 1 ) , J  X 1 , ∆ X ( k )  with J 1 ( X 1 ) = 1 2    Λ − 1 2 2  ˜ Y ( k ) 2 − L 2 X 1 R 2     2 F + 1 2    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F + λ   X 1 − ¯ X 1   2 F . (12) The double forw ard model ( 11 ), as well as the optimization problem ( 12 ), underly the estimation of an image X 1 from an observ ed image Y 1 and a pseudo-observ ed image ˜ Y 2 . V arious instances of this pixel-le vel fusion problem have been widely considered in the literature [ 24 , 32 , 17 , 25 ]. For instance, [ 21 ] and [ 42 ] hav e addressed the problem of single mono-band image superresolution from a single observed image Y 1 , i.e., with L 1 = I m λ 1 and m λ 1 = n λ = 1 . The problem of fusing se veral degraded mono- band images to recover a common high resolution latent image has been considered in [ 11 ]. Similarly , the model ( 11 ) generalizes the con ventional observational model widely adopted by the remote sensing community to conduct multi-band image fusion [ 19 , 9 , 41 , 40 , 23 , 30 , 36 , 37 , 38 ]. W ithin this speciﬁc scenario, a high spatial and high spectral resolution latent image X 1 is estimated from two observ ed images, one of low spatial and high spectral resolutions (i.e., L 1 = I m λ 1 ) and the other of high spatial and low spectral resolutions (i.e., R 2 = I n 2 ). In this context, the CD task considered in this paper can be cast as a so-called r ob ust fusion problem since the multi-band image fusion model ( 11 ) implicitly depends on the 12 (unknown) change image ∆ X . More precisely , since the two latent images X 1 and X 2 are related to the same scene observed at two time instants, they are expected to share a high lev el of similarity , i.e., the change image ∆ X is expected to be spatially sparse. Thus, this additional unkno wn change image ∆ X to be inferred can be considered as an outlier term, akin to those encountered in several robust factorizing models such as robust principal component analysis (RPCA) [ 6 ] and rob ust nonnegati ve factorization [ 16 ]. A particular instance of this strategy has been successfully adopted in [ 14 ] to detect changes between two complementary multi-band images, i.e., in the particular case of L 1 = I m λ 1 and R 2 = I n 2 . In this work, we propose to follow a similar route while generalizing the approach to the much more generic model ( 3 ) to handle all practical scenarios of CD. These different scenarios are discussed in the next paragraph. 3.2. Corr ection step Giv en the current state X 1 of the latent image, the predicted image that would be observed by the sensor S 2 at time t 1 can be deﬁned as ˇ Y ( k ) 2 = L 2 X ( k ) 1 R 2 (13) leading to the predicted change image ∆ ˇ Y 2 = Y 2 − ˇ Y 2 . (14) Then, the correction step in Algo. 1 consists in solving ∆ ˆ X ( k +1) = argmin ∆ X J 2 (∆ X ) , J  X ( k ) 1 , ∆ X  (15) with J 2 (∆ X ) =    Λ − 1 2 2  ∆ ˇ Y ( k ) 2 − L 2 ∆ XR 2     2 F + γ k ∆ X k 2 , 1 . (16) This correction can be interpreted as a joint spatial and spectral deblurring of the pre- dicted change image ∆ ˇ Y ( k ) 2 . Note that this ill-posed inv erse problem is regularized through an ` 2 , 1 -norm penalization, which promotes the spatial sparsity of the change image ∆ X . It is worth noting that the difﬁculty of conducting the two steps of the AM algorithm detailed above is highly related to the spatial and/or spectral degradations operated on 13 Forward model ] 1 Forward model ] 2 Comments Spectral Spatial Spectral Spatial degradation degradation degradation degradation S 1 − − − − Con ventional CD framew ork – Y 1 and Y 2 of same spatial and spectral resolutions S 2 L 1 − − − Y 1 of lower spectral resolution Y 1 and Y 2 of same spatial resolutions S 3 − R 1 − − Y 1 of lower spatial resolution Y 1 and Y 2 of same spectral resolutions S 4 − R 1 L 2 − Y 1 and Y 2 of complementary resolutions S 5 L 1 R 1 − − Y 1 of low spatial and spectral resolutions S 6 − R 1 − R 2 Generalization of S 3 with non-integer relativ e spatial downsampling factor S 7 L 1 R 1 − R 2 Generalization of S 4 with non-integer relativ e spatial downsampling factor S 8 L 1 − L 2 − Generalization of S 2 with some non-ov erlapping spectral bands S 9 L 1 R 1 L 2 − Generalization of S 4 with some non-ov erlapping spectral bands S 10 L 1 R 1 L 2 R 2 Generalization of S 4 with some non-ov erlapping spectral bands and non-integer relati ve spatial downsampling f actor T able 1: Overvie ws of the spectral and spatial degradations w .r .t. to experimental scenarios. The symbol − stands for “no degradation”. 14 the two latent images, according to applicati ve scenarios which are detailed in the ne xt section. Interestingly , the following section will also show that these steps generally reduce to ubiquitous (multi-band) image processing tasks, namely denoising, spectral deblurring or spatial super -resolution from a single or multiple images, for which efﬁ- cient and reliable strategies ha ve been already proposed in the literature. 4. Algorithmic implementations for applicative scenarios The general model presented in ( 3 ) and the AM algorithm proposed in Section 1 can be implemented to handle all scenarios deriv ed from two multi-band optical images. These scenarios dif fer by the corresponding spatial and spectral degradations relating the pair of observ ed images { Y 1 , Y 2 } and the pair of latent images { X 1 , X 2 } . T able 1 summarizes the 10 distinct scenarios (denoted S 1 to S 10 ) according to the de gradations operated on the two latent images X 1 and X 2 . The speciﬁcities of these scenarios are also discussed in what follows. S 1 is dev oted to a pair of observed images sharing the same spatial and spectral reso- lutions. In this case, CD can be conducted by pixel-wise comparisons, as classi- cally addressed in the literature, e.g., in [ 31 ] and [ 3 ]. S 2 consists in conducting CD between two images with the same spatial resolution b ut different spectral resolution, considered in [ 28 ] and [ 27 ]. S 3 consists in conducting CD between two images with the same spectral resolution but dif ferent spatial resolution. S 4 relies on two complementary images: the ﬁrst image with high spectral and low spatial resolutions, the second image with low spectral and high spatial resolu- tions. This is the CD scenario considered in [ 13 , 15 , 14 ]. When the two observed images hav e been acquired at the same time instants ( t i = t j ), this scenario cor- responds to the multi-band image fusion task considered in numerous w orks, e.g., in [ 37 ], [ 40 ] and [ 30 ]. 15 S 5 represents an ev en less fav orable instance of S 2 and S 3 where one image is of high spatial and spectral resolutions while the other is of lo w spatial and spectral resolutions. S 6 generalizes S 3 . As for S 3 , both observed images ha ve the same spectral resolutions and dif ferent spatial resolutions. Howev er, contrary to S 3 , the relativ e do wnsam- pling factor between images is non-integer , which precludes the use of a unique spatial degradation matrix R = BS . As a consequence, the latent images X 1 and X 2 are characterized by a common spatial resolution which is higher than those of both observ ed images. The choice of this virtual downsampling factor is based on the greatest common divisor between spatial resolutions. S 7 generalizes S 4 with a non-integer relati ve downsampling f actor (as for S 6 ). S 8 generalizes S 2 where the two observ ed image share the same spatial resolution b ut hav e distinct spectral resolutions. Howe ver , contrary to S 2 , this dif ference in spectral resolutions cannot be expressed using a unique spectral degradation ma- trix. This may happen when the two spectral ranges of observed images contain non-ov erlapping bands. S 9 generalizes S 4 , b ut the dif ference in spectral resolutions cannot be expressed using a single degradation matrix (as for S 8 ). S 10 generalizes S 4 , but the dif ference in spatial resolutions cannot be expressed using a unique spatial degradation matrix (as for S 6 ) and the dif ference in spectral resolutions cannot be expressed using a single spectral degradation matrix (as for S 8 ). Finally , the following paragraphs instantiate the AM algorithm for each scenario. These speciﬁc instantiations will relate the fusion and correction steps with ubiquitous image processing tasks that can be performed ef ﬁciently thanks to recent contrib utions proposed in the image processing literature. T able 2 summarizes these implementations w .r .t. the discussed scenarios. 16 Fusion Step Correction Step Algorithm Operation Algorithm Operation S 1 Least squares Denoising ` 2 , 1 -prox. mapping Denoising S 2 Least squares Spectral deblurring ` 2 , 1 -prox. mapping Denoising S 3 [ 42 ] Spatial super-resolution ` 2 , 1 -prox. mapping Denoising S 4 [ 38 ] Multi-band image fusion Forward-backw ard Spectral deblurring S 5 ADMM Least squares Spectral deblurring ` 2 , 1 -prox. mapping Denoising [ 42 ] Spatial super-resolution S 6 ADMM [ 42 ] Spatial super-resolution ADMM [ 42 ] Spatial super-resolution [ 42 ] Spatial super-resolution ` 2 , 1 -prox. mapping Denoising S 7 ADMM [ 38 ] Multi-band image fusion ADMM [ 42 ] Spatial super-resolution [ 42 ] Spatial super-resolution ` 2 , 1 -prox. mapping Denoising S 8 Least squares Spectral deblurring Forward-backw ard Spectral deblurring S 9 ADMM Least squares Spectral deblurring Forward-backw ard Spectral deblurring [ 42 ] Multi-band image fusion S 10 ADMM [ 38 ] Multi-band image fusion ADMM ` 2 , 1 -prox. mapping Denoising [ 42 ] Spatial super-resolution [ 42 ] Spatial super-resolution Least squares Spectral deblurring Least squares Spectral deblurring T able 2: Overvie w of the steps of the AM algorithm w .r .t. applicative scenarios. 17 4.1. Scenario S 1 Considering the de gradation matrices speciﬁed in T able 1 for this scenario, the forward model ( 9 ) can be re written as Y 1 = X 1 + N 1 (17a) Y 2 = ( X 1 + ∆ X ) + N 2 (17b) As e xpected, for this scenario, the observ ed, latent and change images share the same spatial and spectral resolutions. The resulting objectiv e function, initially in ( 7 ), is simpliﬁed as J S 1 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − ( X 1 + ∆ X ))    2 F + 1 2    Λ − 1 2 1 ( Y 1 − X 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . (18) The two steps of the AM algorithm are detailed belo w . 4.1.1. Fusion: optimization w .r .t. X 1 At the k th iteration of the AM algorithm, let assume that the current value of the change image is denoted by ∆ X ( k ) . As suggested in Section 3.1 , a corrected image ˜ Y ( k ) 2 that would be observ ed at time t 1 by the sensor S 2 giv en the image Y 2 observed at time t 2 and the change image ∆ X ( k ) can be introduced as ˜ Y ( k ) 2 = Y 2 − ∆ X ( k ) . (19) Updating the latent image X 1 consists in minimizing w .r .t. X 1 the partial function J S 1 , 1 ( X 1 ) , J S 1  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − X 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1     2 F + λ   X 1 − ¯ X 1   2 F . This formulation sho ws that recovering X 1 in Scenario S 1 reduces to a denoising prob- lem from an observed image Y 1 and a pseudo-observed image ˜ Y ( k ) 2 . A closed-form solution of this ` 2 -penalized least-square problem can be easily and efﬁciently com- puted. 18 4.1.2. Corr ection: optimization w .r .t. ∆ X Follo wing the same strate gy proposed in [ 14 ], let ˇ Y ( k ) 2 denote the pr edicted image that would be observed by the sensor S 2 at time t 1 giv en the current state of the latent image X ( k ) 1 . Since the two sensors share the same spatial and spectral characteristics, one has ˇ Y ( k ) 2 = X ( k ) 1 . (20) Similarly to ( 5 ), the pr edicted change image can thus be deﬁned as ∆ ˇ Y ( k ) 2 = Y 2 − ˇ Y ( k ) 2 . (21) The objectiv e function ( 18 ) w .r .t. ∆ X is then rewritten by combining ( 20 ) and ( 21 ) with ( 18 ), leading to J S 1 , 2 (∆ X ) , J S 1 ( X ( k ) 1 , ∆ X ) =    Λ − 1 2 2  ∆ ˇ Y ( k ) 2 − ∆ X     2 F + γ k ∆ X k 2 , 1 . (22) Again, since the observ ed, latent and change images share the same spatial and spec- tral resolutions, this correction step reduces to a denoising task of the predicted change image ∆ ˇ Y ( k ) 2 . W ith the particular CD-driv en choice of φ 2 ( · ) in ( 8 ), minimizing J S 1 , 2 (∆ X ) is an ` 2 , 1 -penalized least square problem. Minimizing ( 22 ) also deﬁnes the proximal operator associated with the ` 2 , 1 -norm and can be directly achieved by applying a group-soft thresholding on the predicted change image ∆ ˇ Y ( k ) 2 . 4.2. Scenario S 2 In this scenario, the tw o observ ed images are of same spatial resolution (as for scenario S 1 ) but with different optical spectral information, which preclude a simple comparison between pixels. For this scenario, the joint forward observation model deriv ed from ( 9 ) can be written as Y 1 = L 1 X 1 + N 1 , (23a) Y 2 = ( X 1 + ∆ X ) + N 2 , (23b) 19 which results in the objectiv e function J S 2 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − ( X 1 + ∆ X ))    2 F + 1 2    Λ − 1 2 1 ( Y 1 − L 1 X 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . W ithin an AM algorithmic schemes, the two sub-problems of interest are detailed be- low . 4.2.1. Fusion: optimization w .r .t. X 1 The same strategy as for scenario S 1 in paragraph 4.1.1 is adopted. As model ( 23b ) is the the same as model ( 17b ), the corrected image ˜ Y ( k ) 2 is deﬁned following ( 19 ). Then, updating the latent image X 1 consists in minimizing the partial objectiv e function J S 2 , 1 ( X 1 ) , J S 2  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − L 1 X 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1     2 F + λ   X 1 − ¯ X 1   2 F . (24) This problem can be interpreted as a spectral deblurring of the observ ed image Y 1 where the corrected image ˜ Y ( k ) 2 plays the role of prior information. Minimizing ( 24 ) can be easily conducted by computing the standard least square solution. 4.2.2. Corr ection: optimization w .r .t. ∆ X As both models ( 23b ) and ( 17b ) are the same, optimizing w .r .t ∆ X can be con- ducted follo wing the procedure detailed in paragraph 4.1.2 (i.e., denoising of the pre- dicted change image). 4.3. Scenario S 3 In this scenario, the two observed images share the same spectral resolution but differ by their spatial resolutions. These spatial resolutions are related by an integer relativ e do wnsampling factor , which allows a unique spatial degradation matrix R 1 20 to be used 2 . The joint forward observ ation model deriv ed from ( 9 ) using the speciﬁc degradation matrices presented in T able 1 can be written as Y 1 = X 1 R 1 + N 1 . (25a) Y 2 = ( X 1 + ∆ X ) + N 2 . (25b) with the objectiv e function J S 3 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − ( X 1 + ∆ X ))    2 F + 1 2    Λ − 1 2 1 ( Y 1 − X 1 R 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . 4.3.1. Fusion: optimization w .r .t. X 1 The same strategy as for previous scenarios is adopted here. As model ( 25b ) is the same as model ( 17b ), the corrected image ˜ Y ( k ) 2 is deﬁned following ( 19 ). Then, updating the latent image consists in minimizing w .r .t. X 1 the partial function J S 3 , 1 ( X 1 ) , J S 3  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − X 1 R 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1     2 F + λ   X 1 − ¯ X 1   2 F . This fusion task can be interpreted as a set of n λ super-resolution problems associated with each band of the observed image Y 1 , where the corrected image ˜ Y ( k ) 2 acts here as a prior information. Closed-form expressions of these n λ solutions are giv en by [ 42 ]. 4.3.2. Corr ection: optimization w .r .t. ∆ X As the model ( 25b ) is the same as model ( 17b ) of scenarios S 1 and S 2 , optimizing w .r .t. ∆ X can be conducted follo wing the procedure detailed in paragraph 4.1.2 (i.e., denoising of the predicted change image). 2 The case of observed images with non-integer relative spatial downsampling factor is discussed in sce- nario S 6 . 21 4.4. Scenario S 4 Scenario S 4 is speciﬁcally addressed in [ 14 ] with the joint forward model Y 1 = X 1 R 1 + N 1 , Y 2 = L 2 ( X 1 + ∆ X ) + N 2 . The two observed images hav e complementary information since Y 1 and Y 2 are of high spectral and spatial resolutions, respecti vely . The resulting objecti ve function writes J S 4 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − L 2 ( X 1 + ∆ X ))    2 F + 1 2    Λ − 1 2 1 ( Y 1 − X 1 R 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . (26) When these images hav e been acquired at the same time instant, the change image is ∆ X = 0 and this conﬁguration boils down to a multiband image fusion problem addressed in [ 38 ]. Thus, minimizing ( 26 ) can be conducted following the AM strategy by combining a multiband image fusion step [ 38 ] and a spectral deblurring step of the predicted change image. The interested reader is in vited to consult the work in [ 14 ] for a comprehensiv e description of the resolution. 4.5. Scenario S 5 Under this scenario, the observed image Y 2 is of higher spatial and spectral reso- lutions than the observed image Y 1 . W ithin a con ventional fusion conte xt, one would probably discard Y 1 since it would not bring additional information to the one provided by Y 2 . Con versely , within a CD context, both observed images are of interest and can be exploited. More precisely , here, the joint forward observ ation model derived from ( 9 ) is speciﬁcally written Y 1 = L 1 X 1 R 1 + N 1 , (27a) Y 2 = ( X 1 + ∆ X ) + N 2 , (27b) 22 with the resulting objectiv e function J S 5 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − ( X 1 + ∆ X ))    2 F + 1 2    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . Its minimization relies on the two steps detailed belo w . 4.5.1. Fusion: optimization w .r .t. X 1 The same strategy as for pre vious scenarios is adopted here. After deﬁning the corrected image ˜ Y ( k ) 2 by ( 19 ), updating the the latent image X 1 consists in minimizing J S 5 , 1 ( X 1 ) , J S 5  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1     2 F + λ   X 1 − ¯ X 1   2 F . (28) Minimizing ( 28 ) can be interpreted as a simultaneous spatial super-resolution and spec- tral deblurring of the multiband image Y 1 , with prior information brought by ˜ Y ( k ) 2 . This minimization is a much more challenging task than the fusion steps encountered for scenarios S 1 – S 4 . Indeed, the simultaneous spatial and spectral degradations applied to X 1 prev ents a closed-form solution to be efﬁciently computed. Thus, one proposes to resort to an iterativ e algorithm, namely the alternating direction method of multipli- ers (ADMM). It consists in introducing the splitting v ariable U ∈ R m λ 1 × n = L 1 X 1 . The resulting scaled augmented Lagrangian for the problem is expressed as L µ ( X 1 , U , V ) =    Λ − 1 2 1 ( Y 1 − UR 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1     2 F + λ   X 1 − ¯ X 1   2 F + µ 2 k L 1 X 1 − U + V k 2 F . (29) The ADMM iterativ ely minimizes L µ w .r .t. U and X 1 and updates the dual variable V . By comparing the partial objectiv e function ( 28 ) and its augmented counterpart ( 29 ), it clearly appears that the splitting strategy allows the spectral and spatial degra- dations to be decoupled. Thus, each of these steps can be easily conducted. More 23 precisely , optimizing w .r .t. U consists in conducting a super-resolution step achie ved as for scenario S 3 by resorting to the algorithm proposed in [ 42 ]. Con versely , optimiz- ing w .r .t. X 1 consists in solving a least-square problem whose closed-form solution can be computed (akin to scenario S 2 ). 4.5.2. Corr ection: optimization w .r .t. ∆ X Again, as the forward model ( 27b ) is the same as ( 17b ) of Scenario S 1 , optimizing w .r .t. ∆ X can be conducted follo wing the procedure detailed in paragraph 4.1.2 (i.e., denoising of the predicted change image). 4.6. Scenario S 6 As for scenario S 3 , scenario S 6 considers two observed images of same spectral resolutions b ut with distinct spatial resolutions. Howe ver , contrary to scenario S 3 , this difference in spatial resolutions cannot be expressed thanks to a unique spatial degradation matrix R 1 due to a non-integer relativ e downsampling factor . Thus the forward model is written Y 1 = X 1 R 1 + N 1 . (30a) Y 2 = ( X 1 + ∆ X ) R 2 + N 2 . (30b) with the following objecti ve function J S 6 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − ( X 1 + ∆ X ) R 2 )    2 F + 1 2    Λ − 1 2 1 ( Y 1 − X 1 R 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . (31) In ( 30 ), both latent images are supposed to suffer from spatial degradations. Thus, choosing which spatial degradation affects the change image ∆ X results in a particular spatial resolution for this change map. T o deriv e a CD map at a high spatial resolution, the spatial degradation applied to ∆ X should be chosen as the one with the lowest virtual downsampling factor . The minimization of ( 31 ) according to the AM strategy is addressed in the following paragraphs. 24 4.6.1. Fusion: optimization w .r .t. X 1 For this scenario, the corrected image in ( 10 ) is deﬁned as ˜ Y ( k ) 2 = Y 2 − ∆ X ( k ) R 2 . Then, updating the latent image X 1 consists in minimizing w .r .t. X 1 the partial func- tion J S 6 , 1 ( X 1 ) , J S 6  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − X 1 R 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1 R 2     2 F + λ   X 1 − ¯ X 1   2 F . (32) As for scenario S 3 , minimizing ( 32 ) can be interpreted as recov ering a spatially super- resolved image X 1 from the observed image Y 1 and the corrected image ˜ Y ( k ) 2 . How- ev er, contrary to scenario S 3 , here, ˜ Y ( k ) 2 rather deﬁnes an additional data-ﬁtting term instead of a prior information [ 11 ]. Moreov er , this sub-problem cannot be solved di- rectly since no closed-form solution can be efﬁciently derived, mainly due to the simul- taneous presence of the two spatial degradation operators. Thus, as for scenario S 5 , one resorts to the ADMM scheme by introducing the splitting v ariable U ∈ R n λ × n = X 1 . The resulting scaled augmented Lagrangian can be written as L µ ( X 1 , U , V ) =    Λ − 1 2 1 ( Y 1 − UR 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1 R 2     2 F + λ   X 1 − ¯ X 1   2 F + µ 2 k X 1 − U + V k 2 F . (33) Both minimizations of ( 33 ) w .r .t. U and X 1 can be conducted band-by-band follo wing the strate gy proposed in [ 42 ], which provides closed-form solutions of the underlying single-image super-resolution problems and also ensures the con ver gence of the AM algorithm. 4.6.2. Corr ection: optimization w .r .t. ∆ X For this scenario, a predicted image that would be observed by the sensor S 2 at time t 1 can be deﬁned as ˇ Y ( k ) 2 = X ( k ) 1 R 2 (34) 25 with the resulting predicted change image ∆ ˇ Y ( k ) 2 = Y 2 − ˇ Y ( k ) 2 . (35) The objecti ve function ( 7 ) w .r .t. ∆ X is then re written by combining ( 34 ) and ( 35 ) with ( 7 ), leading to J S 6 , 2 (∆ X ) , J S 6 ( X ( k ) 1 , ∆ X ) =    Λ − 1 2 2  ∆ ˇ Y ( k ) 2 − ∆ XR 2     2 F + γ k ∆ X k 2 , 1 . (36) The minimization of ( 36 ) can be interpreted as a super-resolution problem. Ev en if a forward-backward algorithm could be used to iterati vely minimize this objecti ve func- tion, the size of the spatial degradation matrix R 2 suggests to resort to an ADMM. By introducing the splitting variable W ∈ R n λ × m 2 = ∆ XR 2 , the resulting scaled augmented Lagrangian for the problem is expressed as L µ (∆ X , W , V ) =    Λ − 1 2 2  ∆ ˇ Y ( k ) 2 − W     2 F + λ k ∆ X k 2 , 1 + µ 2 k ∆ XR 1 − W + V k 2 F . (37) Closed-form expressions of the minimizers of ( 37 ) w .r .t. ∆ X and W can be deri ved, following a group soft-thresholding operation and the technique proposed in [ 42 ], re- spectiv ely . 4.7. Scenario S 7 Scenario S 7 generalizes scenario S 4 with the speciﬁc case of a non-inte ger relative spatial downsampling factor , which precludes the use of a unique spatial degradation matrix. The resulting joint observation model is Y 1 = L 1 X 1 R 1 + N 1 . (38a) Y 2 = ( X 1 + ∆ X ) R 2 + N 2 (38b) which leads to the following objecti ve function J S 7 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − ( X 1 + ∆ X ) R 2 )    2 F + 1 2    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . 26 The choice of assuming that the image acquired by the sensor S 2 does not suffers from spectral degradation is moti vated by an easier and more accurate estimation of the change image ∆ X by av oiding additional spectral deblurring steps. The two sub- problems underlying the AM algorithm are detailed below . 4.7.1. Fusion: optimization w .r .t. X 1 By deﬁning the corrected image as for Scenario S 6 , i.e., ˜ Y ( k ) 2 = Y 2 − ∆ X ( k ) R 2 , updating the latent image X 1 consists in minimizing the partial function J S 7 , 1 ( X 1 ) , J S 7  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1 R 2     2 F + λ   X 1 − ¯ X 1   2 F . (39) Unfortunately , it is not possible to deri ve a closed-form solution of the minimizer ( 39 ). As for Scenarios S 5 and S 6 , capitalizing on the conv exity of the objective function, an ADMM strategy is followed. By deﬁning the splitting variable U ∈ R m λ 1 × n = L 1 X 1 . The scaled augmented Lagrangian can be written L µ ( X 1 , U , V ) =    Λ − 1 2 1 ( Y 1 − UR 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − X 1 R 2     2 F + λ   X 1 − ¯ X 1   2 F + µ 2 k L 1 X 1 − U + V k 2 F . (40) Iterativ e minimizations of ( 40 ) w .r .t. both U and X 1 can be conducted efﬁciently . More precisely , optimizing w .r .t. U consists in solving a set of super-resolution problems whose closed-form solutions are gi ven band-by-band in [ 42 ]. Regarding the minimiza- tion w .r .t. X 1 , it consists in solving a ` 2 -penalized super-resolution problem, whose closed-form solution is giv en in [ 38 ]. 4.7.2. Corr ection: optimization w .r .t. ∆ X Since the observation model ( 38b ) related to ∆ X is the same as the one of Scenario S 6 (see ( 30b )), optimizing w .r .t. ∆ X can be achieved thanks to ADMM, as described in paragraph 4.6.2 (spatial super-resolution of the predicted change image). 27 4.8. Scenario S 8 This scenario is similar to the Scenario S 2 described in paragraph 4.2 . It relies on tw o images of same spatial resolution b ut of distinct spectral resolution. Howe ver , contrary to Scenario S 2 , this difference in spectral resolutions cannot be e xpressed with a unique spectral degradation matrix, e.g., due to respecti ve spectral ranges with non-ov erlapping bands. In this case the joint forward observation model is Y 1 = L 1 X 1 + N 1 . (41a) Y 2 = L 2 ( X 1 + ∆ X ) + N 2 . (41b) with the resulting objectiv e function J S 8 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − L 2 ( X 1 + ∆ X ))    2 F + 1 2    Λ − 1 2 1 ( Y 1 − L 1 X 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . The choice of which degradation matrices applies to the change image ∆ X is driven by considering the matrix with larger number of bands, which results in a change image of higher spectral resolution. The associated sub-problems are described in what follo ws. 4.8.1. Fusion: optimization w .r .t. X 1 Similarly to Scenario S 4 , by deﬁning the corrected image as ˜ Y ( k ) 2 = Y 2 − L 2 ∆ X ( t ) , updating the latent image X 1 consists in minimizing J S 8 , 1 ( X 1 ) , J S 8  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − L 1 X 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − L 2 X 1     2 F + λ   X 1 − ¯ X 1   2 F . (42) Minimizing ( 42 ) formulates a joint spectral deblurring problem from an observed im- age Y 1 and a pseudo-observed image ˜ Y ( k ) 2 . Thanks to its quadratic form, this least- square problem can be easily solved. 28 4.8.2. Corr ection: optimization w .r .t. ∆ X The predicted image that would be observed by sensor S 2 at time t 1 can be deﬁned as ˇ Y ( k ) 2 = L 2 X ( k ) 1 with the resulting predicted change image ∆ ˇ Y ( k ) 2 = Y 2 − ˇ Y ( k ) 2 . The objecti ve function ( 7 ) w .r .t. ∆ X is then re written by combining ( 34 ) and ( 35 ) with ( 7 ), leading to J S 8 , 2 (∆ X ) , J S 8 ( X ( k ) 1 , ∆ X ) =    Λ − 1 2 2  ∆ ˇ Y ( k ) 2 − L 2 ∆ X     2 F + γ k ∆ X k 2 , 1 . (43) As for scenario S 4 , minimizing ( 43 ) is a spectral deblurring of the predicted change image ∆ ˇ Y ( k ) 2 , which can be achieved using a forward-backward algorithm as proposed in [ 14 ]. 4.9. Scenario S 9 This scenario generalizes scenario S 4 , but with relativ e spectral responses in volving non-ov erlapping bands. The joint forward observation model is then Y 1 = L 1 X 1 R 1 + N 1 . (44a) Y 2 = L 2 ( X 1 + ∆ X ) + N 2 . (44b) which yields the objectiv e function J S 9 ( X 1 , ∆ X ) = 1 2    Λ − 1 2 2 ( Y 2 − L 2 ( X 1 + ∆ X ))    2 F + 1 2    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F + λ   X 1 − ¯ X 1   2 F + γ k ∆ X k 2 , 1 . Note that the estimated latent and change images are deﬁned at the highest spatial resolution while beneﬁting from the spectral resolutions of both observed images. The choice of assuming that the image acquired by sensor S 2 does not suf fer from spatial 29 degradation has been motiv ated by an easier and accurate estimation of the change image ∆ X by avoiding additional spatial super-resolution steps. The resulting sub- problems in volved in the AM algorithm are detailed belo w . 4.9.1. Fusion: optimization w .r .t. X 1 As for scenarios S 4 and S 8 , the corrected image ˜ Y ( k ) 2 can be deﬁned as ˜ Y ( k ) 2 = Y 2 − L 2 ∆ X ( k ) . Thus, updating the latent image X 1 consists in minimizing J S 9 , 1 ( X 1 ) , J S 9  X 1 , ∆ X ( k )  =    Λ − 1 2 1 ( Y 1 − L 1 X 1 R 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − L 2 X 1     2 F + λ   X 1 − ¯ X 1   2 F . (45) Minimizing ( 45 ) is challenging mainly due to the simultaneous presence of spatial and spectral degradation matrices R 1 and L 2 with an additional spatial de gradation L 1 . Therefore, there is no closed-form solution for this problem, which can be eventu- ally solved thanks to ADMM. By introducing the splitting v ariable U ∈ R m λ × m 1 = X 1 R 1 . The resulting scaled augmented Lagrangian is L µ ( X 1 , U , V ) =    Λ − 1 2 1 ( Y 1 − L 1 U )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − L 2 X 1     2 F + λ   X 1 − ¯ X 1   2 F + µ 2 k X 1 R 1 − U + V k 2 F . (46) Closed-form expression of the minimizers of ( 46 ) w .r .t. X 1 and U can be deri ved, following [ 38 ] and a least-square formulation, respecti vely . 4.9.2. Corr ection: optimization w .r .t. ∆ X As both models ( 44b ) and ( 41b ) are the same, optimizing w .r .t. ∆ X can be achiev ed following the strategy detailed in paragraph 4.8.2 , i.e., by spectrally deblurring a pre- dicted change image ∆ ˇ Y ( k ) 2 thanks to the forward-backward algorithm proposed in [ 14 ]. 4.10. Scenario S 10 This scenario generalizes all the pre vious scenario with the particular difﬁculties of non-ov erlapping bands in the spectral responses and non-integer relati ve spatial do wn- 30 sampling factor of the respective spatial degradations. The joint forward observ ation model is given by ( 9 ), which results in the objectiv e function J S 10 in ( 7 ). Again, as for scenarios S 7 and S 9 , the choice of the spatial and spectral degradations applied to the change image ∆ X should be moti vated by reaching the highest spatial and spectral resolutions of this change image. The optimization sub-problems are ﬁnally discussed below . 4.10.1. Fusion: optimization w .r .t. X 1 For this scenario, the corrected image ˜ Y ( k ) 2 is gi ven by ( 10 ), leading to an updating rule of the X 1 consists in minimizing ( 12 ). This minimization cannot be conducted in a straightforward manner , since it requires to conduct a spectral deblurring and a spatial super-resolution simultaneously . Howe ver , the optimal solution can be reached by resorting to a ADMM with two splitting variables U 1 = L 1 X 1 ∈ R m λ 1 × n and U 2 = X 1 R 2 ∈ R n λ × m 2 . The resulting scaled augmented Lagrangian for the problem is expressed as L µ ( X 1 , U 1 , U 2 , V 1 , V 2 ) =    Λ − 1 2 1 ( Y 1 − U 1 R 1 )    2 F +    Λ − 1 2 2  ˜ Y ( k ) 2 − L 2 U 2     2 F + λ   X 1 − ¯ X 1   2 F + µ 2 k L 1 X 1 − U 1 + V 1 k 2 F + µ 2 k X 1 R 2 − U 2 + V 2 k 2 F . (47) Closed-form expressions of the minimizers of ( 47 ) w .r .t. X 1 , U 1 and U 2 can be deriv ed as proposed in [ 38 ], [ 42 ] and following a least-square formulation, respecti vely . 4.10.2. Corr ection: optimization w .r .t. ∆ X For this scenario, gi ven the current state X ( k ) 1 of the latent image, the predicted im- age that would be observed by the sensor S 2 at time t 1 can be deﬁned as in ( 13 ) leading to the predicted change image ( 14 ). Then, the correction step consists in minimizing the objectiv e function J S 10 , 1 (∆ X ) in ( 16 ). It consists in conducting a spectral deblur- ring and spatial super-resolution jointly . This problem has no closed-form solution. Therefore, the objective function is iteratively minimized using an ADMM with two splitting variables W 1 ∈ R m λ 1 × n = L 1 ∆ X and W 2 ∈ R n λ × n = ∆ X . The resulting 31 scaled augmented Lagrangian for the problem is expressed as L µ (∆ X , W 1 , W 2 , V 1 , V 2 ) =    Λ − 1 2 2  ∆ ˇ Y ( k ) 2 − W 1 R 2     2 F + γ k W 2 k 2 , 1 + µ 2 k L 1 ∆ X − W 1 + V 1 k 2 F + µ 2 k ∆ X − W 2 + V 2 k 2 F . (48) Closed-form expression of the minimizers of ( 48 ) w .r .t. ∆ X , W 1 and W 2 can be deriv ed, following a least-square formulation, the computation proposed in [ 42 ] and a group soft-thresholding, respectiv ely . 5. Experiments 5.1. Refer ence images T o illustrate the performance of the proposed algorithmic framework using real multi-band optical data on each speciﬁc scenario discussed in paragraph 4 , observed images from 4 largely studied open access multi-band sensors have been chosen, namely Landsat-8 from [ 35 ], Sentinel-2 from [ 12 ], Earth observing-1 Advanced Land Imager (EO-1 ALI) [ 34 ] and Airborne V isible Infrared Imaging Spectrometer (A VIRIS) from [ 20 ]. These images hav e been acquired over the same geographical location, i.e., the Mud Lake region in Lake T ahoe, CA, USA between June 8th, 2011 and October 29th, 2016. Unfortunately , no ground truth information is av ailable for the chosen image pairs, as e xperienced in numerous experimental situations [ 4 ]. Howe ver , this re gion is characterized by interesting natural meteorological changes, e.g., drought of the Mud Lake, snow falls and vegetation gro wth, occurring along the seasons which help to vi- sually infer the major changes between two dates and to assess the relev ance of the detected changes. All considered images have been manually geographically and geo- metrically aligned to fulﬁll the requirements imposed by the considered CD setup. In addition to the data pro vided by these sensors, complementary images have been synthetically generated by considering so-called virtual sensors derived from the real ones. The speciﬁcations of these virtual sensors, summarized in Figure 1 , are chosen such that all applicativ e scenarios pre viously discussed can be div ersely represented. 32 They are met by selecting a subset of the initial spectral bands or by artiﬁcially degrad- ing the spatial resolution of the real sensors. 300 500 700 900 1100 1300 1500 1700 1900 2100 2300 2500 Wavelenght (nm) HS-29 AVIRIS HS-224 MS-3 MS-6 Sentinel-2 MS-4 MS-3 MS-9 EO-1 ALI PAN MS-3 MS-5 MS-8 Landsat 8 PAN 15m 15m 10m 20m 10m 30m 10m 30m 30m 30m 15m 30m Spatial Resolution Figure 1: Spectral and spatial characteristics of real (green) and virtual (red) sensors. 5.1.1. Landsat-8 images Landsat-8 is the eighth Earth observ ation satellite series of the US LANDSA T Program [ 35 ], launched on February 11th, 2013 with 16-days revisiting period. It is equipped with the Operational Land Imager (OLI) and the Thermal InfraRed Sen- sor (TIRS). In the conducted experiments, 3 sets of real images acquired at the dates 10/18/2013, 04/15/2015 and 09/22/2015 hav e been considered. For each acquisition, Landsat-8 provides • one panchromatic image ov er the spectral range 0 . 503 – 0 . 676 µ m (band ] 8 ) at a 15 m spatial resolution (denoted P AN), • one multispectral image of 8 spectral bands (bands ] 1 – ] 7 and ] 9 ) at a 30 m reso- lution (denoted MS- 8 ). For experimental purpose, as explained above, these real images are complemented with the following virtually acquired images 33 • one multispectral image of 5 spectral bands (bands ] 1 – ] 4 and ] 7 ) at a 30 m spatial resolution (denoted MS- 5 ), • one red-green-blue (RGB) multispectral image of 3 spectral bands (bands ] 2 – ] 4 ) at a 30 m spatial resolution (denoted MS- 3 ). 5.1.2. Sentinel-2 images Sentinel-2 is a series of two identical satellites for Earth observation missions de- veloped by ESA [ 12 ] as part of the Copernicus Program launched in 2015 and 2017 with 5-days revisiting period. The multi-spectral instrument embedded on each plat- form is composed of two different sensors for acquisition in the visible and infrared spectral domains, respectively . The actual dataset used in the experiments is composed of tw o images acquired on 04/12/2016 and 10/29/2016 and, for each real scene, among all av ailable spectral bands, one considers • one multispectral image of 4 visible/near infrared (VNIR) spectral bands (bands ] 2 – ] 4 and ] 8 ) at a 10 m spatial resolution (denoted MS-4) • one multispectral image of 6 short wav e infrared spectral range (SWIR) spectral bands (bands ] 5 – ] 8 a and ] 11 – ] 12 ) at a 20 m spatial resolution (denotes MS- 6 ) and one additional virtually image, namely , • one RGB multispectral image of 3 spectral bands (bands ] 2 – ] 4 ) at a 10 m spatial resolution (denoted MS- 3 ). 5.1.3. EO-1 ALI images Operated by N ASA, EO-1 ALI is a Earth observation satellite part of the Ne w Mil- lennium Program launched in 2000 with 16-days repeat cycle and decommissioned in 2017 [ 34 ]. The main embedded sensor Adv anced Land Imager (ALI) is comple- mented with the Hyperion spectrometer and the Linear Etalon Imaging Spectrometer Array (LEISA) for atmospheric correction. The considered dataset corresponds to 2 acquisition dates, 06/08/2011 and 08/04/2011, for 34 • one panchromatic image over the spectral range 0 . 48 – 0 . 69 µ m (band ] 1 ) at a 10 m spatial resolution (denoted P AN), • one multispectral image of 9 spectral bands (bands ] 2 – ] 10 ) at a 30 m resolution (denoted MS- 9 ), in addition to the virtual acquisition of • one RGB multispectral image of 3 spectral bands (bands ] 3 – ] 5 ) at a 30 m spatial resolution (denoted MS- 3 ). 5.1.4. A VIRIS images A VIRIS is the second aircraft embedding an image spectrometer de veloped by Jet Propulsion Laboratory (JPL) for Earth remote sensing [ 20 ]. It deli vers calibrated im- ages in 224 contiguous 10 nm-width spectral channels ranging from 0 . 4 µ m to 2 . 5 µ m. Since it is an airborne-dependent system, the spatial resolution is not a priori ﬁxed and is designed for each individual acquisition. The dataset considered in the conducted experiments is composed by two real images acquired on 04/10/2014 and 09/19/2014. For each scene, one considers • the original h yperspectral image of 224 spectral bands at a 15 m spatial resolution (denoted HS- 224 ) • one virtual hyperspectral image of 29 spectral bands (corresponding to the RGB domain) at a 15 m spatial resolution (denoted HS- 29 ) 5.2. Design of the spatial and spectral degr adations The proposed model requires the prior knowledge of spectral and spatial degrada- tion matrices L and R = BS , respecti vely . Reg arding the spectral degradation matri- ces required in each simulation scenario, they can be easily deriv ed from the intrinsic sensor characteristics freely a vailable by av eraging the spectral bands corresponding to the prescribed response. Conv ersely , the spatial de gradation is not a sensor speciﬁca- tion. It depends not only on the considered systems as well as external factors b ut also on the targeted resolution of the fused image. This work relies on commonly adopted 35 assumptions by considering R as a Gaussian blur and by adjusting the downsampling factor in S as an inte ger value corresponding to the relative ratio between spatial reso- lution of both observed images. 5.3. Compar ed methods As previously exposed, the proposed rob ust fusion-based CD framework (referred to as RF) is able to deal with all combinations of mono- and multi-band optical images of different spatial and spectral resolutions. Howe ver , up to author’ s kno wledge, there is no technique in the literature with such a versatility , i.e., able to address all these scenarios. For this reason, the technique referred to as the worst-case (WC) and also considered in [ 14 ] has been used as a baseline and state-of-the-art CD technique. This WC method consists in preprocessing the observed images by spatially and/or spec- trally degrading them in order to reach a set of observed images of the same spectral and spatial resolutions. Then, when handling images of same resolutions, classical CD technique, e.g., CV A proposed in [ 22 ], can be easily conducted to build a lo w spatial resolution change mask denoted ˆ d WC . 5.4. Results The follo wing paragraphs discuss the CD performance of the proposed RF method and of the WC approach for each applicativ e scenario detailed in paragraph 4 (see also T able 1 ). Depending on the considered scenario, pairs of real and/or virtual images described in paragraph 5.1 are selected, beneﬁting from different acquisition times but common acquisition location. T able 3 summarizes the pair of observ ed images provided by the real and/or virtual sensors used in each scenario. Note that several combinations of images can be made for Scenarios S 1 – S 5 . 5.4.1. Scenario S 1 In the ﬁrst scenario, CD is conducted on a pair of images of same spatial and spec- tral resolutions, which corresponds to the most fa vorable and commonly considered CD frame work. Figures 2 to 4 present the CD binary masks reco vered by the proposed RF-based CD method as well as by the WC method for three pairs of panchromatic, 36 Image ] 1 Image ] 2 Sensor Spatial resol. Spectral resol. Sensor Spatial resol. Spectral resol. S 1 Landsat-8 15 P AN Landsat-8 15 P AN Landsat-8 30 MS- 3 Landsat-8 30 MS- 3 A VIRIS 15 HS- 224 A VIRIS 15 HS- 224 S 2 EO-1 ALI 10 P AN Sentinel-2 10 MS- 3 Landsat-8 15 P AN A VIRIS 15 HS- 29 S 3 Sentinel-2 10 MS- 3 EO-1 ALI 30 MS- 3 Sentinel-2 10 MS- 3 Landsat-8 30 MS- 3 S 4 Landsat-8 15 P AN Landsat-8 30 MS- 3 EO-1 ALI 10 P AN Landsat-8 30 MS- 3 Landsat-8 15 P AN EO-1 ALI 30 MS- 3 S 5 EO-1 ALI 30 MS- 3 A VIRIS 15 HS- 29 Landsat-8 30 MS- 3 A VIRIS 15 HS- 29 S 6 EO-1 ALI 10 P AN Landsat-8 15 P AN S 7 Sentinel-2 10 MS- 3 Landsat-8 15 P AN S 8 Landsat-8 30 MS- 8 EO-1 ALI 30 MS- 9 S 9 Landsat-8 30 MS- 5 Sentinel-2 10 MS-4 S 10 Sentinel-2 20 MS- 6 EO-1 ALI 30 MS- 9 T able 3: Pairs of real and/or virtual images, and their spatial and spectral characteristics, used for each applicativ e scenario. multispectral and hyperspectral images, respectiv ely . Note that, in this scenario, the WC boils down to conduct CV A directly on the observed images since they already share the same spatial and spectral resolutions and, thus, do not require to be degraded to be able to conduct pixel wise comparison. These CD maps sho w that both CD meth- ods detect the most signiﬁcant changes, in particular the draught of the lak e. Howe ver , for all conﬁgurations, the proposed method visually present CD maps with better detec- tion/false alarm rates when compared with the WC method. This can be explained by the f act that the proposed method denoises the observed image while jointly estimating 37 the change image ∆ X . Con versely , the WC method directly uses the observ ed images to deriv e the change image, which may suffer from noise introducing false alarms and misdetections. This is particularly visible in Fig. 4 depicting the results obtained from a hyperspectral image, kno wn to be of lo wer signal-to-noise ratio. (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 2: Scenario S 1 : (a) Landsat-8 15 m P AN observed image Y 1 acquired on 04/15/2015, (b) Landsat-8 15 m P AN observ ed image Y 2 acquired on 09/22/2015, (c) change mask ˆ d WC estimated by the WC method and (d) change mask ˆ d RF estimated by the proposed approach. 5.4.2. Scenario S 2 This CD scenario deals with observ ed images of same spatial resolution b ut dif fer- ent spectral resolutions. Figures 5 and 6 illustrate two possible situations and sho w the CD results of the proposed RF-based CD method compared with the WC method. In 38 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 3: Scenario S 1 : (a) Landsat-8 30 m MS- 3 observed image Y 1 acquired on 04/15/2015, (b) Landsat-8 30 m MS- 3 observed image Y 2 acquired on 09/22/2015, (c) change mask ˆ d WC estimated by the WC method and (d) change mask ˆ d RF estimated by the proposed approach. this scenario, similarly to scenario S 1 , both estimated CD maps have the same spatial resolution as the observed image pair , which means that there is not any loss of spatial resolution. On the other hand, the proposed method deliv ers a change map estimated from ∆ X of same spectral resolution as the highest spectral resolution among the two observed images. Con versely , the WC method conducts CV A on a pair of images after spectral degradation to reach the lowest spectral resolution, which possibly results in loss of signiﬁcant information. The consequent impact on the change/no-change deci- sion is the visual reduction of false alarm rate for the proposed RF method, e ven if both 39 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 4: Scenario S 1 : (a) A VIRIS 15 m HS- 224 observed image Y 1 acquired on 04/10/2014, (b) A VIRIS 15 m HS- 224 observed image Y 2 acquired on 09/19/2014, (c) change mask ˆ d WC estimated by the WC method and (d) change mask ˆ d RF estimated by the proposed approach. CD maps hav e the same spatial resolution. 5.4.3. Scenario S 3 In scenario S 3 , corresponding to the re verse situation encountered in scenario S 2 , observed images share the same spectral resolution but with different spatial resolution. Figures 7 and 8 present the results obtained for two possible real situations. Note that CD maps obtained by the proposed RF method are of higher spatial resolutions than the ones estimated by the WC approach. Thus, this scenario is the ﬁrst to illustrate the most important dif ferences between both approaches, i.e., the difference in spatial 40 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 5: Scenario S 2 : (a) EO-1 ALI 10 m P AN observed image Y 1 acquired on 06/08/2011, (b) Sentinel- 2 10 m MS- 3 observed image Y 2 acquired on 04/12/2016, (c) change mask ˆ d WC estimated by the WC approach from a pair of 10 m P AN degraded images and (d) change mask ˆ d RF estimated by the proposed approach from a 10 m MS- 3 change image ∆ ˆ X . resolutions of the CD maps. In scenario S 2 , the results have already shown that the loss of spectral information inherent to the WC approach leads to an increase of false alarms and misdetections. Here, the loss of spatial information when conducting the WC approach results in an inaccurate localization of the possible changes. 5.4.4. Scenario S 4 This scenario has been deeply inv estigated in [ 14 ] who conducted a comprehen- siv e analysis of the performance of the proposed RF-based CD method. This scenario 41 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 6: Scenario S 2 : (a) Landsat-8 15 m P AN observed image Y 1 acquired on 09/22/2015, (b) A VIRIS 15 m HS- 29 observed image Y 2 acquired on 04/10/2014, (c) change mask ˆ d WC estimated by the WC approach from a pair of 15 m P AN degraded images and (d) change mask ˆ d RF estimated by the proposed approach from a 15 m HS- 29 change image ∆ ˆ X . corresponds to a more difﬁcult CD in vestigation than all previous ones since the pair of observed images ha ve not the same spatial neither spectral resolutions. As a conse- quence, the conv entional WC approach is constrained to compare a spatially de graded version of one observed image with a spectrally degraded v ersion of the other observed image. Irredeemably , these degradations result in a loss of spectral information, es- sential to assess the presence of change, and a loss of spatial information, required to accurately localize the possible changes. On the contrary , the proposed method is able 42 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 7: Scenario S 3 : (a) Sentinel-2 10 m MS- 3 observed image Y 1 acquired on 10/29/2016, (b) EO-1 ALI 30 m MS- 3 observed image Y 2 acquired on 08/04/2011, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m MS- 3 de graded images and (d) change mask ˆ d RF estimated by the proposed approach from a 10 m MS- 3 change image ∆ ˆ X . to deriv e the CD mask from a change image characterized by the best of the spectral and spatial resolution of the observed images. Figures 9 to 11 depict the CD results obtained for three common conﬁgurations and illustrate the superiority of the proposed RF-based CD method. 5.4.5. Scenario S 5 As in the previous case, this scenario handles images which do not share the same spatial neither spectral resolutions. Howe ver , contrary to scenario S 4 , this scenario 43 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 8: Scenario S 3 : (a) Sentinel-2 10 m MS- 3 observed image Y 1 acquired on 04/12/2016, (b) Landsat- 8 30 m MS- 3 observed image Y 2 acquired on 09/22/2015, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m MS- 3 de graded images and (d) change mask ˆ d RF estimated by the proposed approach from a 10 m MS- 3 change image ∆ ˆ X . considers one of the two images of higher spatial and spectral resolution. Again, the WC is expected to be less reliable (in terms of decision and localization) due to the loss of spectral and spatial information consecutive to the de gradations before conducting CV A. Figures 12 and 13 present the results obtained from tw o possible real conﬁgura- tions. As e xpected the proposed RF-based CD method pro vides visually more satisfac- tory results. In particular , as shown in Fig. 13 , the WC method is unable to accurately localize the change due to lake draught from the pair of multispectral and h yperspectral 44 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 9: Scenario S 4 : (a) Landsat-8 15 m P AN observed image Y 1 acquired on 09/22/2015, (b) Landsat- 8 30 m MS- 3 observed image Y 2 acquired on 04/15/2015, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m P AN degraded images and (d) change mask ˆ d RF estimated by the proposed approach from a 15 m MS- 3 change image ∆ ˆ X . images. 5.4.6. Scenario S 6 This scenario represents a particular instance of scenario S 3 , i.e., with two observed images of dif ferent spatial resolution b ut same spectral resolution. Nev ertheless, here, the two spatial resolutions are related by a non-inte ger do wnsampling ratio which pre- cludes the use of a unique spatial degradation matrix in the proposed RF-based CD method. As detailed in paragraph 4.6 , super-resolutions are conducted during the fu- 45 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 10: Scenario S 4 : (a) EO-1 ALI 10 m P AN observed image Y 1 acquired on 06/08/2011, (b) Landsat-8 30 m MS- 3 observed image Y 2 acquired on 09/22/2015, (c) change mask ˆ d WC estimated by the WC ap- proach from a pair of 30 m P AN degraded images (d) change mask ˆ d RF estimated by the proposed approach from a 10 m MS- 3 change image ∆ ˆ X . sion and correction steps of the AM algorithm, which leads to a change image ∆ ˆ X with a spatial resolution higher than the ones of the two observed images (deﬁned as the greatest common divisor of the resolutions). For instance, Fig. 14 illustrates one possible conﬁguration for which the observed images Y 1 and Y 2 , depicted in Fig. 14(a) and 14(b) , are of 15 m and 10 m spatial resolutions, respecti vely . Thus the change image ∆ ˆ X and change mask ˆ d RF estimated by the proposed method are at a 5 m reso- lution. Con versely , the WC method pro vides a CD map at a spatial resolution based on 46 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 11: Scenario S 4 : (a) Landsat-8 15 m P AN observed image Y 1 acquired on 09/22/2015, (b) EO- 1 ALI 30 m MS- 3 observed image Y 2 acquired on 06/08/2011, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m P AN degraded images (d) change mask ˆ d RF estimated by the proposed approach from a 15 m MS- 3 change image ∆ ˆ X . the least common multiple, which is, in this case, 30 m. The signiﬁcantly higher spatial resolution of the change map is clear in Fig. 14(d) . 5.4.7. Scenario S 7 This scenario consists in a more challenging context than scenario S 6 since, in ad- dition to the non-integer relative do wnsampling factor , the tw o observed images do not share the same spectral resolution. As before, the change image ∆ ˆ X and the bi- nary change mask ˆ d RF estimated by the proposed RF-based CD method are deﬁned 47 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 12: Scenario S 5 : (a) EO-1 ALI 30 m MS- 3 observed image Y 1 acquired on 08/04/2011, (b) A VIRIS 15 m HS- 29 observed image Y 2 acquired on 04/10/2014, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m MS- 3 degraded images (d) change mask ˆ d RF estimated by the proposed approach from a 15 m HS- 29 change image ∆ ˆ X . at a higher spatial resolution than both resolutions of the observed image. Figure 15 presents one example of this scenario. As e xpected, the proposed method beneﬁts from the estimated highest spatial and spectral resolution change image ∆ X to localize the changes, contrary to the WC method which can only exploit a pair of spatially and spectrally degrade images. This important dual resolution gap contributes a lot in the observed dif ferences on the false alarm and good detection rates. 48 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 13: Scenario S 5 : (a) Landsat-8 30 m MS- 3 observ ed image Y 1 acquired on 04/15/2015, (b) A VIRIS 15 m HS- 29 observed image Y 2 acquired on 09/19/2014, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m MS- 3 de graded images and (d) change mask ˆ d RF estimated by the proposed approach from a 15 m HS- 29 change image ∆ ˆ X . 5.4.8. Scenario S 8 Scenario S 8 generalizes scenario S 2 , with the particular case of presence of non- ov erlapping bands in the two sensor spectral responses, which requires the simultane- ous use of two spectral degradation matrices in the proposed RF method. Figure 16 provides one instance of this scenario. Due to the presence of non-ov erlapping bands, before conducting CV A, the WC requires to ignore the spectral bands which are not commonly shared by the two observed images. Con versely , by fully exploiting the 49 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 14: Scenario S 6 : (a) Landsat-8 15 m P AN observed image Y 1 acquired on 10/18/2013, (b) EO- 1 ALI 10 m P AN observed image Y 2 acquired on 08/04/2011, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m P AN degraded images (d) change mask ˆ d RF estimated by the proposed approach from 5 m P AN change image ∆ ˆ X . whole av ailable spectral information, the proposed method combines the overlapped bands and the non-overlapping bands to estimate a change image ∆ ˆ X of higher spec- tral resolution than the two observed images. This higher amount of information leads to visually more consistent results in Fig. 16(d) . 5.4.9. Scenario S 9 This scenario correspond to a modiﬁed instance of scenario S 4 (images of different spatial and spectral resolutions) with some non-ov erlapping bands (as for the pre vious 50 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 15: Scenario S 7 : (a) Sentinel-2 10 m MS- 3 observed image Y 1 acquired on 04/12/2016, (b) Landsat- 8 15 m P AN observed image Y 2 acquired on 09/22/2015, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m P AN degraded images and (d) change mask ˆ d RF estimated by the proposed approach from 5 m MS- 3 change image ∆ ˆ X . scenario). The results obtained for one conﬁguration are depicted in Figure 17 . In this case, the change image ∆ ˆ X is characterized by a spatial resolution deﬁned by the highest one among the observ ed image with a spectral resolution higher than both observed images. Once again, the results show the accuracy of the proposed method in terms of detection and spatial resolution of the estimated change map. 51 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 16: Scenario S 8 : (a) Landsat-8 30 m MS- 8 observed image Y 1 acquired on 04/15/2015, (b) EO- 1 ALI 30 m MS- 9 observed image Y 2 acquired on 06/08/2011, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m MS- 7 degraded images (d) change mask ˆ d RF estimated by the proposed approach from 30 m MS- 10 change image ∆ ˆ X . 5.4.10. Scenario S 10 The last scenario combines of the difﬁculties previously encountered: images of different spatial and spectral resolution, characterized by a non-inte ger relati ve do wn- sampling factor and non-ov erlapping spectral bands. As for scenarios S 6 and S 7 , the change image ∆ ˆ X and change mask ˆ d RF recov ered by the proposed RF-based CD method is of higher spatial resolution than the tw o observ ed images. In addition, as for scenario S 8 and S 9 , the change image is also deﬁned at a higher spectral resolution. 52 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 17: Scenario S 9 : (a) Landsat-8 30 m MS- 5 observed image Y 1 acquired on 09/22/2015, (b) Sentinel- 2 10 m MS- 4 observed image Y 2 acquired on 04/12/2016, (c) change mask ˆ d WC estimated by the WC approach from a pair of 30 m MS- 3 de graded images and (d) change mask ˆ d RF estimated by the proposed approach from a 10 m MS- 6 change image ∆ ˆ X . Con versely , the WC approach deri ves a change image of lower spatial and spectral res- olutions before conducting CV A. Figure 18 depicts the results obtained by both meth- ods. On this particularly challenging scenario, the proposed approach demonstrates its superiority in recov ering rele vant changes and in localizing them accurately . 53 (a) Y 1 (b) Y 2 (c) ˆ d WC (d) ˆ d RF Figure 18: Scenario S 10 : (a) Sentinel-2 20 m MS- 6 observed image Y 1 acquired on 04/12/2016, (b) EO-1 ALI 30 m MS- 9 observed image Y 2 acquired on 06/08/2011, (c) change mask ˆ d WC estimated by the WC approach from a pair of 60 m MS- 4 de graded images and (d) change mask ˆ d RF estimated by the proposed approach from a 10 m MS- 11 change image ∆ ˆ X . 6. Conclusion This paper deri ved a robust fusion frame work to perform change detection between optical images of dif ferent spatial and spectral resolutions. The versatility of the pro- posed approach allowed all possible real scenarios to be handled ef ﬁciently . The tech- nique was based on the deﬁnition of two high spatial and spectral resolution latent images related to the observed images via a double physically-inspired forward model. The difference between these two latent images was assumed to be spatially sparse, 54 implicitly locating the changes at a high resolution scale. Inferring these two latent images w as formulated as an in verse problem which was solved within a 2 -step itera- tiv e scheme. Depending on the considered scenario, these 2 steps can be interpreted as ubiquitous signal & image processing problems (namely spatial super -resolution, spec- tral deblurring, denoising or multi-band image fusion) for which closed-form solutions or efﬁcient algorithms had been recently proposed in the literature. Real images ac- quired by four different sensors were used to illustrate the accurac y and the ﬂexibility of the proposed method, as well as its superiority with respect to the state-of-the-art change detection methods. Future works will assess the robustness of the proposed technique w .r .t. nonlinear effects (e.g., due to atmospheric effects, geometric and ra- diometric distortions). Detecting changes between optical and non-optical data is also under in vestigation. Appendix A. Matrix normal distribution The probability density function p ( X | M , Σ r , Σ r ) of a matrix normal distribution MN r,c ( M , Σ r , Σ c ) is giv en in [ 18 ] p ( X | M , Σ r , Σ r ) = exp ( − 1 2 tr [ Σ − 1 c ( X − M ) T Σ − 1 r ( X − M ) ]) (2 π ) rc/ 2 | Σ c | r/ 2 | Σ r | c/ 2 where M ∈ R r × c is the mean matrix, Σ r ∈ R r × r is the row cov ariance matrix and Σ c ∈ R c × c is the column cov ariance matrix. Appendix B. Acknowledgments Part of this w ork has been supported by Coordenação de Aperfeiçoamento de En- sino Superior (CAPES), Brazil, and EU FP7 through the ERANETMED JC-W A TER program [MapIn vPlnt Project ANR-15-NMED- 0002-02]. References [1] Bell, T . E., 1995. Remote sensing. IEEE Spectrum 32 (3), 24–31. 55 [2] Bov olo, F ., Bruzzone, L., Jan. 2007. A theoretical framework for unsupervised change detection based on change vector analysis in the polar domain. IEEE T rans. Geosci. Remote Sens. 45 (1), 218–236. [3] Bov olo, F ., Bruzzone, L., Sept. 2015. The time variable in data fusion: A change detection perspectiv e. IEEE Geosci. Remote Sens. Mag. 3 (3), 8–26. [4] Bov olo, F ., Marchesi, S., Bruzzone, L., June 2012. A frame work for automatic and unsupervised detection of multiple changes in multitemporal images. IEEE T rans. Geosci. Remote Sens. 50 (6), 2196–2212. [5] Campbell, J. B., W ynne, R. H., 2011. Introduction to remote sensing, 5th Edition. Guilford Press, New Y ork. [6] Candés, E. J., Li, X., Ma, Y ., Wright, J., 2011. Robust principal component anal- ysis? Journal of the A CM (J A CM) 58 (3), 11. [7] Dalla Mura, M., Prasad, S., Paciﬁci, F ., Gamba, P ., Chanussot, J., Benediktsson, J. A., Sept. 2015. Challenges and Opportunities of Multimodality and Data Fusion in Remote Sensing. Proc. IEEE 103 (9), 1585–1601. [8] Du, P ., Liu, S., Xia, J., Zhao, Y ., 2013. Information fusion techniques for change detection from multi-temporal remote sensing images. Information Fusion 14 (1), 19–27. [9] Eismann, M. T ., Hardie, R. C., March 2005. Hyperspectral resolution enhance- ment using high-resolution multispectral imagery with arbitrary response func- tions. IEEE T rans. Image Process. 43 (3), 455–465. [10] Elachi, C., V an Zyl, J., 2006. Introduction to the physics and techniques of re- mote sensing, 2nd Edition. Wile y series in remote sensing. W iley-Interscience, Hoboken, N.J. [11] Elad, M., Feuer, A., 1997. Restoration of a single superresolution image from sev eral blurred, noisy , and undersampled measured images. IEEE Trans. Image Process. 6 (12), 1646–1658. 56 [12] European Space Agency, 2017. Sentinel-2. URL https://sentinel.esa.int/web/sentinel/missions/ sentinel- 2 [13] Ferraris, V ., Dobigeon, N., W ei, Q., Chabert, M., March 2017. Change detection between multi-band images using a robust fusion-based approach. In: 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). pp. 3346–3350. [14] Ferraris, V ., Dobigeon, N., W ei, Q., Chabert, M., June 2017. Robust fusion of multiband images with dif ferent spatial and spectral resolutions for change detec- tion. IEEE T rans. Computational Imaging 3 (2), 175–186. [15] Ferraris, V ., Dobigeon, N., W ei, Q., Chabert, M., March 2018. Detecting changes between optical images of different spatial and spectral resolutions: A fusion- based approach. IEEE T rans. Geosci. Remote Sens. 56 (3), 1566–1578. [16] Févotte, C., Dobigeon, N., 2015. Nonlinear hyperspectral unmixing with robust nonnegati ve matrix factorization. IEEE Trans. Image Process. 24 (12), 4810– 4819. [17] Ghassemian, H., 2016. A re view of remote sensing image fusion methods. Infor - mation Fusion 32 (Part A), 75–89. [18] Gupta, A. K., Nagar , D. K., 1999. Matrix V ariate Distribution. No. 104 in Mono- graphs and Surve ys in Pure and Applied Mathematics. Chapman and Hall. [19] Hardie, R. C., Eismann, M. T ., Wilson, G. L., Sept. 2004. MAP estimation for hy- perspectral image resolution enhancement using an auxiliary sensor . IEEE T rans. Image Process. 13 (9), 1174–1184. [20] Jet Propulsion Laboratory, 2017. Airborne visible / infrared imaging spectrometer (A VIRIS). URL https://aviris.jpl.nasa.gov 57 [21] Jianchao Y ang, Wright, J., Huang, T . S., Y i Ma, Nov . 2010. Image super- resolution via sparse representation. IEEE Trans. Image Process. 19 (11), 2861– 2873. [22] Johnson, R. D., Kasischke, E. S., Jan. 1998. Change vector analysis: A technique for the multispectral monitoring of land cov er and condition. Int. J. Remote Sens. 19 (3), 411–426. [23] K otwal, K., Chaudhuri, S., 2013. A bayesian approach to visualization-oriented hyperspectral image fusion. Information Fusion 14 (4), 349–360. [24] K otwal, K., Chaudhuri, S., 2013. A nov el approach to quantitati ve ev aluation of hyperspectral image fusion techniques. Information Fusion 14 (1), 5–18. [25] Li, S., Kang, X., Fang, L., Hu, J., Y in, H., 2017. Pixel-lev el image fusion: A surve y of the state of the art. Information Fusion 33 (Supplement C), 100–112. [26] Loncan, L., de Almeida, L. B., Bioucas-Dias, J. M., Briottet, X., Chanussot, J., Dobigeon, N., Fabre, S., Liao, W ., Licciardi, G. A., Simoes, M., T ourneret, J.-Y ., V e ganzones, M. A., V ivone, G., W ei, Q., Y okoya, N., Sept. 2015. Hyperspectral pansharpening: A revie w . IEEE Geosci. Remote Sens. Mag. 3 (3), 27–46. [27] Nielsen, A. A., Feb . 2007. The Re gularized Iterativ ely Re weighted MAD Method for Change Detection in Multi- and Hyperspectral Data. IEEE Trans. Image Pro- cess. 16 (2), 463–478. [28] Nielsen, A. A., Conradsen, K., Simpson, J. J., 1998. Multiv ariate alteration de- tection (MAD) and MAF postprocessing in multispectral, bitemporal image data: New approaches to change detection studies. Remote Sens. En vironment 64 (1), 1–19. [29] Richards, J. A., Jia, X., 2006. Remote sensing digital image analysis: an intro- duction, 4th Edition. Springer , Berlin. [30] Simões, M., Bioucas Dias, J., Almeida, L., Chanussot, J., June 2015. A conv ex formulation for hyperspectral image superresolution via subspace-based re gular- ization. IEEE T rans. Geosci. Remote Sens. 6 (53), 3373–3388. 58 [31] Singh, A., June 1989. Revie w Article Digital change detection techniques using remotely-sensed data. Int. J. Remote Sens. 10 (6), 989–1003. [32] Song, H., Huang, B., Zhang, K., Zhang, H., 2014. Spatio-spectral fusion of satel- lite images based on dictionary-pair learning. Information Fusion 18 (Supplement C), 148–160. [33] Union of Concerned Scientists, 2017. UCS Satellite Database. URL https://www.ucsusa.org/nuclear- weapons/ space- weapons/satellite- database#.WigLykqnFPY [34] United States Geological Surve y, 2017. EO-1 Adv anced Land Imager (ALI). URL https://eo1.usgs.gov/sensors/ali [35] United States Geological Surve y, 2017. Landsat-8. URL https://landsat.usgs.gov/landsat- 8 [36] W ei, Q., Bioucas-Dias, J., Dobigeon, N., T ourneret, J.-Y ., 2015. Hyperspectral and multispectral image fusion based on a sparse representation. IEEE T rans. Geosci. Remote Sens. 53 (7), 3658–3668. [37] W ei, Q., Dobigeon, N., T ourneret, J.-Y ., Sept. 2015. Bayesian Fusion of Multi- Band Images. IEEE J. Sel. T opics Signal Process. 9 (6), 1117–1127. [38] W ei, Q., Dobigeon, N., T ourneret, J.-Y ., Nov . 2015. Fast Fusion of Multi-Band Images Based on Solving a Sylvester Equation. IEEE Trans. Image Process. 24 (11), 4109–4121. [39] Y oko ya, N., Mayumi, N., Iw asaki, A., April 2013. Cross-Calibration for Data Fusion of EO-1/Hyperion and T erra/ASTER. IEEE J. Sel. T opics Appl. Earth Observations Remote Sens. 6 (2), 419–426. [40] Y oko ya, N., Y airi, T ., Iwasaki, A., Feb. 2012. Coupled nonnegati ve matrix fac- torization unmixing for hyperspectral and multispectral data fusion. IEEE T rans. Geosci. Remote Sens. 50 (2), 528–537. 59 [41] Zhang, Y ., De Backer , S., Scheunders, P ., Nov . 2009. Noise-resistant w av elet- based Bayesian fusion of multispectral and hyperspectral images. IEEE Trans. Geosci. Remote Sens. 47 (11), 3834–3843. [42] Zhao, N., W ei, Q., Basarab, A., Dobigeon, N., K ouame, D., T ourneret, J.-Y ., Aug. 2016. Fast Single Image Super-Resolution Using a New Analytical Solution for ` 2 – ` 2 Problems. IEEE T rans. Image Process. 25 (8), 3683–3697. 60

Robust fusion algorithms for unsupervised change detection between multi-band optical images - A comprehensive case study

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment