Hyperspectral unmixing with spectral variability using adaptive bundles and double sparsity

Hyperspectral unmixing with spectral variability using adaptive bundles
and double sparsity

Tatsumi Uezato, Mathieu Fauvel, 
and Nicolas Dobigeon, 
Part of this work has been funded by EU FP7 through the ERANETMED JC-WATER program, MapInvPlnt Project ANR-15-NMED- 0002-02 and by the MUESLI IDEX ATS project, Toulouse INP.T. Uezato and N. Dobigeon are with University of Toulouse, IRIT/INP-ENSEEIHT, CNRS, 2 rue Charles Camichel, BP 7122, 31071 Toulouse Cedex 7, France (e-mail: {Tatsumi.Uezato, Nicolas.Dobigeon}@enseeiht.fr).M. Fauvel is with University of Toulouse, INRA, DYNAFOR, BP 32607, Auzeville-Tolosane, 31326 Castanet Tolosan, France (e-mail: mathieu.fauvel@ensat.fr).
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 state-of-the-art 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.

Hyperspectral imaging, spectral unmixing, sparse unmixing, endmember variability.

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 so-called 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, state-of-the-art 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 sparsity-based 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 sparsity-based 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

Ii-a Conventional (variability-free) 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 non-negativity constraint (ANC) and the abundance sum-to-one 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].

Ii-B 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 multi-task 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 element-wise 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 element-wise 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 sparsity-based methods still suffer from the following limitations:

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

  2. 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

Iii-a 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 so-called 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 block-diagonal 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 non-negativity constraint (ANC) and the abundance sum-to-one 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.

Iii-B 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.

Iii-C Relationships between MEMM and existing models

Iii-C1 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 III-B 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.

Iii-C2 MEMM and pixel-wise 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 pre-defined endmember bundles do not completely capture this variability within each class.

Iii-C3 MEMM and sparsity-based 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 bi-layer 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 bi-layer hierarchy with a twofold structured, physically-motivated 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 MEMM-based unmixing algorithm

Unmixing according to the proposed MEMM can be formulated as the minimization problem

(18)

This minimization problem is similar to the double sparsity-inducing method proposed in [24]. Using an alternative formulation, the minimization problem can be written as the following non-convex 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.

Iv-a 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 .

Iv-B 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].

1:
2:Initialization: and .
3:Set using an unmixing method (e.g. FCLS).
4:
5:
6:Main procedure:
7:while  the stopping criterion is not satisfied do
8:     
9:     
10:     
11:end while
12:
Algorithm 1 Algorithm for MEMM-based unmixing

Iv-C 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 state-of-the-art LMM-based 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 element-wise 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.

Fig. 1: Synthetically generated endmember bundles.

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.

V-a 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.

V-A1 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 signal-to-noise ratios of dB, dB and dB were considered to the mixed spectra.

V-A2 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.

V-B 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.

V-C 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 signal-to-reconstruction 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
TABLE I: SRE per endmember class ().
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
TABLE II: SRE per endmember spectrum ().

V-D 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 MESMA-based 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 sparsity-methods 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
TABLE III: Sparsity level per endmember class (). Reference: in SIM1 and in SIM2.)
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
TABLE IV: Sparsity level per endmember class (). Reference: in SIM1 and in SIM2.

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
TABLE V: Distance between actual and estimated supports per endmember class ().
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
TABLE VI: Distance between actual and estimated supports per endmember spectrum ().

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
TABLE VII: Computational time for unmixing pixels using endmember bundles of endmember spectra.

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.

(a)
(b)
Fig. 2: Real hyperspectral images: (a) MUESLI image and (b) AVIRIS image.

Vi Experiments using real data

Vi-a 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 Saint-André, 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 n-Dimensional 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.

Fig. 3: First row: endmember bundles used for unmixing the MUESLI image. Second row: endmember bundles generated by MEMM.

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).

Fig. 4: First column: endmember bundles used for unmixing the AVIRIS image. Second column: endmember bundles generated by MEMM.

Vi-B 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.

Fig. 5: MUESLI image: estimated abundance maps. From top to bottom: building, road, shrub, crop land 1, crop land 2 and grass.

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: AVIRIS image: abundance maps. From top to bottom: vegetation, soil and water.

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 sparsity-based methods to unmix this well-studied 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 bundle-based method and data driven-based 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 state-of-the-art 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. Bioucas-Dias, 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. Camps-Valls, 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. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, D. Qian, P. Gader, and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based 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 data-aided 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 field-based 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 image-based 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. Dalla-Mura, 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. Torres-Madronero, M. Velez-Reyes, 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 low-coherence 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. 2-3, 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. Bioucas-Dias, 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 multiple-endmember 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. 1-2, 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. Bioucas-Dias 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. Bioucas-Dias, 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/INP-ENSEEIHT, Tech. Rep., April 2018. [Online]. Available: {http://dobigeon.perso.enseeiht.fr/papers/Uezato_TechReport_2018b.pdf}
  • [44] M.-D. Iordache, J. M. Bioucas-Dias, 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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
169206
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description