# Backtracking-Based Simultaneous Orthogonal Matching Pursuit for Sparse Unmixing of Hyperspectral Data

Authors: Fanqiang Kong; Wenjun Guo; Yunsong Li; Qiu Shen; Xin Liu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Sparse unmixing is a promising approach in a semisupervised fashion by assuming that the observed signatures of a hyperspectral image can be expressed in the form of linear combination of only a few spectral signatures (endmembers) in an available spectral library. Simultaneous orthogonal matching pursuit (SOMP) algorithm is a typical simultaneous greedy algorithm for sparse unmixing, which involves finding the optimal subset of signatures for the observed data from a spectral library. But the numbers of endmembers selected by SOMP are still larger than the actual number, and the nonexisting endmembers will have a negative effect on the estimation of the abundances corresponding to the actual endmembers. This paper presents a variant of SOMP, termed backtracking-based SOMP (BSOMP), for sparse unmixing of hyperspectral data. As an extension of SOMP, BSOMP incorporates a backtracking technique to detect the previous chosen endmembers’ reliability and then deletes the unreliable endmembers. Through this modification, BSOMP can select the true endmembers more accurately than SOMP. Experimental results on both simulated and real data demonstrate the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

With the rapid development of space technology, hyperspectral remote sensing image has gained more and more attention in many application domains, such as environmental monitoring, target detection, material identification, mineral exploration, and military surveillance. However, due to the low spatial resolution of the hyperspectral imaging sensor, each pixel in the hyperspectral image often contains a mixture of several different materials. In order to deal with the problem of spectral mixing, hyperspectral unmixing is used to decompose each pixel’s spectrum to identify and quantify the fractional abundances of the pure spectral signatures or endmembers in each mixed pixel [1, 2]. There are two basic models of hyperspectral unmixing used to analyze the mixed pixel problem: the linear mixture model and the nonlinear mixture model. The linear mixture model assumes that each mixed pixel can be expressed as a linear combination of endmembers weighted by their corresponding abundances [3]. The linear mixture model has been widely applied for spectral unmixing, due to its computational tractability and flexibility. Under the linear mixture model, the traditional linear spectral unmixing algorithms based on geometry [4–8], statistics [5, 9], and nonnegative matrix factorization [10–12] have been proposed. However, some of these methods [10–12] are unsupervised and could extract virtual endmembers with no physical meaning. In [5], the presence in the data of at least one pure pixel per endmember is assumed. If the pure pixel assumption is not fulfilled because of the inadequate spatial resolution and the microscopic mixture of distinct materials, the unmixing results will not be accurate and the unmixing process is a rather challenging task.Sparse unmixing, as a semisupervised method, has been proposed to overcome this challenge. It assumes that the observed image can be expressed as a linear combinations of spectral signatures from a spectral library that is known in advance [3, 13]. But the number of spectral signatures in the spectral library is much larger than the number of endmembers in the hyperspectral image; the sparse unmixing model is combinatorial and difficult to find a unique, stable, and optimal solution. Fortunately, sparse linear regression techniques can be used [14, 15] to solve it.Several sparse regression techniques, such as greedy algorithms (GAs) [16–19] and convex relaxation methods [20–23], are usually adopted to solve the sparse unmixing problem. Convex relaxation methods, such as SUnSAL [20], SUnSAL-TV [21], and CLSUnSAL [22], use the alternating direction method of multipliers to efficiently solve the l1 norm sparse regression problem which can decompose a complex problem into several simpler ones. Convex relaxation methods can obtain the global sparse optimization solution and are more sophisticated than the Gas; however they are far more complicated than the GAs. The GAs, such as OMP [16], SP [17], and CGP [18, 19], adopt one or more potential endmembers from the spectral library in each iteration that explains the largest correlation between the current residual of one input signal and the spectral library. The GAs can get an approximate solution for the l0 norm problem without smoothing the penalty function and have low computational complexity. However the endmembers selection criterion of GAs is not optimal in the sense of minimizing the residual of the new approximation; it means that a nonexisting endmember once selected into the supporting set will never be deleted. So the GAs tend to be trapped into the local optimum and are likely to miss some of the actual endmembers. To solve the local optimal solutions problem of GAs, several representative simultaneous greedy algorithms (SGAs) are presented, including simultaneous subspace pursuit (SSP) [24], subspace matching pursuit (SMP) [24], and simultaneous orthogonal matching pursuit (SOMP) [25]. These SGAs divide the whole hyperspectral image into several blocks and pick some potential endmembers from the spectral library in each block. Then the endmembers picked in each block are associated as the endmembers sets of the whole hyperspectral data. Finally, the abundances are estimated using the whole hyperspectral data with the obtained endmembers sets. The SGAs have the same low computational complexity as the GAs and can find the actual endmembers far more accurately than the GAs. The SGAs adopt a block-processing strategy to efficiently solve the local optimal problem, but the endmembers picked in all the blocks are not all actual endmembers and the nonexisting endmembers will affect the estimation of the abundances corresponding to the actual endmembers. To solve the drawback of the SGAs, RSFoBa [26] combinates a forward greedy step and a backward greedy step, which can select the actual endmembers more accurately than the SGAs.Inspired by the existing SGA methods, we propose a sparse unmixing algorithm termed backtracking-based simultaneous orthogonal matching pursuit (BSOMP) in this paper. Similar to SOMP and SMP, it uses a block-processing strategy to select some potential endmembers and adds them to the estimated endmembers set, which divides the whole hyperspectral image into several blocks and picks some potential endmembers from the spectral library in each block. Furthermore, BSOMP incorporates a backtracking process to detect the previous chosen endmembers’ reliability and then deletes the unreliable endmember from the estimated endmembers set in each iteration. Through this modification, BSOMP can identify the true endmembers set more accurately than the other considered SGA methods.The remainder of the paper is organized as follows. Section2 introduces the simultaneous sparse unmixing model. In Section 3, we present the proposed BSOMP algorithm and give out the theoretical analysis for the algorithm. Section 4 provides a quantitative comparison between BSOMP and previous sparse unmixing algorithms, using both simulated hyperspectral and real hyperspectral data. Finally, we conclude in Section 5.


## 2. Simultaneous Sparse Unmixing Model

The linear sparse unmixing model assumes that the observed spectrum of a mixed pixel is a linear combination of a few spectral signatures presented in a known spectral library. Lety∈RL denote the measured spectrum vector of a mixed pixel with L bands, A∈RL×m, where m is the number of signatures in the library A; denote a L×m spectral library; then the linear sparse unmixing model can be expressed as follows [23]:(1)y=Ax+n,where x∈Rm denotes the fractional abundance vector with regard to the library A and n∈RL is the noise. Considering physical constraints, abundance nonnegativity constraint (ANC) and abundance sum-to-one constraint (ASC) are imposed on the linear sparse model as follows: (2)x≥0,∑i=1mxi=1,where xi is the ith element of x.The sparse unmixing problem can be expressed as follows:(3)min⁡xx0s.t.y-Ax22≤δ,x≥0,where x0 (called the l0 norm) denotes the number of nonzero atoms in x and δ≥0 is the tolerated error due to the noise and model error. It is worth mentioning that we do not explicitly add the ASC in (3), because the hyperspectral libraries generally contain only nonnegative components, and the nonnegativity of the sources automatically imposes a generalized ASC [23].The simultaneous sparse unmixing model assumes that several input signals can be expressed in the form of different linear combinations of the same elementary signals. This means that all the pixels in the hyperspectral image are constrained to share the same subset of endmembers selected from the spectral library. Then we can use SGA methods for sparse unmixing; the sparse unmixing model in (1) becomes(4)Y=AX+N,where Y∈RL×K denote the hyperspectral data matrix with L bands and K mixed pixels, A∈RL×m denote a L×m spectral library, X∈Rm×K denotes the fractional abundance matrix, each column of which corresponds with the abundance fractions of a mixed pixel, and N∈RL×K is the noise matrix.Under the simultaneous sparse unmixing model, the simultaneous sparse unmixing problem can be expressed as follows:(5)min⁡XXrow-0s.t.Y-AXF2≤δ,X≥0,where XF denotes the Frobenius norm of X, Xrow-0 is the number of nonzero rows in matrix X [25], and δ is the tolerated error due to the noise and model error.It should be noted that the model in (5) is reasonable because there should be only a few rows with nonzero entries in the abundance matrix in light of only a small number of endmembers in the hyperspectral image, compared with the dimensionality of the spectral library [22].


## 3. Backtracking-Based SOMP

In this section, we first present our new algorithm, BSOMP, for sparse unmixing of hyperspectral data. Then, a theoretical analysis of the algorithm will be given. Then, we give a convergence theorem for the proposed algorithm.


### 3.1. Statement of Algorithm

The whole process of using BSOMP for sparse unmixing of hyperspectral data is summarized in Algorithm1. The algorithm includes three main parts: SOMP for endmember selection, backtracking processing, and abundance estimation. In the first part, we adopt a block-processing strategy for SOMP to efficiently select endmembers. This strategy divides the whole hyperspectral image into several blocks. Then, in each block, SOMP will pick several potential endmembers from the spectral library and add them to the estimated endmembers set. In the second part, we utilize a backtracking strategy to remove some endmembers chosen wrongly from the estimated endmembers set in previous processing and identify the true endmembers set more accurately. The backtracking processing is halted when the maximum total correlation between an endmember in the spectral library and the residual drops below threshold δ. Finally, the abundances are estimated using the obtained endmembers set under the constraint of nonnegativity.Algorithm 1:BSOMP for hyperspectral sparse unmixing.
Part 1 (SOMP):
(1) Initialize hyperspectral dataY and spectral library A
(2) Divide hyperspectral dataY into several blocks: Y=[Y1,Y2,…,YB], and initialize the index set SSOMP=⌀
(3)Foreach block do
(4)  Set index setSb=⌀ and iteration counter k=1. Initialize the residual data of block b: R0=Yb;
(5)Whilestopping criterion 1 has not been met do
(6)Compute the index of the best correlated member of A to the actual residual: j=arg⁡maxi⁡(Rl-1)TAi2, where Ai is the
ith column of A
(7)Update support set: Sb=Sb∪j
(8)Compute Xl=ASb∗ASb-1ASb∗Yb  {ASb is the matrix containing the columns of A having the indexes from Sb}
(9)Update residual:  Rl=Yb-ASbXl
(10)l←l+1
(11)End while
(12)  SetSSOMP=SSOMP∪Sb
(13)End for
Part 2 (Backtracking process):
(14) Initialize hyperspectral dataY, the index set St=SSOMP and At=ASSOMP  {At is the matrix containing the columns of A having
the indexes from St}
(15)While stopping criterion 2 has not been met do
(16)Compute solution: Xt=At∗At-1At∗Y
(17)Compute the member of At having the lowest abundance: index←minj⁡(Xtj)
(18)Remove the member having the lowest fractional abundance from the index set St←St∖Sindex and the endmember set:
At←At∖Aindex
(19)t←t+1
(20)End while
Part 3 (Abundance estimation):
(21) Estimate abundances using the original hyperspectral data matrix and the endmember set under the constraint of
nonnegativity.
X←arg⁡minX⁡AtX-YF, subject to X≥0.In the following paragraphs, we provide some definitions and then give two stopping conditions for the BSOMP algorithm.Definition 1.
Ifx is a vector, where xi is the ith element of x, the lp norm of x is defined as(6)xp=∑ixip1/pfor1≤p<∞,x∞=max⁡i⁡xi.Definition 2.
The(p,q) operator norms of the matrix A is defined as (7)Ap,q=max⁡x≠0⁡Axqxp.Several of the(p,q) operator norms can be computed easily using the following lemma [25, 27].Lemma 3.
(
1
) The (1,q) operator norm is the maximum lp norm of any column of A.
(
2
) The (2,2) operator norm is the maximum singular value of A.
(
3
) The (p,∞) operator norm is the maximum lp norm of any row of A.We give a stopping condition for the first part of the BSOMP algorithm:(8)RlF-Rl-1FRl-1F≤σ,where RlF is the Frobenius norm of the residual corresponding to the lth iteration and σ is a threshold.For the second part of the BSOMP algorithm, we give a stopping condition for deciding when to halt the iteration:(9)A∗Rt∞,∞-A∗Rt-1∞,∞A∗Rt∞,∞≤δ,where Yt is the residual corresponding to the tth iteration, A∗ denotes the conjugate transpose of A, A∗Rt∞,∞ is the maximum total correlation between an endmember in the spectral library A and the residual Rt,  and δ is a threshold.


### 3.2. Theoretical Analysis

In this section, we will introduce some definitions and then give a theoretical analysis for the BSOMP algorithm.Definition 4.
The coherenceμ of a spectral library A equals the maximum correlation between two distinct spectral signatures [28]:(10)μ=max⁡k≠h⁡Ak,Ah,where Ak is the kth column of A. If the coherence parameter is small, each pair of spectral signatures is nearly orthogonal.Definition 5.
The cumulative coherence [29, 30] is defined as (11)μk=max⁡Λ≤k⁡max⁡j∉Λ⁡∑i∈ΛAi,Aj,where the index set Λ∈Sall, Sall is the support of all the spectral signatures in the spectral library A. The cumulative coherence measures the maximum total correlation between a fixed spectral signature and t distinct spectral signatures. Particularly, μ(0)=0.Definition 6.
Given hyperspectral image data matrixY∈RL×K and available spectral library matrix A∈RL×m, Sopt is an index set which lists at most T endmembers presented in the hyperspectral image scene;Aopt is an optimal submatrix of the spectral library A which contains the T columns indexed in Sopt;Xopt is the optimal abundance matrix corresponding to Aopt;Yopt is a T-term best approximation of Y, and Yopt=AoptXopt;St is an index set which lists k endmembers after t iterations and assumes Sopt∈St;At is a submatrix of the spectral library A which contains the k columns indexed in St;Xt is the abundance matrix corresponding to At;Yt is the reconstruction of Y by BSOMP algorithm after t iterations (assume that BSOMP algorithm selects a nonexisting endmember to remove from index set St in each iteration), and Yt=AtXt;Rt is the residual error of Y after t iterations, and Rt=Y-Yt.Theorem 7.
Suppose that the best approximationYopt of the hyperspectral image data matrix Y over the index set Sopt satisfies the error bound: A∗(Y-Yopt)∞,∞≤τ, the first part of BSOMP has selected k potential endmembers from the spectral library to constitute the initial index set St, and Sopt∈St. After iteration t of the second part of BSOMP, halt the second part of the algorithm if (12)A∗Rt∞,∞≤1-μT1-μT-μk-Tτ,where k≥T. If this inequality fails, then BSOMP must remove another nonexisting endmember in iteration (t+1).The proof of Theorem7 is shown in the Appendix. If the second part of the BSOMP algorithm terminates at the end of iteration t, we may conclude that the BSOMP algorithm removes t nonexisting endmembers from St and has chosen k optimal endmembers.


## 3.1. Statement of Algorithm

The whole process of using BSOMP for sparse unmixing of hyperspectral data is summarized in Algorithm1. The algorithm includes three main parts: SOMP for endmember selection, backtracking processing, and abundance estimation. In the first part, we adopt a block-processing strategy for SOMP to efficiently select endmembers. This strategy divides the whole hyperspectral image into several blocks. Then, in each block, SOMP will pick several potential endmembers from the spectral library and add them to the estimated endmembers set. In the second part, we utilize a backtracking strategy to remove some endmembers chosen wrongly from the estimated endmembers set in previous processing and identify the true endmembers set more accurately. The backtracking processing is halted when the maximum total correlation between an endmember in the spectral library and the residual drops below threshold δ. Finally, the abundances are estimated using the obtained endmembers set under the constraint of nonnegativity.Algorithm 1:BSOMP for hyperspectral sparse unmixing.
Part 1 (SOMP):
(1) Initialize hyperspectral dataY and spectral library A
(2) Divide hyperspectral dataY into several blocks: Y=[Y1,Y2,…,YB], and initialize the index set SSOMP=⌀
(3)Foreach block do
(4)  Set index setSb=⌀ and iteration counter k=1. Initialize the residual data of block b: R0=Yb;
(5)Whilestopping criterion 1 has not been met do
(6)Compute the index of the best correlated member of A to the actual residual: j=arg⁡maxi⁡(Rl-1)TAi2, where Ai is the
ith column of A
(7)Update support set: Sb=Sb∪j
(8)Compute Xl=ASb∗ASb-1ASb∗Yb  {ASb is the matrix containing the columns of A having the indexes from Sb}
(9)Update residual:  Rl=Yb-ASbXl
(10)l←l+1
(11)End while
(12)  SetSSOMP=SSOMP∪Sb
(13)End for
Part 2 (Backtracking process):
(14) Initialize hyperspectral dataY, the index set St=SSOMP and At=ASSOMP  {At is the matrix containing the columns of A having
the indexes from St}
(15)While stopping criterion 2 has not been met do
(16)Compute solution: Xt=At∗At-1At∗Y
(17)Compute the member of At having the lowest abundance: index←minj⁡(Xtj)
(18)Remove the member having the lowest fractional abundance from the index set St←St∖Sindex and the endmember set:
At←At∖Aindex
(19)t←t+1
(20)End while
Part 3 (Abundance estimation):
(21) Estimate abundances using the original hyperspectral data matrix and the endmember set under the constraint of
nonnegativity.
X←arg⁡minX⁡AtX-YF, subject to X≥0.In the following paragraphs, we provide some definitions and then give two stopping conditions for the BSOMP algorithm.Definition 1.
Ifx is a vector, where xi is the ith element of x, the lp norm of x is defined as(6)xp=∑ixip1/pfor1≤p<∞,x∞=max⁡i⁡xi.Definition 2.
The(p,q) operator norms of the matrix A is defined as (7)Ap,q=max⁡x≠0⁡Axqxp.Several of the(p,q) operator norms can be computed easily using the following lemma [25, 27].Lemma 3.
(
1
) The (1,q) operator norm is the maximum lp norm of any column of A.
(
2
) The (2,2) operator norm is the maximum singular value of A.
(
3
) The (p,∞) operator norm is the maximum lp norm of any row of A.We give a stopping condition for the first part of the BSOMP algorithm:(8)RlF-Rl-1FRl-1F≤σ,where RlF is the Frobenius norm of the residual corresponding to the lth iteration and σ is a threshold.For the second part of the BSOMP algorithm, we give a stopping condition for deciding when to halt the iteration:(9)A∗Rt∞,∞-A∗Rt-1∞,∞A∗Rt∞,∞≤δ,where Yt is the residual corresponding to the tth iteration, A∗ denotes the conjugate transpose of A, A∗Rt∞,∞ is the maximum total correlation between an endmember in the spectral library A and the residual Rt,  and δ is a threshold.


## 3.2. Theoretical Analysis

In this section, we will introduce some definitions and then give a theoretical analysis for the BSOMP algorithm.Definition 4.
The coherenceμ of a spectral library A equals the maximum correlation between two distinct spectral signatures [28]:(10)μ=max⁡k≠h⁡Ak,Ah,where Ak is the kth column of A. If the coherence parameter is small, each pair of spectral signatures is nearly orthogonal.Definition 5.
The cumulative coherence [29, 30] is defined as (11)μk=max⁡Λ≤k⁡max⁡j∉Λ⁡∑i∈ΛAi,Aj,where the index set Λ∈Sall, Sall is the support of all the spectral signatures in the spectral library A. The cumulative coherence measures the maximum total correlation between a fixed spectral signature and t distinct spectral signatures. Particularly, μ(0)=0.Definition 6.
Given hyperspectral image data matrixY∈RL×K and available spectral library matrix A∈RL×m, Sopt is an index set which lists at most T endmembers presented in the hyperspectral image scene;Aopt is an optimal submatrix of the spectral library A which contains the T columns indexed in Sopt;Xopt is the optimal abundance matrix corresponding to Aopt;Yopt is a T-term best approximation of Y, and Yopt=AoptXopt;St is an index set which lists k endmembers after t iterations and assumes Sopt∈St;At is a submatrix of the spectral library A which contains the k columns indexed in St;Xt is the abundance matrix corresponding to At;Yt is the reconstruction of Y by BSOMP algorithm after t iterations (assume that BSOMP algorithm selects a nonexisting endmember to remove from index set St in each iteration), and Yt=AtXt;Rt is the residual error of Y after t iterations, and Rt=Y-Yt.Theorem 7.
Suppose that the best approximationYopt of the hyperspectral image data matrix Y over the index set Sopt satisfies the error bound: A∗(Y-Yopt)∞,∞≤τ, the first part of BSOMP has selected k potential endmembers from the spectral library to constitute the initial index set St, and Sopt∈St. After iteration t of the second part of BSOMP, halt the second part of the algorithm if (12)A∗Rt∞,∞≤1-μT1-μT-μk-Tτ,where k≥T. If this inequality fails, then BSOMP must remove another nonexisting endmember in iteration (t+1).The proof of Theorem7 is shown in the Appendix. If the second part of the BSOMP algorithm terminates at the end of iteration t, we may conclude that the BSOMP algorithm removes t nonexisting endmembers from St and has chosen k optimal endmembers.


## 4. Experimental Results

In this section, we use both simulated hyperspectral data and real hyperspectral data to demonstrate the performance of the proposed algorithm. The BSOMP algorithm is compared with three SGAs (SOMP [25], SMP [24], and RSFoBa [26]) and three convex relaxation methods (i.e., SUnSAL [20], SUnSAL-TV [21], and CLSUnSAL [22]). All the considered algorithms have taken into account the abundance nonnegativity constraint. The TV regularizer used in SUnSAL-TV is a nonisotropic one. For the RSFoBa algorithm, the l∞ norm is considered in our experiment. As mentioned above, when the spectral library is overcomplete, the GAs behave worse than the SGAs and the convex relaxation methods; thus, we do not include the GAs in our experiment.In the simulated data experiments, we evaluated the performances of the sparse unmixing algorithms in situations of different SNR level of noise; the signal-to-reconstruction error (SRE) [21] is used to measure the quality of the reconstruction abundances of spectral endmembers: (13)SRE=EX22EX-X^22,SREdB=10log⁡10⁡SRE,where X denotes the true abundances of spectral endmembers and X^ denotes the reconstruction abundances of spectral endmembers.


### 4.1. Evaluation with Simulated Data

In the simulated experiments, the United States Geological Survey (USGS) digital spectral library (splib06a) [31] is used to build the spectral library A. The reflectance values of 498 spectral signatures are measured for 224 spectral bands distributed uniformly in the interval of 0.4–2.5 μm (A∈R224×498). Fifteen spectral signatures are chosen from the spectral library to generate our simulated data. Figure 1 shows five spectral signatures used for all the simulated data experiments. The other ten spectral signatures that are not displayed in the figure include Monazite HS255.3B, Zoisite HS347.3B, Rhodochrosite HS67, Neodymium Oxide GDS34, Pigeonite HS199.3B, Meionite WS700.Hlsep, Spodumene HS210.3B, Labradorite HS17.3B, Grossular WS484, and Wollastonite HS348.3B.Figure 1
Five spectral signatures from the USGS used in our simulated data experiments. The title of each subimage denotes the mineral corresponding to the signature.In Simulated Data 1 experiment, we generated seven datacubes of50×50 pixels and 224 bands per pixel, each containing a different number of endmembers: 3, 5, 7, 9, 11, 13, and 15. In each simulated pixel, the fractional abundances of the endmembers were randomly generated following a Dirichlet distribution [23]. Note that there was no pure pixel in Simulated Data 1, and the fractional abundances of the endmembers were less than 0.8. After Simulated Data 1 was generated, Gaussian white noise or correlated noise was added to Simulated Data 1, having different levels of the signal-to-noise ratio (SNR), that is, 20, 25, 30 35, 40, 45, and 50 dB. The Gaussian white noise was generated using the awgn function in MATLAB. The correlated noise was obtained by low-pass filtering i.i.d. Gaussian noise.Simulated Data 2 was generated using nine randomly selected spectral signatures fromA1, with 100×100 pixels and 224 bands per pixel, provided by Dr. M. D. Iordache and Prof. J. M. Bioucas-Dias. The fractional abundances satisfy the ANC and the ASC and are piecewise smoothed. The true abundances of three endmembers are illustrated in Figure 2. After Simulated Data 2 was generated, Gaussian white noise or correlated noise with seven levels of SNR (20, 25, 30, 35, 40, 45, and 50 dB) was added to Simulated Data 2 as well.Figure 2
Results on Simulated Data 1 with 30 dB white noise: SRE as function of endmember number.In the two simulated data experiments, all the SGA methods have adopted the block-processing strategy to get better performances. The block sizes of the SGA methods were empirically set to 10, and the stopping threshold parameterσ of SOMP and SMP was set to 1·10−2. The stopping threshold parameters σ and δ of BSOMP were set to 1·10−2 and 0.1. RSFoBa was tested using different values of the parameters: λ and σ: 0 and 5·10−3 for Simulated Data 1 experiment; 1 and 1·10−2 for Simulated Data 2 experiment. SUnSAL-TV was tested using different values of the parameters λ and λTV: 1·10−2 and 1·10−2, 1·10−2 and 1·10−2, 1·10−2 and 5·10−3, 5·10−3 and 1·10−3, 1·10−3 and 5·10−4, 5·10−4 and 1·10−4, and 5·10−4 and 1·10−4 for SNR = 20, 25, 30, 35, 40, 45, and 50 dB. SUnSAL and CLSUnSAL was tested using different values of the parameter λ: 1·10−2, 5·10−3, 1·10−3, 5·10−4, 5·10−4, 5·10−4, and 5·10−4 for SNR = 20, 25, 30, 35, 40, 45, and 50 dB.Figure2 shows the SRE (dB) results obtained on Simulated Data 1 with white noise as a function of endmember number, using the different tested methods. The SREs of all the algorithms decrease as the endmember number increases. We can find that BSOMP, SOMP, SMP, RSFoBa, and SUnSAL-TV perform better than SUnSAL and CLSUnSAL. All the SGA methods have comparable performances with SUnSAL-TV; BSOMP can always obtain the best estimation of the abundances when the number of endmembers is less than 15. Figure 3 shows the SRE (dB) results obtained on Simulated Data 1 with white noise as a function of SNR when the endmember number is eleven. The SREs of all the algorithms decrease as the SNR decreases. We can find that all the SGA methods perform better than SUnSAL-TV, SUnSAL, and CLSUnSAL. This result indicates that all the SGA methods adopt the block-processing strategy to effectively pick up all the actual endmembers. Among all the SGA methods, BSOMP and RSFoBa behave better than SOMP and SMP, and BSOMP can always obtain the best estimation of the abundances. Figures 4 and 5 show the SRE results obtained on Simulated Data 1 with correlated noise as a function of endmember number and SNR, respectively. As shown in Figures 4 and 5, BSOMP gets an overall better performance than the other algorithms.Figure 3
Results on Simulated Data 1 with white noise when the endmember number is eleven: SRE as function of SNR.Figure 4
Results on Simulated Data 1 with 30 dB correlated noise: SRE as function of endmember number.Figure 5
Results on Simulated Data 1 with correlated noise when the endmember number is eleven: SRE as function of SNR.Figure6 shows the number of the potential endmembers obtained from the spectral library by all the SGA methods. The SGA methods adopt the block-processing strategy to effectively pick up all the actual endmembers. It can be observed that BSOMP and RSFoBa can select the actual endmembers more accurately than SOMP and SMP. It is worth mentioning that many endmembers in the estimated endmembers set of SOMP are nonexisting and affect the performance of SOMP. However, BSOMP incorporates a backtracking process to remove some endmembers chosen wrongly from the estimated endmembers set, which can select the actual endmembers more accurately than SOMP. Differing from RSFoBa, BSOMP uses a backtracking process in the whole hyperspectral data to remove some endmembers chosen wrongly from the estimated endmembers set in previous block-processing, which can identify the true endmembers set more accurately than RSFoBa. It indicates that BSOMP behaves better than the other SGA methods.Figure 6
Results obtained by the four SGAs on Simulated Data 1 with 30 dB white noise or correlated noise: number of retained endmembers as function of endmember number.Figures7 and 8 show the SRE results obtained on Simulated Data 2 with white noise and correlated noise as a function of SNR, respectively. We can find that BSOMP, SOMP, SMP, RSFoBa, and SUnSAL-TV perform better than SUnSAL and CLSUnSAL. This result indicates that all the SGA methods adopt the block-processing strategy to effectively pick up all the actual endmembers, and SUnSAL-TV combinates spatial-contextual coherence regularizer, which can get better performance over SUnSAL and CLSUnSAL. All the SGA methods have comparable performances with SUnSAL-TV; BSOMP and RSFoBa can always obtain the better estimation of the abundances. RSFoBa incorporates the spatial-contextual coherence regularizer to select the actual endmembers more accurately than SOMP and SMP, and BSOMP incorporates a backtracking process in the whole hyperspectral data to identify the true endmembers set more accurately than SOMP and SMP. It can be observed that BSOMP and RSFoBa can select the actual endmembers more accurately than SOMP and SMP.Figure 7
Results on Simulated Data 2 with white noise: SRE as function of SNR.Figure 8
Results on Simulated Data 2 with correlated noise: SRE as function of SNR.For visual comparison, Figure9 shows the true abundance maps and the abundance maps estimated by the different algorithms on Simulated Data 2 with 20 dB white noise. In Figure 9, the abundance maps of endmembers, which were obtained by SUnSAL-TV, contain fewer noise points. However, the abundance maps estimated by SUnSAL-TV may exhibit an oversmooth visual effect around some pixels; some of the details and edges are not well preserved. Compared with SUnSAL-TV, the estimated abundances by BSOMP may contain more noise points but preserve more edge regions and are closer to the true abundances. Moreover, it is much easier to tell the different materials apart. Compared with SOMP and SMP, the estimated abundances by BSOMP contain fewer noise points and preserve more edge regions, which are more accurate and have a better visual effect.Figure 9
Comparison of abundance maps on Synthetic Data 2 with 20 dB white noise. From left column to right column are the abundance maps corresponding to Endmember 1, Endmember 4, and Endmember 9, respectively. From top row to bottom row are true abundance maps or abundance maps obtained by SUnSAL, CLSUnSAL, SUnSAL-TV, SOMP, SMP, RSFoBa, and BSOMP, respectively.Table1 shows the processing time measured after applying the tested algorithms to the two simulated data sets with 30 dB white noise. The algorithms were implemented in MATLAB 2009a and executed on a desktop PC equipped with an Intel Core 2 Duo CPU (2.93 GHz) and 2 GB of RAM memory. It can be observed that the SUnSAL-TV algorithm is far more time consuming compared with the other algorithms. We can find that BSOMP incorporating a backtracking technique is more time consuming than SOMP.Table 1
Processing times(s) measured after applying the tested methods to the two considered simulated data sets with 30 dB white noise.
Datacube
SUnSAL
CLSUnSAL
SUnSAL-TV
SOMP
SMP
RSFoBa
BSOMP
Simulated Data 1 (11 endmembers)
67.67
78.41
363.27
2.083
21.289
3.938
6.418
Simulated Data 2
254.71
279.02
1271.12
9.618
35.444
28.277
106.524In order to make the differences among the BSOMP algorithm and the other sparse unmixing algorithms clearer, we display their comparison concerned with several properties in Table2.Table 2
Comparison of the sparse unmixing algorithms with whether the backward step, joint sparsity, and contextual information are considered.
Property
Backward step
Joint sparsity
Contextual information
SMP
·
SOMP
·
RSFoBa
·
·
·
BSOMP
·
·
SUnSAL
CLSUnSAL
·
SUnSAL-TV
·SUnSAL is a convex relaxation method for sparse unmixing which solves thel1 norm sparse regression problem. CLSUnSAL is convex relaxation method for solving the l2,1 optimization problem. The main difference between SUnSAL and CLSUnSAL is that the former employs pixelwise independent regressions, while the latter enforces joint sparsity among all the pixels. SUnSAL-TV is the special case of SUnSAL which incorporates the spatial-contextual coherence regularizer into the SUnSAL. As observed in Figures 7 and 8, SUnSAL-TV behaves better than SUnSAL and CLSUnSAL, which indicates that the spatial-contextual coherence regularizer can make SUnSAL-TV more effective.SOMP and SMP are the forward SGAs which select one or more potential endmembers from the spectral library in each iteration. RSFoBa is a forward-backward SGA which incorporates the backward step, spatial-contextual information, and joint sparsity into the greedy algorithm. BSOMP is a forward-backward SGA which incorporates the backward step and joint sparsity to the greedy algorithm. As observed in Figure6, the numbers of endmembers selected by the SGAs are still larger than the actual number. BSOMP and RSFoBa behave better than SMP and SOMP, which indicates that the backward step can make BSOMP and RSFoBa more accurate to select the actual endmembers. Differing from RSFoBa, which adopts a forward step and a backward step to pick some potential endmembers in each block, BSOMP adopts a forward step to pick some potential endmembers in each block and then adopts a backward step in the whole hyperspectral data to remove some endmembers chosen wrongly from the estimated endmembers set. Figure 10 shows the SRE (dB) results obtained on Simulated Data 1 with white noise as a function of block size. We can find that BSOMP performs better than the other SGAs. Figure 11 shows the number of the potential endmembers obtained from the spectral library by all the SGA methods. It can be observed that BSOMP can select the actual endmembers more accurately than RSFoBa. This result indicates that a backward step in the whole hyperspectral data can make BSOMP more effective.Figure 10
Results on Simulated Data 1 with white noise when the endmember number is eleven: SRE as function of block size.Figure 11
Results on Simulated Data 1 with 30 dB white noise when the endmember number is eleven: number of retained endmembers as function of block size.


### 4.2. Evaluation with Real Data

The hyperspectral data set used in our real data experiment is the well-known AVIRIS Cuprite data set (http://aviris.jpl.nasa.gov/html/aviris.freedata.html). The size of the hyperspectral scene is a 250×191 pixels subset, with 188 spectral bands (water absorption bands and low SNR bands are removed). The spectral library for this data set is the USGS spectral library (containing 498 pure signatures) with the corresponding bands removed. The minerals map (http://speclab.cr.usgs.gov/cuprite95.tgif.2.2um_map.gif) which was produced by a Tricorder 3.3 software (http://speclab.cr.usgs.gov/PAPERS/tetracorder/) is shown in Figure 12. Note that the Tricorder map is only suitable for the hyperspectral data collected in 1995; the AVIRIS Cuprite data was collected was in 1997. Thus we only use the mineral map as a reference indicator for qualitative analysis of the fractional abundance maps produced by different sparse unmixing methods.Figure 12
USGS map showing the distribution of different minerals in the Cuprite mining district in Nevada.In this section, for the real data experiment, the sparsity of solution and the reconstruction error of the hyperspectral image are used to evaluate the performances of the sparse unmixing algorithms, since the true abundance maps of the real data are not available. The sparsity is obtained by counting the average number of nonzero abundances estimated by each algorithm. The same as in [22], we define the fractional abundances larger than 0.001 as nonzero abundances to avoid counting negligible values. The reconstruction error is calculated by the original hyperspectral image and the one reconstructed by the actual endmembers and their corresponding abundances [32]. The root mean square error (RMSE) is used to evaluate the reconstruction error between the original hyperspectral image and the reconstructed one: (14)RMSE=1K∑j=1KYj-Y^j2,where Y denotes the original hyperspectral image and Y^ denotes the reconstruction hyperspectral image. If the RMSE of a sparse unmixing algorithm is smaller, it also reflects the unmixing performance of this algorithm.The regularization parameter used for CLSUnSAL was empirically set to 0.01, while the one corresponding to SUnSAL was set to 0.001. The parametersλ and λTV of SUnSAL-TV were set to 0.001 and 0.001. The block sizes of the SGA methods were empirically set to 10. The stopping threshold parameter σ of SOMP and SMP was set to 0.01. The stopping threshold parameters σ and δ of BSOMP were set to 0.01 and 0.1. The parameters λ and σ of RSFoBa were set to 1 and 0.01. Table 2 shows the sparsity of the abundance vectors and the reconstruction error obtained by the different algorithm. In Table 3, we can find that BSOMP, SOMP, SMP, and SUnSAL perform better than SUnSAL-TV in the sparsest solutions, although SUnSAL-TV has the lowest reconstruction errors. BSOMP can achieve nearly as few reconstruction errors as SUnSAL-TV but obtain much sparser solutions. In general, BSOMP can obtain the sparsest solutions and achieve the lower reconstruction errors, which indicates the effectiveness of the proposed algorithm.Table 3
Sparsity of the abundance vectors and the reconstruction errors.
Algorithms
SUnSAL
CLSUnSAL
SUnSAL-TV
SOMP
SMP
RSFoBa
BSOMP
Sparsity
17.5629
20.7037
20.472
13.684
15.102
12.702
11.810
RMSE
0.0051
0.0035
0.0028
0.0040
0.0034
0.0032
0.0033For visual comparison, Figure13 shows the estimated fractional abundance maps by applying the proposed BSOMP, SMP, RSFoBa, SOMP, SUnSAL, and SUnSAL-TV algorithms to the AVIRIS Cuprite data for three different minerals (Alunite, Buddingtonite, and Chalcedony). For comparison, the classification maps of three materials extracted from the Tricorder 3.3 software are also presented in Figure 13. It can be observed that the abundances estimated by BSOMP are comparable or higher in the regions assigned to the respective minerals in comparison to the other considered algorithms, with regard to the classification maps produced by the Tricorder 3.3 software. We also analyzed the sparsity of the abundances estimated by BSOMP; the average number of endmembers with abundances higher than 0.001 estimated by BSOMP is 11.810 (per pixel). This means that the algorithm uses a minimum number of members to explain the data. Overall, the qualitative results and visual results reported in this section indicate that the BSOMP algorithm is more effective for sparse unmixing than the other algorithms.Figure 13
Abundance maps estimated by SUnSAL, CLSUnSAL, SUnSAL-TV, SOMP, SMP, RSFoBa, and BSOMP and the distribution maps produced by Tricorder 3.3 software for the AVIRIS Cuprite scene. From top row to bottom row are the maps produced or estimated by Tricorder software, SUnSAL, CLSUnSAL, SUnSAL-TV, SOMP, SMP, RSFoBa, and BSOMP, respectively. From left column to right column are the maps corresponding to Alunite, Buddingtonite, and Chalcedony, respectively.


## 4.1. Evaluation with Simulated Data

In the simulated experiments, the United States Geological Survey (USGS) digital spectral library (splib06a) [31] is used to build the spectral library A. The reflectance values of 498 spectral signatures are measured for 224 spectral bands distributed uniformly in the interval of 0.4–2.5 μm (A∈R224×498). Fifteen spectral signatures are chosen from the spectral library to generate our simulated data. Figure 1 shows five spectral signatures used for all the simulated data experiments. The other ten spectral signatures that are not displayed in the figure include Monazite HS255.3B, Zoisite HS347.3B, Rhodochrosite HS67, Neodymium Oxide GDS34, Pigeonite HS199.3B, Meionite WS700.Hlsep, Spodumene HS210.3B, Labradorite HS17.3B, Grossular WS484, and Wollastonite HS348.3B.Figure 1
Five spectral signatures from the USGS used in our simulated data experiments. The title of each subimage denotes the mineral corresponding to the signature.In Simulated Data 1 experiment, we generated seven datacubes of50×50 pixels and 224 bands per pixel, each containing a different number of endmembers: 3, 5, 7, 9, 11, 13, and 15. In each simulated pixel, the fractional abundances of the endmembers were randomly generated following a Dirichlet distribution [23]. Note that there was no pure pixel in Simulated Data 1, and the fractional abundances of the endmembers were less than 0.8. After Simulated Data 1 was generated, Gaussian white noise or correlated noise was added to Simulated Data 1, having different levels of the signal-to-noise ratio (SNR), that is, 20, 25, 30 35, 40, 45, and 50 dB. The Gaussian white noise was generated using the awgn function in MATLAB. The correlated noise was obtained by low-pass filtering i.i.d. Gaussian noise.Simulated Data 2 was generated using nine randomly selected spectral signatures fromA1, with 100×100 pixels and 224 bands per pixel, provided by Dr. M. D. Iordache and Prof. J. M. Bioucas-Dias. The fractional abundances satisfy the ANC and the ASC and are piecewise smoothed. The true abundances of three endmembers are illustrated in Figure 2. After Simulated Data 2 was generated, Gaussian white noise or correlated noise with seven levels of SNR (20, 25, 30, 35, 40, 45, and 50 dB) was added to Simulated Data 2 as well.Figure 2
Results on Simulated Data 1 with 30 dB white noise: SRE as function of endmember number.In the two simulated data experiments, all the SGA methods have adopted the block-processing strategy to get better performances. The block sizes of the SGA methods were empirically set to 10, and the stopping threshold parameterσ of SOMP and SMP was set to 1·10−2. The stopping threshold parameters σ and δ of BSOMP were set to 1·10−2 and 0.1. RSFoBa was tested using different values of the parameters: λ and σ: 0 and 5·10−3 for Simulated Data 1 experiment; 1 and 1·10−2 for Simulated Data 2 experiment. SUnSAL-TV was tested using different values of the parameters λ and λTV: 1·10−2 and 1·10−2, 1·10−2 and 1·10−2, 1·10−2 and 5·10−3, 5·10−3 and 1·10−3, 1·10−3 and 5·10−4, 5·10−4 and 1·10−4, and 5·10−4 and 1·10−4 for SNR = 20, 25, 30, 35, 40, 45, and 50 dB. SUnSAL and CLSUnSAL was tested using different values of the parameter λ: 1·10−2, 5·10−3, 1·10−3, 5·10−4, 5·10−4, 5·10−4, and 5·10−4 for SNR = 20, 25, 30, 35, 40, 45, and 50 dB.Figure2 shows the SRE (dB) results obtained on Simulated Data 1 with white noise as a function of endmember number, using the different tested methods. The SREs of all the algorithms decrease as the endmember number increases. We can find that BSOMP, SOMP, SMP, RSFoBa, and SUnSAL-TV perform better than SUnSAL and CLSUnSAL. All the SGA methods have comparable performances with SUnSAL-TV; BSOMP can always obtain the best estimation of the abundances when the number of endmembers is less than 15. Figure 3 shows the SRE (dB) results obtained on Simulated Data 1 with white noise as a function of SNR when the endmember number is eleven. The SREs of all the algorithms decrease as the SNR decreases. We can find that all the SGA methods perform better than SUnSAL-TV, SUnSAL, and CLSUnSAL. This result indicates that all the SGA methods adopt the block-processing strategy to effectively pick up all the actual endmembers. Among all the SGA methods, BSOMP and RSFoBa behave better than SOMP and SMP, and BSOMP can always obtain the best estimation of the abundances. Figures 4 and 5 show the SRE results obtained on Simulated Data 1 with correlated noise as a function of endmember number and SNR, respectively. As shown in Figures 4 and 5, BSOMP gets an overall better performance than the other algorithms.Figure 3
Results on Simulated Data 1 with white noise when the endmember number is eleven: SRE as function of SNR.Figure 4
Results on Simulated Data 1 with 30 dB correlated noise: SRE as function of endmember number.Figure 5
Results on Simulated Data 1 with correlated noise when the endmember number is eleven: SRE as function of SNR.Figure6 shows the number of the potential endmembers obtained from the spectral library by all the SGA methods. The SGA methods adopt the block-processing strategy to effectively pick up all the actual endmembers. It can be observed that BSOMP and RSFoBa can select the actual endmembers more accurately than SOMP and SMP. It is worth mentioning that many endmembers in the estimated endmembers set of SOMP are nonexisting and affect the performance of SOMP. However, BSOMP incorporates a backtracking process to remove some endmembers chosen wrongly from the estimated endmembers set, which can select the actual endmembers more accurately than SOMP. Differing from RSFoBa, BSOMP uses a backtracking process in the whole hyperspectral data to remove some endmembers chosen wrongly from the estimated endmembers set in previous block-processing, which can identify the true endmembers set more accurately than RSFoBa. It indicates that BSOMP behaves better than the other SGA methods.Figure 6
Results obtained by the four SGAs on Simulated Data 1 with 30 dB white noise or correlated noise: number of retained endmembers as function of endmember number.Figures7 and 8 show the SRE results obtained on Simulated Data 2 with white noise and correlated noise as a function of SNR, respectively. We can find that BSOMP, SOMP, SMP, RSFoBa, and SUnSAL-TV perform better than SUnSAL and CLSUnSAL. This result indicates that all the SGA methods adopt the block-processing strategy to effectively pick up all the actual endmembers, and SUnSAL-TV combinates spatial-contextual coherence regularizer, which can get better performance over SUnSAL and CLSUnSAL. All the SGA methods have comparable performances with SUnSAL-TV; BSOMP and RSFoBa can always obtain the better estimation of the abundances. RSFoBa incorporates the spatial-contextual coherence regularizer to select the actual endmembers more accurately than SOMP and SMP, and BSOMP incorporates a backtracking process in the whole hyperspectral data to identify the true endmembers set more accurately than SOMP and SMP. It can be observed that BSOMP and RSFoBa can select the actual endmembers more accurately than SOMP and SMP.Figure 7
Results on Simulated Data 2 with white noise: SRE as function of SNR.Figure 8
Results on Simulated Data 2 with correlated noise: SRE as function of SNR.For visual comparison, Figure9 shows the true abundance maps and the abundance maps estimated by the different algorithms on Simulated Data 2 with 20 dB white noise. In Figure 9, the abundance maps of endmembers, which were obtained by SUnSAL-TV, contain fewer noise points. However, the abundance maps estimated by SUnSAL-TV may exhibit an oversmooth visual effect around some pixels; some of the details and edges are not well preserved. Compared with SUnSAL-TV, the estimated abundances by BSOMP may contain more noise points but preserve more edge regions and are closer to the true abundances. Moreover, it is much easier to tell the different materials apart. Compared with SOMP and SMP, the estimated abundances by BSOMP contain fewer noise points and preserve more edge regions, which are more accurate and have a better visual effect.Figure 9
Comparison of abundance maps on Synthetic Data 2 with 20 dB white noise. From left column to right column are the abundance maps corresponding to Endmember 1, Endmember 4, and Endmember 9, respectively. From top row to bottom row are true abundance maps or abundance maps obtained by SUnSAL, CLSUnSAL, SUnSAL-TV, SOMP, SMP, RSFoBa, and BSOMP, respectively.Table1 shows the processing time measured after applying the tested algorithms to the two simulated data sets with 30 dB white noise. The algorithms were implemented in MATLAB 2009a and executed on a desktop PC equipped with an Intel Core 2 Duo CPU (2.93 GHz) and 2 GB of RAM memory. It can be observed that the SUnSAL-TV algorithm is far more time consuming compared with the other algorithms. We can find that BSOMP incorporating a backtracking technique is more time consuming than SOMP.Table 1
Processing times(s) measured after applying the tested methods to the two considered simulated data sets with 30 dB white noise.
Datacube
SUnSAL
CLSUnSAL
SUnSAL-TV
SOMP
SMP
RSFoBa
BSOMP
Simulated Data 1 (11 endmembers)
67.67
78.41
363.27
2.083
21.289
3.938
6.418
Simulated Data 2
254.71
279.02
1271.12
9.618
35.444
28.277
106.524In order to make the differences among the BSOMP algorithm and the other sparse unmixing algorithms clearer, we display their comparison concerned with several properties in Table2.Table 2
Comparison of the sparse unmixing algorithms with whether the backward step, joint sparsity, and contextual information are considered.
Property
Backward step
Joint sparsity
Contextual information
SMP
·
SOMP
·
RSFoBa
·
·
·
BSOMP
·
·
SUnSAL
CLSUnSAL
·
SUnSAL-TV
·SUnSAL is a convex relaxation method for sparse unmixing which solves thel1 norm sparse regression problem. CLSUnSAL is convex relaxation method for solving the l2,1 optimization problem. The main difference between SUnSAL and CLSUnSAL is that the former employs pixelwise independent regressions, while the latter enforces joint sparsity among all the pixels. SUnSAL-TV is the special case of SUnSAL which incorporates the spatial-contextual coherence regularizer into the SUnSAL. As observed in Figures 7 and 8, SUnSAL-TV behaves better than SUnSAL and CLSUnSAL, which indicates that the spatial-contextual coherence regularizer can make SUnSAL-TV more effective.SOMP and SMP are the forward SGAs which select one or more potential endmembers from the spectral library in each iteration. RSFoBa is a forward-backward SGA which incorporates the backward step, spatial-contextual information, and joint sparsity into the greedy algorithm. BSOMP is a forward-backward SGA which incorporates the backward step and joint sparsity to the greedy algorithm. As observed in Figure6, the numbers of endmembers selected by the SGAs are still larger than the actual number. BSOMP and RSFoBa behave better than SMP and SOMP, which indicates that the backward step can make BSOMP and RSFoBa more accurate to select the actual endmembers. Differing from RSFoBa, which adopts a forward step and a backward step to pick some potential endmembers in each block, BSOMP adopts a forward step to pick some potential endmembers in each block and then adopts a backward step in the whole hyperspectral data to remove some endmembers chosen wrongly from the estimated endmembers set. Figure 10 shows the SRE (dB) results obtained on Simulated Data 1 with white noise as a function of block size. We can find that BSOMP performs better than the other SGAs. Figure 11 shows the number of the potential endmembers obtained from the spectral library by all the SGA methods. It can be observed that BSOMP can select the actual endmembers more accurately than RSFoBa. This result indicates that a backward step in the whole hyperspectral data can make BSOMP more effective.Figure 10
Results on Simulated Data 1 with white noise when the endmember number is eleven: SRE as function of block size.Figure 11
Results on Simulated Data 1 with 30 dB white noise when the endmember number is eleven: number of retained endmembers as function of block size.


## 4.2. Evaluation with Real Data

The hyperspectral data set used in our real data experiment is the well-known AVIRIS Cuprite data set (http://aviris.jpl.nasa.gov/html/aviris.freedata.html). The size of the hyperspectral scene is a 250×191 pixels subset, with 188 spectral bands (water absorption bands and low SNR bands are removed). The spectral library for this data set is the USGS spectral library (containing 498 pure signatures) with the corresponding bands removed. The minerals map (http://speclab.cr.usgs.gov/cuprite95.tgif.2.2um_map.gif) which was produced by a Tricorder 3.3 software (http://speclab.cr.usgs.gov/PAPERS/tetracorder/) is shown in Figure 12. Note that the Tricorder map is only suitable for the hyperspectral data collected in 1995; the AVIRIS Cuprite data was collected was in 1997. Thus we only use the mineral map as a reference indicator for qualitative analysis of the fractional abundance maps produced by different sparse unmixing methods.Figure 12
USGS map showing the distribution of different minerals in the Cuprite mining district in Nevada.In this section, for the real data experiment, the sparsity of solution and the reconstruction error of the hyperspectral image are used to evaluate the performances of the sparse unmixing algorithms, since the true abundance maps of the real data are not available. The sparsity is obtained by counting the average number of nonzero abundances estimated by each algorithm. The same as in [22], we define the fractional abundances larger than 0.001 as nonzero abundances to avoid counting negligible values. The reconstruction error is calculated by the original hyperspectral image and the one reconstructed by the actual endmembers and their corresponding abundances [32]. The root mean square error (RMSE) is used to evaluate the reconstruction error between the original hyperspectral image and the reconstructed one: (14)RMSE=1K∑j=1KYj-Y^j2,where Y denotes the original hyperspectral image and Y^ denotes the reconstruction hyperspectral image. If the RMSE of a sparse unmixing algorithm is smaller, it also reflects the unmixing performance of this algorithm.The regularization parameter used for CLSUnSAL was empirically set to 0.01, while the one corresponding to SUnSAL was set to 0.001. The parametersλ and λTV of SUnSAL-TV were set to 0.001 and 0.001. The block sizes of the SGA methods were empirically set to 10. The stopping threshold parameter σ of SOMP and SMP was set to 0.01. The stopping threshold parameters σ and δ of BSOMP were set to 0.01 and 0.1. The parameters λ and σ of RSFoBa were set to 1 and 0.01. Table 2 shows the sparsity of the abundance vectors and the reconstruction error obtained by the different algorithm. In Table 3, we can find that BSOMP, SOMP, SMP, and SUnSAL perform better than SUnSAL-TV in the sparsest solutions, although SUnSAL-TV has the lowest reconstruction errors. BSOMP can achieve nearly as few reconstruction errors as SUnSAL-TV but obtain much sparser solutions. In general, BSOMP can obtain the sparsest solutions and achieve the lower reconstruction errors, which indicates the effectiveness of the proposed algorithm.Table 3
Sparsity of the abundance vectors and the reconstruction errors.
Algorithms
SUnSAL
CLSUnSAL
SUnSAL-TV
SOMP
SMP
RSFoBa
BSOMP
Sparsity
17.5629
20.7037
20.472
13.684
15.102
12.702
11.810
RMSE
0.0051
0.0035
0.0028
0.0040
0.0034
0.0032
0.0033For visual comparison, Figure13 shows the estimated fractional abundance maps by applying the proposed BSOMP, SMP, RSFoBa, SOMP, SUnSAL, and SUnSAL-TV algorithms to the AVIRIS Cuprite data for three different minerals (Alunite, Buddingtonite, and Chalcedony). For comparison, the classification maps of three materials extracted from the Tricorder 3.3 software are also presented in Figure 13. It can be observed that the abundances estimated by BSOMP are comparable or higher in the regions assigned to the respective minerals in comparison to the other considered algorithms, with regard to the classification maps produced by the Tricorder 3.3 software. We also analyzed the sparsity of the abundances estimated by BSOMP; the average number of endmembers with abundances higher than 0.001 estimated by BSOMP is 11.810 (per pixel). This means that the algorithm uses a minimum number of members to explain the data. Overall, the qualitative results and visual results reported in this section indicate that the BSOMP algorithm is more effective for sparse unmixing than the other algorithms.Figure 13
Abundance maps estimated by SUnSAL, CLSUnSAL, SUnSAL-TV, SOMP, SMP, RSFoBa, and BSOMP and the distribution maps produced by Tricorder 3.3 software for the AVIRIS Cuprite scene. From top row to bottom row are the maps produced or estimated by Tricorder software, SUnSAL, CLSUnSAL, SUnSAL-TV, SOMP, SMP, RSFoBa, and BSOMP, respectively. From left column to right column are the maps corresponding to Alunite, Buddingtonite, and Chalcedony, respectively.


## 5. Conclusions

In this paper, we present a new SGA termed BSOMP for sparse unmixing of hyperspectral data. The main idea of the proposed algorithm is that BSOMP incorporates a backtrack step to detect the previous chosen endmembers’ reliability and then deletes the unreliable endmembers. It is proved that BSOMP can recover the optimal endmembers from the spectral library under certain conditions. Experiments on both synthetic data and real data also demonstrate that the BSOMP algorithm is more effective than the other SGAs for sparse unmixing.The topic we want to discuss is why BSOMP is effective for sparse unmixing. Endmember selection and abundance estimation are the two main parts of the SGAs. The number of the potential endmembers selected by the SGAs is far smaller than the number of spectral signatures in the spectral library, and then the abundance estimation method makes the abundance vectors sparse further with the relatively small endmembers subset. The main difference among BSOMP, SOMP, SMP, and RSFoBa is the endmember selection. The SGAs adopt a block-processing strategy to efficiently select endmembers, but the numbers of endmembers selected are still larger than the actual number. However the nonexisting endmembers will have a negative effect on the estimation of the abundances corresponding to the actual endmembers.Compared with SOMP and SMP, BSOMP and RSFoBa adopt a backward step to detect the previous chosen endmembers’ reliability and delete the nonexisting endmembers. So BSOMP and RSFoBa can select the actual endmembers more accurately than SOMP and SMP. The main difference between BSOMP and RSFoBa is that the former adopts a backward step in the whole hyperspectral data, while the latter adopts a backward step in each block. As shown in Figures10 and 11, we can observe that the block size will significantly affect the number of the endmembers obtained by RSFoBa, while BSOMP is not affected by the block size. In conclusion, BSOMP can select the actual endmembers more accurately than the other SGAs and be beneficial to improve estimation effectiveness.



---

*Source: 842017-2015-05-14.xml*