Sparse Nonnegative Matrix Factorization for Hyperspectral Unmixing Based on Endmember Independence and Spatial Weighted Abundance

Abstract

Hyperspectral image unmixing is an important task for remote sensing image processing. It aims at decomposing the mixed pixel of the image to identify a set of constituent materials called endmembers and to obtain their proportions named abundances. Recently, number of algorithms based on sparse nonnegative matrix factorization (NMF) have been widely used in hyperspectral unmixing with good performance. However, these sparse NMF algorithms only consider the correlation characteristics of abundance and usually just take the Euclidean structure of data into account, which can make the extracted endmembers become inaccurate. Therefore, with the aim of addressing this problem, we present a sparse NMF algorithm based on endmember independence and spatial weighted abundance in this paper. Firstly, it is assumed that the extracted endmembers should be independent from each other. Thus, by utilizing the autocorrelation matrix of endmembers, the constraint based on endmember independence is to be constructed in the model. In addition, two spatial weights for abundance by neighborhood pixels and correlation coefficient are proposed to make the estimated abundance smoother so as to further explore the underlying structure of hyperspectral data. The proposed algorithm not only considers the relevant characteristics of endmembers and abundances simultaneously, but also makes full use of the spatial-spectral information in the image, achieving a more desired unmixing performance. The experiment results on several data sets further verify the effectiveness of the proposed algorithm.

1. Introduction

With the continuous improvement of science and technology, remote sensing images have been developed by leaps and bounds. Hyperspectral image (HSI), a kind of remote sensing image, has attracted the attention of many researchers due to its rich spectral information [1,2,3]. HSI contains dozens or even hundreds of continuous bands, and each pixel can be extracted to a complete spectral curve that reflects the characteristics of ground objects. Thus, it has been successfully applied to many aspects, such as agriculture, meteorology, exploration and so on. However, there is a variety of materials mixed in a pixel, i.e., the phenomenon of spectral mixing. The spectral mixing will seriously affect the subsequent processing of HSI, such as classification [4], detection [5,6], etc. Therefore, the decomposition of mixed pixels for HSI becomes more and more crucial.

Mixed pixel decomposition of HSI, referred to as hyperspectral unmixing (HU), is to decompose the mixed pixel into several materials (endmembers) and to obtain their proportions (abundances). The models for HU are mainly divided into linear mixing model (LMM) and nonlinear mixing model [7]. The LMM assumes that each photon acts on only one material with no interactions among them. It is easy to be solved and can meet the basic needs of research, which has been widely used in HU. This paper is also based on the LMM for unmixing research. The process of HU generally contains three steps: endmember number estimation, endmember extraction and abundance estimation. There are many traditional algorithms for endmember number estimation, such as hyperspectral signal identification by minimum error (Hysime) [8], virtual dimensionality (VD) [9], minimum noise fraction [10]. Based on different assumptions, the endmember extraction algorithms can be grouped into the methods based on pure pixel, minimum volume, statistics, etc. The well-known algorithms for endmember extraction contain pixel purity index (PPI) [11], N-FINDR [12], vertex component analysis (VCA) [13] and so on. Besides, two kinds of constraints for abundances exist in HU, i.e., the abundance nonnegative constraint (ANC) and the abundance sum to one constraint (ASC). The former requires that the abundance value is nonnegative, and the latter asks the sum of abundance value of each pixel to be one. Nonnegative constrained least squares algorithm and fully constrained least squares (FCLS) algorithm [14] are two unmixing methods based on LMM combined with different constraints of abundances. The algorithms of VCA and FCLS are often adopted as the initial method for endmember extraction and abundance estimation in the experiment. In addition, numerous different algorithms have been proposed, including the geometric analysis method [15], filtering method [16], deep learning [17], etc. Based on the characteristics of hyperspectral images and some prior information, researchers present a series of abundance estimation algorithms.

The sparse unmixing (SU) algorithm provides an important direction for the research. Studies have stated that as for HSI, not all endmembers of the image participate in the mixing of each pixel, but only a few [18]. Correspondingly, the abundance is sparse. In addition, due to the widely available spectral library, the SU algorithm is able to measure the situation when the pure pixels do not exist in the image [19]. Lots of classical sparse unmixing algorithms appear in view of the different priors and understandings for HSI. The commonly used norms for sparsity constraint are L₁ regularization, L_1/2 regularization and L_2,1 regularization. The solution obtained by L_1/2 regularization is sparser than that by L₁ regularization [20]. The sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL) [19] algorithm explores the L₁ norm as the constraint for abundance based on the alternating direction method of multipliers to solve the sparse regression problem. However, it just analyzes the hyperspectral data and does not incorporate the spatial information. In order to make use of the spatial information in the image, the sparse unmixing via variable splitting augmented Lagrangian and total variation (SUnSAL-TV) algorithm adds the total variation of regularization to the model of SUnSAL algorithm [21]. Furthermore, the collaborative SUnSAL (CLSUnSAL) [22] considers the pixels in HSI sharing the same active set of atoms in the library and employs the L_2,1 regularization for collaborative sparse regression. Local collaborative sparse regression unmixing algorithm imposes the collaborative sparsity among neighboring pixels, which assumes that the neighboring pixels share the same active set of endmembers [23]. The spatial discontinuity-weighted sparse unmixing [24] adopts a spatial discontinuity weight for SUnSAL to preserve the spatial details of abundance. The joint local block grouping with the noise-adjusted principal component analysis sparse method [25] utilizes the local block grouping to get spatial information and draws the representative spatial correlations obtained by the noise-adjusted principal component analysis for unmixing. Besides, with the intention to take advantage of the spatial information, spectral and spatial weighting factors are exploited by the spectral-spatial weighted sparse unmixing framework, imposing the sparsity of solution [26]. Although the SU algorithm has made some achievements, the spectra collected in the spectral library and the image show there are still many differences that raise doubts.

The nonnegative matrix factorization (NMF)-based unmixing approach is another highlighted branch for HSI, which has attracted the attentions of researchers due to its successful applications in many fields. The main task of NMF [27] is decomposing a nonnegative matrix into the product of two nonnegative matrices to reduce the dimension of high-dimensional data. Since its goal is similar to that of spectral unmixing, the NMF model has been widely employed in unmixing. However, the NMF model has an ill-conditioned problem, which tends to fall into a local optimal solution. Therefore, it is necessary to add some specific constraints of endmembers and abundances to the NMF model based on the characteristic of HU. Casalino et al. added both sparsity constraint and spatial information to a new nonnegative matrix under approximation model [28]. Scholars put forward lots of excellent NMF-based algorithms that achieve good effects for HU [29,30,31,32]. Qian et al. introduced the L_1/2 sparsity NMF for HU through the L_1/2 regularization raised in [20] to make the solution sparser and more accurate [18]. Miao et al. presented the minimum-volume constraint NMF (MVCNMF) by a geometric constraint of endmembers [33]. Li et al. performed three steps of HU together in the robust collaborative NMF (CoNMF) [34]. Inspired by the manifold learning, Lu et al. added the graph regularized constraint in NMF (GLNMF) to fully exploit the latent manifold structure of HSI [35]. Wang et al. divided the pixels of HSI into groups based on their correlation and used a spatial group sparsity regularization term for abundance to unmix [36]. Under the self-learning semi-supervised framework, Wang et al. integrated the prior information into NMF as the constraints of endmembers and abundances in the unmixing process [37]. Xiong et al. brought a nonconvex non-separable sparse NMF approach via a generalized minimax concave sparse regularization preserving the convexity of NMF for each variable [38].

Inspired by the advantage of NMF model, we develop a sparse NMF unmixing algorithm based on endmember independence and spatial weighted abundance for HSI (EASNMF). The purpose of the proposed algorithm is to make the extracted endmembers independent of each other and obtain smooth abundances. For the endmembers, it is considered that the more independent the endmembers are, the better they can characterize the HSI. Thus, the constraint of endmembers via autocorrelation matrix is added to the NMF model. In addition, only a subset of endmembers participates in the mixing of pixels, which leads to the sparsity of abundances. Therefore, we adopt the sparse constraint for abundances and introduce a weight in view of spatial information to make the abundances smoother. Furthermore, in order to exploit the latent manifold structure of the HSI data, manifold regularization is also employed in our model. The results on both the simulated data set and the real data set demonstrate the effectiveness of the proposed EASNMF algorithm with the flowchart shown in Figure 1. In general, the EASNMF algorithm not only puts forward the appropriate constraints based on the characteristics of endmembers and abundances, but also fully integrates the spatial-spectral information for HU.

The rest of this paper is arranged as follows. Section 2 is the related work, Section 3 introduces the proposed EASNMF algorithm in detail, followed by Section 4 of the experiment, and finally, the conclusion is in Section 5.

2. Related Work

2.1. LMM

The unmixing algorithms often rely on the establishment of the mixing models, and the LMM is an important mixing model. Let $Y\in {ℝ}^{L\times P}$ represent the HSI observation matrix of L bands and P pixels, $E\in {ℝ}^{L\times K}$ indicate the endmember matrix with K endmembers, $A\in {ℝ}^{K\times P}$ mean the abundance matrix and $N\in {ℝ}^{L\times P}$ refer to the noise matrix, thus the LMM can be formed as follows:

$$Y=EA+N$$

(1)

Two constraints of abundances including the ANC and ASC are below:

$$\mathrm{ANC}:a_{\mathrm{ij}}\ge 0,\forall i,j$$

(2)

$$\mathrm{ASC}:{\displaystyle \sum _{i=1}^K{a}_{ij}}=1,j=1:P$$

(3)

where $a_{ij}$ is the abundance value of i-th endmember at j-th pixel of the HSI.

2.2. NMF

NMF, a powerful tool for statistical analysis, is one commonly used model for HU due to its significant advantages. The standard form of the NMF model based on the cost function of Euclidean distance is as follows:

$$\underset{E,A}{\min }\frac{1}{2}{\Vert Y-EA\Vert }_F^2,\hspace{1em}\mathrm{s}.\mathrm{t}.E\ge 0,A\ge 0$$

(4)

where ${\Vert \cdot \Vert }_F$ denotes the Frobenius norm. The purpose of NMF is to seek two nonnegative matrices decomposed from the HSI data.

To optimize the function with respect to $E$ and $A$ in Equation (4), the updated rules of the iterative algorithm proposed in [27] are below:

$$E\leftarrow E\odot (Y{A}^T)\oslash (EA{A}^T)$$

(5)

$$A\leftarrow A\odot (E^{T}Y)\oslash (E^{T}EA)$$

(6)

where ${(\cdot )}^T$ refers to the transpose of matrix, $\odot$ and $\oslash$ are the elementwise multiplication and division, called Hadamard product and Hadamard division, respectively.

However, due to the nonconvex objective function of the NMF model in Equation (4), it suffers from the problem of nonunique solution. Therefore, to reduce the feasible solution set, some constraints based on the characteristics of endmembers and abundances are introduced to the NMF model. There are various constraints to solve this problem, such as manifold constraint [39,40], sparseness constraint [41], low-rank constraint [42], smooth constraint [43] and so on. These NMF-based approaches are all named constrained NMF, with the formulation as follows:

$$\underset{E,A}{\min }\frac{1}{2}{\Vert Y-EA\Vert }_F^2+\lambda f(E)+\mu \phi (A)\hspace{1em}\mathrm{s}.\mathrm{t}.A\ge 0,1_K^{T}A=1_P^T$$

(7)

where $f(E)$ and $\phi (A)$ are the constraints of endmembers and abundances, and the two parameters λ and μ separately adjust the effects of the corresponding regularization term in Equation (7).

3. Sparse NMF for Hyperspectral Unmixing Based on Endmember Independence and Spatial Weighted Abundance

In this section, the proposed sparse unmixing algorithm based on endmember independence and spatial weighted abundance with manifold regularization is introduced in detail. The raised EASNMF algorithm can get the independent endmembers and the smooth abundances, which fully exploits the spatial-spectral information and the intrinsic geometrical characteristics of HSI data.

3.1. Endmember Independence Constraint

As we know, the solution space of the NMF model is very large. In addition, the endmembers are very important to unmixing research, which will affect the effect of HU. Therefore, we can utilize this characteristic of endmembers as the prior knowledge added to NMF model. This way, the accurate endmembers can be received to further improve the unmixing effect. The HSI data is formed by different endmembers with a certain proportion, and it is easy to find that the endmembers should be independent of each other. For the independence, the autocorrelation matrix can be adopted to constrain the endmembers. If the endmembers are independent from each other, their autocorrelation matrix should be a diagonal matrix. That is, the off-diagonal elements of its autocorrelation matrix should be as close to 0 as possible. Therefore, the NMF model with endmember independence constraint is as follows:

$$\begin{array}{l}\underset{E,A}{\min }\frac{1}{2}{\Vert Y-EA\Vert }_F^2+\alpha ({\Vert E^{T}E\Vert }_1-{\Vert E\Vert }_F^2)\\ \mathrm{s}.\mathrm{t}.A\ge 0,1_K^{T}A=1_P^T\end{array}$$

(8)

where α is the parameter to balance the data fidelity and endmember independence term. The second term refers to the sum of the off-diagonal elements of the autocorrelation matrix for endmembers, i.e., the difference between the sum of all the elements (the first sub-term) and the sum of the diagonal elements (the second sub-term). The purpose for the second term in Equation (8) is to make the endmembers independent of each other as much as possible; that is, the correlation between different endmembers should be as small as possible.

3.2. Abundance Sparse and Spatial Weighted Constraint

Studies have shown that most of the mixed pixels are mixtures with only a few endmembers on the scene [41]. That is to say, the mixed pixel is likely to be the superposition of only a few endmembers, not all endmembers. Thus, the corresponding abundance is sparsity, which can be considered an intrinsic property of HU. Therefore, the sparsity constraint as an effective tool has been introduced to HU. As mentioned before, the L_1/2 regularizer proposed by [20] is proved to provide a sparse and accurate result. Taking the sparsity of abundance into consideration, we add the sparse constraint of abundance into the model, which is formed as follows:

$$\begin{array}{l}\underset{E,A}{\min }\frac{1}{2}{\Vert Y-EA\Vert }_F^2+\alpha ({\Vert E^{T}E\Vert }_1-{\Vert E\Vert }_F^2)+\beta {\Vert A\Vert }_{\frac{1}{2}}\\ \mathrm{s}.\mathrm{t}.A\ge 0,1_K^{T}A=1_P^T\end{array}$$

(9)

where β is the weight parameter to adjust the effect of the last term in Equation (9), and ${\Vert A\Vert }_{1/2}={\displaystyle \sum _{i,j}^{K,P}{(a_{ij})}^{1/2}}$.

Moreover, the neighboring pixels are more likely to have similar fractional abundance values, which is considered spatial structure information. This information can be constructed as a weight matrix for abundances to make full use of. Suppose the pixel $y_j$, whose corresponding abundance value is $a_j$, is one neighbor of the pixel $y_i$, and there are m neighbors for pixel $y_i$. For each iteration of abundance, the abundance average of the neighborhood for each pixel is calculated to construct a weight matrix $W=[w_1,w_2,\dots ,w_i]\in {ℝ}^{K\times P}$ for next iteration. The element in weight matrix W is computed as follows:

$$w_i^{(k+1)}=\frac{1}{{(\frac{1}{m}{\displaystyle \sum _{j=1}^m{a}_j})}^{(k)}+eps}\hspace{1em}{y}_j\in \mathbb{N}(y_i)$$

(10)

where eps is a predetermined positive constant. Here the Euclidean distance is adopted to calculate the similarity of pixels in the image, and then m pixels with smallest values are chosen as the neighbors to obtain the element of weight matrix W in Equation (10). It is hoped that if the spectral signatures of pixels are similar, their abundance values should be similar. The model with the weight matrix $W$ is expressed as below:

$$\begin{array}{l}\underset{E,A}{\min }\frac{1}{2}{\Vert Y-EA\Vert }_F^2+\alpha ({\Vert E^{T}E\Vert }_1-{\Vert E\Vert }_F^2)+\beta {\Vert W\odot A\Vert }_{\frac{1}{2}}\\ \mathrm{s}.\mathrm{t}.A\ge 0,1_K^{T}A=1_P^T\end{array}$$

(11)

where $\odot$ is the term-by-term Hadamard product.

In this part, the priors of sparseness and spatial information are integrated into the NMF model to shrink the solution space and further promote the unmixing performance. However, it just considers the sparse characteristic for unmixing and neglects the intrinsic geometrical structure of HSI. Therefore, it is necessary to further explore the potential characteristic of HSI data for unmixing.

3.3. Manifold Regularization Constraint

As is well known, HSI is a kind of high-dimensional data. Recently, researchers showed that the hyperspectral data vary smoothly along the geodesics of the data manifold and tend to lie on a low-dimensional subspace embedded in the high-dimensional data space [35]. Moreover, the manifold learning finds the representation in low-dimensional manifold space for high-dimensional space data. It can dig into the essence of data and discover its inherent laws. In Equation (11), only the sparse characteristic and the Euclidean structure of hyperspectral data are taken into account as we have posted before. Therefore, it is necessary to introduce the intrinsic manifold structure into the proposed model to render better performance of HU.

There are P pixels in HSI and each pixel can be considered a data point. Thus, the nearest neighbor graph is constructed by each pixel as its vertices and its weight matrix is denoted as $W_g$. The weight between two pixels $y_i$ and $y_j$ is defined as follows:

$$w_{g_{ij}}=\{\begin{array}{cc}\hfill corrcoef(y_i,y_j)& y_i\in \mathbb{N}(y_j)or{y}_j\in \mathbb{N}(y_i)\hfill \\ \hfill 0& otherwise\hfill \end{array}$$

(12)

Here $corrcoef(\cdot )$ means the correlation coefficient and it is calculated by $\frac{cov(y_i,y_j)}{\sqrt{var(y_i)var(y_j)}}$ where $cov(\cdot )$ and $var(\cdot )$ separately mean the covariance and variance. That is, if the pixel $y_j$ is a neighbor of pixel $y_i$, the weight between these two pixels is obtained by computing their correlation coefficient. The correlation coefficient is usually used to describe the degree of correlation between two variables in statistics, whose absolute value is between 0 and 1. Generally speaking, the closer its absolute value is to 1, the greater the correlation between two variables is.

Furthermore, based on the analysis before, if two pixels $y_i$ and $y_j$ are close in the original space, their representations $a_i$ and $a_j$ in the new space should also be close [35]. For this purpose, the manifold constraint is proposed as below:

$$\begin{array}{ll}\frac{1}{2}{\displaystyle \sum _{i,j=1}^P{\Vert a_i-a_j\Vert }^2{w}_{g_{ij}}}& ={\displaystyle \sum _{i=1}^P{a}_i^T{a}_i{d}_{ii}}-{\displaystyle \sum _{i,j=1}^P{a}_i^T{a}_j{w}_g{}_{{}_{ij}}}\hfill \\ & =Tr(AD{A}^T)-Tr(A{W}_g{A}^T)=Tr(AL{A}^T)\hfill \end{array}$$

(13)

where $Tr(\cdot )$ indicates the trace of the matrix, $d_{ii}={\displaystyle \sum _{j=1}^P{w}_{g_{ij}}}$ and $L=D-W_g$. Then, incorporating the manifold regularization into the model, the finial objective function is exhibited as follows:

$$\begin{array}{l}\underset{E,A}{\min }\frac{1}{2}{\Vert Y-EA\Vert }_F^2+\alpha ({\Vert E^{T}E\Vert }_1-{\Vert E\Vert }_F^2)+\beta {\Vert W\odot A\Vert }_{\frac{1}{2}}+\frac{\gamma }{2}Tr(AL{A}^T)\\ \mathrm{s}.\mathrm{t}.A\ge 0,1_K^{T}A=1_P^T\end{array}$$

(14)

where γ acts as the penalty parameter to control the manifold regularization term.

According the updated rule in [20], the iterative solution of Equation (14) is presented as follow.

$$E\leftarrow E\odot (Y{A}^T)\oslash (EA{A}^T+2\alpha E{I}_1-2\alpha E)$$

(15)

$$A\leftarrow A\odot (E^{T}Y+\gamma A{W}_g)\oslash (E^{T}EA+\frac{\beta }{2}{W}^{\frac{1}{2}}\odot A^{-\frac{1}{2}}+\gamma AD)$$

(16)

where I₁ is the matrix with all 1 elements. Considering the ASC, a simple and effective technique in [35,41] is employed. When updating the abundance A by Equation (16), the matrices Y and E will be replaced by $Y_f$ and $E_f$ by adding a row as the inputs to achieve the ASC, which are defined as below:

$$Y_f=\left[\begin{array}{c}Y\\ \epsilon 1_P^T\end{array}\right]\hspace{1em}{E}_f=\left[\begin{array}{c}E\\ \epsilon 1_K^T\end{array}\right]$$

(17)

where the parameter ε controls the influence of ASC, and in our experiment, it is set to be 15, which will be mentioned later. Then taking the ASC into consideration, the iterative criterion for abundance is as follows:

$$A\leftarrow A\odot (E_f^T{Y}_f+\gamma A{W}_g)\oslash (E_f^T{E}_{f}A+\frac{\beta }{2}{W}^{\frac{1}{2}}\odot A^{-\frac{1}{2}}+\gamma AD)$$

(18)

The whole algorithm has been described in detail. Our algorithm not only proposes the appropriate constraints based on the characteristics of endmembers and abundances simultaneously, but also makes full use of the spatial-spectral information in the image, achieving a desired unmixing performance. Algorithm 1 briefly presents the solution to Equation (14) and summarizes the aforementioned description; the values of parameters α, β and γ will be discussed in detail later.

Algorithm 1 Sparse NMF for HU Based on Endmember Independence and Spatial Weighted Abundance

1. Input: The hyperspectral image Y, the number of endmember K, the parameters α, β and γ.

2. Output: Endmember matrix E and abundance matrix A.

3. Initialize E and A by VCA-FCLS algorithm, W by Equation (10), Wg by Equation (12), and D.

4. Repeat:

5. Update E by Equation (15).

6. Augment Y and A separately to get Yf and Af.

7. Update A by Equation (18).

8. Update W by Equation (10).

9. Until stopping criterion is satisfied.

4. Experiments Results

This section mainly describes a series of experiments designed to evaluate the effectiveness of the proposed EASNMF method. We first introduce the evaluation metrics and the data sets including the simulated data set and the real data set. Then the experimental setting is explained. Finally, the results of the EASNMF algorithm and the comparisons composed of MVCNMF, L_1/2-NMF, GLNMF and CoNMF on both the simulated data set and the real data set are displayed and analyzed.

4.1. Performance Evaluation Criteria

In the experiment of this paper, two widely adopted evaluation metrics are used to measure the accuracy separately for endmembers and abundances. First is the spectral angle distance (SAD), which can qualify the similarity of the extracted endmember and its real endmember by calculating their spectral angle. When the SAD value is smaller, the performance for endmember extraction is better. Besides, the SAD is not affected by spectral scale either, whose definition is as below:

$$SAD=\arccos (\frac{E^T\overline{E}}{\Vert E\Vert \Vert \overline{E}\Vert })$$

(19)

where $E$ and $\overline{E}$ are the real endmember and the extracted endmember.

The error between the abundance and its real abundance is computed by the root-mean-square error (RMSE) in the experiment, which is formed as follows:

$$RMSE={(\frac{1}{P}{\Vert A-\overline{A}\Vert }^2)}^{\frac{1}{2}}$$

(20)

where $A$ represents the real abundance and $\overline{A}$ denotes the estimated abundance. When the estimated abundance is close to the real abundance, the error is small corresponding to the good performance for abundance estimation.

4.2. Data Sets

There are three data sets employed in the experiment to evaluate the effectiveness of the EASNMF algorithm, which contains two simulated data sets and one real data set called Cuprite.

• Simulated data set 1:

The first simulated data set provided by Dr. M. D. Iordache and Prof. J. M. Bioucas-Dias is generated with 100 × 100 pixels and nine spectra randomly selected from the USGS spectral library [44]. Its abundance maps are shown in Figure 2 for illustrative purposes. This data set has 224 bands and its abundance follows a Dirichlet distribution. Owing to its good spatial homogeneity, it becomes the data set widely used in HU [25,38]. Finally, the Gaussian noise with 30 dB was added.

• Simulated data set 2:

The second data set is provided in the HyperMix tool [45] with 100 × 100 pixels and 221 bands for testing the spectral unmixing algorithms. There are nine endmembers randomly selected from the USGS library after removing certain bands for this data set. The fractional abundance maps associated with each endmember are displayed in Figure 3. Similar to the simulated data set 1, the Gaussian noise with 30 dB was included in the experiment.

• Cuprite data set

The scene adopted in the real data experiment is named Cuprite data set, which was captured by airborne visible infra-red imaging spectrometer (AVIRIS) in 1997. Since the Cuprite data set contains rich minerals that are usually highly mixed, it is a popular data set for researchers to verify the effectiveness of the HU algorithm [23,37,41]. A sub-image with 250 × 191 pixels is selected from the scene containing 224 spectral bands ranging from 400 to 2500 nm. Figure 4 shows the real data set (left) and the reference maps (right) produced by Tricorder 3.3 software product in 1995, which maps different minerals presented in the mining district. Although it is inappropriate to compare the distribution map directly with Cuprite data set, the reference map can still be used in qualitative analysis of abundance map evaluation. Besides, its resolutions of spectral and spatial are approximately 10 nm and 20 m. The bands 1–2, 105–115, 150–170, and 223–224 affected by water vapor and atmospheric were removed remaining 188 bands. The agreement for the endmember number is not available, a frequently used and widely recognized number is twelve, including alunite, andradite, buddingtonite, dumortierite, kaolinite1, kaolinite2, montmorillonite, muscovite, nontronite, pyrope, sphene and chalcedony.

4.3. Compared Algorithms

In our experiment, four unmixing algorithms listed as follows are selected as the comparisons for the proposed EASNMF algorithm:

• L_1/2-NMF algorithm: it extends the NMF method by incorporating the L_1/2 sparsity constraint, which provides a more sparser and accurate results [18].
• GLNMF algorithm: it incorporates the manifold regularization into sparsity NMF, which can preserve the intrinsic geometrical characteristic of HSI data during the unmixing process [35].
• MVCNMF algorithm: it adds the minimum volume constraint into the NMF model and extracts the endmember from highly mixed image data [33].
• CoNMF algorithm: it performs all stages involved in HU process including the endmember number estimation, endmember estimation and abundance estimation [34].

4.4. Initializations and Parameter Settings

There are several important issues that need to be addressed in advance. The details of these issues are discussed below.

• Initialization: the initialization of endmember and abundance is the first issue. In our experiment, we choose the VCA-FCLS algorithm, one basic method for endmember extraction and abundance estimation, as our initialization method to speed up the optimization. VCA algorithm [13] exploits two facts to extract the endmembers: the endmembers are the vertices of a simplex and the affine transformation of a simplex is also a simplex. FCLS algorithm, a quadratic programming technique, is developed to address the fully constrained linear mixing problems, which uses the efficient algorithm to simultaneously implement both the ASC and ANC [14].
• Stopping criterion: it is another important issue and two stopping criteria are adopted for the optimization, i.e., error tolerance and maximum iteration number. When any stop condition is reached, the algorithm stops. When the error is successively within the limits of tolerance, a predefined value, the iteration is stopped. The error tolerance is set as 1.0 × 10⁻⁴ for a simulated data set and 1.0 × 10⁻³ for the real data set in our experiment. The times of iteration meet the maximum iteration number, the optimization ends. The maximum iteration number is set as 1.0 × 10⁶ in experiment.
• ANC and ASC: for the abundance, its initial value obtained by VCA-FCLS algorithm is generally nonnegative. Thus, according to the update rule recorded in Equations (15) and (16), the E and A are obviously nonnegative. Besides, considering the ASC, the A adopted by Equation (18) also satisfies the constraint. Moreover, the parameter ε in Equation (17) controls the convergence rate of ASC. When its value is large, it will lead to an accurate result but with lower convergence rate. As in many papers [35,41], the parameter ε is set as 15 in the experiments for desired tradeoff.
• Parameter setting: there are three parameters in the proposed model, i.e., α, β, γ. They separately control the independence constraint of the endmember, abundance sparse constraint, and the manifold constraint, which will be analyzed in detail in next part of the experiment.
• Endmember number: the endmember number is one of the crucial processes in HU, which is another independent topic. In our experiment, it is considered a topic that does not have much relation to this paper and it is assumed to be known. In fact, the algorithms of HySime [8] and VD [9] could be adopted to estimate the number of endmembers. Hysime algorithm [8] is a new minimum mean square error-based approach to infer the signal subspace in hyperspectral imagery. In the experiment, we can also analyze the number of endmembers around the number estimated by Hysime algorithm via the reconstruction error.
• Computational complexity: here, we analyze the computational complexity of the proposed EASNMF algorithm. It is noticeable that the matrix W_g is sparse and there are m nonzero elements in each row. Therefore, the floating-point addition and multiplication for AW_g in Equation (16) cost mPK times. Additionally, the computing cost of A^−1/2 is (PK)². Except for these costs, the other three floating-point calculation times for each iteration are listed in

  Table 1. The floating-point calculation times for each iteration in EASNMF algorithm.

  	Update E	Update A	Total
  Addition	LPK + (2L + P)K2 + 2LK	LPK + (L + P)K2 + (4 + m)PK	2LPK + (3L + 2P)K2 + 2LK + (4 + m)PK
  Multiplication	LPK + (2L + P)K2 + LK	LPK + (L + P)K2 + (3 + m)PK	2LPK + (3L + 2P)K2 + LK + (3 + m)PK
  Division	LK	PK	(L + P)K

  .

4.5. Experiment on Simulated Data Set 1

In this section, we evaluate the proposed EASNMF method by the simulated data set 1 to investigate it precisely. The three parameters, including α, β and γ, need to be determined in advance. As mentioned earlier, the parameter α controls the endmember independence term, the parameter β adjusts the effect of abundance sparse constraint and the parameter γ is the penalty parameter for manifold regularization. Figure 5 shows the curves of these three parameters with respect to SAD and RMSE. From Figure 5a,b, it can be easily found that both of the SAD and RMSE curves are not sensitive to the parameter α. Besides, the curves of SAD and RMSE generally rise with the increasing of parameter β. Moreover, Figure 5b demonstrates that when the parameter γ is around 1, the values of SAD and RMSE are small. It corresponds to the good unmixing effect for endmember extraction and abundance estimation. Therefore, the parameters α, β and γ are separately set as 0.01, 0.001 and 1.

After determining the parameters of the proposed model, we perform our algorithm with simulated data set 1. Figure 6 shows the reference endmember curves with a red solid line and the estimated endmember curves obtained by EASNMF method with a blue dotted line. Through the observation and analysis of Figure 6, we can see that most estimated endmembers are very close to the reference, and there are some small differences between the references and estimations for endmembers 3 and 9. In general, the endmembers obtained by the proposed method are in good accordance with the referenced ones, which demonstrates the satisfactory endmember estimations provided by EANMF algorithm.

At the same time, we also exhibit the abundance maps of the proposed method and compare it with some related algorithms to illustrate its effectiveness for unmixing. The comparison algorithms include the methods of L_1/2-NMF, GLNMF, MVCNMF and CoNMF. The results of the EASNMF algorithm and the comparisons are displayed in Figure 7. Due to the limited space of the paper, here the abundance maps of only three representative endmembers are exhibited, and they are endmembers 2, 6, and 9, respectively. Comparing to the real abundance maps in Figure 2, it can be observed that the abundance map obtained by the EASNMF algorithm is smoother than that of L_1/2-NMF algorithm, especially in the homogeneous part for endmember 9. In addition, the result of GLNMF algorithm is satisfactory with some details missing, such as the texture of homogeneous region in Figure 7a. Since there is no constraint on abundance, the background of multiple abundance maps extracted by MVCNMF algorithm is messier than the other methods. Although the CoNMF algorithm can handle the steps of endmember number estimation, endmember extraction and abundance estimation together, it does not put forward some specific constraints, making the extracted abundance not ideal. In general, the performance of the proposed EASNMF algorithm is satisfied, which illustrates its effectiveness for HU.

Furthermore, in order to quantitatively evaluate the algorithms, we also calculate the values of SAD and RMSE by Equations (19) and (20). For the purpose of comparing, the SAD value and RMSE value of the EASNMF algorithm and comparisons are listed in

Table 2. The values of RMSE and SAD for EASNMF algorithm and the comparisons on simulated data set 1.

		L1/2-NMF	GLNMF	MVCNMF	CoNMF	EASNMF
RMSE	Average	0.0257	0.0264	0.0391	0.0973	0.0252
	Endmember 1	0.0173	0.0183	0.0273	0.0883	0.0168
	Endmember 2	0.0165	0.0174	0.0255	0.0839	0.0163
	Endmember 3	0.0335	0.0368	0.0527	0.1104	0.0341
	Endmember 4	0.0203	0.0248	0.0299	0.0945	0.0199
	Endmember 5	0.0202	0.0257	0.0259	0.0934	0.0204
	Endmember 6	0.0382	0.0357	0.0598	0.1124	0.0368
	Endmember 7	0.0129	0.0191	0.0301	0.0631	0.0137
	Endmember 8	0.0195	0.0210	0.0371	0.0837	0.0189
	Endmember 9	0.0527	0.0385	0.0638	0.1459	0.0500
SAD	Average	0.0218	0.0318	0.0444	0.2215	0.0188
	Endmember 1	0.0141	0.0346	0.0203	0.2233	0.0138
	Endmember 2	0.0083	0.0151	0.0200	0.1250	0.0076
	Endmember 3	0.0396	0.0319	0.0974	0.3549	0.0333
	Endmember 4	0.0060	0.0101	0.0096	0.0821	0.0054
	Endmember 5	0.0151	0.0190	0.0192	0.1102	0.0132
	Endmember 6	0.0540	0.1077	0.0713	1.4470	0.0411
	Endmember 7	0.0099	0.0074	0.0240	0.0661	0.0100
	Endmember 8	0.0075	0.0153	0.0235	0.0476	0.0075
	Endmember 9	0.0415	0.0454	0.1141	0.8180	0.0378

. Simultaneously, the SAD value and the RMSE value for each endmember are also recorded in Table 2. Based on the values in Table 2, the results of MVCNMF and CoNMF are worse than the other comparisons, whether for the endmembers or the abundances. In addition, it can be noted that compared to the listed comparisons, the best and second-best results are mostly obtained by EASNMF algorithm. Owing to the appropriate constraints based on the endmember independence and spatial weight, the performance of the proposed method is slightly higher than the listed comparison algorithms for unmixing.

4.6. Experiment on Simulated Data Set 2

In this section, we perform the proposed method on simulated data set 2, whose results will be shown and analyzed in detail. Similarly, it needs to analyze the parameters in model before the experiment. Figure 8 presents the curves of SAD and RMSE values with different values of three parameters α, β and γ. From Figure 8, the overall values of SAD and RMSE are relatively low, demonstrating the effectiveness of endmember extraction and abundance estimation. Besides, it is not difficult to find that the parameter α has a small effect on the values of SAD and RMSE in the local interval. When parameter β is small, the corresponding SAD and RMSE values are relatively small, indicating the good performance of unmixing. With the increase of parameter γ, the values of SAD and RMSE gradually decrease, and their curves tend to be stable around γ = 1. Thus, the parameters α, β and γ are separately set as 0.1, 0.01, 1.

In our paper, the endmember number is another important issue for HU, which is assumed to be known. Here, we can also use the reconstruction error defined by ${\Vert Y-\overline{E}\overline{A}\Vert }^2$ to analyze the endmember numbers in the experiment [34]. Figure 9 exhibits the curve of reconstruction error with respect to different endmember numbers. As the endmember number increases, the error decreases and tends to be stable. It can be seen from the curve in Figure 9, when the endmember number is 9, the error is the smallest.

In order to indicate clearly the difference between the experimental results of different methods and the references on simulated data set 2, Figure 10 displays the error maps of abundances obtained by EASNMF algorithm and comparisons. In the error maps, the closer the color is to blue, the smaller the error is. Due to the limitation of the space, we only present the error maps of abundances corresponding to three typical endmembers, i.e., endmembers 3, 5 and 9. It can be seen from Figure 10 that the result of CoNMF algorithm is the worst and there are some scattered points on the error map of L_1/2-NMF algorithm. Since the MVCNMF algorithm does not have any abundance constraints in its model for unmixing, the error distributes in the whole image without any spatial structure information. Although the error maps of the GLNMF algorithm and EASNMF algorithm are somewhat similar, the overall color of the error maps obtained by EASNMF method is darker than that of GLNMF algorithm, indicating a small error in general. In addition, owing to the smoothness constraint on abundance, the error distribution of the proposed algorithm is relatively smoother than the comparisons, especially for endmembers 3 and 9. Due to the complexity of the data set scene, the advantage of the proposed algorithm is not obvious which is the limitation for the proposed algorithm. In general, the proposed algorithm is effective for unmixing.

Furthermore,

Table 3. The values of RMSE and SAD for the proposed EASNMF method and comparisons on simulated data set 2.

		L1/2-NMF	GLNMF	MVCNMF	CoNMF	EASNMF
RMSE	Average	0.0820	0.0812	0.0863	0.1149	0.0783
	Endmember 1	0.1824	0.1596	0.2311	0.1359	0.1567
	Endmember 2	0.0410	0.0479	0.0414	0.1060	0.0442
	Endmember 3	0.0837	0.0755	0.0839	0.1266	0.0743
	Endmember 4	0.0785	0.0406	0.0517	0.1194	0.0458
	Endmember 5	0.0570	0.0544	0.0685	0.1418	0.0496
	Endmember 6	0.2066	0.1857	0.2109	0.1028	0.1852
	Endmember 7	0.0305	0.0459	0.0397	0.0907	0.0377
	Endmember 8	0.0514	0.0716	0.0656	0.0979	0.0630
	Endmember 9	0.0402	0.0498	0.0800	0.1129	0.0483
SAD	Average	0.0164	0.0195	0.0184	0.1274	0.0149
	Endmember 1	0.0416	0.0231	0.0466	0.0371	0.0255
	Endmember 2	0.0061	0.0089	0.0115	0.1087	0.0065
	Endmember 3	0.0088	0.0244	0.0118	0.6776	0.0051
	Endmember 4	0.0068	0.0092	0.0066	0.0747	0.0070
	Endmember 5	0.0199	0.0245	0.0251	1.5555	0.0161
	Endmember 6	0.3723	0.7477	0.3692	0.0464	0.6028
	Endmember 7	0.0070	0.0258	0.0064	0.1516	0.0156
	Endmember 8	0.0208	0.0101	0.0209	0.0593	0.0124
	Endmember 9	0.0047	0.0075	0.0120	0.0857	0.0052

records the SAD value and RMSE value of the EASNMF algorithm and comparisons, including the average value and the value corresponding to each endmember. Due to the fact that the spatial distribution of this data set is trivial, the advantage of spatial constraints in unmixing model is not obvious. That is the reason that for some endmembers, the performance of the L_1/2-NMF is better than the other algorithms. Through comparison and analysis, we can find that the algorithms of GLNMF and MVCNMF is similar and the performance of CoNMF algorithm is the worst. Since we utilize the characteristic of endmembers and abundances and exploit the latent structure of data, the proposed method makes full used of the spatial-spectral information in the image. On the whole, in terms of the average value of SAD and RMSE values, the EASNMF algorithm obtains the smallest values. Although the improvement is not obvious to average values, it still demonstrates its effectiveness for HU.

4.7. Experiment on Cuprite Data Set

We turn our attention to the real data set. The real data set adopted in the experiment is Cuprite data set, which is commonly used for HU [39]. Firstly, the parameters in the proposed method need to be determined. Since there are no real abundances for Cuprite data set, the method to determine the value of parameter is usually by the experience. Here the metric of SAD is integrated to compare the extracted endmembers with the spectra in spectral library to further adjust the parameters. The purpose is to get good endmembers that are considered the base for unmixing by researchers. Figure 11 shows the performance of EASNMF algorithm for endmember extraction on Cuprite data set with the parameters α, β and γ. As shown in Figure 11a, the curve increases with the increase of parameters α and β. The difference between the minimum and maximum of the curve is small in the local interval of parameters value. From Figure 11b, it can be seen that the curve increases with the increase of parameter γ. In addition, from the parameter analysis of simulated data sets, we can see that the values of parameters α and β are proportional, which corresponds to 10. Therefore, based on the above analysis, the parameters are finally set as 0.1, 0.01, 0.1.

For the issue of endmember number supposed to be known, we can also use the reconstruction error to analyze the effect of endmember number. Figure 12 plots the curve of reconstruction error with different endmember number. Overall, the curve first drops, then stabilizes and finally rises. When the number of the endmember is 12, the smallest value of reconstruction error is achieved. In addition, according to the analysis in many articles [39,41], the estimated number of endmembers in Cuprite image is 12 for unmixing due to the tiny differences between some spectra of the same mineral.

The comparison between the USGS library spectra (red solid line) and the endmember signatures (blue dotted line) obtained by EASNMF algorithm on Cuprite data set are displayed in Figure 13. It can be found that most endmember signatures are similar to the spectra in spectral library. Moreover, for quantitative analysis,

Table 4. The SAD values for EASNMF algorithm and the comparisons on Cuprite data set.

		L1/2-NMF	GLNMF	MVCNMF	CoNMF	EASNMF
SAD	Average	0.0772	0.0782	0.0804	0.1428	0.0769
	Alunite	0.1137	0.1190	0.1140	0.4430	0.1136
	Andradite	0.0700	0.0709	0.0708	0.1510	0.0697
	Buddingtonite	0.0743	0.0731	0.0771	0.6229	0.0700
	Dumortierite	0.0848	0.0840	0.0866	0.2227	0.0825
	Kaolinite1	0.0984	0.1005	0.1039	0.2934	0.1002
	Kaolinite2	0.0742	0.0685	0.0746	0.4582	0.0748
	Muscovite	0.0892	0.0856	0.0897	0.3318	0.0878
	Montmorillonite	0.0594	0.0607	0.0643	0.1357	0.0607
	Nontronite	0.0710	0.0746	0.0778	0.2425	0.0739
	Pyrope	0.0596	0.0644	0.0602	0.1416	0.0588
	Sphene	0.0571	0.0674	0.0621	1.4085	0.0584
	Chalcedony	0.0866	0.0810	0.0878	0.0830	0.0883

lists the SAD value between each endmember and its corresponding spectral in spectral library. From the table it can again be found that the EASNMF method achieves the low average SAD value. However, the advantage of the proposed method is not so obvious over the comparisons. On the one hand, due to the fragmentary abundance maps of some endmember signatures, it leads the influence of the manifold and smoothness constraints in the model weakly. On the other hand, the parameter value of the proposed algorithm is not optimal, which takes into account the parameter analysis of simulated data. The grayscale abundance maps obtained by the proposed EASNMF algorithm are exhibited in Figure 14. Based on the analyses mentioned, it can be concluded that the proposed algorithm is effective for unmixing.

5. Conclusions

In this paper, we present a sparse NMF algorithm based on endmember independence and spatial weighted abundance for hyperspectral image unmixing. The proposed method not only considers the characteristic of endmembers and their abundances at the same time, but also makes full use of the spatial-spectral information in the image. First, we add the endmember independence constraint to the NMF model based on the assumption that the extracted endmembers should be independent from each other. Then, a weight matrix is constructed by the neighborhood pixel for abundance to make it smooth. In addition, inspired by the manifold learning, we construct the connection weight between two pixels by the correlation coefficient to further explore the structure of HSI data. The experiment results on three data sets including the simulated data set and the real data set demonstrate the effect of the proposed EASNMF algorithm.