Hyperspectral Unmixing via Low-Rank Representation with Space Consistency Constraint and Spectral Library Pruning

Spectral unmixing is a popular technique for hyperspectral data interpretation. It focuses on estimating the abundance of pure spectral signature (called as endmembers) in each observed image signature. However, the identification of the endmembers in the original hyperspectral data becomes a challenge due to the lack of pure pixels in the scenes and the difficulty in estimating the number of endmembers in a given scene. To deal with these problems, the sparsity-based unmixing algorithms, which regard a large standard spectral library as endmembers, have recently been proposed. However, the high mutual coherence of spectral libraries always affects the performance of sparse unmixing. In addition, the hyperspectral image has the special characteristics of space. In this paper, a new unmixing algorithm via low-rank representation (LRR) based on space consistency constraint and spectral library pruning is proposed. The algorithm includes the spatial information on the LRR model by means of the spatial consistency regularizer which is based on the assumption that: it is very likely that two neighbouring pixels have similar fractional abundances for the same endmembers. The pruning strategy is based on the assumption that, if the abundance map of one material does not contain any large values, it is not a real endmember and will be removed from the spectral library. The algorithm not only can better capture the spatial structure of data but also can identify a subset of the spectral library. Thus, the algorithm can achieve a better unmixing result and improve the spectral unmixing accuracy significantly. Experimental results on both simulated and real hyperspectral datasets demonstrate the effectiveness of the proposed algorithm.


Introduction
Hyperspectral imaging (HSI) has gained more attention in the past two decades.Spectral unmixing is an important task for remotely sensed hyperspectral data exploitation [1].HSI collects and processes electromagnetic spectrum information from the ground via many narrow bands so that it has a high spectral resolution.The characteristics of HSI make it possible to identify ground objects based on their unique spectral signature.Thus, HSI has been successfully applied to various fields such as target detection, mineralogy, agriculture and environment monitoring [2].However, with the low spatial resolution of HSI, multiple pure materials will jointly occupy a single pixel and become difficult to distinguish.Thus, as an important technique for hyperspectral data exploitation, the spectral unmixing problem [2,3] has recently attracted more attention from researchers.Usually, the spectral unmixing includes two major tasks: identifying the pure materials (endmembers), and estimating their corresponding proportions (abundances) presented in the mixed pixel.
There are two basic models to analyse the mixed pixel: the linear mixture model (LMM) [4] and the nonlinear mixture model [5][6][7].Compared with the nonlinear mixture model, the LMM has been widely applied for abundance estimation due to its computational tractability and flexibility.This model assumes that the observed (usually mixed) spectral signal is a linear combination of all the pure spectral signatures presented in that pixel.Ahmed et al. [8] proposed a method by using the artificial neural network to switch the linear and nonlinear for unmixing.
The traditional spectral unmixing algorithms with LMM focus on the study of the endmember extraction [9,10] and the abundance estimation [11].A variety of endmember identification algorithms based on geometrical and statistical approaches have been exploited [3].These algorithms assume that there are at least one pure pixel existing for each endmember in the scene, such as the pixel purity index [12], N-FINDR [13] and vertex component analysis (VCA) [14].Some researchers also have focused on the endmember generation algorithms, such as iterative error analysis [15] and iterative constrained endmembers (ICE) [16], which are based on the assumption that there are no pure pixels existing.However, these algorithms are very likely to fail in highly mixed or noisy hyperspectral data [17].Recently, Li et al. developed a new algorithm, called robust collaborative nonnegative matrix factorization (R-CoNMF), that can perform the whole steps of the hyperspectral unmixing chain [18].
In many cases, the identification of the endmembers in the original hyperspectral data becomes a challenge due to the lack of pure pixels in the scenes and the difficulty in estimating the number of endmembers in a given scene.To deal with these problems, the sparsity-based unmixing algorithms [17], which regard a large standard spectral library as endmembers, have recently been proposed.Compared with the size of the spectral library, the number of endmembers in a given scene is usually much smaller so that the abundances of the spectral library will be sparse.In practice, this is a combinatorial issue that requires efficient sparse regression techniques based on sparsity-inducing regularizers.The sparse unmixing algorithm via variable splitting and augmented Lagrangian (SUnSAL) [19] as one of the representative algorithms was proposed early.
However, it is not always true that the sparsity-based unmixing algorithms can obtain better spectral unmixing results than traditional unmixing techniques.As we know, in a spectral section, two different objects may have the same spectral curves, or the same object in different states shows different spectral characteristics which will lead to many similar spectral curves in the library.That is due to the high mutual coherence [20][21][22] of the library signatures which limits the success of these algorithms.To mitigate this drawback, some types of structured sparsity in the unmixing solutions via suitable regularization terms have been proposed recently [23][24][25][26].For instance, the sparse unmixing via variable splitting augmented Lagrangian and total variation (SUnSAL-TV) [23] was developed to utilize both spectral information and spatial information.SUnSAL-TV accounts for the spectral homogeneity of each pixel and its neighbours in the sparse unmixing formulation by means of a TV regularizer.Iordache et al. proposed the collaborative sparse unmixing via variable splitting and augmented Lagrangian (CLSUnSAL) algorithm [26].CLSUnSAL exploits the fact that a hyperspectral image always contains a small number of endmembers so that the fractional abundances matrix of the spectral library signatures contains only a few lines with nonzero entries.Therefore, both SUnSAL-TV and CLSUnSAL have performed better than the traditional sparse-based unmixing methods.Li et al. [27] proposed a new method for sparse unmixing of hyperspectral data with noise level estimation.Moreover, learning structured LRR has also been introduced [28], where a discriminative low-rank dictionary learning algorithm for sparse representation was proposed.The algorithm separates the sparse noise from the signals while simultaneously optimizing the dictionary atoms to reconstruct the denoised signals.Williams et al. [29] proposed a method of validation of abundance map reference data for spectral unmixing.Recently, the method of tensor is a new proposed method.Qian et al. [30] adopted the method called matrix-vector nonnegative tensor factorization for blind unmixing of hyperspectral imagery.Fan et al. [31] proposed a method of hyperspectral image restoration using low-rank tensor recovery.
However, it is often the fact that the size of endmember dictionary is much small than the spectral library.In other words, the endmember dictionary is determined by VCA algorithm instead of using spectral library.Recently, to avoid identifying the endmembers in the original hyperspectral data, some researchers use the spectral library as the endmembers and utilize some models based on LRR to estimate the abundance at the same time [32,33].However, the high mutual coherence of the library signatures makes this model perform badly.It does not distinguish well between the true endmembers and the others.To mitigate the drawback, in this paper, we propose a new hyperspectral unmixing method based on LRR model with dictionary pruning strategy.
In addition, most of the traditional models which are based on LRR only take spectral information of hyperspectral data into consideration.The hyperspectral image has special characteristics of space and the information of the space is very important for hyperspectral unmixing.The use of spatial information is concerned with the relationship between each pixel and its neighbours, considering the spatial consistency of the image.Although using spectral information unmixing technology can already obtain very good unmixing results, more previous algorithms using spatial information, such as pre-treatment [34] or unmixing algorithms [35,36], all show the positive impact on the improvement of the unmixing accuracy.
To mitigate the drawback, a method which not only considers the global correlation of the data via LRR with space consistency constraint, but also avoids the accurate estimation of the endmembers by pruning the spectral library is proposed in this paper.The proposed algorithm can be described as an iterative process.First, we use the model of the LRR with the space consistency constraint to estimate the abundance while the whole spectral library acting as the endmembers dictionary.Because the number of the real endmembers is much smaller than the whole spectral library, there exists several endmembers whose abundance values are very low.Secondly, we prune the spectral library based on the abundance matrix obtained in the previous step.The pruning strategy follows such an assumption that if the abundance map of one material does not contain any large values, it is not a real endmember and will be removed from the spectral library.After this step, we can get a pruned spectral library which is a subset of the original one.Then, the pruned spectral library is the new endmembers dictionary and the same process in the first step is performed again, and so on.When meeting the stop conditions, the algorithm terminates.Based on the dictionary pruning strategy, the proposed algorithm can use the spectral library as endmembers and avoid extracting them from the hyperspectral image, compared with the algorithm in [36].In addition, the algorithm can also better capture the spatial structure of data than the sparsity-based unmixing algorithms.The whole process of this method is summarized in Figure 1.
The main contributions of our method can be summarized as following.(1) The important characteristic of HSI, space consistency, is integrated in the low-rank representation model for hyperspectral unmixing, which will make the abundance estimation more accurate.(2) A strategy of spectral library pruning is designed by using the coefficients of the low-rank representation, which not only make the result more stable, but also avoids the extraction of endmembers.In addition, we will demonstrate the two advantages in the experiments.
The remainder of this paper is organized as follows.In Section 2, we review LMM and introduce the approach.Section 3 describes the proposed methodology.Section 4 analyses the performance of the proposed approach with simulated data and real hyperspectral data.Conclusions are drawn in Section 5.

Linear Spectral Unmixing
The widely used LMM assumes that the spectral response of a pixel is a linear combination of all the pure spectral signatures (endmembers) present in the pixel.In this case, the model can be described as: where y is an L × 1 column vector (the mixed pixel), L is the number of bands, A is an L × m matrix containing m spectral signatures (endmembers), x is a m × 1 vector containing the fractional abundances of the endmembers, and n is an L × 1 vector denoting the noise and model error.
Considering the ground truth, there are two constraints imposed on the LMM: the abundance non-negativity constraint (ANC) and the abundance sum-to-one constraint (ASC), respectively: x 1 be the observed image and N is the total number of mixed pixels in the image.Let R m N × ∈ X be the abundances matrix.Then, Equation ( 1) can be rewritten as:

Linear Spectral Unmixing
The widely used LMM assumes that the spectral response of a pixel is a linear combination of all the pure spectral signatures (endmembers) present in the pixel.In this case, the model can be described as: where y is an L × 1 column vector (the mixed pixel), L is the number of bands, A is an L × m matrix containing m spectral signatures (endmembers), x is a m × 1 vector containing the fractional abundances of the endmembers, and n is an L × 1 vector denoting the noise and model error.
Considering the ground truth, there are two constraints imposed on the LMM: the abundance non-negativity constraint (ANC) and the abundance sum-to-one constraint (ASC), respectively: ASC : Let Y ∈ R L×N be the observed image and N is the total number of mixed pixels in the image.Let X ∈ R m×N be the abundances matrix.Then, Equation (1) can be rewritten as: where E ∈ R L×N denoting the noise and model error.
In this paper, we will focus on the LMM, due to its computational tractability and flexibility in different applications.

Abundance Estimation via LRR
Abundance estimation via LRR is based on the following theorem [37].
) and rank(A) = m, then we have: In [38], the dictionary A is extracted from the hyperspectral image itself.Matrix A usually satisfies the full column rank property as the spectra of the extracted pure endmembers are distinct from each other and the number of bands L is much larger than the number of endmembers m.Normally, the columns of Y are highly correlated and it means that the matrix Y is a low-rank matrix.Thus, we can infer that the corresponding representation matrix X is also low rank.
To use this property, researchers employ the LRR model [39,40] to spectral unmixing problem.Starting with the simple LMM, this method can be described as the following optimization problem: where Y ∈ R L×N is the observed data matrix, A ∈ R L×m is the endmembers matrix and X ∈ R m×N is the fractional abundances matrix.X* is the lowest rank solution and the terms X ≥ 0, 1 T X = 1 T represent the ANC and ASC constraints.
Because the rank of the matrix has a discrete type, Equation ( 6) is not easy to solve.Thus, we can use the nuclear norm of the matrix to replace the rank of the matrix as follows: min where X ∈ R m×N is the lowest rank representation of the observed signal Y ∈ R L×N , A ∈ R L×m is a linear data space of the dictionary, E ∈ R L×N is the matrix of noise error, the parameter λ is coordination value, • * is the nuclear norm and • 2,1 represents the mixed l 2,1 norm of the matrix.
N is the number of signals, L is the signal dimension, m is the number of dictionary atoms.

Low-Rank Representation of Coefficient Constraints
Low-rank representation is a compressed sensing technology that can effectively restore the original real data, even if the data exist in several sub-spaces.Low-rank representation not only can make full use of the global structure information of the data, but is also robust to noise which has received wide attention since it was put forward.Low-rank representation seeks the lowest rank representation of the observed signal in relative to a suitable dictionary.The basic LRR model with noise is shown in Equation (7).
Under normal circumstances, when the dictionary is certain, the similar observation signal has a similar low-rank coefficient.Then, using this property, the regular term is added to the objective function of the basic LRR model.It restricts that similar observation signals have a similar low-rank coefficient.Then, Equation ( 7) can be modified as follows: min where I ∈ R n×n is the unit matrix, and • F is the Frobenius norm of the matrix.The size of matrix B is n × n, and it can be calculated by the method that, if y j is similar to y i , then the matrix B equals to ω ij , or equals to zero, where the ω ij is the similarity of the observation vector y i and y j , according to the different situations having different construction methods, generally in the range between 0 and 1, more commonly using a heat kernel function and Gauss kernel function.
The regular term added in Equation ( 8) is to make the similar observation vectors have similar LRR coefficients, and the coordination parameter β controls the proportion of the regularization in the whole optimization objective function.The optimization problem of the formula can not only make full use of the high correlation among the different bands of the observed data, but also protect the spatial local structure of the data.

Space Consistency Constraint
In the previous section, we introduce the method of applying LRR technique to the unmixing problem, where the VCA is used to extract the endmember dictionary A from the hyperspectral image itself.Qu et al. [38] proposed this concept, and used the LRR to solve the problem of unmixing for the first time.In addition, they put forward an improved bilinear model, which applied the LRR to obtain the result of unmixing.
In this paper, we propose an algorithm based on LRR with spatial consistency constraint for abundance estimation.This algorithm not only makes use of the spectral information of the data, but also considers the spatial information of the data, that is, the relationship between the pixel and its neighbours.It shows the core idea of the algorithm, that is, how we use spatial information in the process of unmixing in Figure 2. is n n × , and it can be calculated by the method that, if j y is similar to i y , then the matrix B equals to ij ω , or equals to zero, where the ij ω is the similarity of the observation vector i y and j y , according to the different situations having different construction methods, generally in the range between 0 and 1, more commonly using a heat kernel function and Gauss kernel function.

Remote
The regular term added in Equation ( 8) is to make the similar observation vectors have similar LRR coefficients, and the coordination parameter β controls the proportion of the regularization in the whole optimization objective function.The optimization problem of the formula can not only make full use of the high correlation among the different bands of the observed data, but also protect the spatial local structure of the data.

Space Consistency Constraint
In the previous section, we introduce the method of applying LRR technique to the unmixing problem, where the VCA is used to extract the endmember dictionary A from the hyperspectral image itself.Qu et al. [38] proposed this concept, and used the LRR to solve the problem of unmixing for the first time.In addition, they put forward an improved bilinear model, which applied the LRR to obtain the result of unmixing.
In this paper, we propose an algorithm based on LRR with spatial consistency constraint for abundance estimation.This algorithm not only makes use of the spectral information of the data, but also considers the spatial information of the data, that is, the relationship between the pixel and its neighbours.It shows the core idea of the algorithm, that is, how we use spatial information in the process of unmixing in Figure 2. As can be seen from Figure 2, if the current pixel is the i-th observation signal, that is, i y , selecting its n neighbourhoods as its nearest neighbour (n = 3).Typically, the current pixel is similar to its spatial neighbour, but this requires a prerequisite, that is, the uniform region.If the current pixel is at the border, then the pixels in the red box of Figure 2 may belong to two or more substances.Therefore, in addition to the spatial neighbour as the prior information, we also need to measure the distance between the current pixel and its nearest neighbour to distinguish the pixels which are adjacent to each other but belong to different kinds.As shown in Figure 2, the 3 3 × neighbourhood of the current pixel has a total of 8 pixels, the pixel i is at the boundary, and only the j-th, k-th and m-th endmembers are the same kind of the pixel i.At this point, it is needed to judge by As can be seen from Figure 2, if the current pixel is the i-th observation signal, that is, y i , selecting its n neighbourhoods as its nearest neighbour (n = 3).Typically, the current pixel is similar to its spatial neighbour, but this requires a prerequisite, that is, the uniform region.If the current pixel is at the border, then the pixels in the red box of Figure 2 may belong to two or more substances.Therefore, in addition to the spatial neighbour as the prior information, we also need to measure the distance between the current pixel and its nearest neighbour to distinguish the pixels which are adjacent to each other but belong to different kinds.As shown in Figure 2, the 3 × 3 neighbourhood of the current pixel has a total of 8 pixels, the pixel i is at the boundary, and only the j-th, k-th and m-th endmembers are the same kind of the pixel i.At this point, it is needed to judge by the distance measure to avoid the condition that the boundary is misidentified as the homogeneous region.The detail of measurement process is: first, calculate the Euclidean distance between the current pixel and its nearest neighbour, and then, select the nearest p pixels as the real neighbour of the current pixel.
In summary, the so-called spatial consistency constraint reflects the following facts: (1) The pixels of the near neighbourhood in the image have similar endmembers and their corresponding abundance.(2) The pixels of the near spectral distance have similar endmembers and their corresponding abundance.
(3) At the same time, combined with the above two points, selecting p pixels with the nearest distance to the current pixels in the n neighbourhood as its real neighbourhood, constrain its abundances to be similar.
The above thoughts are shown to be of great advantage through experiments.If we only consider the first point, we cannot avoid the boundary problem which will cause the error increase.Because of the presence of noise and the interference of the atmosphere, even though the spectral information of the pixels is very similar, but they are also likely to belong to the different materials.Thus, it is not appropriate to only consider the second point.The combination of the first point and the second point at the same time is the key for the algorithm to improve the accuracy of unmixing.
Then, we propose an optimization model based on the LRR with the spatial consistency constraint for the estimation of the abundance as follows: min where Y ∈ R L×N is the hyperspectral data, X ∈ R L×m is the abundance matrix, E ∈ R L×N is the noise matrix, and ∈ R N×N is the unit matrix.The λ and β are coordinate parameters, m is the number of the endmember and p refers to the p pixels of the nearest distance to the current pixel in n × n neighbourhood of the current pixel as its real neighbour.The size of matrix B is n × n, and it can be calculated by the method that if y j is similar to y i , then the matrix B equals to ω ij , or equals to zero, where ω ij is the similarity of the observation vector y i and y j .We set ω ij to be 1.
The optimization problem shown in Equation ( 9) can be solved by the method of Augmented Lagrange Multiplier (ALM) [41].The ANC is added in each step of the iterative process, while the ASC is embodied in the normalization process of the last step of the iteration.Because the ANC has obvious physical meaning, that is, the value of abundance cannot be negative, so the general situation is to add it.The ASC need to choose to add or not based on the actual situation.
The whole solving process of the optimization problem is as follows.Firstly, we remove the ANC and the ASC of Equation ( 9), and translate it into the following equivalent problem: Then, Equation ( 10) can be converted to the ALM problem as follows: min where Y 1 , Y 2 and Y 3 are the Lagrangian multipliers, and µ is the penalty parameter.For efficiency, we choose the inexact ALM, which we outline in Algorithm 1.

Algorithm 1: Solve the SCC-LRR problem by ALM
Input: library A, data matrix Y, near neighbourhood consistency constraint matrix H, and regularization parameter λ, β.

Spectral Library Pruning
After solving optimization (Equation ( 9)), we can obtain the fractional abundances matrix X corresponding to the spectral library A. Compared with the size of the spectral library, the number of endmembers in a given scene is usually much smaller.Thus, sparsity exists among the lines of the obtained abundances matrix X and each line of X denotes the abundance map corresponds to one endmember.To utilize the sparsity, we propose a dictionary pruning strategy to prune the abundances matrix X and its corresponding spectral signatures in A.
This strategy can be described as two steps: (1) Compute the number of the pixels whose abundance value corresponding to one endmember is smaller than a preset threshold denoted by ε. (2) If the number is equal to the total number of pixels in the scene, we will get rid of the spectral signature from the endmember matrix.
It means that if one endmember is rarely distributed over the scene, it is not a real endmember.In general, if one spectral signature is a real endmember, its abundance values corresponding to some areas will be high in some areas in the scene.That is our motivation of using the strategy to prune the spectral library A. In the experiments on simulated data, we set ε to be 0.02 and then if one endmember is distributed under 2% over the scene, it will be removed from the library.The value of ε should be set to a small value to avoid removing the true endmembers from the library.Thus, we set ε to 0.01 in the experiments of real data because real data are more complex than simulated data.In each iteration, the size of spectral library will decrease and the abundances of the retained endmembers will increase.Thus, in each iteration, the value of ε will be updated by ε × iterations.
Then, we can get a pruned spectral library A as the endmember matrix in which the number of spectral signatures is smaller than the original one.By using this pruned spectral library A, we solve the optimization (Equation ( 9)) again and obtain a new fractional abundances matrix Y corresponding to the pruned spectral library A. According to the dictionary pruning strategy, we can prune the spectral library A repeatedly until the number of spectral signatures retained meets the stopping condition.
A threshold denoted by T acting as the stop condition is defined to control this iterative process.The stop condition can be described like that when the difference between the number of spectral signatures retained and the estimated number of endmembers in a scene is less than T, the algorithm stops.Thus, a higher value of T means that it retains more endmembers from the spectral library.In this paper, we use the hyperspectral subspace identification by minimum error (HySime) algorithm [42] to estimate the number of endmembers.In the experiments with simulated data, the range of values for T is 1 to 10 and in the experiments with real data, we set T to be 25.We retain more endmembers for real data because they are more complex than simulated data.When the situation is complex (e.g., low SNR, big endmember number or a very large library), we set ε to be a small value and T be a large value to avoid removing the true endmembers from the library.After several iterations, we can get a good result of the optimization problem.The whole process of this method has been summarized in Figure 1.

Experiments with Simulated Data and Real Data
To validate the two main advantages of our method, the space consistency and the strategy of spectral pruning, we first compare the method of LRR with our proposed method to validate the effectiveness of the space consistency constraint for abundance estimation, and we also show the spectral library pruned by our method to demonstrate the effectiveness of the spectral library pruning strategy.
In Section 4.1, we use three simulated datasets to analyse the performance of the proposed approach.Firstly, we introduce the three datasets, and then we analysis the performances of them separately including the parameter settings, error analysis, etc.In Section 4.2, we use the real dataset and analyse its performance with other traditional methods.Finally, we analyse the parameters of the experiments in Section 4.3.The results of the experiments show the good performance of the proposed method.

Simulated Datasets
The spectral library we used in our simulated image experiments is generated from a random selection of 240 materials from the USGS library, denoted splib061 and released in September 2007.It comprises spectral signatures with reflectance values given in 224 spectral bands and distributed uniformly in the interval 0.4-2.5 µm.Our library, denoted by A1, is used to generate three different simulated hyperspectral data cubes.
(1) Simulated Data Cube 1 (DC1): This simulated data cube is generated following the methodology of [26], using five randomly selected spectral signatures from A1. DC1 has 75 × 75 pixels and each simulated pixel was generated using a LMM, with the five endmembers and imposing the ASC in it.In the resulting simulated image, shown in Figure 3a, there are pure regions as well as mixed regions constructed using mixtures ranging between two and five endmembers, distributed spatially in the form of distinct square regions.Figure 3b-f, respectively, shows the true fractional abundances for each of the five endmembers.The background pixels are made up of mixtures of the same five endmembers, but their respective fractional abundances values were randomly fixed as 0.1149, 0.0741, 0.2003, 0.2055 and 0.4051, respectively.The obtained data cube was then contaminated with white noise, having different levels of the signal-to-noise ratio (SNR): 20, 30 and 40 dB.(2) Simulated Data Cube 2 (DC2): Using the library A1, we generated a data cube of 48 × 48 pixels and it contains six endmembers.The endmembers were randomly selected from library A1.In each simulated pixel, the fractional abundances of the endmembers follow a Dirichlet distribution [14].As DC1, the scene was again contaminated with white noise using the same SNR value adopted for DC1.(3) Simulated Data Cube 3 (DC3): Using the library A1, we generate various data cubes of 75 × 75 pixels, each containing five endmembers.The simulated method is similar to the first simulated data.There are pure regions as well as mixed regions constructed using mixtures ranging between two and five endmembers, distributed spatially in the form of distinct square regions.The background pixels are made up of mixtures of the same five endmembers, but their respective fractional abundances values were randomly fixed as 0.1149, 0.0741, 0.2003, 0.2055 and 0.4051, respectively.The obtained data cube was then contaminated with white noise pixels, each containing five endmembers.The simulated method is similar to the first simulated data.There are pure regions as well as mixed regions constructed using mixtures ranging between two and five endmembers, distributed spatially in the form of distinct square regions.
The background pixels are made up of mixtures of the same five endmembers, but their respective fractional abundances values were randomly fixed as 0.1149, 0.0741, 0.2003, 0.2055 and 0.4051, respectively.The obtained data cube was then contaminated with white noise In this work, the proposed algorithm is tested on three groups of simulated datasets and a real hyperspectral image.To verify the superiority of the proposed algorithm, the performance of it is compared with both the traditional unmixing algorithms (nonnegative constrained least squares (NCLS), full constrained least squares (FCLS)) [11] and the sparse unmixing algorithms (SUnSAL [17], SUnSAL-TV [23] and CLSUnSAL [26]).The parameter threshold ε plays an important role in the proposed algorithm and we set it to 0.02 in our simulated image experiments.
The quality metric adopted in our experiments to assess the unmixing results is the signal to reconstruction error (SRE) [17], which can be defined as follows: where X is the true fractional abundance vector and x is the estimated fractional abundance vector.
( ) E ⋅ stands for mean value.The higher the SRE is, the better the quality of the unmixing is.As it can give more information regarding the power of the signal in relation with the power of the error, we use this measure instead of the classical root mean squared error (RMSE) [34].We also compute the abundance angle distance (AAD) as another performance discriminator adopted in this paper.AAD can express as follows: In this work, the proposed algorithm is tested on three groups of simulated datasets and a real hyperspectral image.To verify the superiority of the proposed algorithm, the performance of it is compared with both the traditional unmixing algorithms (nonnegative constrained least squares (NCLS), full constrained least squares (FCLS)) [11] and the sparse unmixing algorithms (SUnSAL [17], SUnSAL-TV [23] and CLSUnSAL [26]).The parameter threshold ε plays an important role in the proposed algorithm and we set it to 0.02 in our simulated image experiments.
The quality metric adopted in our experiments to assess the unmixing results is the signal to reconstruction error (SRE) [17], which can be defined as follows: ), where X is the true fractional abundance vector and is the estimated fractional abundance vector.E(•) stands for mean value.The higher the SRE is, the better the quality of the unmixing is.As it can give more information regarding the power of the signal in relation with the power of the error, we use this measure instead of the classical root mean squared error (RMSE) [34].We also compute the abundance angle distance (AAD) as another performance discriminator adopted in this paper.AAD can express as follows: where X i is the true fractional abundance vector of the i-th endmember and xi is the estimated fractional abundance vector of the i-th endmember.m stands for the number of the endmembers.Figure 4 shows the abundance maps obtained by different unmixing methods for endmember #5 in DC1.Throughout the estimated abundance results, the abundance maps obtained by NCLS and SUnSAL are full of noise points and it is hard to make out the endmember signature from the mixed spectral.The reason for this is that the spectral library A1 is generated from a random selection of 240 materials from the USGS library.The spectral angle distance is quite small, which makes it difficult to separate the endmembers from noise.Since DC1 has distinct spatial information, SUnSAL-TV, which considers spatial information, achieves a better result than NCLS and SUnSAL.The abundance maps obtained by SUnSAL-TV have smoother spatial changes in the homogeneous regions, and contain fewer noise points.However, in some cases, the abundance images obtained by SUnSAL-TV may exhibit an over-smooth visual effect in some regions.
, where Xi is the true fractional abundance vector of the i-th endmember and ˆi x is the estimated fractional abundance vector of the i-th endmember.m stands for the number of the endmembers.
Figure 4 shows the abundance maps obtained by different unmixing methods for endmember #5 in DC1.Throughout the estimated abundance results, the abundance maps obtained by NCLS and SUnSAL are full of noise points and it is hard to make out the endmember signature from the mixed spectral.The reason for this is that the spectral library A1 is generated from a random selection of 240 materials from the USGS library.The spectral angle distance is quite small, which makes it difficult to separate the endmembers from noise.Since DC1 has distinct spatial information, SUnSAL-TV, which considers spatial information, achieves a better result than NCLS and SUnSAL.The abundance maps obtained by SUnSAL-TV have smoother spatial changes in the homogeneous regions, and contain fewer noise points.However, in some cases, the abundance images obtained by SUnSAL-TV may exhibit an over-smooth visual effect in some regions.In Table 1, we can also find that LRR model performs about similar to SUnSAL in the unmixing problem.LRR model can capture global structures of the abundance matrix and is also robust to noise.However, due to the large library and its high mutual coherences, LRR does not perform much better than SUnSAL.In addition, the method based on the LRR with the space consistency constraint (SCC-LRR) also achieves a better result than the other methods besides the proposed method.The reason is that there are some homogeneous regions in DC1 and the method of the SCC-LRR takes the space consistency into consideration.Due to consideration of the space consistency constraint and the dictionary pruning strategy, the abundance maps obtained by the proposed algorithm are closer to the true abundance maps than the ones estimated by SUnSAL-TV.As shown in Table 1, the accuracy of the results is clearly improved when the proposed algorithm is used.Therefore, the proposed algorithm significantly outperforms the other methods.According to the comparative analysis of the experiment, the selections of n = 3 and p = 3 are more appropriate.When n = 5, the experiment shows that the effect is not improved and there are more calculations.In addition, it will lead to the advantage of space consistency not not being easily In Table 1, we can also find that LRR model performs about similar to SUnSAL in the unmixing problem.LRR model can capture global structures of the abundance matrix and is also robust to noise.However, due to the large library and its high mutual coherences, LRR does not perform much better than SUnSAL.In addition, the method based on the LRR with the space consistency constraint (SCC-LRR) also achieves a better result than the other methods besides the proposed method.The reason is that there are some homogeneous regions in DC1 and the method of the SCC-LRR takes the space consistency into consideration.Due to consideration of the space consistency constraint and the dictionary pruning strategy, the abundance maps obtained by the proposed algorithm are closer to the true abundance maps than the ones estimated by SUnSAL-TV.As shown in Table 1, the accuracy of the results is clearly improved when the proposed algorithm is used.Therefore, the proposed algorithm significantly outperforms the other methods.
According to the comparative analysis of the experiment, the selections of n = 3 and p = 3 are more appropriate.When n = 5, the experiment shows that the effect is not improved and there are more calculations.In addition, it will lead to the advantage of space consistency not not being easily reflected when the value of p is too small.On the contrary, if the value of p is too large, it will cause a greater error.In addition, the optimal setting and running times are shown in Table 2.
To validate the proposed algorithm extensively, we use DC2 to demonstrate the effectiveness.Figure 5 shows the ground-truth and the estimated abundances obtained by different unmixing methods in the A1 library.The datasets contain six true endmembers and the positions of the true endmembers in the original library are 5, 6, 26, 34, 67 and 177, as shown in Figure 5a.In Figure 5, it can be seen graphically that the lines (denoting the abundance of a certain endmember in all pixels of the image) estimated by CLSUnSAL are more sparse than the ones estimated by SUnSAL.After applying SUnSAL, as shown in Figure 5b, there are many low abundance values estimated for endmembers which are not actually present in the image.Due to the high mutual coherence of the library signatures, SUnSAL does not perform well.Unlike SUnSAL, which employs pixelwise independent regression, CLSUnSAL enforces joint sparsity among all the pixels.Thus, it improves significantly the accuracy of the unmixing solutions over those provided by SUnSAL.Although the abundances obtain by CLSUnSAL are similar to those in the ground-truth, there are also several lines that are not true abundances in Figure 5c. Figure 5d shows the estimated abundances obtained by the proposed algorithm.After applying this method, the remained subset of the spectral library contains six endmembers.However, almost every line of the image is much brighter than the one in the true abundance map.This is because that the image in DC2 contains very few homogeneous regions, which leads to the failure of spatial consistency constraints.However, it improves the accuracy in the whole, which is shown in Table 3.Using DC2, we also demonstrate the effectiveness of the space consistency and the dictionary pruning method separately.In Table 3, we can find that, if using the whole library A1 as the endmembers matrix, CLSUnSAL can obtain the highest SRE (dB) values.This is because CLSUnSAL constrains the pixels to share the same set of endmembers.AAD is defined as the angle between estimated abundances and the true abundances.The AAD values achieved by LRR are lower than  Using DC2, we also demonstrate the effectiveness of the space consistency and the dictionary pruning method separately.In Table 3, we can find that, if using the whole library A1 as the Then, we also test the performance of the proposed method on the DC3, as shown in Figure 6. Figure 6 shows the simulated data and the true fractional abundances for each of the five endmembers.There are both homogeneous region and heterogeneous region, as shown in Figure 6a. Figure 6b-f shows the true abundance maps of the five endmembers.Figure 7a-e shows the estimated abundance maps obtained by the proposed method, which contains five endmembers.Considering the sparse consistency constraint, the pixel of the heterogeneous region is not similar to its neighbour pixels which will makes the effect of the unmixing become worse, as shown in Figure 7. On the contrary, in the homogeneous region, the effects of the unmixing can be obtained.In Figure 7, the proposed method can obtain better results overall.We also compare the results with the other methods, as shown in Table 5.In addition, the parameter setting and the time of the iterations are shown in Table 5.In Table 5, we can find that our method has better performance than other traditional methods.However, the time of running is higher, which we need to improve in the future.Then, we also test the performance of the proposed method on the DC3, as shown in Figure 6. Figure 6 shows the simulated data and the true fractional abundances for each of the five endmembers.There are both homogeneous region and heterogeneous region, as shown in Figure 6a. Figure 6b-f shows the true abundance maps of the five endmembers.Figure 7a-e shows the estimated abundance maps obtained by the proposed method, which contains five endmembers.Considering the sparse consistency constraint, the pixel of the heterogeneous region is not similar to its neighbour pixels which will makes the effect of the unmixing become worse, as shown in Figure 7. On the contrary, in the homogeneous region, the effects of the unmixing can be obtained.In Figure 7, the proposed method can obtain better results overall.We also compare the results with the other methods, as shown in Table 5.In addition, the parameter setting and the time of the iterations are shown in Table 5.In Table 5, we can find that our method has better performance than other traditional methods.However, the time of running is higher, which we need to improve in the future.In summary, from the qualitative and quantitative comparison between the results obtained by these methods, we can conclude that the proposed algorithm offers more promising results than the other unmixing methods.However, further experiments should be conducted with real hyperspectral scenes to fully substantiate our findings in real scene.

Real Datasets
The scene used in our real data experiments is the well-known AVIRIS Cuprite Dataset in west-central Nevada in 1997.This scene has been widely used to validate the performance of unmixing algorithms.The portion used in experiments corresponds to a 250 × 191-pixel subset of the data.Due to water absorption and low SNR, the number of spectral bands is reduced from 224 to 188.For illustrative purposes, Figure 8 shows a mineral map produced in 1995 by USGS.Although the map is available for hyperspectral data collected in 1995 and the publicly available AVIRIS Cuprite data were collected in 1997, the USGS map serves as a good indicator for qualitative assessment of the fractional abundance maps produced by the different unmixing algorithms.
The standard spectral library for these data, denoted by A, is the USGS library containing 498 pure endmember signatures.Before unmixing the real hyperspectral data, essential calibration was undertaken to mitigate the mismatches between the hyperspectral image and the signatures in the library [17].
Figure 9 shows a comparison between the fractional abundance maps by applying the SUnSAL, LRR, CLSUnSAL, and the proposed algorithm to the AVIRIS Cuprite scene using the library A. As shown in Figure 9, these endmembers are: alunite, buddingtonite and chalcedony, which are known to be present (in prominent fashion) in the Cuprite mining district [26].The parameters are set as follows: λ = 1 × 10 −3 for SUnSAL, λ = 10 for LRR and λ = 5 × 10 −3 for CLSUnSAL.For the proposed algorithm, we set threshold ε to be 0.01, λ = 10 and T = 25.In addition, the ASC constraint is not applied in our models in this section.In summary, from the qualitative and quantitative comparison between the results obtained by these methods, we can conclude that the proposed algorithm offers more promising results than the other unmixing methods.However, further experiments should be conducted with real hyperspectral scenes to fully substantiate our findings in real scene.

Real Datasets
The scene used in our real data experiments is the well-known AVIRIS Cuprite Dataset in west-central Nevada in 1997.This scene has been widely used to validate the performance of unmixing algorithms.The portion used in experiments corresponds to a 250 × 191-pixel subset of the data.Due to water absorption and low SNR, the number of spectral bands is reduced from 224 to 188.For illustrative purposes, Figure 8 shows a mineral map produced in 1995 by USGS.Although the map is available for hyperspectral data collected in 1995 and the publicly available AVIRIS Cuprite data were collected in 1997, the USGS map serves as a good indicator for qualitative assessment of the fractional abundance maps produced by the different unmixing algorithms.
The standard spectral library for these data, denoted by A, is the USGS library containing 498 pure endmember signatures.Before unmixing the real hyperspectral data, essential calibration was undertaken to mitigate the mismatches between the hyperspectral image and the signatures in the library [17].
Figure 9 shows a comparison between the fractional abundance maps by applying the SUnSAL, LRR, CLSUnSAL, and the proposed algorithm to the AVIRIS Cuprite scene using the library A. As shown in Figure 9, these endmembers are: alunite, buddingtonite and chalcedony, which are known to be present (in prominent fashion) in the Cuprite mining district [26].The parameters are set as follows: λ = 1 × 10 −3 for SUnSAL, λ = 10 for LRR and λ = 5 × 10 −3 for CLSUnSAL.For the proposed algorithm, we set threshold ε to be 0.01, λ = 10 and T = 25.In addition, the ASC constraint is not applied in our models in this section.Figure 9 shows the abundance maps estimated for the minerals alunite, buddingtonite, and chalcedony by applying the SUnSAL, LRR, CLSUnSAL and SLP-LRR algorithms to the AVIRIS Cuprite scene using library A. Because each endmember has several variations, the spectral library A which contains 498 members, includes about 229 kinds of materials.Thus, the mutual coherence of the library is very close to one.In Figure 9, it can be observed that the fractional abundances estimated by CLSUnSAL are generally higher in the regions assigned to the respective materials than SUnSAL.Due to the high correlation of the spectral library, the traditional unmixing algorithm SUnSAL may not obtain good performance.CLSUnSAL constrains the pixels to share the same set of endmembers.Thus, it obtains better results than SUnSAL.It is also worth noting that the abundance maps inferred by LRR do not exhibit good spatial consistency of minerals of interest compared with our proposed method.This is because the proposed method can better capture the global structure of the abundance matrix by using LRR model based on space consistency constraint as well as identify a subset of the spectral library by the spectral dictionary pruning, which only contains 42 pure endmembers.Compared with the size of the original spectral library, the number of endmembers of the subset is much smaller.Thus, using the subset as the endmember matrix improves the accuracy of the unmixing solutions significantly.Overall, the qualitative results reported in this section indicate that our proposed method which combines LRR model based on space consistency constraint and the dictionary pruning strategy, can obtain more accurate results than other methods.Further experiments with additional hyperspectral scenes and quantitative comparisons should be conducted in future work to fully objectify our findings.Figure 9 shows the abundance maps estimated for the minerals alunite, buddingtonite, and chalcedony by applying the SUnSAL, LRR, CLSUnSAL and SLP-LRR algorithms to the AVIRIS Cuprite scene using library A. Because each endmember has several variations, the spectral library A which contains 498 members, includes about 229 kinds of materials.Thus, the mutual coherence of the library is very close to one.In Figure 9, it can be observed that the fractional abundances estimated by CLSUnSAL are generally higher in the regions assigned to the respective materials than SUnSAL.Due to the high correlation of the spectral library, the traditional unmixing algorithm SUnSAL may not obtain good performance.CLSUnSAL constrains the pixels to share the same set of endmembers.Thus, it obtains better results than SUnSAL.It is also worth noting that the abundance maps inferred by LRR do not exhibit good spatial consistency of minerals of interest compared with our proposed method.This is because the proposed method can better capture the global structure of the abundance matrix by using LRR model based on space consistency constraint as well as identify a subset of the spectral library by the spectral dictionary pruning, which only contains 42 pure endmembers.Compared with the size of the original spectral library, the number of endmembers of the subset is much smaller.Thus, using the subset as the endmember matrix improves the accuracy of the unmixing solutions significantly.Overall, the qualitative results reported in this section indicate that our proposed method which combines LRR model based on space consistency constraint and the dictionary pruning strategy, can obtain more accurate results than other methods.Further experiments with additional hyperspectral scenes and quantitative comparisons should be conducted in future work to fully objectify our findings.

Discussion of the Parameters Setting and Time Complexity
We discuss the parameters setting and the time complexity.The relationship between the value of λ and the unmixing accuracy (SRE) with β = 0.01 and T = 1 is presented in Figure 10a.As the λ value increases, the SRE value increases and then gets stable.λ controls the relative weight of the noise of the solution.An appropriate value of λ makes a more accurate intermediate, i.e., J, X, and E can be achieved, and the unmixing results are very close to the true abundances.
The relationship between the value of β and the unmixing accuracy (SRE) with λ = 36 and T = 1 is presented in Figure 10b.When the β value is low, the SRE value is high.As the value of the β increases, the value of SRE decreases.The reason is that the data in DC2 have low space consistency, so the value of β needs to be low to get the high SRE value.The relationship between the T value and the unmixing accuracy (SRE) with λ = 36 and β = 0.01 is presented in Figure 10c.The threshold T acting as the stop condition controls the size of retained subset of the spectral library.As the T value increases, the number of retained endmembers increases, which leads to a lower SRE value.However, when the dataset is complex, i.e., the noise is high or the number of true endmembers is

Discussion of the Parameters Setting and Time Complexity
We discuss the parameters setting and the time complexity.The relationship between the value of λ and the unmixing accuracy (SRE) with β = 0.01 and T = 1 is presented in Figure 10a.As the λ value increases, the SRE value increases and then gets stable.λ controls the relative weight of the noise of the solution.An appropriate value of λ makes a more accurate intermediate, i.e., J, X, and E can be achieved, and the unmixing results are very close to the true abundances.
The relationship between the value of β and the unmixing accuracy (SRE) with λ = 36 and T = 1 is presented in Figure 10b.When the β value is low, the SRE value is high.As the value of the β increases, the value of SRE decreases.The reason is that the data in DC2 have low space consistency, so the value of β needs to be low to get the high SRE value.The relationship between the T value and the unmixing accuracy (SRE) with λ = 36 and β = 0.01 is presented in Figure 10c.The threshold T acting as the stop condition controls the size of retained subset of the spectral library.As the T value increases, the number of retained endmembers increases, which leads to a lower SRE value.However, when the dataset is complex, i.e., the noise is high or the number of true endmembers is large, the threshold should be set larger to avoid cutting off the true endmembers.We prefer to retain a larger size of subset of the spectral library rather than cut off any one of the true endmembers.large, the threshold should be set larger to avoid cutting off the true endmembers.We prefer to retain a larger size of subset of the spectral library rather than cut off any one of the true endmembers.
(a) (b) (c) Tables 2 and 4 show the computation time costs for the different algorithms on DC1 and DC2.In comparison to the other algorithms, our proposed algorithm costs more time due to the iterations.Thus, future works may focus on integrating the dictionary pruning into the optimization models to make them more efficient.

Conclusions
To improve the accuracy of spectral unmixing, this paper proposed a LRR algorithm based on the space consistency constraint and spectral dictionary pruning.Compared with the existing sparsity-based unmixing algorithms, the proposed algorithm, which takes the space consistency into account, can capture the global structure.Moreover, the high mutual coherence of spectral libraries always affects the performance of sparse unmixing.The proposed algorithm, which contains dictionary pruning steps, can identify a subset of the spectral library.The size of the subset, which contains the true endmembers in the spectral library, is much smaller than the original one.Thus, this strategy can not only mitigate the aforementioned effect but also improve the accuracy of unmixing significantly.Besides, combining our proposed dictionary pruning method with other unmixing algorithms can further boost the performances of the unmixing algorithms.
Both simulated hyperspectral datasets and real hyperspectral remote sensing image were used to test the performance of the proposed algorithm.Future work will focus on the method to prune the spectral library more efficiently, as it is quite significant for the unmixing methods.Our future work will focus on the study of tensor based low-rank representation model for hyperspectral unmixing.Tensor-based models have shown good performance in the field of computer vision and HSI classification.Having 3D representation capability, tensor-based models can maintain the spectral and spatial information of HSI data naturally.Therefore, we will design a tensor-based low-rank method to improve the hyperspectral unmixing further.Tables 2 and 4 show the computation time costs for the different algorithms on DC1 and DC2.In comparison to the other algorithms, our proposed algorithm costs more time due to the iterations.Thus, future works may focus on integrating the dictionary pruning into the optimization models to make them more efficient.

Conclusions
To improve the accuracy of spectral unmixing, this paper proposed a LRR algorithm based on the space consistency constraint and spectral dictionary pruning.Compared with the existing sparsity-based unmixing algorithms, the proposed algorithm, which takes the space consistency into account, can capture the global structure.Moreover, the high mutual coherence of spectral libraries always affects the performance of sparse unmixing.The proposed algorithm, which contains dictionary pruning steps, can identify a subset of the spectral library.The size of the subset, which contains the true endmembers in the spectral library, is much smaller than the original one.Thus, this strategy can not only mitigate the aforementioned effect but also improve the accuracy of unmixing significantly.Besides, combining our proposed dictionary pruning method with other unmixing algorithms can further boost the performances of the unmixing algorithms.
Both simulated hyperspectral datasets and real hyperspectral remote sensing image were used to test the performance of the proposed algorithm.Future work will focus on the method to prune the spectral library more efficiently, as it is quite significant for the unmixing methods.Our future work will focus on the study of tensor based low-rank representation model for hyperspectral unmixing.Tensor-based models have shown good performance in the field of computer vision and HSI classification.Having 3D representation capability, tensor-based models can maintain the spectral and spatial information of HSI data naturally.Therefore, we will design a tensor-based low-rank method to improve the hyperspectral unmixing further.

Figure 1 .
Figure 1.The flow chart of the proposed method.

Figure 1 .
Figure 1.The flow chart of the proposed method.

Figure 2 .
Figure 2. The schematic for space consistency of hyperspectral image.

Figure 2 .
Figure 2. The schematic for space consistency of hyperspectral image.

Figure 4 .
Figure 4. Abundance maps obtained by different unmixing methods for endmember #5 in DC1 and, from top to bottom, SNR is 20 dB, 30 dB, and 40 dB.

Figure 4 .
Figure 4. Abundance maps obtained by different unmixing methods for endmember #5 in DC1 and, from top to bottom, SNR is 20 dB, 30 dB, and 40 dB.

Figure 5 .
Figure 5. Ground-truth and estimated abundances obtained by different unmixing methods in the scene DC2, with SNR = 40 dB: (a) ground-truth abundance; (b) estimated abundances obtained by SUnSAL; (c) estimated abundances obtained by CLSUnSAL; and (d) estimated abundances obtained by SLP-SCC-LRR.

Figure 5 .
Figure 5. Ground-truth and estimated abundances obtained by different unmixing methods in the scene DC2, with SNR = 40 dB: (a) ground-truth abundance; (b) estimated abundances obtained by SUnSAL; (c) estimated abundances obtained by CLSUnSAL; and (d) estimated abundances obtained by SLP-SCC-LRR.

Figure 7 .
Figure 7. Abundance maps obtained by the proposed unmixing method in DC3 and the SNR is 40 dB: (a) the true abundance of endmember 1; (b) the true abundance of endmember 2; (c) the true abundance of endmember 3; (d) the true abundance of endmember 4; and (e) the true abundance of endmember 5.

Figure 7 .
Figure 7. Abundance maps obtained by the proposed unmixing method in DC3 and the SNR is 40 dB: (a) the true abundance of endmember 1; (b) the true abundance of endmember 2; (c) the true abundance of endmember 3; (d) the true abundance of endmember 4; and (e) the true abundance of endmember 5.

Figure 8 .
Figure 8. USGS map showing the location of different minerals in the Cuprite mining district in Nevada.The map is available online at http://speclab.cr.usgs.gov/cuprite95.tgif.2.2um_map.gif.

Figure 8 .
Figure 8. USGS map showing the location of different minerals in the Cuprite mining district in Nevada.The map is available online at http://speclab.cr.usgs.gov/cuprite95.tgif.2.2um_map.gif.

Figure 9 .
Figure 9. Abundance maps estimated for the minerals alunite, buddingtonite, and chalcedony by applying the SUnSAL, CLSUnSAL, SCC-LRR, and SLP-SCC-LRR algorithms to the AVIRIS Cuprite scene using the library A.

Table 1 .
SRE (dB)and AAD values achieved after applying different unmixing methods to Simulated DATA 1.

Table 2 .
The optimal parameters settings and running times(s) with the Simulated DATA 1.

Table 3 .
SRE (dB)and AAD values achieved after applying different unmixing methods to Simulated DATA 2.

Table 3 .
SRE (dB)and AAD values achieved after applying different unmixing methods to Simulated DATA 2.

Table 5 .
SRE (dB)and AAD values to Simulated Dataset 3 with the parameter settings and running times.

Table 5 .
SRE (dB)and AAD values to Simulated Dataset 3 with the parameter settings and running times.