Hyperspectral Unmixing with Robust Collaborative Sparse Regression

Chang Li 1, Yong Ma 2,∗, Xiaoguang Mei 2, Chengyin Liu 1 and Jiayi Ma 2 1 School of Electronic Information and Communications, Huazhong University of Science and Technology, Wuhan 430074, China; lichang@hust.edu.cn (C.L.); Liuchengyin@hust.edu.cn (C.L.) 2 Electronic Information School, Wuhan University, Wuhan 430072, China; meixiaoguang@gmail.com (X.M.); jiayima@whu.edu.cn (J.M.) * Correspondence: mayong@whu.edu.cn; Tel.: +86-186-0713-8001


Introduction
Over the last few decades, hyperspectral imaging (HSI) has been receiving considerable attention in different remote sensing applications such as spectral unmixing, object classification and matching [1][2][3][4][5].Due to insufficient spatial resolution of the imaging sensor and mixing effects of the ground surface, mixed pixels are widespread in hyperspectral images, which leads to difficulties for conventional pixel-level applications [6][7][8].Therefore, spectral unmixing is an essential step for the deep exploitation of hyperspectral image, which decomposes mixed pixels into a collection of pure spectra signatures, called endmembers, and their corresponding proportions in each pixel, called abundances [9,10].
When considering the problem of unmixing hyperspectral images, most of the literature in the geoscience and remote sensing areas adopts the widely used linear mixing model (LMM) due to the relative simplicity and straightforward interpretation.If the spectral endmembers are selected from a library containing a large number of spectral samples available a priori [11,12], then finding the optimal subset of signatures to best model the mixed pixel in the scene leads to a sparse solution [13].
Sparse unmixing (SU) assumes that the observed image can be formulated as finding the optimal subset of pure spectral signatures from a prior large spectral library, and it can typically be formulated as a linear sparse regression problem.To solve this problem, Bioucas et al. proposed sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL) [14], which ignores the spatial information.Iordache et al. proposed SUnSAL and total variation (SUnSAL-TV) [15] to exploit the spatial information for SU, which can obtain better unmixing performance than SUnSAL.In addition, some greedy algorithms have been developed for SU, such as the orthogonal matching pursuit (OMP) [16] and subspace matching pursuit (SMP) [17].Moreover, Iordache et al. proposed collaborative SUnSAL (CLSUnSAL) [18], which improves the unmixing results by adopting the collaborative (also called "multitask" or "simultaneous") sparse regression framework.The above-mentioned unmixing algorithms are all based on the commonly admitted linear mixing model (LMM).However, the LMM may be not valid in many situations, for example, when there are multi-scattering effects or intimate interactions, and nonlinear mixing models (NLMMs) provide an alternative to overcoming the inherent limitations of the LMM [19,20].
NLMMs have been proposed in the hyperspectral image processing and can be divided into two main classes [21].The first class of NLMMs consists of physical models based on the nature of the environment.These models include the bidirectional reflectance based model [22], Fan bilinear model (FM) [23], generalized bilinear model (GBM) [24], modified GBM (MGBM) [25] and multilinear mixing (MLM) model [26].The second class of NLMMs allows for more flexible models for other approximating physics-based models.These flexible models include the neural network model [27], kernel model [28,29], post-nonlinear model [30] and so on.However, one major drawback of these NLMMs is that they require a specific form of nonlinearity, which makes them limited in practice [31].Févotte et al. proposed the robust LMM (rLMM) [31] to overcome the above mentioned problems, which does not require for specification of any analytical form of the nonlinearity.Instead, nonlinearities are merely treated as outliers.
In this paper, to make the SU more flexible for all kinds of HSI unmixing in practice, we propose a new SU method called robust collaborative sparse regression (RCSR) based on rLMM.The RCSR simultaneously takes the collaborative sparse property of the abundance and the sparsely distributed additive property of the outlier into consideration, which can be formed as a robust joint sparse regression problem.The RCSR can be solved by the inexact augmented Lagrangian method (IALM) [32].
The main contribution of this work lies in that we propose a new SU method named RCSR based on the rLMM, which can be solved by the IALM.Experiments on both synthetic datasets and real hyperspectral images demonstrate that the proposed RCSR is efficient for solving the SU problem compared with the other four state-of-the-art algorithms.
The remainder of this paper is organized as follows.In Section 2, we briefly introduce the related work rLMM and describe the proposed RCSR for SU.In Section 3, we evaluate the performances of the proposed RSUs and the other four state-of-the-art algorithms on the synthetic datasets and real HSI, and conclude this paper in Section 4.

Robust Collaborative Sparse Regression
The LMM assumes that the spectral response of a pixel in any given spectral band is a linear combination of all of the endmembers presented in the pixels at the respective spectral band.The LMM can be written as follows: where y denotes a D × 1 vector of observed pixel spectra in a hyperspectral image, with D denoting the number of bands, A = [a 1 , ..., a D ] ∈ R D×M denoting the endmember signatures, with M denoting the number of endmembers, x ∈ R M×1 the abundance vector, and n the additive noise.The matrix formulation of LMM can be written as follows: where Y ∈ R D×B denotes the collected mixtures matrix, with B denoting the number of pixels, X ∈ R M×B denotes the abundance matrix, and N ∈ R D×B the collected additive noise.The abundances have to obey two constraints, namely, abundance nonnegativity constraint (ANC) and the abundance sum-to-one constraint (ASC), i.e., However, for real HSI applications, the ASC constraint does not always hold true in practice, since signature variability is usually intense in HSI [16,33].Therefore, the ASC constraint is not taken into consideration.
In [34], it has been proved that the probability of recovery failure decays exponentially in the number of channels, which demonstrates that multichannel sparse recovery is better than single channel methods.In addition, the probability bounds still hold true even for a small number of signals.In other words, for a real HSI, the number of endmembers is often much smaller than the number of pixels, which makes the SU have more chances to succeed.
In SU, hyperspectral vectors are approximated by a linear combination of a "small" number of spectral signatures in the library, and the number of columns are equal to the number of pixels, thus the nonzero abundance lines should appear in only a few lines [35], which implies sparsity along the pixels of an HSI.Since the collaborative (also called "simultaneous" or "multitask") sparse regression approach has shown advantages over the noncollaborative ones, i.e., the mutual coherence has a weaker impact on the unmixing [18,34,36].The "collaborative" means to impose sparsity among the endmembers simultaneously for all pixels.The CLSUnSAL (also called collaborative hierarchical lasso) imposes sparsity both at the group and individual level, which leads to a structured solution as the matrix of fractional abundances contains only a few nonzero lines [18].It is assumed that the abundance has the underlying collaborative sparse property, which is characterized by the 2,1 norm, and 2,1 -norm is defined as follows: The 2,1 norm imposes sparsity among the lines of X, which can promote a small number of nonzero lines of X.Since the pixels all share the same support, it is reasonable to enforce joint sparsity among all the pixels, which can be characterized by the 2,1 norm.Therefore, mathematically, the CLSUnSAL [18] can be written as follows: However, the CLSUnSAL is based on the LMM, and the LMM may be not valid when there are multi-scattering effects, and NLMMs provide an alternative to overcoming the inherent limitations of the LMM [19].Févotte et al. proposed the rLMM to make the blind unmixing of HSI more flexible to analyze a large variety of remotely sensed scenes, which takes the possible nonlinear effects into consideration, and the nonlinearities are merely treated as outliers [31].As for the blind unmixing of HSI, the endmember signatures A and the abundance matrix X are both unknown.However, for SU of HSI, the endmember signatures A are selected from a library containing a large number of spectral samples available a priori, and only the abundance matrix X is unknown.Until now, the rLMM has been not yet used for SU of HSI.In addition, the superiority of rLMM over LMM has been demonstrated in [31].
The rLMM assumes that the spectral response of a pixel in any given spectral band is approximated by a linear combination of all of the endmembers present in the pixel at the respective spectral band and the additive outlier [31]: where e denotes the outlier term (accounting for nonlinearity).The matrix formulation of rLMM can be written as follows: where Y ∈ R D×B denotes the collected mixtures matrix, with B denoting the number of pixels, X ∈ R M×B denotes the abundance matrix, E ∈ R D×B denotes the collected outlier, and N ∈ R D×B the collected additive noise.Févotte et al. [31] proposed a blind nonlinear hyperspectral unmixing method named robust nonnegative matrix factorization (rNMF) based on the rLMM, and the outliers are treated as sparsely distributed, which can also be characterized by the 2,1 norm.Mathematically, the rNMF [31] can be written as follows: min where ) is used to measure the dissimilarity, and d(x|y) is either the squared Euclidean distance or the Kullback-Leibler divergence.Blind unmixing of HSI aims at obtaining the endmembers and corresponding fractional abundances, knowing only the collected mixing spectral data.Thus, the endmember signatures A and the abundance matrix X are both unknown.
For sparse unmixing of HSI, the endmember signatures A are selected from a library containing a large number of spectral samples available a priori, and only the abundance matrix X is unknown.
The main difference lies in the endmember matrix A; for blind unmixing, the endmember matrix A is generally assumed to represent the pure materials present in the HSI.However, for sparse unmixing, the endmember matrix A relies on the existence of spectral libraries usually acquired in the laboratory, some of the endmembers in the spectral library are pure materials not present in the HSI.Mathematically, when we just specify A as the same spectral library and put a sparse constraint on the abundance matrix X, the sparse unmixing problem is a simpler version of the blind sparse unmixing problem.
To better pursue the outlier in rLMM, which has the underling sparsely distributed additive property, we also adopt the 2,1 -norm to impose group sparsity, which has the advantage of rotation invariant compared with the 1 norm [37,38].Therefore, mathematically, the proposed RCSR based on the rLMM can be written as follows: where .F represents the matrix Frobenius norm, and λ and α are two regularization parameters.Since the rLMM is a generalization of LMM, thus the proposed RCSR is a natural extension of CLSUnSAL with an additional outlier term, which makes the collaborative SU of HSI more robust for outliers.
The optimization problem Equation ( 9) can be solved by the IALM [32].By adding the auxiliary matrix P ∈ R M×N , the problem in Equation ( 9) can be reformulated as follows: Thus, the augmented Lagrangian function can be formed as follows: and then we apply the alternating minimization scheme to update the seven variables P, X, E, Λ, µ, i.e., update one of the five variables with the others fixed.To update P, we solve To update X, we solve whose solution is the well-known vect-soft threshold [39], applied independently to each row r of the update variable as follows: where ζ = P k+1 − Λ k /µ, and vect-soft(b, τ) denotes the row-wise application of the vect-soft-threshold function To update E, we solve which can be also solved by the well-known vect-soft threshold [39] E k+1 (r, :) = max(vect-soft(γ(r, :), α 2 ), 0), (16) where γ = Y − AP k+1 .To sum up, the detailed procedure for solving the proposed RCSR is listed in Algorithm 1.
Algorithm 1: Solving (9) with IALM Input: Y, A; Output: X, Ê; 1 Initialize X 0 , E 0 , P 0 Λ 0 , λ, α, ρ µ, µ max = 10 6 , k=0; Compute X k+1 by ( 14); 5 Compute E k+1 by ( 16); The IALM is a variation of the exact augmented Lagrangian method, and its convergence has been studied for at most two blocks (i.e., unknown matrix variables) [40].For our problem Equation ( 8), there is no guarantee for the convergence in theory, and ε is an error threshold.Furthermore, the IALM is known to generally perform well in reality [40].In practice, when we choose the parameters appropriately, it can be observed that the proposed RCSR convergences before the maximum iteration is reached.

Experiments
In this section, we first carry out simulated experiments to demonstrate the advantages of the proposed RCSRs compared with four algorithms based on the LMM, i.e., SUnSAL [14], CLSUnSAL [18], OMP [16] and SMP [17].To evaluate the performance of different HSI SU algorithms, the signal-to-reconstruction error (SRE) [16] is adopted to measure the power between the signal and error, which is defined as follows: where X and X are the actual and estimated abundance, respectively.Generally speaking, larger SRE means better hyperspectral sparse unmixing performance.

Experimental Results with Synthetic Data
We use the spectral library randomly selected from the United States Geological Survey (USGS) digital spectral library (Available at: http://speclab.cr.usgs.gov/spectral-lib.html), which has 224 spectral bands uniformly ranging from 0.4 µm to 2.5 µm, and contains 498 spectral signatures of endmembers.We generate the synthetic HSI based on the LMM [16], FM [23], GBM [24] and MGBM [25], and the latter three models are nonlinear unmixing models.
We tune the compared SUnSAL and CLSUnSAL to their best performances by using different regularization parameters: 10 −5 , 10 −4 , 10 −3 , 10 −2 , 10 −1 and 1.The maximum number of iterations and error tolerances of SUnSAL and CLSUnSAL are set to 1000 and 10 −6 , respectively.For OMP, we set the correct number of endmembers as the input parameter.For SMP, we set the given threshold δ = 10 −3 .For the proposed RCSR, the performance is tuned to the best by setting λ and α to the following parameters: 10 −5 , 10 −4 , 10 −3 , 10 −2 , 10 −1 and 1.The maximum number of iterations and error tolerances of RCSR is set to the same as those of SUnSAL and CLSUnSAL.To avoid unnecessary deviation, we repeat the simulations 10 times to obtain the mean SREs.
The synthetic HSIs all have 100 × 100 pixels using endmembers randomly chosen from the USGS library, and all of the abundance fractions are generated following the Dirichlet distribution, which satisfy the ANC.The obtained datacubes are then contaminated by Gaussian white noise and correlated noise with different signal-to-noise ratio SNR = 10 log 10 ( Y 2 F / N 2 F ). Figure 1 shows the performance of SRE as a function of the number of endmember under Gaussian white noise when the SNR is 10 with the LMM, FM, GBM and MGBM, respectively.It can be easily seen from Figure 1 that the proposed RCSR generally obtains the best SRE.In addition, the performances of most algorithms tend to get worse as the number of endmembers increase, which is due to the fact that the spectral signatures in the spectral library is usually highly correlated.To test the performances of different algorithms at high SNR, Figure 2 shows the performance of SRE as a function of the number of endmembers under Gaussian white noise when the SNR is 50 with the LMM, FM, GBM and MGBM, respectively.It is shown that the RCSRs all have the best SRE, and the RCSR generally performs better than CLSUnSAL when the SNR is 50 with the three nonlinear models, which demonstrates that the RCSR is more robust for outliers than CLSUnSAL.Compared with Figure 1, the SREs in Figure 2 are generally higher than those in Figure 1, which indicates that the level of noise has a large impact on the final performance of SU.In addition, we study the level of noise on the performance of SU of HSI. Figure 3 shows the performance of SRE as a function of SNR under Gaussian white noise when the number of endmembers is four with the LMM, FM, GBM and MGBM, respectively.It can be seen that the RCSR obtains the best SRE at some time, and SMP and CLSUnSAL sometimes obtain the best SREs.However, the performances of RCSR are much more stable than those of SMP and CLSUnSAL.Moreover, Figure 4 shows the performance of SRE as a function of SNR under Gaussian white noise when the number of endmembers is 20 with the LMM, FM, GBM and MGBM, respectively.From Figure 4, it can be observed that the RCSR has comparably good SREs with CLSUnSAL.In addition, when the SNR is 10, the performances of RCSR is obvious better than CLSUnSAL, which demonstrates that the RCSR is more robust for noise than CLSUnSAL.Since it is hard to calibrate the hyperspectral data obtained from an airborne or spaceborne sensor, the noise and the spectra in real hyperspectral imaging applications are usually low-pass type, which makes the noise highly correlated [16].
Thus, it is very necessary to conduct experiments when the obtained datacubes are contaminated by correlated noise.
The synthetic HSI has 100 × 100 pixels, and all of the abundance fractions are generated following the Dirichlet distribution, which satisfy the ASC.The obtained datacubes are then contaminated with correlated noise, and we generate the correlated noise with the correlated noise function available at: http://www.mathworks.com/matlabcentral/fileexchange/21156-correlated-Gaussian-noise/content/correlatedGaussianNoise.m.Figures 5 and 6 show the SU results contaminated by correlated noise when the SNR is 10 and 50, respectively.Figures 7 and 8 show the SU results contaminated by correlated noise when the number of endmembers is 4 and 20, respectively.From Figures 5-8, it can be deduced that the RCSR can obtain the best SRE in most circumstances.In addition, the RCSR has comparably good results with CLSUnSAL at high SNRs, and the RCSR can better handle the outliers than CLSUnSAL.Moreover, SMP sometimes obtains the best SREs.However, the performances of RCSR are much more stable than those of SMP.From Figures 1-8, it can be observed that when fixing the SNR and the number of endmembers, the variation of SRE is small with regard to four different unmixing models, which is due to the fact that the energy of the linear mixing part is much bigger than that of the outlier.

Experimental Results with Real Data
The real dataset used in our experiment is the most benchmarked dataset for hyperspectral unmixing, which was captured by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) over a Cuprite mining district in June 1997 in the state of Nevada.The mineral map of the cuprite image is available at http://speclab.cr.usgs.gov/cuprite95.tgif.2.2um_map.gif.Some spectral bands (1-2, 104-113, 148-167 and 221-224) have been removed due to noise corruption and atmospheric absorption, leading to P = 188 spectral bands ranging from 0.4 µm to 2.5 µm with a nominal bandwidth of 10 nm.The false color image is shown in Figure 9, which is the size 250 × 191.Since the minerals of the Cuprite are all included in the USGS library, we adopt a subset of the USGS library as the spectral library for the SU of Cuprite, which contains 240 members.We set the regularization parameters of SUnSAL and CLSUnSAL to 10 −3 and 10 −1 , respectively.The maximum number of iterations and error tolerances of SUnSAL and CLSUnSAL are set to 1000 and 10 −5 , respectively.For OMP, we set the correct number of endmembers as the input parameter.For SMP, we set the given threshold δ = 10 −3 .For the proposed RCSR, the regularization parameters λ and α are set to 10 −1 and 10 −1 , respectively.The maximum number of iterations and error tolerances of RCSR are set to the same as those of SUnSAL and CLSUnSAL.Figure 10 shows the fractional abundance maps estimated by RCSR and the other four compared SU methods using the subimages of the AVIRIS Cuprite dataset.It can be clearly seen from Figure 10 that the visual performances of fractional abundance using the RCSR and the other four compared SU methods are very consistent with each other.

Conclusions
In this paper, we propose an RCSR for SU of HSI, which is based on the rLMM.The RCSR takes the possible nonlinear effects (i.e., outliers) into consideration, which exploits the collaborative sparse property of the abundance and sparsely distributed additive properties of the outliers.The RCSR can be formed as a robust joint sparse regression problem, which can be solved by the IALM.Experiments on both synthetic datasets and real hyperspectral images demonstrate that the proposed RCSR is efficient for solving the hyperspectral unmixing problem compared with the other four state-of-the-art algorithms.
Super resolution based SU is a recently developed spectral unmixing approach, in the future, we will consider how to apply the proposed SU model to the hyperspectral face image to super-resolve a high-resolution face [41].

Figure 1 .
Figure 1.Performance of SRE as a function of the number of endmember under Gaussian white noise when the SNR is 10 with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 2 .
Figure 2. Performance of SRE as a function of the number of endmembers under Gaussian white noise when the SNR is 50 with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 3 .
Figure 3. Performance of SRE as a function of SNR under Gaussian white noise when the number of endmembers is four with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 4 .
Figure 4. Performance of SRE as a function of SNR under Gaussian white noise when the number of endmembers is 20 with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 5 .
Figure 5. Performance of SRE as a function of the number of endmembers under correlated noise when the SNR is 10 with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 6 .
Figure 6.Performance of SRE as a function of the number of endmembers under correlated noise when the SNR is 50 with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 7 .
Figure 7. Performance of SRE as a function of SNR under correlated noise when the number of endmembers is 4 with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 8 .
Figure 8. Performance of SRE as a function of SNR under correlated noise when the number of endmembers is 20 with the (a) LMM; (b) FM; (c) GBM and (d) MGBM.

Figure 10 .
Figure 10.Fractional abundance maps estimated by different unmixing methods using the subimage of the AVIRIS Cuprite dataset.