Hyperspectral unmixing with spectral variability using adaptive bundles
and double sparsity
Abstract
Spectral variability is one of the major issue when conducting hyperspectral unmixing. Within a given image composed of some elementary materials (herein referred to as endmember classes), the spectral signature characterizing these classes may spatially vary due to intrinsic component fluctuations or external factors (illumination). These redundant multiple endmember spectra within each class adversely affect the performance of unmixing methods. This paper proposes a mixing model that explicitly incorporates a hierarchical structure of redundant multiple spectra representing each class. The proposed method is designed to promote sparsity on the selection of both spectra and classes within each pixel. The resulting unmixing algorithm is able to adaptively recover several bundles of endmember spectra associated with each class and robustly estimate abundances. In addition, its flexibility allows a variable number of classes to be present within each pixel of the hyperspectral image to be unmixed. The proposed method is compared with other stateoftheart unmixing methods that incorporate sparsity using both simulated and real hyperspectral data. The results show that the proposed method can successfully determine the variable number of classes present within each class and estimate the corresponding class abundances.
I Introduction
Hyperspectral analysis has received an increasing attention because of its high spectral resolution, which enables a variety of objects to be identified and classified [1]. Mixed pixels caused by the presence of multiple objects within a single pixel adversely affect the performance of hyperspectral analysis [2]. To address this problem, a wide variety of spectral unmixing methods have been developed over the last decades [3, 4, 5, 6]. Spectral unmixing methods aim at decomposing a mixed spectrum into a collection of reference spectra (known as endmembers) characterizing the macroscopic materials present in the scene, and their respective proportions (known as abundances) in each image pixel [2]. Despite the large number of developed spectral unmixing methods, there are still major challenges for accurate estimates of endmember signatures and abundances [3]. Among these challenges, endmember variability may lead to large amounts of errors in abundance estimates [7]. It results from the fact that each endmember can rarely be represented by a unique spectral signature. Conversely, it is subject to socalled spectral variability, e.g., caused by variations in the acquisition process, the intensity of illumination or other physical characteristics of the materials [8, 9]. Taking this endmember variability into account during the spectral unmixing process is one of keys for successful application of spectral unmixing [10].
The methods that incorporate endmember variability can be categorized into two main approaches (see [10, 7, 11] for recent overviews). The first approach relies on the definition of a set of multiple spectral signatures, referred to as endmember bundles, to characterize each endmember class. Endmember bundles can be collected from field campaign or can be extracted from data itself using endmember bundle extraction methods [12, 13, 9]. The advantage of this approach is that the method can use a priori information representing endmember variability. Endmember bundles can be validated by experts in order to provide accurate representation about endmember variability [12]. Although traditional methods incorporating endmember bundles (e.g. [14]) are known to be computationally expensive, more efficient methods have been recently developed and have shown great potential [15, 16, 17, 18]. However, as pointed out in [19], it is unlikely that the endmember bundles completely represent endmember variability present in an image. Such incomplete endmember bundles may lead to poor estimates of abundances. The second approach uses physical or statistical descriptions of the endmember variability. More precisely, these methods describes the endmember variability thanks to a statistical distribution [20] or by incorporating additional variability terms in the mixing model [21, 22, 23]. The advantage of this approach results from the adaptive learning of the endmember variability. Indeed, stateoftheart methods such as those recently introduced in [21, 22, 23] enable endmember spectra to spatially vary within each pixel in order to describe endmember variability. This is important since the endmember spectra to be used for the abundance estimation can be different between pixels. However, estimating endmember variability without a priori knowledge is a challenging task, especially when large amounts of endmember variability are present in an image. In addition, the statistical distribution or additional terms used in these methods may be overly simplified to represent endmember variability.
Both approaches demonstrate benefits and drawbacks. A natural question arises: is it possible to combine the strong advantages of both approaches to robustly represent endmember variability? This paper addresses the question and introduces a novel spectral unmixing method that bridges the gap between the aforementioned two approaches with the help of a double sparsitybased method inspired by [24]. Specifically, the proposed method aims at adaptively recovering endmember spectra within each pixel to describe endmember variability while incorporating available a priori information. The proposed method is closely related to the existing methods. Thus, the main contributions of this paper are threefold: 1) to propose a novel spectral unmixing method that incorporates endmember bundles and generates adaptive endmember spectra within each pixel, 2) to give a systematic review of related work and show the relationship between the proposed method and existing methods, 3) to provide comparison between the proposed method and other sparsitybased methods.
This paper is organized as follows. Section II describes related works and existing methods, while highlighting their inherent drawbacks. In Section III, a novel mixing model that incorporates endmember variability is proposed and its relationships with existing methods are discussed. Section IV introduces an associated unmixing algorithm designed to recover the endmember classes, adaptive bundles and abundances. Section V and Section VI show experimental results obtained from simulated data and real hyperspectral images. Finally, conclusions are drawn in Section VII.
Ii Related works and issues raised by existing methods
Iia Conventional (variabilityfree) linear mixing model
Let denote the spectrum measured at the th pixel of an hyperspectral image. According to the linear mixing model (LMM), the observed spectrum of the th pixel is approximated by a weighted linear combination of endmember spectra and abundance fractions
(1) 
where is the matrix of the spectral signatures associated with the endmember classes, is the abundance fractions of the pixel and represents noise and modeling error. LMM is generally accompanied by abundance nonnegativity constraint (ANC) and the abundance sumtoone constraint (ASC)
(2) 
where is the abundance fraction of the th class in the th pixel. This model implicitly relies on two assumptions: i) each endmember class is described by a unique spectrum and ii) the endmember matrix is fixed and commonly used to unmix all the pixels of a given image. In other words, LMM does not account for spectral variability. However, as discussed previously, this is likely unrealistic because spectral variability is naturally observed in hyperspectral images, e.g., because of variations in illumination or physical intrinsic characteristics of materials [19].
IiB Linear mixing models incorporating endmember bundles
Multiple endmember spectral mixture analysis (MESMA) [14] allows the variability of the endmember spectrum representative of each class and a varying number of endmember classes present within each pixel. Although MESMA has been widely used for a variety of applications [14, 25], it owns several major limitations: i) it is highly computationally expensive because MESMA needs to test a large number of combinations of endmember spectra [26], ii) MESMA tends to select an overestimated number of endmember classes because it uses the reconstruction error to select the appropriate combination of endmember spectra [27] and iii) the performance of MESMA may significantly decrease when endmember spectra (or bundles) within each class do not completely represent the spectral variability [19].
To overcome these limitations, recent works have proposed a new class of methods that incorporate all endmember bundles defined as [16, 15]
(3) 
where represents a set of endmember spectra (aka bundle) characterizing the th class, is the number of endmember spectra in the th class and is the total number of endmember spectra of all classes with . Generalizing LMM in (1), those methods firstly model a given observed pixel spectrum with respect to all spectra in endmember bundles and corresponding multiple abundances
(4) 
where is multiple abundance fractions corresponding to each spectrum of the endmember bundles . As for LMM, ASC or ANC can also be imposed to . As a second step, multiple abundance fractions are summed within each class to generate a single abundance fraction for each class
(5) 
with
(6) 
where is a column vector of ones and represent a dimensional vector whose components are zeros. While these two steps are conducted separately in [16, 15], they can be also considered jointly within a multitask Gaussian process framework [19, 17]. Even if these methods have been shown to be effective, a large number of endmember spectra within each class may be redundant. In such case, following a model selection inspiration, Veganzones et al. introduce a complementary sparsity regularization on the multiple abundance vectors [15]
(7) 
where represents the elementwise comparison, is the norm, is the norm which is known to promote sparsity. Once the multiple abundance vector has been estimated, it is normalized in order to reduce the effects of multiplicative factors and satisfy ASC. Following the same approach, further sparsity can be imposed using norm [28] or reweighted approaches [29, 30]. Overall, this sparsity property allows the selection of a smaller number of endmember spectra. However, it may not lead to the selection of a smaller number of endmember classes. Conversely, to promote sparsity on the number of endmember classes, one strategy consists in formulating the unmixing problem through a sparse group lasso [31]
(8) 
where is the th column of , is the elementwise product and thus extracts the elements in belonging to the th class. This approach has the great advantage of promoting sparsity in both the number of endmember spectra and the number of endmember classes. Another strategy relies on the concept of “social sparsity” that can exploit the structure of endmember bundles more explicitly [18]. The method assumes that can be partitioned into groups representing each endmember class, leading to the optimization problem
(9) 
where is the norm and is the th multiple abundance fraction of the th pixel. Finally, abundances associated each endmember class can be obtained by summing the multiple abundances within each class as in (5). This method can be considered as a generalized model since, by adjusting the values of , it boils down to the group lasso, the elitist lasso or the fractional case [18].
However, all aforementioned sparsitybased methods still suffer from the following limitations:

Physically unrealistic abundance fractions: They explicitly generate unrealistic multiple abundances corresponding to each spectrum in endmember bundles.

Lack of adaptability to describe endmember variability: ASC imposed on does not allow a consistant description of the endmembers within each pixel. In addition, it cannot capture adaptive and hierarchical structure of endmember spectra for each pixel.
To overcome these two shorcomings, the present paper capitalizes on this abundant literature to design a new multiple endember mixing model introduced in the next section.
Iii Multiple endmember mixing models
Iiia Memm
The proposed model relies on 3 main ingredients, namely endmember bundles, bundling coefficients and abundances. According to this model, each endmember bundle is mixed to provide a suitable and adaptive endmember spectrum used to unmix a given pixel. The proposed multiple endmember mixing model (MEMM) is defined as
(10) 
where gathers socalled bundling coefficients of the th pixel which decompose the endmember signatures according to the endmember bundles for the considered pixel. To enforce the bundle structure, the bundling coefficients associated with the pixel is defined as the following blockdiagonal matrix
(11) 
where is the bundling coefficients for the th class at the th pixel. Each bundling coefficient must be nonnegative and the bundling vector is expected to be sparse. Indeed multiple endmember spectra within each class are usually redundant and only a few endmember spectra within each class should be enough to unmix a pixel. This property can be induced by considering the following bundling constraints
(12) 
where is the norm that counts the number of nonzero elements and is the maximum number of nonzero elements in , i.e., the maximum number of endmembers to be used within each class to describe the pixel. The abundance nonnegativity constraint (ANC) and the abundance sumtoone constraint (ASC) are usually imposed. In addition, in this work, complementary sparsity is imposed on each abundance vector, i.e.,
(13) 
where is the number of endmember classes to be used to decompose the image pixel.
IiiB MEMMs
The sparsity constraint (12) applied to can be slightly modified to obtain another meaningful set of constraints
(14) 
The resulting model, referred to as MEMM in what follows, is designed to generate at most one scaled endmember spectrum for each class.
IiiC Relationships between MEMM and existing models
IiiC1 Memm and MESMA
When the sparsity constraint on the abundances in (13) is not considered, the optimization problem associated with the MEMM model described in paragraph IIIB is equivalent to MESMA and sparse MESMA [32]. Unlike MESMA that considers the reconstruction error to determine the optimal combination of endmember classes within each pixel, MEMM incorporates the sparsity constraint to select the optimal combination. This prevents a larger number of endmember classes to be selected for each pixel.
IiiC2 MEMM and pixelwise endmember variability models
By denoting the equivalent endmember matrix associated with the th pixel, MEMM models the observed pixel spectra as
(15) 
where can be interpreted as a set of spatially varying endmember spectra. This approach has been also adpoted in recent works to incorporate endmember variability as additive factors [22], multiplicative factors [21, 33] or a combination of additive and multiplicative factors [23]. In particular, when in (3), the endmember bundles are reduced to unique endmember spectra characterizing each class. The associated bundling coefficient matrix is diagonal where each coefficient scales the corresponding endmember spectrum (). Thus, MEMM generalizes the recently introduced extended linear mixing model [21].
However, MEMM is different from the aforementioned methods since it resorts to a priori information (i.e., endmember bundles) to model the endmember variability. More precisely, MEMM describes the admissible variability within an endmember class as the convex cone spanned by the corresponding bundles. As a consequence, per se, MEMM offers an adaptive description of the spectral variability even when predefined endmember bundles do not completely capture this variability within each class.
IiiC3 MEMM and sparsitybased unmixing methods
By setting , MEMM in (10) can be rewritten as
(16) 
similarly to the existing models discussed in Section II. The main difference is that MEMM enables the multiple abundances to be decomposed into bundling coefficients within each class and abundances , resulting into a bilayer description of the abundances. This hierarchical decomposition has been also adopted by unmixing methods based on multilayer nonnegative matrix factorization (MLNMF) [34]. However, each layer induced by MEMM (i.e., bundling matrix and abundance vector) has a clear and meaningful role. Moreover, when the existing methods impose ASC onto , they assume that mixed spectra belong to the simplex spanned by the endmember bundles. However, this is a limited assumption since observed spectra may be outside the simplex, e.g., when affected by variations in illumination. Conversely, MEMM imposes ASC only onto the abundances of endmember classes and enables the bundling coefficients to scale the endmember signature, e.g., to capture variability induced by varying illumination (see experiments in Section VII). In addition, MEMM complements this bilayer hierarchy with a twofold structured, physicallymotivated sparsity imposed on the multiple abundance vector . This bilevel sparsity has the significant advantage of reducing overfitting and one may expect a significant improvement of stability and interpretability of the abundance estimates.
Finally, the intrinsic structure of MEMM is also similar to the recently developed methods based on robust constrained matrix factorization [35] and kernel archetypoid analysis [36]. These methods model a set of pixel spectra as
(17) 
where is a matrix gathering a set of coefficients, is the abundance matrix and is the error and noise matrix. Sparsity (induced by ASC or the use of pseudonorm) and nonnegativity constraints are imposed onto each column of and can be interpreted as synthetic endmember spectra. These methods use the subset of whole image pixels to generate synthetic endmember spectra that are fixed within each image. On the other hand, MEMM uses the subset of endmember bundles to generate synthetic endmember spectra that may be different for each pixel.
Iv MEMMbased unmixing algorithm
Unmixing according to the proposed MEMM can be formulated as the minimization problem
(18) 
This minimization problem is similar to the double sparsityinducing method proposed in [24]. Using an alternative formulation, the minimization problem can be written as the following nonconvex minimization problem:
(19) 
with
(20)  
(21)  
(22) 
where and are parameters which control the balance between the data fitting term and the sparse regularizations, is the indicator function on the set (i.e., when whereas when ), and is the simplex defined by the ASC and ANC. Solving this optimization problem is challenging since the regularization functions and are nonconvex and nonsmooth. However, it can be tackled thanks to the proximal alternating linearized minimization (PALM) [37]. With guarantees to converge to a critical point, PALM iteratively updates the parameters and by alternatively minimizing the objective function with respect to (w.r.t.) these parameters, i.e., by solving the following proximal problems
(23) 
The pseudocode for MEMM is shown in Algorithm 1 and these two steps are described in what follows.
Iva Optimization w.r.t.
To optimize only w.r.t. the diagonal entries in , the objective function can be rewritten with the following decomposition
(24) 
where
This leads to the following updating rule
where . Using similar computations as in [37], this can be conducted as
(25) 
where represents a step size for each iteration. The proximal operator associated with can be computed using the approach [37]. Finally, the bundling matrix can be reconstructed as where generates the block diagonal matrix from the vector .
IvB Optimization with respect to
To optimize w.r.t. , the objective function can be rewritten using the decomposition
(26)  
(27) 
where . Thus, updating the abundance vector can be formulated as
where . Using the proximal operator, this can be written as
(28) 
where represents a step size for each iteration. Moreover the proximal mapping associated with can be performed using the method developed in [38].
IvC Initialization and stopping rule
MEMM requires initial estimates and of the abundance vector and bundling matrix, respectively. To do so, first, MEMM estimates a multiple abundance vector using a stateoftheart LMMbased unmixing method (e.g., FCLS, see line 3). Then an initial estimate of the single abundance vector is computed according to (5) (see line 4). Finally, the bundling matrix is arbitrarily initialized as the corresponding scaling factor (see line 5, where stands for the elementwise division). Once initial estimates have been obtained, and are iteratively updated in lines 8–10. The algorithm stops when the difference between updated and previous values of the objective function is smaller than a predetermined threshold.
V Experiments using simulated data
First, the relevance of the proposed MEMM and its variant MEMM has been evaluated thanks to experiments conducted on simulated datasets.
Va Generation of bundles
First, spectra were selected from the USGS spectral library. The spectra were chosen so that the minimum angle between any two spectra was larger than . This pruning prevented synthetic endmember bundles to overlap each other. Second, endmember bundles () were designed by randomly generated endmember spectra for each endmember bundle using the approach proposed in [39]. These bundles, depicted in Fig. 1, were used for generating the following two simulated data.
VA1 Simulated dataset 1 (SIM1)
The first dataset, referred to as SIM1 in what follows, was generated using MESMA. A mixed spectrum was generated using the following 5 steps. First, the number of endmember classes was randomly determined in the set (). Second, a random combination of endmember classes was selected. Third, one spectrum within each selected endmember class was randomly chosen. Forth, the abundances of the selected endmember classes were randomly generated using a Dirichlet distribution to jointly ensure ANC and ASC. Finally, a mixed spectrum was generated by a linear combination of endmember spectra of the selected endmember classes and the randomly generated abundances. A set of mixed spectra were generated in this study. Different amounts of additive Gaussian noise with corresponding signaltonoise ratios of dB, dB and dB were considered to the mixed spectra.
VA2 Simulated dataset 2 (SIM2)
SIM2 was generated using MEMM. A mixed spectrum was generated similarly as done in SIM1. The main difference is the bundling coefficients in MEMM. In order to generate the bundling coefficients, the number of spectra was randomly chosen in the set (). The bundling coefficients of the randomly selected spectra were generated from a Dirichlet distribution. A mixed spectrum was generated by a linear combination of spectra, the bundling coefficients and the abundances. As for SIM1, mixed spectra were generated and different amounts of Gaussian noise were also added to SIM2.
VB Compared methods
MEMM and MEMM were compared with other 5 methods that incorporate endmember bundles and promote sparsity: fully constrained least squares (FCLS) [40], sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL) [41], alternating angle minimization (AAM) [42] and methods based on group lasso and elitist lasso [18]. Note that FCLS and SUnSAL can incorporate endmember bundles by considering the model in (4). FCLS and SUnSAL were chosen for comparing the proposed methods with most widely used unmixing methods that promote sparsity. AAM was selected for comparison because it is a latest variant of MESMA and it is more computationally efficient than MESMA while achieving good performance. Group lasso and elitist lasso were also included for comparison because they can incorporate the structure of endmember bundles as well as MEMM. MEMM, MEMM and the methods based on group lasso and elitist lasso require initial estimates of abundances. In order to fairly compare methods, this study used abundances estimated by FCLS for the initial estimates required for the methods. SUnSAL, the methods based on group lasso or elitist lasso and MEMM require a parameter or controlling sparsity regularization. MEMM requires two parameters and . In order to fairly compare the methods, these parameters were empirically determined in the set for each simulated data so that the selected values produced the highest SRE. Parameter sensitivity of the proposed method is reported in the supplementary document [43]. Finally, computational time were also discussed.
VC Performance criteria
The main objective of this experiment was to assess the algorithm performance while selecting a combination of endmember classes and spectra, and estimating abundances corresponding to endmember classes and spectra. Three criteria were chosen for quantitative validation of the methods. To evaluate the quality of the reconstruction, one defines the signaltoreconstruction error (SRE) per endmember class and per endmember spectrum as [44]
(29)  
where and are the actual and estimated abundance vectors of all pixels, and are the actual and estimated multiple abundance vectors of all pixels.
Second, the number of nonzero abundances were used to evaluate the sparsity level per endmember class and per endmember spectrum recovered by the methods, i.e., [44, 32]
(30)  
As in [44], abundances smaller than were considered as zero abundances.
Finally, to validate the performance in selecting a relevant combination of endmember classes or spectra, one defines the distance between the two actual and estimated supports (DIST) [45, 32]
(31)  
where and are true and estimated support sets (i.e., indexes of nonzero values), represents the total number of elements in the set and stands for the intersection operator. The figures of merit DIST and DIST evaluated the distance between two supports of endmember classes and the distance between two supports of endmember spectra, respectively.
SNR  FCLS  AAM  SUnSAL  Group lasso  Elitist lasso  MEMM  MEMM  

SIM1  30dB  18.0507  18.7656  18.9285  19.4215  19.2333  14.5739  18.6904 
40dB  22.9074  25.8112  23.0496  23.589  23.303  14.8229  23.3532  
50dB  27.358  33.8279  27.5013  27.4184  26.8979  14.1579  28.0017  
SIM2  30dB  16.1865  15.9501  16.5182  18.0736  16.4315  12.1614  16.4785 
40dB  21.6942  19.4906  21.6671  22.0263  21.3809  13.6682  22.1381  
50dB  25.6617  22.9952  25.5168  25.9215  24.7702  14.9405  26.5442 
SNR  FCLS  AAM  SUnSAL  Group lasso  Elitist lasso  MEMM  MEMM  

SIM1  30dB  1.1752  2.2035  1.7234  2.0398  1.5188  0.4570  0.7692 
40dB  2.8611  4.6569  3.0262  3.0357  2.9957  0.5171  2.7653  
50dB  3.6312  12.7796  3.7692  3.7591  3.7576  0.1700  3.5721  
SIM2  30dB  1.0084  0.2324  1.4405  1.8935  1.3043  2.0312  0.4313 
40dB  2.3323  1.4226  2.4976  2.5741  2.48  2.1441  2.2462  
50dB  3.2175  3.3104  3.3273  3.3266  3.3666  3.0151  3.2304 
VD Results
SRE per class was calculated for each method and reported in Table I. For SIM1 with dB, AAM performed best among all methods. The performance of AAM, however, was degraded as SNR became lower. For data with dB, the results derived from AAM were worse than those derived from SUnSAL, elitist lasso and group lasso. In addition, AAM performed poorly compared with other methods in SIM2. This showed that the MESMAbased approach (AAM) was less effective when given endmember bundles did not completely represent endmember variability present in the data and SNR of the data was low (dB). This finding was also observed in [19]. MEMM produced better results for SIM2 with dB than the other sparsitymethods and produced comparable results with dB and dB. MEMM performed poorly, compared with all methods. SRE per spectrum was also calculated from each method and reported in Table II. Compared with SRE per class, SRE per spectrum was very low for all methods. This showed that the exact recovery of multiple abundances was challenging under conditions where a large number of endmember spectra were present within each class.
SNR  FCLS  AAM  SUNSAL  Group lasso  Elitist lasso  MEMM  MEMM  

SIM1  30dB  4.8  4.25  3.32  4.59  6.83  3.05  1.84 
40dB  4.62  4.29  3.82  4.54  5.64  3.64  2.14  
50dB  4.34  4.05  4.24  4.58  5  3.55  1.99  
SIM2  30dB  5.02  4.45  3.56  4.72  6.63  4.27  2.03 
40dB  5.21  4.59  4.16  5.02  5.6  4.25  2.5  
50dB  4.97  4.49  4.87  4.72  5.54  3.9  2.28 
SNR  FCLS  AAM  SUNSAL  Group lasso  Elitist lasso  MEMM  MEMM  

SIM1  30dB  10.83  4.25  10.52  43.12  11.28  3.04  5.07 
40dB  14.25  4.29  12.39  29.88  15.81  3.55  8.72  
50dB  18.66  4.05  18.07  23.57  21.31  3.45  11.13  
SIM2  30dB  11.98  4.45  13.08  66.64  13.16  4.22  5.26 
40dB  16.79  4.59  15.6  45.99  19.67  4.15  10.88  
50dB  25.34  4.49  24.84  65.97  29.76  3.78  19.38 
The average number of nonzero abundances in and were shown in Table III and Table IV, respectively. For the number of nonzero abundances per class, MEMM performed best among all methods for both SIM1 and SIM2. This showed that the sparse constraint used in MEMM successfully led to the selection of smaller numbers of endmember classes. For the number of nonzero abundances per spectrum, AAM and MEMM produced the smaller number of nonzero abundances than other methods. This was because these methods selected at most one spectrum within each class and enforced greater sparsity when selecting endmember spectra.
SNR  FCLS  AAM  SUNSAL  Group lasso  Elitist lasso  MEMM  MEMM  

SIM1  30dB  0.5893  0.5489  0.3939  0.5866  0.6967  0.3491  0.1195 
40dB  0.5689  0.5163  0.4616  0.5594  0.6368  0.4498  0.1265  
50dB  0.5137  0.4828  0.5035  0.5479  0.5795  0.4180  0.0758  
SIM2  30dB  0.5725  0.5156  0.3775  0.5487  0.6538  0.5025  0.1843 
40dB  0.5789  0.4843  0.4521  0.5534  0.5952  0.4938  0.162  
50dB  0.5345  0.4696  0.5247  0.5185  0.5825  0.4179  0.1033 
SNR  FCLS  AAM  SUNSAL  Group lasso  Elitist lasso  MEMM  MEMM  

SIM1  30dB  0.9167  0.767  0.853  0.94  0.8963  0.7493  0.8115 
40dB  0.8914  0.6707  0.8587  0.9137  0.8881  0.7351  0.774  
50dB  0.8649  0.5622  0.8578  0.8971  0.8911  0.6586  0.7431  
SIM2  30dB  0.9039  0.8426  0.8476  0.9332  0.8915  0.8729  0.8142 
40dB  0.8791  0.7672  0.8525  0.9167  0.883  0.8175  0.7993  
50dB  0.8738  0.6949  0.8696  0.9172  0.8816  0.7277  0.8095 
The errors in the support sets were shown in Table V and Table VI. When estimating the support set per class, MEMM outperformed other methods for both SIM1 and SIM2. This showed that the double sparsity imposed by MEMM also led to the appropriate selection of the combination of endmember classes for each pixel. The performance of MEMM, however, was degraded when estimating the support set per spectrum. Other methods also performed poorly for the support set per spectrum. This shows that when multiple highly correlated endmember spectra are present within each class, the existing methods experience difficulty to select an optimal combination of endmember spectra.
Data  FCLS  AAM  SUNSAL  Group lasso  Elitist lasso  MEMM  MEMM 

SIM1  0.3006  517.7361  0.155  0.9517  1.2694  0.4846  2.6625 
SIM2  0.1966  549.1244  0.2616  0.8426  1.2122  2.2311  2.5705 
Finally, the computational times of all methods were shown in Table VII. MEMM was more computationally expensive than FCLS, SUnSAL, Group lasso and Elitist lasso. The proposed methods, however, were computationally cheaper than AAM because they did not need to test a large number of the combinations of endmember spectra.
Vi Experiments using real data
Via Description of real hyperspectral images
The methods were finally compared on two real hyperspectral images. The first hyperspectral image was acquired in June 2016, over the city of SaintAndré, France, during the MUESLI airborne acquisition campaign. The image was composed of spectral bands. The spectral bands affected by noise (between m and m) were removed, leading to spectral bands. In the image scene, spatially discrete objects were present. In this study, each spatially discrete region was assumed to be composed of a single endmember class. Endmember bundles were extracted from each region. As a consequence, large amounts of endmember variability were expected to be present within each spatially discrete region associated with a particular class and mixed pixels were expected to be located in the boundary of these spatially discrete regions. Moreover, the scene of interest was composed of two flight lines under significantly different illumination conditions, as shown in Fig. 1(a). Thus, this image (referred to as MUESLI image) was used to evaluated whether the methods could accurately estimate abundances when large amounts of endmember variability were present. From this image, endmember bundles composed of a total of spectral signatures representing spatially discrete objects were extracted using the nDimensional Visualizer provided by the ENVI software. These bundles are represented in Fig. 3(top). Some endmember classes present in the studied area were affected by different illumination conditions. Unlike simulated data, ground truth was not available. Estimated abundances were qualitatively validated by visual inspection of the abundance maps.
The second image was acquired over Moffett Field, CA, USA, by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS). The image, depicted in Fig. 1(b), initially comprised spectral bands. After the noisy spectral bands were removed, spectral bands remained. The area of interest, composed of a lake and a vegetated coastal area, was considered in many previous studies, e.g., to assess the performance of unmixing methods. Thus, this second image (referred to as AVIRIS image) was used to test whether the proposed method could perform at least as well as the existing methods to analyze this widely used image scene. As for the AVIRIS image, endmember bundles were extracted and qualitative validation was conducted because of the lack of ground truth. The extracted bundles are depicted in Fig. 4(left).
ViB Results
For both images, the proposed method was compared with FCLS, SUnSAL, AAM, the methods based on group lasso and elitist lasso. The parameters (, and ) required for the methods were empirically determined by qualitatively evaluating abundances derived from different values of the parameters. Finally, synthetic endmember bundles generated by MEMM were compared with endmember bundles initially used for unmixing.
Abundance maps estimated by the methods on the MUESLI image are depicted in Fig. 5. These maps show that the abundances estimated by MEMM were more consistent at the boundary affected by the different illumination conditions. This showed that MEMM was more robust to the different illumination conditions than concurrent methods. The abundances estimated by MEMM were also high for each endmember class and showed less noisy. This suggested that MEMM also promoted more sparsity than other methods.
Fig. 6 shows the abundance maps estimated for the AVIRIS image. All methods except elitist lasso generated similar abundances. Abundances estimated by elitist lasso were different because it was designed to use a larger number endmember classes for unmixing each pixel. MEMM produced similar abundances when compared with FCLS, group lasso. This showed that MEMM could perform at least as well as other sparsitybased methods to unmix this wellstudied test site. MEMM, however, generated more noisy abundance maps. This showed that the initial estimates of abundances used in MEMM did not lead to an optimal combination of endmember spectra and the optimal abundances.
Finally, the endmember bundles recovered by MEMM were compared with endmember bundles initially used to unmix both hyperspectral images, see Fig. 3 and Fig. 4. In Fig. 3, the synthetic endmember bundles filled the gaps that were present in the original endmember bundles. The extended endmember bundles showed more detailed spectral variability within each class in terms of both spectrum amplitudes and shapes also in Fig. 4. This enabled MEMM to generate adaptive endmember spectra within each pixel and estimate more accurate abundances even when initial endmember bundles did not completely represent endmember variability.
Vii Conclusion
This paper proposed a multiple endmember mixing model that bridges the gap between endmember bundlebased method and data drivenbased methods. MEMM appeared to be superior to the existing methods as follows: i) it incorporated endmember bundles to generate adaptive endmember spectra for each pixel, ii) it had explicit physical meaning and generated hierarchical structure of endmember spectra, iii) it imposed double sparsity for the selection of both endmember classes and endmember spectra. MEMM were tested and compared to the stateoftheart methods using simulated data and real hyperspectral images. MEMM showed comparable results for estimating abundances while it outperformed other methods in terms of selecting a set of endmember classes within each pixel. This paper deeply focused on sparsity constraints for both bundling coefficients and abundances. However, other constraints (e.g., spatial constraints) can be easily incorporated in the proposed unmixing framework. Future work will consist of reducing the computational complexity of the proposed method.
Acknowledgment
The authors would like to thank Prof. Jose M. BioucasDias, Dr. Rob Heylen and Dr. Travis R. Meyer for sharing the MATLAB codes of the SUnSAL, AAM and goup/elitist lasso unmixing algorithms.
References
 [1] A. Plaza, J. A. Benediktsson, J. W. Boardman, J. Brazile, L. Bruzzone, G. CampsValls, J. Chanussot, M. Fauvel, P. Gamba, A. Gualtieri, M. Marconcini, J. C. Tilton, and G. Trianni, “Recent advances in techniques for hyperspectral image processing,” Remote Sens. Environment, vol. 113, pp. 110–122, 2009.
 [2] N. Keshava and J. F. Mustard, “Spectral unmixing,” IEEE Signal Process. Mag., vol. 19, no. 1, pp. 44–57, 2002.
 [3] J. M. BioucasDias, A. Plaza, N. Dobigeon, M. Parente, D. Qian, P. Gader, and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regressionbased approaches,” IEEE J. Sel. Topics Appl. Earth Observations Remote Sens., vol. 5, no. 2, pp. 354–379, 2012.
 [4] N. Dobigeon, J. Y. Tourneret, C. Richard, J. C. M. Bermudez, S. McLaughlin, and A. O. Hero, “Nonlinear unmixing of hyperspectral images: Models and algorithms,” IEEE Signal Process. Mag., vol. 31, no. 1, pp. 82–94, 2014.
 [5] R. Heylen, M. Parente, and P. Gader, “A review of nonlinear hyperspectral unmixing methods,” IEEE J. Sel. Topics Appl. Earth Observations Remote Sens., vol. 7, no. 6, pp. 1844–1868, 2014.
 [6] T. Uezato, M. Fauvel, and N. Dobigeon, “Hyperspectral image unmixing with LiDAR dataaided spatial regularization,” IEEE Trans. Geosci. Remote Sens., 2018.
 [7] A. Zare and K. C. Ho, “Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing,” IEEE Signal Process. Mag., vol. 31, no. 1, pp. 95–104, 2014.
 [8] R. J. Murphy, S. T. Monteiro, and S. Schneider, “Evaluating classification techniques for mapping vertical geology using fieldbased hyperspectral sensors,” IEEE Trans. Geosci. Remote Sens., vol. 50, no. 8, pp. 3066–3080, 2012.
 [9] T. Uezato, R. J. Murphy, A. Melkumyan, and A. Chlingaryan, “A novel endmember bundle extraction and clustering approach for capturing spectral variability within endmember classes,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 11, pp. 6712–6731, 2016.
 [10] B. Somers, G. P. Asner, L. Tits, and P. Coppin, “Endmember variability in spectral mixture analysis: A review,” Remote Sens. Environment, vol. 115, no. 7, pp. 1603–1616, 2011.
 [11] L. Drumetz, J. Chanussot, and C. Jutten, “Endmember variability in spectral unmixing: recent advances,” in Proceedings of IEEE Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), 2016, pp. 1–4.
 [12] C. A. Bateson, G. P. Asner, and C. A. Wessman, “Endmember bundles: a new approach to incorporating endmember variability into spectral mixture analysis,” IEEE Trans. Geosci. Remote Sens., vol. 38, no. 2, pp. 1083–1094, 2000.
 [13] B. Somers, M. Zortea, A. Plaza, and G. P. Asner, “Automated extraction of imagebased endmember bundles for improved spectral unmixing,” IEEE J. Sel. Topics Appl. Earth Observations Remote Sens., vol. 5, no. 2, pp. 396–408, 2012.
 [14] D. A. Roberts, M. Gardner, R. Church, S. Ustin, G. Scheer, and R. O. Green, “Mapping chaparral in the Santa Monica mountains using multiple endmember spectral mixture models,” Remote Sens. Environment, vol. 65, no. 3, pp. 267–279, 1998.
 [15] M. A. Veganzones, L. Drumetz, G. Tochon, M. DallaMura, A. J. Plaza, J. M. Bioucas Dias, and J. Chanussot, “A new extended linear mixing model to address spectral variability,” in Proc. IEEE GRSS Workshop Hyperspectral Image SIgnal Process.: Evolution in Remote Sens. (WHISPERS), 2014, pp. 1–4.
 [16] M. A. Goenaga, M. C. TorresMadronero, M. VelezReyes, S. J. Van Bloem, and J. D. Chinea, “Unmixing analysis of a time series of Hyperion images over the Guanica dry forest in Puerto Rico,” IEEE J. Sel. Topics Appl. Earth Observations Remote Sens., vol. 6, no. 2, pp. 329–338, 2013.
 [17] T. Uezato, R. J. Murphy, A. Melkumyan, and A. Chlingaryan, “Incorporating spatial information and endmember variability into unmixing analyses to improve abundance estimates,” IEEE Trans. Image Process., vol. 25, no. 12, pp. 5563–5575, 2016.
 [18] T. R. Meyer, L. Drumetz, J. Chanussot, A. L. Bertozzi, and C. Jutten, “Hyperspectral unmixing with material variability using social sparsity,” in Proc. IEEE Int. Conf. Image Process. (ICIP), 2016, pp. 2187–2191.
 [19] T. Uezato, R. J. Murphy, A. Melkumyan, and A. Chlingaryan, “A novel spectral unmixing method incorporating spectral variability within endmember classes,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 5, pp. 2812–2831, 2016.
 [20] A. Halimi, N. Dobigeon, and J. Tourneret, “Unsupervised unmixing of hyperspectral images accounting for endmember variability,” IEEE Trans. Image Process., vol. 24, no. 12, pp. 4904–4917, 2015.
 [21] L. Drumetz, M. A. Veganzones, S. Henrot, R. Phlypo, J. Chanussot, and C. Jutten, “Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability,” IEEE Trans. Image Process., vol. 25, no. 8, pp. 3890–3905, 2016.
 [22] P. A. Thouvenin, N. Dobigeon, and J. Y. Tourneret, “Hyperspectral unmixing with spectral variability using a perturbed linear mixing model,” IEEE Trans. Signal Process., vol. 64, no. 2, pp. 525–538, 2016.
 [23] D. Hong, N. Yokoya, J. Chanussot, and X. Zhu, “Learning a lowcoherence dictionary to address spectral variability for hyperspectral unmixing,” in Proc. IEEE Int. Conf. Image Process. (ICIP), 2017, pp. 1–5.
 [24] R. Rubinstein, M. Zibulevsky, and M. Elad, “Double sparsity: Learning sparse dictionaries for sparse signal approximation,” IEEE Trans. Signal Process., vol. 58, no. 3, pp. 1553–1564, 2010.
 [25] P. E. Dennison and D. A. Roberts, “Endmember selection for multiple endmember spectral mixture analysis using endmember average RMSE,” Remote Sens. Environment, vol. 87, no. 23, pp. 123–135, 2003.
 [26] L. Tits, B. Somers, R. Heylen, P. Scheunders, and P. Coppin, “Improving the efficiency of MESMA through geometric unmixing principles,” in SPIE Remote Sensing, vol. 8892, 2013.
 [27] L. Demarchi, F. Canters, J. C. Chan, and T. Van de Voorde, “Multiple endmember unmixing of chris/proba imagery for mapping impervious surfaces in urban and suburban environments,” IEEE Trans. Geosci. Remote Sens., vol. 50, no. 9, pp. 3409–3424, 2012.
 [28] J. Sigurdsson, M. O. Ulfarsson, and J. R. Sveinsson, “Hyperspectral unmixing with regularization,” IEEE Trans. Geosci. Remote Sens., vol. 52, no. 11, 2014.
 [29] C. Y. Zheng, H. Li, Q. Wang, and C. P. Chen, “Reweighted sparse regression for hyperspectral unmixing,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 1, pp. 479–488, 2016.
 [30] W. He, H. Zhang, and L. Zhang, “Total variation regularized reweighted sparse nonnegative matrix factorization for hyperspectral unmixing,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 7, pp. 3909–3921, 2017.
 [31] M. Iordache, J. BioucasDias, and A. Plaza, “Hyperspectral unmixing with sparse group lasso,” in Proc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS), 2011, pp. 3586–3589.
 [32] F. Chen, K. Wang, and T. F. Tang, “Spectral unmixing using a sparse multipleendmember spectral mixture model,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 10, pp. 5846–5861, 2016.
 [33] T. Imbiriba, R. A. Borsoi, and J. C. M. Bermudez, “Generalized linear mixing model accounting for endmember variability,” arXiv preprint arXiv:1710.07723, 2017.
 [34] R. Rajabi and H. Ghassemian, “Spectral unmixing of hyperspectral imagery using multilayer NMF,” IEEE Geosci. Remote Sens. Lett., 2015.
 [35] N. Akhtar and A. Mian, “Rcmf: Robust constrained matrix factorization for hyperspectral unmixing,” IEEE Trans. Geosci. Remote Sens., 2017.
 [36] C. Zhao, G. Zhao, and X. Jia, “Hyperspectral image unmixing based on fast kernel archetypal analysis,” IEEE J. Sel. Topics Appl. Earth Observations Remote Sens., 2017.
 [37] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized minimization for nonconvex and nonsmooth problems,” Mathematical Programming, vol. 146, no. 12, pp. 459–494, 2014.
 [38] K. Anastasios, B. Stephen, and C. Volkan, “Sparse projections onto the simplex,” in Proc. Int. Cong. Machine Learning, 2013.
 [39] P. A. Thouvenin, N. Dobigeon, and J. Y. Tourneret, “Online unmixing of multitemporal hyperspectral images accounting for spectral variability,” IEEE Trans. Image Process., vol. 25, no. 9, pp. 3979–3990, 2016.
 [40] J. M. BioucasDias and M. A. T. Figueiredo, “Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing,” in Proc. IEEE GRSS Workshop Hyperspectral Image SIgnal Process.: Evolution in Remote Sens. (WHISPERS), 2010, pp. 1–4.
 [41] M. D. Iordache, J. M. BioucasDias, and A. Plaza, “Sparse unmixing of hyperspectral data,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 6, pp. 2014–2039, 2011.
 [42] R. Heylen, A. Zare, P. Gader, and P. Scheunders, “Hyperspectral unmixing with endmember variability via alternating angle minimization,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 8, 2016.
 [43] T. Uezato, M. Fauvel, and N. Dobigeon, “Hyperspectral unmixing with spectral variability using adaptive bundles and double sparsity – complementary results,” University of Toulouse, IRIT/INPENSEEIHT, Tech. Rep., April 2018. [Online]. Available: {http://dobigeon.perso.enseeiht.fr/papers/Uezato_TechReport_2018b.pdf}
 [44] M.D. Iordache, J. M. BioucasDias, and A. Plaza, “Collaborative sparse regression for hyperspectral unmixing,” IEEE Trans. Geosci. Remote Sens., vol. 52, no. 1, pp. 341–354, 2014.
 [45] M. Elad, Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. Springer, 2010.