Robust Fusion of Multi-Band Images with Different Spatial and Spectral Resolutions for Change Detection
11 Robust fusion of multi-band images withdifferent spatial and spectral resolutions forchange detection
Vinicius Ferraris, Nicolas Dobigeon,
Senior Member, IEEE ,Qi Wei,
Member, IEEE , and Marie Chabert
Abstract
Archetypal scenarios for change detection generally consider two images acquired through sen-sors of the same modality. However, in some specific cases such as emergency situations, the onlyimages available may be those acquired through different kinds of sensors. More precisely, this paperaddresses the problem of detecting changes between two multi-band optical images characterized bydifferent spatial and spectral resolutions. This sensor dissimilarity introduces additional issues inthe context of operational change detection. To alleviate these issues, classical change detectionmethods are applied after independent preprocessing steps (e.g., resampling) used to get the samespatial and spectral resolutions for the pair of observed images. Nevertheless, these preprocessingsteps tend to throw away relevant information. Conversely, in this paper, we propose a method thatmore effectively uses the available information by modeling the two observed images as spatial andspectral versions of two (unobserved) latent images characterized by the same high spatial and highspectral resolutions. As they cover the same scene, these latent images are expected to be globallysimilar except for possible changes in sparse spatial locations. Thus, the change detection task isenvisioned through a robust multi-band image fusion method which enforces the differences betweenthe estimated latent images to be spatially sparse. This robust fusion problem is formulated as aninverse problem which is iteratively solved using an efficient block-coordinate descent algorithm.The proposed method is applied to real panchormatic/multispectral and hyperspectral images with
Part of this work has been submitted to the IEEE Int. Conf. Acoust., Speech and Signal Process. (ICASSP), 2017 [1].Part of this work has been supported by Coordenação de Aperfeiçoamento de Ensino Superior (CAPES), Brazil, and EUFP7 through the ERANETMED JC-WATER Program, MapInvPlnt Project ANR-15-NMED-0002-02.V. Ferraris, N. Dobigeon and M. Chabert are with University of Toulouse, IRIT/INP-ENSEEIHT, France (email:{vinicius.ferraris, nicolas.dobigeon, marie.chabert}@enseeiht.fr).Q. Wei is with Department of Engineering, University of Cambridge, CB2 1PZ, Cambridge, UK (email:[email protected]).
September 21, 2016 DRAFT a r X i v : . [ c s . C V ] S e p simulated realistic changes. A comparison with state-of-the-art change detection methods evidencesthe accuracy of the proposed strategy. Index Terms
Change detection, image fusion, different resolutions, hyperspectral imagery, multispectral im-agery.
I. I
NTRODUCTION
Remote sensing is a reliable technique for Earth surface monitoring and observation [2], [3]. Oneof the most important applications using remotely sensed data is the so-called change detection (CD)problem. CD has many definitions and it is generally considered as the ability of analyzing two ormore multi-date (i.e., acquired at different time instants) and possibly multi-source (i.e., acquired bydifferent sensors) images of the same scene to detect areas where potential changes have occurred[4], [5]. Because of the increasing number of satellites and of new policies for data distribution, moremulti-temporal data becomes available. While it increases the amount of information on the presentscene, it highlights some additional issues when designing operational change detection techniques.Each remotely sensed observation image is intimately connected to the acquisition modality pro-viding a particular excerpt of the observed scene according to the sensor specifications. For instance,optical images are generally well suited to map horizontal structures, e.g., land-cover type at largescales [6]. More particularly, remote sensing images acquired by multi-band optical sensors can beclassified according to their spectral and spatial resolutions. The spectral resolution is related tothe capability in sensing the electromagnetic spectrum. This term can also refer to the number ofspectral bands [3], [7], which generally leads to a commonly adopted classification of these images: panchromatic (PAN) images, characterized by a low spectral resolution, multispectral (MS) and hyperspectral (HS) images which sense part of the spectrum with higher precision. Alternatively,multi-band optical images can be classified with respect to (w.r.t.) their spatial resolution [3], [6].The concept of spatial resolution should be understand as the capability of representing the smallestobject that can be resolved up to a specific pixel size. Images having small resolution size and finerdetails are generally identified as of (high resolution (HR) in contrast to low resolution (LR) imageswhere only coarse features are observable. Because of the physical limitations of optical passivesensors, multi-band optical images suffer from a trade-off between spectral and spatial resolution [2],[8]. To ensure that any sensor has sufficient amount of energy to guarantee a proper acquisition (interms of, e.g., signal-to-noise ratio), one of the resolutions must be decreased allowing the other tobe increased. For this reason, PAN images are generally characterized by higher spatial resolutionand lower spectral resolution than MS or a HS images.
September 21, 2016 DRAFT
Optical images have been the most studied remote sensing modality for CD since the widelyadmitted additive Gaussian modeling of optical optical sensor noises allows CD techniques to beimplemented through a simple operation of image differencing [4], [5]. Originally designed for single-band images, CD differencing methods have been adapted to handle multi-band images by consideringspectral change vectors [9], [10] and transform analysis [11], [12]. The possibility of detecting changesby exploiting both spatial and spectral information is one of the greatest advantages of these multi-band images. Nevertheless, images of same modality are not always available. In some specificscenarios, for instance consecutive to natural disasters, the availability of data imposes observationimages acquired through different kind of sensors. Such disadvantageous emergency situations yetrequire fast, flexible and accurate methods able to handle also the incompatibilities introduced bythe each sensor modality [13]–[16]. Most of the CD classical methods do not support differences inresolutions. Generally, each observed image is independently preprocessed in order to get the sameresolution and then classical CD techniques are applied. However, independent resampling operationsdo not take into account the pair of observed images and even throw away important information.Recently, a general CD framework has been proposed in [17] to deal with multi-band images withdifferent spatial and spectral resolutions based on a -step procedure (fusion, prediction, detection).Instead of independently preprocessing each observed image, this approach consists in recoveringa latent (i.e., unobserved) HR-HS image containing changed and unchanged regions by fusing bothobserved images. Then, it predicts pseudo-observed images by artificially degrading the estimated HR-HS latent image using the same forward models underlying the actually observed images. As the pairsof predicted and observed observations have the same spatial and spectral resolutions, any classicalmulti-band CD method can be finally applied to build a change map. Albeit significantly improvingdetection performance when compared to crude methods relying on independent preprocessing, the -step sequential formulation appears to be non-optimal for the following twofold reasons: i) anyinaccuracies in the fusion step are propagated throughout the subsequent degradation and detectionsteps, ii) relevant information regarding the change may be lost during the prediction steps, since itconsists in spatially or spectrally degrading the latent images to estimate the pseudo-observed images.Thus, significant improvements in terms of change detection performance may be expected providedone is able to overcome both limitations.In this paper, capitalizing on the general framework developed in [17], we show that the CD taskcan be formulated as a particular instance of the multi-band image fusion problem. However, contraryto the -step procedure in [17], the proposed approach jointly estimates a couple of distinct HR-HSlatent images corresponding to the two acquisition times as well as the change image. Since the twoHR-HS latent images are supposed to represent the same scene, they are expected to share a highlevel of similarity or, equivalently, to differ only in a few spatial locations. Thus, akin to numerous September 21, 2016 DRAFT robust factorizing models such as robust principal component analysis [18] and robust nonnegativematrix factorization [19], the two observed images are jointly approximated by a standard lineardecomposition model complemented with an HR-HS outlier term corresponding to the change image.This so-called CD-driven robust fusion of multi-band images is formulated as an inverse problemwhere, in particular, the outlier term is characterized by a spatial sparsity-inducing regularization. Theresulting objective function is solved through the use of a block coordinate descent (BCD) algorithm,which iteratively optimizes w.r.t. one latent image and the change image. Remarkably, optimizingw.r.t. the latent image boils down to a classical multi-band image fusion step and can be efficientlyconducted following the algorithmic solutions proposed in [20]. The CD map can be finally generatedfrom the recovered HR-HS change image.The paper is organized as follows. Section II formulates the change detection problem for multi-band optical image. Section III presents the solution for the formulated problem based on robustfusion. The simulation strategy as well as the results and considerations are present in Section IV.II. F
ROM CHANGE DETECTION TO ROBUST FUSION
A. Generic forward model
Let us consider the image formation process as a sequence of transformations, denoted T [ · ] , ofthe original scene into an output image. The output image of a particular sensor is referred to as theobserved image and denoted Y ∈ R n λ × m where m and n λ are the numbers of pixels and spectralbands in the observed image, respectively. It provides a limited version of the original scene withcharacteristics imposed by the image signal processor (ISP) characterizing the sensor. The originalscene can be conveniently represented by an (unknown) latent image of higher spatial and spectralresolutions, X ∈ R m λ × n , where n ≥ m and m λ ≥ n λ are the numbers of pixels and spectral bands,respectively, related to the observed image following Y = T [ X ] . (1)The intrinsic sequence of transformations of the sensor over the latent image X can be typicallyclassified as spectral or spatial degradations. On one hand, spatial degradations are related to thespatial characteristics of the sensor such as sampling scheme and optical transfer function. On theother hand, spectral degradations refer to the wavelength sensitivity and the spectral sampling. Thereare many ways to represent the degradation process. In this paper, is is considered as a sequence oflinear operations leading to the following generic forward model [21]–[23] Y = LXR + N (2)where September 21, 2016 DRAFT • L ∈ R n λ × m λ is the spectral degradation matrix, • R ∈ R n × m is the spatial degradation matrix, • N is the additive term comprising sensor noise and modeling errors.In (2), the left-multiplying matrix L ∈ R n λ × m λ degrades the latent image by combination ofsome spectral bands for each pixel while the right-multiplying matrix R ∈ R n × m degrades thelatent image by linear combination of pixels within the same spectral band. The former degradationcorresponds to a spectral resolution reduction with respect to the latent image X as in [20], [22],[23]. In practice, this degradation models an intrinsic characteristic of the sensor, namely the spectralresponse. It can be either learned by cross-calibration or known a priori [23], [24]. Conversely, thespatial degradation matrix R models the combination of different transformations which are specificof the sensor architecture taking into account external factors including wrap, blurring, translationand decimation [20], [24], [25]. In this work, since geometrical transformations such as wrap andtranslations can be corrected using image co-registration techniques in pre-processing steps, onlya spatially invariant blurring and a decimation (i.e., subsampling) will be considered. A space-invariant blur can be modeled by a symmetric convolution kernel associated with a sparse symmetricToeplitz matrix B ∈ R n × n which operates a cyclic convolution on the each individual band [26].The decimation operation, denoted by the n × m matrix S , corresponds to a uniform downsamplingoperator of factor d = d r × d c with m = n/d ones on the block diagonal and zeros elsewhere, suchthat S T S = I m [20]. To summarize, the overall spatial degradation process corresponds to the matrixcomposition R = BS ∈ R n × m .The noise corrupting multi-band optical images is generally modeled as additive and Gaussian [2],[5], [20], [27]. Thus the noise matrix N in (2) is assumed to be distributed according to the followingmatrix normal distribution N ∼ MN n λ ,m ( n λ × m , Λ , Π ) . (3)The row covariance matrix Λ carries information regarding the between-band spectral correlation.Following [20], in what follows, this covariance matrix Λ will be assumed to be diagonal, whichimplies that the noise is independent from one band to the other and characterized by a specificvariance in each band. Conversely, the column covariance matrix Π models the noise correlation The corresponding operator S T represents an upsampling transformation by zero-interpolation from m to n . The probability density function p ( X | M , Σ r , Σ r ) of a matrix normal distribution MN r,c ( M , Σ r , Σ c ) is given by [28] p ( X | M , Σ r , Σ r ) = exp ( − tr [ Σ − c ( X − M ) T Σ − r ( X − M ) ]) (2 π ) rc/ | Σ c | r/ | Σ r | c/ where M ∈ R r × c is the mean matrix, Σ r ∈ R r × r is the row covariance matrix and Σ c ∈ R c × c is the column covariancematrix. September 21, 2016 DRAFT w.r.t. to the pixel locations. Following a widely admitted hypothesis of the literature, this matrixis assumed to be identity, Π = I m , to reflect the fact the noise is spatially independent. In realapplications, both matrices Λ and Π can be estimated by previous calibrations [24]. B. Problem statement
Let us denote t j and t i the acquisition times of two co-registered multi-band optical images. Itis not assumed any specific information about time ordering, either t i < t j or t i > t j are possiblecases. Hence, without loss of generality, the HR-PAN/MS image acquired at time t i is assumed to bea low spectral resolution (i.e., PAN or MS) image of high spatial resolution denoted Y t i HR ∈ R n λ × n .The image acquired at time t j is a LR-HS image denoted Y t j LR ∈ R m λ × m . The problem addressed inthis paper consists of detecting significant changes between these two images. This is a challengingtask mainly due to the spatial and spectral resolution dissimilarity which prevents any use of simpleyet efficient differencing operation [4], [5]. To alleviate this issue, this work proposes to generalizethe CD framework introduced in [17]. More precisely, following the widely admitted forward modeldescribed in Section II-A and adopting consistent notations, the observed images Y t i HR and Y t j LR canbe related to two HR-HS latent images X t i and X t j , respectively, as follows Y t i HR = LX t i + N HR (4a) Y t j LR = X t j BS + N LR . (4b)Note that (4a) and (4b) are a specific double instance of (2). Indeed, the HR-PAN/MS (resp., LR-HS)image Y t i HR (resp., Y t j LR ) is assumed to be only a spectrally (resp., spatially) degraded version of theHR-HS latent image X t i (resp., X t i ) such that both latent images X t i ∈ R m λ × n and X t j ∈ R m λ × n share the same spectral and spatial resolutions which correspond to the highest resolutions of bothobserved images. Thereby, provided these two latent images can be efficiently inferred, any classicaldifferencing technique can be subsequently implemented on them to detect changes, notably at ahigh spatial resolution. More specifically, it would consist of evaluating an HR-HS change imagedenoted ∆ X = [∆ x , . . . , ∆ x n ] that would gather information related to any change between thetwo observed images ∆ X = X t i − X t j (5)where ∆ x p ∈ R m λ denotes the spectral change vector in the p th pixel ( p = 1 , . . . , n ). This spectralchange image can be exploited by conducting a pixel-wise change vector analysis (CVA) [29] whichexhibits the polar coordinates (i.e., magnitude and direction) of the spectral change vectors. Tospatially locate the changes, a natural approach consists of monitoring the information containedin the magnitude part of this representation [9], [10], [30], by considering the corresponding HR September 21, 2016 DRAFT spectral change energy image e = [ e , . . . , e n ] ∈ R n (6)with e p = (cid:107) ∆ x p (cid:107) , p = 1 , . . . , n. (7)When the CD problem in the p th pixel is formulated as the binary hypothesis testing H ,p : no change occurs in the p th pixel H ,p : a change occurs in the p th pixel (8)the pixel-wise statistical test can be written for a given threshold τ as e p H ,p ≷ H ,p τ. (9)The final binary HR CD map denoted d = [ d , . . . , d n ] ∈ { , } n can be derived as d p = if e p ≥ τ ( H ,p )0 otherwise ( H ,p ) . (10)When complementary information needs to be extracted from the change image ∆ X , e.g., to identifydifferent types of changes, the whole polar representation (i.e., both magnitude and direction) can befully exploited [9], [10]. As a consequence, to solve the multi-band image CD problem, the key issuelies in the joint estimation of the pair of HR-HS latent images (cid:8) X t i , X t j (cid:9) from the forward model(4) or, equivalently, the joint estimation of one of this latent image and the difference image, e.g., (cid:8) X t j , ∆ X (cid:9) . The next paragraph shows that this problem can be formulated as a particular instanceof multi-band image fusion. C. Robust multi-band image fusion
Linear forward models similar to (4) have been extensively investigated in the image processingliterature for various applications. When a unique LR-HS image Y t j LR has been observed at time t j ,recovering the HR-HS latent image X t j from the direct model (4b) can be cast as a superresolutionproblem [31], [32]. Besides, when a complementary HR-PAN/MS image Y t i HR of lower spectralresolution (i.e., PAN or MS) has been simultaneously acquired at time t i = t j under (4a), the twocorresponding latent images are expected to represent exactly the same scene, i.e., ∆ X = or,equivalently, X t i = X t j = X where the time index can be omitted. In such scenario, estimating thecommon HR-HS latent image X from the two observed images Y HR and Y LR is a multi-band imagefusion problem addressed in [20]–[23], [26], [33]–[35], also referred to as MS or HS pansharpening insome specific cases [27]. Whether the problem consists in increasing the resolution of a single imageor fusing multiple images of different spatial and spectral resolutions, the underlying objective consistsin compensating the energy trade-off of optical sensors to get highly spatially and spectrally resolved September 21, 2016 DRAFT images. Those problems are often formulated as an inverse problem, which is generally ill-posed or,at least, ill-conditioned. To overcome this issue, a classical approach consists of penalizing the datafitting terms derived from the linear models (4) and the noise statistics (3) with additional regularizingterms exploiting any prior information on the latent image. Various penalizations have been consideredin the literature, including Tikhonov regularizations expressed in the image domain [21], [36] or a ina transformed (e.g., gradient) domain [37], [38], dictionary- or patch-based regularizations [26], [31],total variation (TV) [23], [39] or regularizations based on sparse wavelet representations [40], [41].In this work, we propose to follow a similar route by addressing, in a first step, the CD problem asa linear inverse problem derived from (4). However, the CD problem addressed here differs from thecomputational imaging problems discussed above by the fact that two distinct HR-HS latent images X t i and X t j need to be inferred, which makes the inverse problem highly ill-posed. However, thisparticular applicative scenario of CD yields a natural reparametrization where relevant prior knowledgecan be conveniently exploited. More precisely, since the two HR-HS latent images are related to thesame scene observed at two time instants, they are expected to share a high level of similarity, i.e.,the change image ∆ X is expected to be spatially sparse. Thus, instead of jointly estimating the pair (cid:8) X t i , X t j (cid:9) of HR-HS latent images, we take benefit from this crucial information to rewrite the jointobservation model (4) as a function of (cid:8) X t j , ∆ X (cid:9) , i.e., Y t i HR = L (cid:0) X t j + ∆ X (cid:1) + N HR (11a) Y t j LR = X t j BS + N LR . (11b)It is worthy to note that this dual observation model parametrized by the new pair (cid:8) X t j , ∆ X (cid:9) ofimages to be inferred can be straightforwardly associated with a particular instance of the multi-bandimage fusion discussed earlier. Indeed, given the HR-HS change image ∆ X and the HR-PAN/MSimage observed at time t i , an HR-PAN/MS corrected image denoted Y t j cHR that would be acquiredby the HR-PAN/MS sensor at time t j can be defined as Y t j cHR = Y t i HR − L ∆ X . (12)In such case, the HR forward model (11a) can be easily rewritten, leading to Y t i cHR = LX t j + N HR (13a) Y t j LR = X t j BS + N LR . (13b)This observation model (13) defines a standard multi-band image fusion problem for the LR-HSobserved image Y t j LR and the corrected HR-PAN/MS image Y t j cHR . Consequently, since the changeimage ∆ X can be considered as an outlier term, akin to those encountered in several robust factorizingmodels such as robust principal component analysis (RPCA) [18] and robust nonnegative factorization September 21, 2016 DRAFT [19] which relies on a similar sparse outlier term, the joint observation model (11) naturally definesa so-called robust fusion scheme whose objective function is detailed in the next paragraph.
D. Robust fusion objective function
Because of the additive nature and the statistical properties of the noise N HR and N LR , bothobserved images Y t i HR and Y t j LR can be assumed matrix normally distributed Y t i HR | X t j , ∆ X ∼ MN n λ ,n (cid:0) L (cid:0) X t j + ∆ X (cid:1) , Λ HR , I n (cid:1) Y t j LR | X t j ∼ MN m λ ,m (cid:0) X t j BS , Λ LR , I m (cid:1) . Besides, since both observations are acquired by different modality sensors, the noise, which issensor-dependent, can be assumed statistically independent. Thus, Y t i HR | X t j , ∆ X and Y t j LR | X t j arealso statistically independent and the joint likelihood function p ( Y t i HR , Y t j LR | X t j , ∆ X ) can be writtenas a simple product of the conditional distributions p ( Y t i HR | X t j , ∆ X ) and p ( Y t j LR | X t j ) .A Bayesian formulation of the robust multi-band image fusion problem allows prior informationto be introduced to regularize the underlying estimation problem [42]. Bayesian estimators can bederived from the joint posterior distribution p ( X t j , ∆ X | Y t i HR , Y t j LR ) ∝ p ( Y t i HR , Y t j LR | X t j , ∆ X ) p ( X t i ) p (∆ X ) (14)where p ( X t i ) and p (∆ X ) correspond to the prior distributions associated with the latent and changeHR-HS images, respectively, assumed to be a priori independent. Under a maximum a posteriori(MAP) paradigm, the joint MAP estimator (cid:110) ˆ X t j MAP , ∆ ˆ X MAP (cid:111) can be derived by minimizing thenegative log-posterior, leading to the following minimization problem (cid:110) ˆ X t i MAP , ∆ ˆ X MAP (cid:111) ∈ Argmin X tj , ∆ X J (cid:0) X t j , ∆ X (cid:1) (15)with J (cid:0) X t j , ∆ X (cid:1) = 12 (cid:13)(cid:13)(cid:13) Λ − HR (cid:0) Y t i HR − L (cid:0) X t j + ∆ X (cid:1)(cid:1)(cid:13)(cid:13)(cid:13) F + 12 (cid:13)(cid:13)(cid:13) Λ − LR (cid:16) Y t j LR − X t j BS (cid:17)(cid:13)(cid:13)(cid:13) F + λφ (cid:0) X t j (cid:1) + γφ (∆ X ) . (16)The regularizing functions φ ( · ) and φ ( · ) can be related to the negative log-prior distributions ofthe HR-HS latent and change images, respectively, and the parameters λ and γ tune the amount ofcorresponding penalizations in the overall objective function J ( X t j , ∆ X ) . These functions should becarefully designed to exploit any prior knowledge regarding the parameters of interest. As discussedin Section II-C, numerous regularizations can be advocated for the HR-HS latent image X t j . In thiswork, a Tikhonov regularization proposed in [21] has been adopted φ (cid:0) X t j (cid:1) = (cid:13)(cid:13) X t j − ¯ X t j (cid:13)(cid:13) F (17) September 21, 2016 DRAFT0 where ¯ X t j refers to a crude estimate of X t j , e.g., resulting from a naive spatial interpolation of theobserved LR-HS image Y t j LR . This choice has been proven to maintain computational efficiency whileproviding accurate results [27]. Additionally, a subspace-based representation can also be adopted toenforce X t j to live in a previously identified subspace, as advocated in [43] and [23].Conversely and more critically, a specific attention should be paid to the regularizing function φ ( · ) . This function should reflect the fact that most of the pixels are expected to remain unchangedin X t i and X t j , i.e., most of the columns of the change image ∆ X are expected to be null vectors.This noticeable property can be easily translated by promoting the sparsity of the spectral changeenergy image e defined by (6). As a consequence, the regularizing function φ ( · ) is chosen as thesparsity-inducing (cid:96) -norm of the change energy image e or, equivalently, as the (cid:96) , -norm of thechange image φ (∆ X ) = (cid:107) ∆ X (cid:107) , = n (cid:88) p =1 (cid:107) ∆ x p (cid:107) . (18)This regularization is a specific instance of the non-overlapping group-lasso penalization [44] whichhas been considered in various applications to promote structured sparsity [19], [45]–[50].The next section describes an iterative algorithm which solves the minimization problem in (15).III. M INIMIZATION ALGORITHM
Computing the joint MAP estimator of the HR-HS latent image X t j at time t j and of the changeimage ∆ X can be achieved by solving the minimization problem in (15). However, no closed-formsolution can be derived for this problem. Thus this section presents a minimization algorithm whichiteratively converges to this solution. It consists in sequentially solving the problem w.r.t. to eachindividual variables X t j and ∆ X . This block coordinate descent algorithm is summarized in Algo.1 whose main steps (fusion and correction) are detailed in what follows. Algorithm 1
BCD algorithm for robust multi-band image fusion
Input: Y t j LR , Y t i HR , L , B , S , Λ HR , Λ LR . Set ∆ X . for k = 1 , . . . , K do X t j k +1 = arg min X tj J ( X t j , ∆ X k ) ∆ X k +1 = arg min ∆ X J ( X t j k +1 , ∆ X ) end forOutput: ˆ X t j MAP (cid:44) X t j K +1 and ∆ ˆ X MAP (cid:44) ∆ ˆ X K +1 September 21, 2016 DRAFT1
A. Fusion: optimization w.r.t X t j At the k th iteration of the BCD algorithm, let assume that the current value of the HR-HS changeimage is denoted ∆ X k . As suggested in Section II-C, an HR-PAN/MS corrected image Y t j cHR ,k thatwould be observed at time t j given the HR-PAN/MS image Y t i HR observed at time t i and the HR-HSchange image ∆ X k can be introduced as Y t j cHR ,k = Y t i HR − L ∆ X k . (19)Updating the current value of the HR-HS latent image consists in minimizing w.r.t. X t j the partialfunction J (cid:0) X t j (cid:1) (cid:44) J (cid:0) X t j , ∆ X k (cid:1) = (cid:13)(cid:13)(cid:13) Λ − LR (cid:16) Y t j LR − X t j BS (cid:17)(cid:13)(cid:13)(cid:13) F + (cid:13)(cid:13)(cid:13) Λ − HR (cid:16) Y t j cHR ,k − LX t j (cid:17)(cid:13)(cid:13)(cid:13) F + λφ (cid:0) X t j (cid:1) . (20)As noticed earlier, this sub-problem boils down to the multi-band image fusion which has receivedconsiderable attention in the recent image processing and remote sensing literature [20], [21], [23],[26], [27], [43]. The two difficulties arising from this formulation lies in the high dimension of theoptimization problem and in the fact that the sub-sampling operator S prevents any fast resolution inthe frequency domain by diagonalization of the spatial degradation matrix R = BS . However, withthe particular choice (17) of the regularization function φ ( · ) adopted in this paper, a closed-formsolution can still be derived and efficiently implemented. It consists in solving a matrix Sylvesterequation [20] of the form C X t j + X t j C = C (21)where the matrices C , C and C depend on the quantities involved in the problem, i.e., the virtualand observed images, the degradation operators, the noise covariance matrices and the spatiallyinterpolated image defined in (17) (see [20] for more details). Note that when a more complexregularization function φ ( · ) is considered (e.g., TV or sparse representation over a dictionary),iterative algorithmic strategies can be adopted to approximate the minimizer of J (cid:0) X t j (cid:1) . B. Correction: optimization w.r.t ∆ X Following the same strategy as in [17], let introduce the predicted
HR-PAN/MS image Y t j pHR ,k = LX t j k (22)that would be observed at time index t j by the HR-PAN/MS sensor given its spectral response L andthe current state of the HR-HS latent image X t j k at the k th iteration of the BCD algorithm. Similarlyto (5), the predicted HR-PAN/MS change image can thus be defined as ∆ Y pHR ,k = Y t i HR − Y t j pHR ,k . (23) September 21, 2016 DRAFT2
The objective function (16) w.r.t ∆ X is then rewritten by combining (22) and (23) with (16), leadingto J (∆ X ) (cid:44) J ( X t j k , ∆ X )= (cid:13)(cid:13)(cid:13) Λ − HR (∆ Y pHR ,k − L ∆ X ) (cid:13)(cid:13)(cid:13) F + γφ (∆ X ) . (24)With the specific CD-driven choice of φ ( · ) in (18), minimizing J (∆ X ) is an (cid:96) , -penalized leastsquare problem. It is characterized by the sum of a convex and differentiable data fitting term with β -Lipschitz continuous gradient ∇ f ( · ) f (∆ X ) (cid:44) (cid:13)(cid:13)(cid:13) Λ − HR (∆ Y pHR ,k − L ∆ X ) (cid:13)(cid:13)(cid:13) F (25)and a convex but non-smooth penalization g (∆ X ) (cid:44) γφ (∆ X ) = γ (cid:107) ∆ X (cid:107) , . (26)Various algorithms have been proposed to solve such convex optimization problems including forward-backward splitting [51], [52], Douglas-Rachford splitting [52], [53] and alternating direction methodof multipliers [54], [55]. Since the proximal operator related to g ( · ) can be efficiently computed(see below), in this work, we propose to resort to an iterative forward-backward algorithm which hasshown to provide the fastest yet reliable results. This algorithmic scheme is summarized in Algo. 2. Itrelies on a forward step which consists in conducting a gradient descent using the data-fitting function f ( · ) in (25), and a backward step relying on the proximal mapping associated with the penalizingfunction g ( · ) in (26). Algorithm 2
Correction step: forward-backward algorithm
Input: ∆ X k , ∆ Y pHR ,k , Λ HR , L , { η j } Jj =1 Set V (cid:44) ∆ X k for j = 1 , . . . , J do % forward step U j +1 = V j − η j ∇ f ( V j ) % backward step V j +1 = prox η j g ( U j +1 ) end forOutput: ∆ X k +1 (cid:44) V J +1 Since the HR-PAN/MS observed image has only a few spectral bands (e.g., n λ ∼ ), the spectraldegradation matrix L ∈ R n λ × m λ is a fat (and generally full-row rank) matrix. Thus, the correspondinggradient operator ∇ f ( · ) defining the forward step (see line 4 of Algo. 2) can be easily and efficiently September 21, 2016 DRAFT3 computed. Conversely, the proximal operator associated with g ( · ) in (26) and required during thebackward step (see line 6 of Algo. 2) is defined asprox ηg ( U ) = arg min Z (cid:18) γ (cid:107) Z (cid:107) , + 12 η (cid:107) Z − U (cid:107) F (cid:19) (27)for some η > . The function g ( U ) in (26) can be split as (cid:80) np =1 g p ( u p ) with, for each column, g p ( · ) = γ (cid:107)·(cid:107) . Based on the separability property of proximal operators [55], the operator (27) canbe decomposed and computed for each pixel location p ( p = 1 , . . . , n ) as (cid:2) prox ηg ( U ) (cid:3) p = prox ηg p ( u p ) (28)where the notations [ · ] p stands for the p th column. Thus, only the proximal operator associated withthe Euclidean distance induced by the (cid:96) -norm is necessary. The Moreau decomposition [55] u p = prox ηg ( u p ) + η prox η − g ∗ p (cid:0) η − u p (cid:1) (29)establishes a relationship between the proximal operators of the function g p ( · ) and its conjugate g ∗ p ( · ) .When the function g ( · ) is a general norm, its conjugate corresponds to the indicator function into theball B defined by its dual norm [48], [55], leading toprox ηg ( u p ) = u p − η P B (cid:18) u p η (cid:19) (30)where P B ( · ) denotes the projection. When g ( · ) is defined by (26), since the (cid:96) -norm is self-dual, thisprojection is P B ( u p ) = γ u p (cid:107) u p (cid:107) if (cid:107) u p (cid:107) > γ u p otherwise. (31)Consequently, replacing (31) in (30), the proximal operator associated with the function g p ( · ) in (28)is prox ηg p ( u p ) = (cid:16) − ηγ (cid:107) u p (cid:107) (cid:17) u p if (cid:107) u p (cid:107) > ηγ otherwise. (32)To conclude, the correction procedure can be interpreted as first a gradient descent step for spectraldeblurring of the HR-HS change image from the HR-PAN/MS predicted change image (forward step),followed by a soft-thresholding of the resulting HR-HS change image to promote spatial sparsity(backward step). IV. E XPERIMENTS
A. Simulation framework
Real dataset for assessing performance of CD algorithms is rarely available. Indeed, this assessmentrequires a couple of images acquired at two different dates, geometrically and radiometrically pre-corrected, presenting changes and, for the scenario considered in this paper, coming from two
September 21, 2016 DRAFT4 (a) d HR (b) X t i (c) X t j (d) Y t i HR (e) Y t j LR Fig. 1: One particular simulation configuration: (a) the HR change mask d HR , (b)-(c) the HR-HSlatent images X t i and X t j , (d)-(e) the spectrally degraded version HR-MS observed image Y t i HR andthe spatially degraded LR-HS observed image Y t j LR . Note that, in this particular configuration, thechange-inducing function ϑ t j ( · , d HR ) is the identity operator (i.e., A t j = A ref ) since it does not applyany change into the corresponding HR-HS latent image X t j while the function ϑ t i ( · , d HR ) includesa triangular region of pixels in X t i affected by changes. Moreover, the HR observed image is herean MS image.different optical sensors. To alleviate this issue, inspired by the well-known Wald’s evaluation protocoldedicated to pansharpening algorithms [56], a framework has been proposed in [17] to assess theperformance of CD algorithms when dealing with optical images of different spatial and spectralresolutions. This framework only requires a single HR-HS reference image X ref and generates apair of latent HR-HS images X t i and X t j resulting from a unmixing-mixing process. This processallows synthetic yet realistic changes to be incorporated within one of these latent images, w.r.t.a pre-defined binary reference HR change mask d HR ∈ R n locating the pixels affected by thesechanges and further used to assess the performance of the CD algorithms. This procedure allowsvarious physically-inspired changes to be considered, e.g., by tuning the relative abundance of a eachendmember or replacing one of them my another. This protocol is briefly described below (see [17]for more details).
1) Reference image:
The HR-HS reference image X ref used in the experiments reported in thispaper is a × × HS image of the Pavia University, Italy, acquired by the reflective opticssystem imaging spectrometer (ROSIS) sensor. This image has undergone a pre-precessing to smooththe atmospheric effects of vapor water absorption by removing some bands. Thus the final HR-HSreference image is of size × × . September 21, 2016 DRAFT5
2) Generating the changes:
Using the same procedure proposed in [17], the HR-HS referenceimage X ref ∈ R m λ × n has been linearly unmixed to define the reference matrix M ref ∈ R m λ × R of R endmember spectral signatures and the corresponding reference abundance matrix A ref ∈ R R × n such that X ref ≈ M ref A ref . The two latent HR-HS images X t i and X t j are then computed as linearmixture of the endmembers in M ref with corresponding abundance matrices A t i and A t j , respectively,derived from the reference abundances A ref and the change mask d HR , i.e., X t i = M ref A t i and X t j = M ref A t j with A t i = ϑ t i (cid:16) A ref , d HR (cid:17) and A t j = ϑ t j (cid:16) A ref , d HR (cid:17) where the two change-inducing functions ϑ t · : R R × n × R n → R R × n are defined to simulate realisticchanges in some pixels of the HR-HS latent images. Three sets of predefined change maskshave been designed according to three specific change rules introduced in [17]. For each simulatedpair (cid:8) X t i , X t j (cid:9) , one of the two functions ϑ t · ( · , d HR ) is defined as a “no-change” operator, i.e., ϑ t · ( A ref , d HR ) = A ref , which leads to an overall set of simulated pairs (cid:8) X t i , X t j (cid:9) of HR-HSlatent images.
3) Generating the observed images:
The HR-PAN/MS observed image Y t i HR is obtained by spec-trally degrading the corresponding HR-HS latent image X t i . Two scenarios are considered. Scenario1 consists in averaging the first bands of the HR-HS latent image to produce an HR-PAN image.Conversely, Scenario 2 considers an HR-MS image by spectrally filtering the HR-HS latent image X t i with a -band LANDSAT-like spectral response. Moreover, to generate a spatially degraded image Y t j LR , the respective latent image X t i has been blurred by a × Gaussian kernel and subsequentlyequally down-sampled in the vertical and horizontal directions with a down-sampling ratio d = 5 . Toillustrate, Fig. 1 shows one of the simulation configurations used during the experiments. B. Compared methods
The proposed robust fusion-based CD technique has been compared to four methods able to dealwith optical images of different spatial and spectral resolutions. The first one has been proposedin [17] and also relies on a fusion-based approach. Up to the authors’ knowledge, it was the firstoperational CD technique able to operate with multi-band optical images of different spatial andspectral images. Contrary to the model (4) proposed in this paper, it consists in recovering a commonlatent image by fusing the two observed images and then predicting an HR-PAN/MS image ˆ Y F ,t i HR from the underlying forward model. An HR-PAN/MS change image ∆ Y F ,t i HR has been then computedas in (5) from the pair of HR-PAN/MS observed and predicted images (cid:110) Y t i HR , ˆ Y F ,t i HR (cid:111) . Finally, as September 21, 2016 DRAFT6 recommended in [17], a spatially-regularized CVA (sCVA) similar to the decision rule detailed inSection II-B has been conducted on ∆ Y F ,t i HR to produce an estimated HR CD mask denoted ˆ d F .The second method aims at producing an HR-PAN/MS predicted image by successive spatialsuperresolution and spectral degradation. More precisely, an HR-HS latent image is first recovered byconducting a band-wise spatial superresolution of the observed LR-HS Y t j LR following the fast methodin [32]. Then this latent image is spectrally degraded according to produce an HR-PAN/MS predictedimage ˆ Y SD ,t j HR . Similarly to the previous fusion-based method, sCVA has been finally conducted onthe pair (cid:110) Y t i HR , ˆ Y SD ,t j HR (cid:111) to produce an HR CD mask denoted ˆ d SD . The third CD method applies thesame procedure with a reverse order of spatial superresolution and spectral degradation, and producesproduces an HR change mask denoted ˆ d DS from the pair of HR-PAN/MS images (cid:110) Y t i HR , ˆ Y DS ,t j HR (cid:111) .The fourth CD method, referred to as the worst-case (WC) as in [17], build a LR change mask ˆ d WC by crudely conducting a sCVA on a spatially degraded version of the HR-PAN/MS image and aspectrally degraded version of the LR-HS image. C. Figures-of-merit
The CD performances of these four methods, as well as the performance of the proposed robustfusion-based method whose HR change mask is denoted ˆ d RF , have been visually assessed fromempirical receiver operating characteristics (ROC), representing the estimated pixel-wise probabilityof detection ( PD ) as a function of the probability of false alarm ( PFA ). Moreover, two quantitativecriteria derived from these ROC curves have been computed, namely, i ) the area under the curve(AUC), corresponding to the integral of the ROC curve and ii ) the distance between the no detectionpoint ( P F A = 1 , PD = 0) and the point at the interception of the ROC curve with the diagonal linedefined by
PFA = 1 − PD . For both metrics, greater the criterion, better the detection.TABLE I: Scenarios 1 & 2: quantitative detection performance (AUC and distance). ˆ d RF ˆ d RF ˆ d WC ˆ d DS ˆ d SD Scenario 1 AUC . . . . . Dist. . . . . . Scenario 2 AUC . . . . . Dist. . . . . . D. Results1) Scenario 1 (HR-PAN vs. LR-HS):
The ROC curves depicted in Fig. 2 with corresponding metricsin Table I (first two rows) correspond to the CD results obtained from a pair of HR-PAN and LR-HS
September 21, 2016 DRAFT7
PFA P D ˆ d RF ˆ d F ˆ d WC ˆ d DS ˆ d SD Fig. 2: Scenario 1 (HR-PAN vs. LR-HS): ROC curves.observed images. Clearly, the proposed robust fusion-based CD technique outperforms the four otherCD techniques. More importantly, it provides almost perfect detections even for very low PFA, i.e.,for very low energy changes. Note that the CD mask d WC estimated by the worst-case method isdefined at an LR.
2) Scenario 2 (HR-MS vs. LR-HS):
Applying the same procedure of Scenario 1 but now consideringan HR-MS observed image instead of the HR-PAN observed image leads to very similar overallperformance. The ROC plot is depicted in Fig. 3 with corresponding metrics in Table I (last tworows). As in Scenario 1, comparing curves in Fig. 2 shows that the proposed method offers a higherprecision even when analyzing a lower spectral resolution HR observed image.As an additional result, for Scenario 2, Fig. 4 compares the abilities of detecting changes of
September 21, 2016 DRAFT8
PFA P D ˆ d RF ˆ d F ˆ d WC ˆ d DS ˆ d SD Fig. 3: Scenario 2 (HR-MS vs. LR-HS): ROC curves.decreasing size of the proposed method against the fusion-based CD method [17] and the worst-caseCD method. Figure 4(a) and 4(b) shows the observed image pair Y t i HR and Y t j LR containing multiplechanges with size varying from × - pixel to × -pixels, with the corresponding change mask d HR presented in Fig. 4(c). Figures 4(d) and 4(e) present the change masks ˆ d F and ˆ d WC recoveredby the two competing methods, respectively, while the CD mask ˆ d RF recovered by the proposedrobust fusion-based method is reported in Fig. 4(f) shows the proposed CD. For each technique,the decision threshold τ required in the CVA in (10) has been tuned to reach the higher distancevalue in the corresponding ROC curves. The first advantage of the proposed method is a significantdecrease of the number of false alarm which are due to propagated errors when implementing the twoother methods. Moreover, these results proves once again that the proposed method achieves a better September 21, 2016 DRAFT9 detection rate with a higher resolution, even when considering extremely localized change regions.Remaining false alarms only occur near edges between change and no-change regions of small sizedue to the difference of spatial resolutions and the width of the blur kernel. Note also that the CDmask estimated by the worst-case method is of coarse scale since based on the comparison of twoLR-MS images. (a) Y t i HR (b) Y t j LR (c) d HR (d) ˆ d F (e) ˆ d WC (f) ˆ d RF Fig. 4: CD precision for Scenario 2 (HR-MS vs. LR-HS): (a) HR-MS observed image Y t i HR , (b)LR-HS observed image Y t j LR , (c) actual change mask d HR , (d) change mask ˆ d F estimated by thefusion-based approach [17], (e) change mask ˆ d WC estimated by the worst-case approach, (f) changemask ˆ d RF estimated by the proposed robust fusion-based approach.V. C ONCLUSION
This paper proposed a robust fusion-based change detection technique to handle two multi-bandoptical observed images of different spatial and spectral resolutions. The technique was based onthe definition of two high resolution hyperspectral latent images related to the observed images viaa double physically-inspired forward model. The difference between these two latent images was
September 21, 2016 DRAFT0 assumed to be spatially sparse, implicitly locating the changes at a high resolution scale. Inferringthese two latent images was formulated as an inverse problem which was solved within a -stepiterative scheme. This algorithmic strategy amounted to solve a standard fusion problem and an (cid:96) , -penalized spectral deblurring step. Contrary to the methods already proposed in the literature, modelingerrors were not anymore propagate in-between steps. A simulation protocol allowed the performanceof the proposed technique in terms of detection and precision to be assessed and compared with theperformance of various algorithms. The detection rate as well as the accuracy of this method wasclearly improved by this robust-fusion based algorithm. Future works include the detection of changebetween optical and non-optical data. R EFERENCES [1] V. Ferraris, N. Dobigeon, Q. Wei, and M. Chabert, “Change detection between multi-band images using a robustfusion-based approach,” in
Proc. IEEE Int. Conf. Acoust., Speech and Signal Process. (ICASSP) , 2016, submitted.[2] C. Elachi and J. Van Zyl,
Introduction to the physics and techniques of remote sensing , 2nd ed., ser. Wiley series inremote sensing. Hoboken, N.J.: Wiley-Interscience, 2006.[3] J. B. Campbell and R. H. Wynne,
Introduction to remote sensing , 5th ed. New York: Guilford Press, 2011.[4] A. Singh, “Review Article Digital change detection techniques using remotely-sensed data,”
Int. J. Remote Sens. ,vol. 10, no. 6, pp. 989–1003, June 1989.[5] F. Bovolo and L. Bruzzone, “The time variable in data fusion: A change detection perspective,”
IEEE Geosci. RemoteSens. Mag. , vol. 3, no. 3, pp. 8–26, Sept. 2015.[6] M. Dalla Mura, S. Prasad, F. Pacifici, P. Gamba, J. Chanussot, and J. A. Benediktsson, “Challenges and Opportunitiesof Multimodality and Data Fusion in Remote Sensing,”
Proc. IEEE , vol. 103, no. 9, pp. 1585–1601, Sept. 2015.[7] D. Landgrebe, “Hyperspectral image data analysis,”
IEEE Signal Process. Mag. , vol. 19, no. 1, pp. 17–28, 2002.[8] J. C. Price, “Spectral band selection for visible-near infrared remote sensing: spectral-spatial resolution tradeoffs,”
IEEE Trans. Geosci. Remote Sens. , vol. 35, no. 5, pp. 1277–1285, 1997.[9] F. Bovolo and L. Bruzzone, “A theoretical framework for unsupervised change detection based on change vectoranalysis in the polar domain,”
IEEE Trans. Geosci. Remote Sens. , vol. 45, no. 1, pp. 218–236, Jan. 2007.[10] F. Bovolo, S. Marchesi, and L. Bruzzone, “A framework for automatic and unsupervised detection of multiple changesin multitemporal images,”
IEEE Trans. Geosci. Remote Sens. , vol. 50, no. 6, pp. 2196–2212, June 2012.[11] A. A. Nielsen, K. Conradsen, and J. J. Simpson, “Multivariate alteration detection (MAD) and MAF postprocessingin multispectral, bitemporal image data: New approaches to change detection studies,”
Remote Sens. Environment ,vol. 64, no. 1, pp. 1–19, 1998.[12] A. A. Nielsen, “The Regularized Iteratively Reweighted MAD Method for Change Detection in Multi- andHyperspectral Data,”
IEEE Trans. Image Process. , vol. 16, no. 2, pp. 463–478, Feb. 2007.[13] J. Inglada, “Similarity measures for multisensor remote sensing images,” in
Proc. IEEE Int. Conf. Geosci. RemoteSens. (IGARSS) , vol. 1. IEEE, 2002, pp. 104–106.[14] V. Alberga, M. Idrissa, V. Lacroix, and J. Inglada, “Performance estimation of similarity measures of multi-sensorimages for change detection applications,” in
Proc. IEEE Int. Workshop Analysis Multitemporal Remote Sensing Images(MultiTemp) . Leuven: IEEE, 2007, pp. 1 – 5.[15] G. Mercier, G. Moser, and S. Serpico, “Conditional copula for change detection on heterogeneous sar data,” in
Proc.IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) . IEEE, 2007, pp. 2394–2397.
September 21, 2016 DRAFT1 [16] J. Prendes, M. Chabert, F. Pascal, A. Giros, and J.-Y. Tourneret, “A new multivariate statistical model for changedetection in images acquired by homogeneous and heterogeneous sensors,”
IEEE Trans. Image Process. , vol. 24,no. 3, pp. 799–812, 2015.[17] V. Ferraris, N. Dobigeon, Q. Wei, and M. Chabert, “Detecting changes between optical images of different spatial andspectral resolutions: a fusion-based approach,” 2016, submitted.[18] E. J. Candés, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis?”
Journal of the ACM (JACM) ,vol. 58, no. 3, p. 11, 2011.[19] C. Févotte and N. Dobigeon, “Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization,”
IEEETrans. Image Process. , vol. 24, no. 12, pp. 4810–4819, 2015.[20] Q. Wei, N. Dobigeon, and J.-Y. Tourneret, “Fast Fusion of Multi-Band Images Based on Solving a Sylvester Equation,”
IEEE Trans. Image Process. , vol. 24, no. 11, pp. 4109–4121, Nov. 2015.[21] ——, “Bayesian Fusion of Multi-Band Images,”
IEEE J. Sel. Topics Signal Process. , vol. 9, no. 6, pp. 1117–1127,Sept. 2015.[22] N. Yokoya, T. Yairi, and A. Iwasaki, “Coupled nonnegative matrix factorization unmixing for hyperspectral andmultispectral data fusion,”
IEEE Trans. Geosci. Remote Sens. , vol. 50, no. 2, pp. 528–537, Feb. 2012.[23] M. Simões, J. Bioucas Dias, L. Almeida, and J. Chanussot, “A convex formulation for hyperspectral imagesuperresolution via subspace-based regularization,”
IEEE Trans. Geosci. Remote Sens. , vol. 6, no. 53, pp. 3373–3388,June 2015.[24] N. Yokoya, N. Mayumi, and A. Iwasaki, “Cross-Calibration for Data Fusion of EO-1/Hyperion and Terra/ASTER,”
IEEE J. Sel. Topics Appl. Earth Observations Remote Sens. , vol. 6, no. 2, pp. 419–426, April 2013.[25] F. Heide, O. Gallo, M. Steinberger, J. Liu, Y.-T. Tsai, W. Heidrich, M. Rouf, K. Egiazarian, D. Pajak, J. Kautz,D. Reddy, and K. Pulli, “FlexISP: A Flexible Camera Image Processing Framework,”
ACM Transactions on Graphics(TOG) - Proceedings of ACM SIGGRAPH Asia 2014 , vol. 33, no. 6, 2014.[26] Q. Wei, J. Bioucas-Dias, N. Dobigeon, and J.-Y. Tourneret, “Hyperspectral and multispectral image fusion based ona sparse representation,”
IEEE Trans. Geosci. Remote Sens. , vol. 53, no. 7, pp. 3658–3668, 2015.[27] L. Loncan, L. B. de Almeida, J. M. Bioucas-Dias, X. Briottet, J. Chanussot, N. Dobigeon, S. Fabre, W. Liao,G. A. Licciardi, M. Simoes, J.-Y. Tourneret, M. A. Veganzones, G. Vivone, Q. Wei, and N. Yokoya, “Hyperspectralpansharpening: A review,”
IEEE Geosci. Remote Sens. Mag. , vol. 3, no. 3, pp. 27–46, Sept. 2015.[28] A. K. Gupta and D. K. Nagar,
Matrix Variate Distribution , ser. Monographs and Surveys in Pure and AppliedMathematics. Chapman and Hall, 1999, no. 104.[29] R. D. Johnson and E. S. Kasischke, “Change vector analysis: A technique for the multispectral monitoring of landcover and condition,”
Int. J. Remote Sens. , vol. 19, no. 3, pp. 411–426, Jan. 1998.[30] A. Singh, “Digital change detection techniques using remotely-sensed data,”
Int. J. Remote Sens. , vol. 10, no. 6, pp.989–1003, 1989.[31] Jianchao Yang, J. Wright, T. S. Huang, and Yi Ma, “Image super-resolution via sparse representation,”
IEEE Trans.Image Process. , vol. 19, no. 11, pp. 2861–2873, Nov. 2010.[32] N. Zhao, Q. Wei, A. Basarab, N. Dobigeon, D. Kouame, and J.-Y. Tourneret, “Fast Single Image Super-ResolutionUsing a New Analytical Solution for (cid:96) – (cid:96) Problems,”
IEEE Trans. Image Process. , vol. 25, no. 8, pp. 3683–3697,Aug. 2016.[33] R. C. Hardie, M. T. Eismann, and G. L. Wilson, “MAP estimation for hyperspectral image resolution enhancementusing an auxiliary sensor,”
IEEE Trans. Image Process. , vol. 13, no. 9, pp. 1174–1184, Sept. 2004.[34] M. T. Eismann and R. C. Hardie, “Hyperspectral resolution enhancement using high-resolution multispectral imagerywith arbitrary response functions,”
IEEE Trans. Image Process. , vol. 43, no. 3, pp. 455–465, March 2005.
September 21, 2016 DRAFT2 [35] Y. Zhang, S. De Backer, and P. Scheunders, “Noise-resistant wavelet-based Bayesian fusion of multispectral andhyperspectral images,”
IEEE Trans. Geosci. Remote Sens. , vol. 47, no. 11, pp. 3834–3843, Nov. 2009.[36] M. Ebrahimi and E. R. Vrscay, “Regularization schemes involving self-similarity in imaging inverse problems,” in
J.Physics: Conf. Series , 2008.[37] Y.-W. Tai, S. Liu, M. S. Brown, and S. Lin, “Super resolution using edge prior and single image detail synthesis,” in
Proc. IEEE Conference on Computer Vision and Pattern Recognition (CVPR) , 2010, pp. 2400 – 2407.[38] J. Sun, J. Sun, Z. Xu, and H.-Y. Shum, “Gradient profile prior and its applications in image super-resolution andenhancement,”
IEEE Trans. Image Process. , vol. 20, no. 6, pp. 1529 – 1542, 2011.[39] H. A. Aly and E. Dubois, “Image up-sampling using total-variation regularization with a new observation model,”
IEEE Trans. Image Process. , vol. 14, no. 10, pp. 1647–1659, 2005.[40] C. V. Jiji, M. V. Joshi, and S. Chaudhuri, “Single-frame image super-resolution using learned wavelet coefficients,”
Int. J. Imaging Syst. Technol. , vol. 14, pp. 105–112, 2004.[41] J. M. Bioucas-Dias, “Bayesian wavelet-based image deconvolution: A GEM algorithm exploiting a class of heavy-tailedpriors,”
IEEE Trans. Image Process. , vol. 15, no. 4, pp. 937–951, 2006.[42] J. Idier,
Bayesian approach to inverse problems , J. Idier, Ed. ISTE ; Wiley, 2008.[43] Q. Wei, N. Dobigeon, and J.-Y. Tourneret, “Bayesian fusion of multispectral and hyperspectral images using a blockcoordinate descent method,” in
Proc. IEEE GRSS Workshop Hyperspectral Image SIgnal Process.: Evolution in RemoteSens. (WHISPERS) , 2015.[44] F. Bach, R. Jenatton, J. Mairal, and G. Obozinski,
Convex Optimization with Sparsity-Inducing Norms , ser. Optimizationfor Machine Learning. MIT Press, 2011.[45] S. Cotter, B. Rao, Kjersti Engan, and K. Kreutz-Delgado, “Sparse solutions to linear inverse problems with multiplemeasurement vectors,”
IEEE Trans. Signal Process. , vol. 53, no. 7, pp. 2477–2488, July 2005.[46] C. Ding, D. Zhou, X. He, and H. Zha, “R1-PCA: rotational invariant L1-norm principal component analysis for robustsubspace factorization,” in
Proc. Int. Conf. Machine Learning (ICML) , 2006, pp. 281–288.[47] J. Liu, S. Ji, and J. Ye, “Multi-task feature learning via efficient L2,1-norm minimization,”
Uncertainty in ArtificialIntelligence , 2009.[48] S. Wright, R. Nowak, and M. Figueiredo, “Sparse Reconstruction by Separable Approximation,”
IEEE Trans. SignalProcess. , vol. 57, no. 7, pp. 2479–2493, July 2009.[49] F. Nie, H. Huang, X. Cai, and C. H. Ding, “Efficient and robust feature selection via joint l2, 1-norms minimization,”in
Advances in neural information processing systems , 2010, pp. 1813–1821.[50] H. Lu, X. Long, and J. Lv, “A fast algorithm for recovery of jointly sparse vectors based on the alternating directionmethods,” in
International Conference on Artificial Intelligence and Statistics , 2011, pp. 461–469.[51] P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,”
Multiscale Modeling &Simulation , vol. 4, no. 4, pp. 1168–1200, 2005.[52] P. L. Combettes and J.-C. Pesquet, “Proximal Splitting Methods in Signal Processing,”
Fixed-Point Algorithm forInverse Problems in Science and Engineering , no. Springer Optimization and Its Applications, 2011.[53] ——, “A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery,”
IEEE J. Sel. TopicsSignal Process. , vol. 1, no. 4, pp. 564–574, 2007.[54] S. Boyd, “Distributed optimization and statistical learning via the alternating direction method of multipliers,”
Foundations and Trends in Machine Learning , vol. 3, no. 1, pp. 1–122, 2010.[55] N. Parikh and S. Boyd, “Proximal Algorithms,”
Foundations and Trends in Optimization , vol. 1, no. 3, pp. 123–231,2013.[56] L. Wald, T. Ranchin, and M. Mangolini, “Fusion of satellite images of different spatial resolutions: assessing thequality of resulting images,”
Photogrammetric engineering and remote sensing , vol. 63, no. 6, pp. 691–699, 1997., vol. 63, no. 6, pp. 691–699, 1997.