UNMIXING MULTITEMPORAL HYPERSPECTRAL IMAGES ACCOUNTING FOR ENDMEMBER VARIABILITY

.


INTRODUCTION
Spectral unmixing (SU) consists of identifying the macroscopic materials present in an hyperspectral image (HI) (called endmembers) and their proportions (called abundances).The spectrum of each material might vary from one pixel to another resulting in the so-called endmember variability (EV) [1,2].This variability appears spatially and thus will be denoted as spatial endmember variability (SEV).In the literature, some statistical methods which address SEV consider the endmembers as random variables.These models include the beta compositional model [3] and the normal compositional model (NCM) [4][5][6][7].In this paper we are interested in applications where successive HIs are acquired for the same scene at different time instants to study the temporal evolution of the physical elements.In this case, a temporal endmember variability (TEV) will also appear depending on the observation season or other weather factors.
The main contribution of this paper is the development of a hierarchical Bayesian model that considers both spatial and temporal EVs for unsupervised hyperspectral unmixing.SEV is introduced by considering the NCM model for each image at successive time instants.To introduce temporal correlation, the endmember means and abundances are assigned a prior enforcing smooth evolution between consecutive images.This temporal prior is defined from the discrete Laplacian of the different parameters and has shown increasing interest for many problems such as image deconvolution [8,9], hyperspectral unmixing [10], medical imaging [11] and altimetry [12].Moreover, the proposed Bayesian model assumes that the endmember variances are different in the different spectral bands, which has shown interesting properties [7,13].
An algorithm is then proposed to estimate the unknown model parameters.However, the minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators cannot be easily computed from the obtained joint posterior.The proposed algorithm alleviates this problem by generating samples distributed according to this posterior using Markov chain Monte Carlo (MCMC) methods.More precisely, we use a Gibbs sampler coupled with a constrained Hamiltonian Monte Carlo (CHMC) algorithm since it has been shown to have good sampling properties for high-dimensional vectors [7,13,14].
The paper is structured as follows.Section 2 introduces our notations and the proposed mixing model.The hierarchical Bayesian model which accounts for the spatial and temporal EVs in HIs is introduced in Section 3. Section 4 analyzes the performance of the proposed algorithm when applied to synthetic images.Conclusions and future works are reported in Section 5.

Notations
Matrix and vectors are denoted with bold upper and lower case letters.A vector is by convention a column vector.When there are many indices, the time index is indicated after a semicolon.For example, m r;t ∈ R L×1 denotes the spectrum of the rth endmember mean for the tth image.The notation 1 : T means the elements from 1 to T .For example, m r,1:T denotes the set of vectors (m r;1 , • • • , m r;T ).

Mixing model and endmember variability
In this paper, we consider T successive HIs observing the same scene at different time instants.For the tth image, and in order to consider SEV, the observation model is based on the NCM defined as [4][5][6][7] a rn;t s rn;t = S n;t a n;t (1) with s rn;t ∼ N (m r;t , Σ t ) which introduces spatial EV since the endmembers vary from one pixel to another, T is the endmember variance vector of the tth image, T is the (R × 1) abundance vector of the nth pixel from the tth image.The abundance vector a n;t contains proportions and thus should satisfy the physical positivity and sum-to-one (PSTO) constraints a rn;t ≥ 0, ∀r ∈ {1, . . ., R} and R r=1 a rn;t = 1 for each image #t.Note that (1) introduces SEV by considering different endmembers for each image pixel.However, it does not consider any correlation between successive temporal images.The next section introduces a hierarchical Bayesian model that takes into account both temporal and spectral EVs.

HIERARCHICAL BAYESIAN MODEL
The unknown parameters of the proposed model include the (L × R × T ) endmember mean matrices M 1:T for all time instants, the (T ×L) matrix containing the endmember variances denoted by σ 2  1:T , the (R × N × T ) abundance matrix A 1:T (whose nth column is A :n;t = a n;t ) for all time instants.

Likelihood
The observation model defined in (1) and the Gaussian properties of the endmembers s rn;t yield the following likelihood Assuming independence between the observed pixels yields

Parameter priors
This section introduces the prior distributions that we have chosen for the parameters of interest A 1:T , M 1:T and σ 2 1:T .

Abundance matrix A
In order to satisfy the PSTO constraints, the abundances belong to the simplex S given by S = a n;t , a rn;t ≥ 0, ∀r, ∀t and R r=1 In [15,16], a uniform distribution on S has been considered for the abundance vector a n;t of a given image.In this paper, we consider a smooth variation of the abundances from one temporal image to another.This correlation is introduced by considering the following truncated Gaussian prior where I A (.) is the indicator function over the set A, D denotes the discrete Laplacian operator and 2 n is a hyperparameter that controls the degree of smoothness for the abundances (it depends on the pixel index because the abundances vary differently from one pixel to another).This prior can also be written as and N (x|μ, Σ) denotes the normal distribution of the variable x with mean μ and covariance matrix Σ.This Gaussian prior distribution constrains a smooth evolution for the abundances.It has been used in different contexts such as image deconvolution [8,9], spectral unmixing of HIs [10] or for medical imaging applications [11].Note that the prior is truncated on the simplex to satisfy the physical PSTO constrains.Note also that we have considered the abundance reparametrization procedure introduced in [13,17] to simplify the sampling procedure.Indeed, this reparametrization expresses the PSTO constraints by only using nonnegativity constraints, easily handled by the sampling procedure as already shown in [7,13]

Prior for the endmember means
To introduce temporal correlation between the endmember means, we assign the following prior for the rth endmember where m r ;t is a fixed spectra (estimated from the data using an endmember extraction algorithm such as VCA [18]) and ψ 2 controls the smoothness of the temporal evolution of m r;1:T for each spectral band.Indeed, the endmember values vary differently from one spectral band to another.Since M contains reflectances, it should satisfy the following constraints 0 < m rl;t < 1, ∀r, ∀l, ∀t, [7,13].Therefore, the Gaussian prior ( 6) has been truncated on the set [0, 1].Assuming prior independence between the endmember means yields f

Prior for the endmember variances
As in [5,7], a non informative Jeffreys prior is chosen for the endmember variances as follows where we have assumed prior independence between the endmember variances.

Hyperparameter priors
As in [19], the hyperparameters L have been fixed empirically by considering the dynamic range of each parameter.

Posterior distribution
The proposed Bayesian model depends on the parameters θ p = A 1:T , M 1:T , σ 2 1:T and the fixed hyperparameters θ h = 2 , ψ 2 .The joint posterior distribution of the unknown parameters can be computed from the following hierarchical structure with T , where we have assumed prior independence between the parameters.The MMSE and MAP estimators associated with the posterior (8) are not easy to determine.These estimators are therefore approximated using samples generated according to (8) by considering an MCMC approach.This can be achieved by using a Gibbs sampler that generates samples according to the conditional distributions of (8) [20].Indeed, the Gibbs algorithm samples sequentially the parameters A, M and σ 2 .The conditional distribution associated with each of the variable σ 2 is an inverse gamma distribution that is easy to sample.However, the conditional distributions of both A and M1 are more complex and require the use of an accept-reject procedure to sample from them.In this paper, we consider the CHMC algorithm that has shown interesting mixing properties for high-dimensional problems [21,Chap. 5] [7,13].Note finally that the sampling algorithm is not described here for brevity (the reader is invited to consult [7,13] for more details about the sampling algorithm).

SIMULATION RESULTS
In this section the performance of the proposed algorithm on synthetic data is assessed.First, the criteria used for the evaluation of the unmixing quality is introduced.Second, the proposed algorithm is compared with state-of-the-art methods considering two sequences of synthetic images.The third part considers a synthetic sequence whose parameters were extracted from an actual image to approximate a real image scenario.

Evaluation criteria
In order to evaluate the quality of the unmixing strategy, we have considered synthetic images with known abundances and endmembers.The unmixing performance can then be measured by using the average root mean square error (aRMSE) and the average spectral angle mapper (aSAM) of the estimates where || • || denotes the standard l 2 norm such that and arccos(•) is the inverse cosine operator.The Earth movers distance (EMD) criterion (based on the Euclidean distance and described in [22]) has also been considered to simultaneously evaluate the estimated endmembers and abundances.Note that the reconstruction error (RE) criterion can also be evaluated for the nth measured and estimated pixel spectra y n;t , ŷn;t as follows

Synthetic images
This section considers two sequences of T = 20 synthetic images.Each image contains 30 × 30 pixels and is generated according to (1) with R = 3 physical elements, (construction concrete, green grass and micaceous loam), corresponding to spectral signatures available in the ENVI software library [23].At each time instant #t, the endmember means are generated by introducing variability on the ENVI-like spectral signatures.These endmember means vary smoothly from one image to the next.Using the spectra obtained, the tth image is generated by considering the NCM model with endmember variances that increase linearly with respect to the spectral bands such that σ2 1 = 3 × 10 (−4) and σ 2 L = 25 × 10 (−4) .A smooth temporal evolution is also considered when generating the abundances under the PSTO constraints.For the first sequence of images (denoted by I 1 ), the abundances are uniformly distributed in the truncated simplex S with (a r < 0.9, ∀r).The constraint a r < 0.9 implies that there is no pure pixel in the image, which makes the problem more challenging.However, in presence of a highly mixed scenario, the abundances can be concentrated in some regions of the simplex [7,24].Thus, in the second sequence (denoted by I 2 ), we have considered a non uniform abundance repartition in the truncated simplex S with (a r < 0.9, ∀r).The abundances have been generated by considering a Dirichlet distribution with parameters (8,8,5) as in [7].The hyperparameters have been fixed to 2 n = ψ 2 = 10 −4 , ∀n, ∀ .These two sequences are processed using different unmixing strategies that are compared to the proposed UsTNCM algorithm (denoted by UsTNCM for unsupervised temporal NCM).More precisely, we have considered the following unmixing algorithms: (i) VCA+FCLS: [18,25] and (ii) UsLMM [15].Table 1 shows the performance for each of the different algorithms.This table shows a reduced performance for VCA+FCLS mainly because of the absence of pure pixels and the variation of the endmember variances with respect to spectral bands.UsLMM and UsTNCM provide good results when processing the first sequence with slightly better results for UsLMM.The interest in the proposed approach is highlighted when processing I 2 where it shows the best results.For this sequence, the pixels are highly mixed and processing the images independently as in UsLMM does not lead to good results.Conversely, the UsTNCM processes the whole sequence jointly and benefits from the temporal correlation between successive images.

Realistic synthetic image
Due to the absence of a sequence of real images, this section considers a 35 × 35 synthetic image that has been generated based on a real Madonna image 2 .This real image was acquired in 2010 by the Hyspex hyperspectral scanner over Villelongue, France (00 03'W and 4257'N).The dataset contains L = 160 spectral bands recorded from the visible to near infrared (400 − 1000nm) with a spatial resolution of 0.5m [26].Three endmember means and abundance maps have been estimated from this real image when considering the VCA+FCLS algorithm.The synthetic tth image was then generated by considering (1) when smoothly varying these estimated parameters with respect to time (with T = 20).The endmember variances vary linearly with respect to spectral bands such that σ 2 1 = 6 × 10 (−5) and σ 2 L = 5 × 10 (−4) .The hyperparameters have been fixed to 2 n = ψ 2 = 10 −4 , ∀n, ∀ .Fig. 1 shows four images (t = 1, 8, 15, 20) of the obtained sequence in true color (bands 72, 33, and 18).These images show three elements: tree, soil and grass that are located in many clusters inside the simplex [7] (which is a challenging scenario as for I 2 ).The performance of the different algorithms is shown in Table 2. UsTNCM shows the best results since it processes the whole sequence jointly.This allows the proposed approach to better capture the temporal variation of both the endmembers and the abundances.

CONCLUSIONS
This paper proposed an unsupervised Bayesian algorithm for unmixing successive hyperspectral images while accounting for temporal and spatial variability.Spatial variability is introduced by considering a normal compositional model.Temporal variability is included as the mean of a Gaussian process that ensures a smooth temporal evolution for the abundances and endmember means between successive images.The parameters are then estimated using an MCMC approach.The algorithm proposed in this paper showed good results when applied to unmixing of synthetic images particularly for high mixed scenarios and in the absence of pure pixels.Future work will address the issue of hyperparameter estimation for the proposed Bayesian model.The validation of the proposed approach on real hyperspectral images would also be very interesting.
containing the endmember means of the tth image and a n;t = [a 1n;t , • • • , a Rn;t ]

Table 1 .
Results on synthetic data.The results should be multiplied by (×10 −2 ) .

Table 2 .
Results on the synthetic Madonna sequence.The results should be multiplied by (×10 −2 ) .