# Low-Rank and Spectral-Spatial Sparse Unmixing for Hyperspectral Remote Sensing Imagery

Authors: Fan Li
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Sparse unmixing is an important technique for hyperspectral data analysis. Most sparse unmixing algorithms underutilize the spatial and spectral information of the hyperspectral image, which is unfavourable for the accuracy of endmember identification and abundance estimation. We propose a new spectral unmixing method based on the low-rank representation model and spatial-weighted collaborative sparsity, aiming to exploit structural information in both the spatial and spectral domains for unmixing. The spatial weights are incorporated into the collaborative sparse regularization term to enhance the spatial continuity of the image. Meanwhile, the global low-rank constraint is employed to maintain the spatial low-dimensional structure of the image. The model is solved by the well-known alternating direction method of multiplier, in which the abundance coefficients and the spatial weights are updated iteratively in the inner and outer loops, respectively. Experimental results obtained from simulation and real data reveal the superior performance of the proposed algorithm on spectral unmixing.


---

## Body

## 1. Introduction

Hyperspectral remote sensing is considered to be a major breakthrough in remote sensing technology due to its high spectral resolutions and simultaneous acquisition of both images and spectra of objects [1, 2]. Hyperspectral remote sensing images provide rich earth observation information, which shows great application potential in mineral resource exploration, precision agriculture, and environmental monitoring [3]. However, the insufficient spatial resolution of hyperspectral images leads to the emergence of a large number of mixed pixels, which influences the fine interpretation of hyperspectral images [4]. Spectral unmixing is an effective way to deal with mixed pixels. It obtains pure spectral signals of ground objects (endmembers) in the imaging area and their corresponding proportions (abundance) [5].The observed spectrum of the mixed pixel is a mixture of endmember spectra, and unmixing is essentially an optimization problem for a given spectral mixture model. The linear mixture model (LMM) assumes that each pixel is formed by linearly combining spectral signatures of endmembers [3]. Many hyperspectral unmixing algorithms based on LMM have been proposed, among which sparse unmixing methodology has received extensive attention [6, 7]. Sparse unmixing is a semisupervised spectral unmixing method, which takes the spectral library as a priori and assumes that the spectra of the endmembers participating in the hyperspectral image can be found in the library [5]. Generally, the number of endmembers in each mixed pixel is very small relative to the number of spectral signatures in the library; that is, the corresponding abundance is sparse. Therefore, the sparse regression algorithm is adopted to select the optimal endmember subset from the spectral library and simultaneously estimate the fractional abundances [8].Sparse unmixing can be regarded as solving the combinatorial optimization problem by a constrained sparse linear regression algorithm, in which sparsity-inducing term promotes the sparsity of the abundance solutions. For example, the sparse unmixing via variable splitting and augmented Lagrangian (SUnSAL) uses theℓ1 norm to characterize the sparsity of each abundance vector [5]. Furthermore, the collaborative SUnSAL (CLSUnSAL) uses the ℓ2,1 norm to enforce the joint sparsity among all pixels [9].The iterative reweightedℓ1 norm minimization is confirmed to enhance the sparsity of solutions [10, 11]. In the work of [12], a weighting factor is introduced in the sparse regularization term to reduce the solution space. The double reweighted sparse unmixing algorithm (DRSU) employs the double weighted ℓ1 norm to simultaneously reduce the nonzero rows in the abundance matrix corresponding to the actual endmembers and the nonzero values in each abundance vector [13].Hyperspectral data presents a three-dimensional cube structure with 2-dimensional spatial information and 1-dimensional spectral information of the target area [14, 15]. Taking spatial information as a priori knowledge greatly improves the unmixing performance, so it has become one of the research hotspots in the field of hyperspectral unmixing [16–18]. For instance, the total variation (TV) regularizer is proposed to integrate the spatial-context information [19, 20], and the TV regularization term is added to the SUnSAL algorithm to construct the SUnSAL-TV algorithm [19]. Moreover, the local collaborative sparse unmixing algorithm (LCSU) focuses on the local collaborative characteristics of the image, which imposes different spatial weights on each local region of the abundance map [21].Low rankness is another inherent characteristic of the abundance matrix, which provides a new perspective for exploring the spatial structure of the data [22, 23]. The high spatial correlation of pixels in the image reflects in the high spectral similarity of these pixels. The similar pixels usually share the same endmembers and similar abundance maps, which means that the corresponding abundance vectors are linearly dependent; that is, the abundance matrix is low rank or approximately low rank. The low-rank constraint helps to obtain a low-rank approximation of the abundance matrix, which captures the global spatial data structure. The alternating direction sparse and low-rank unmixing (ADSpLRU) algorithm first attempts to combine sparse and low-rank constraints, which considers the spatial correlation of the pixels in a 3×3 sliding window [24]. The joint-sparse-blocks and low-rank unmixing (JSpBLRU) algorithm divides the abundance matrix into flexible regular blocks with one row and multiple columns and applies joint sparsity to the blocks, and the low-rank regularizer is merged into the block sparse model synchronously [25].In this paper, a low-rank and spectral-spatial sparse unmixing algorithm is proposed, which provides the low-rank constraint and the collaborative sparse constraint incorporating neighborhood spatial weights. The multiconstraint strategy effectively exploits the spatial and spectral information of the image and improves the accuracy of the abundance estimation.The rest of this paper is structured as follows. Section2 introduces the sparse unmixing model and its variants. Section 3 presents the low-rank and spectral-spatial sparse unmixing model and its optimization algorithm in detail. In Section 4, the effectiveness of the proposed algorithm is demonstrated by experiments on simulated and real hyperspectral data. Section 5 draws conclusions and suggests possible future research.


## 2. Related Works

Sparse unmixing replaces the endmember set used in LMM with a large spectral library and searches the library for the combination of spectral signatures that best represents the hyperspectral images [5]. Let Y∈ℝl×n be a hyperspectral image including n pixels with l bands. The sparse unmixing model can be expressed as
(1)Y=AX+N,where A∈ℝl×m is a spectral library including m spectral signatures, X∈ℝm×n is the abundance matrix corresponding to the spectral library A, and N∈ℝl×n is the noise or model error. The abundance matrix X reflects the proportion of spectral signatures from the library A in each pixel. According to the physical meaning of abundance, two constraints are considered: abundance nonnegativity constraint (ANC) and abundance sum-to-one constraint (ASC), which can be written as X≥0and 1TX=1.In fact, there are only a few entries from the spectral libraryA act on the mixed pixels, which results in many zeros in the relevant fractional abundances. That is, the abundance matrix X is sparse. Based on sparse and redundant representation theory, the unmixing problem can be modeled as follows:
(2)minX12Y−AXF2+λX0s.t.X≥0,where ·F is the Frobenius norm, λ is a nonnegative regularization parameter, and X0 indicates the number of nonzero values in X. Note that the ASC is not strictly conducted in the model (2), which is ascribed to the drawbacks mentioned in [5].The optimization problem related to model (2) is NP-hard. Under the condition of restricted isometric property (RIP), ℓ0 norm can be relaxed to ℓ1 norm [26]. The model (2) is modified as
(3)minX12Y−AXF2+λX1,1s.t.X≥0,where X1,1=∑j=1nxj1, xj is the jth column vector of X. The SUnSAL algorithm efficiently solves model (3) [19]. However, the ℓ1 norm measures the sparsity of each abundance vector independently instead of considering the joint sparsity of all pixels in the scene. The mixed norm ℓ2,1 norm is proposed to improve sparsity among the rows of X collaboratively. The corresponding optimization problem is described as
(4)minX12Y−AXF2+λX2,1s.t.X≥0,where X2,1=∑i=1mxi2, xi is the ith row vector of X. The optimization problem related to (4) is solved by the CLSUnSAL algorithm. As an improvement of the SUnSAL algorithm, it is aimed at finding a structured solution with only a few nonzero lines. In other words, it promotes the row sparsity of the abundance matrix X.The abovementioned unmixing models underutilize the spatial structure features of hyperspectral images. The rich spatial information contained in the image, which has an important impact on the unmixing performance, is not incorporated into these models. To this end, we propose a new sparse unmixing algorithm to address the aforementioned limitation.


## 3. The Proposed Low-Rank and Spectral-Spatial Sparse Unmixing Algorithm

### 3.1. The Proposed LRSSU Model

Inspired by the idea of iterative weighting [27], the spatial weight based on neighborhood is introduced into ℓ2,1 norm sparse regularization term, which integrates the spatial and spectral information. Due to the advantage of the low-rank representation model in maintaining the low-dimensional structure of the image, a low-rank regularization term is incorporated to explore the global spatial correlation of the hyperspectral data. Combining the sparse and low-rank constraints, a new sparse unmixing model with low-rank constraint is proposed as follows:
(5)minX12Y−AXF2+λHspa⊙X2,1+τrankXs.t.X≥0,where Hspa∈ℝm×n is the spatial weight matrix, ⊙ is the Hadamard product, rankX is the rank of the abundance matrix X, and λ and τ are nonnegative regularization parameters. Hspa is updated iteratively as follows:
(6)Hspa,ijt+1=1fNpxpt+ε,where Hspa,ijt represents the element in the ith row and jth column of Hspa at the tth iteration. p represents the pixel corresponding to fractional abundance xij, xp=x1j,x2j,⋯,xmjT is the abundance vector for pixel p, Np represents the neighbor pixel set centered on p in the image, and fN·· denotes a function that explores the spatial correlation of pixels in a neighborhood system N·. In this paper, the 8 neighbor pixels of p constitute the neighborhood system N8p; then, the function f· is defined as
(7)fN8pxp=∑q∈N8pρqxq∑q∈N8pρq,where q is any pixel in the 8 neighbor pixel set of p, xq is the abundance vector corresponding to q, ρq=1/dEp,q denotes the neighborhood weight with regard to q, and dEp,q denotes the Euclidean distance between pixel q and its central pixel p.The rank ofX is closely relevant to the similarity of abundance vectors, and the reduced rank abundance matrix consists of fewer linearly independent vectors accordingly. When ℓ2,1 norm is used to promote the row sparsity of X, the low-rank regularizer helps to determine the endmembers from the spectral library more accurately. Nevertheless, it is difficult to calculate the rank of the matrix directly. The nuclear norm X∗ is the tightest convex relaxation of the rank of X. Therefore, we replace rankX with X∗ and rewrite model (5) as
(8)minX12Y−AXF2+λHspa⊙X2,1+τX∗s.t.X≥0,where X∗=∑iσiX and σiX is the ith singular value of X.The weighted nuclear norm minimization (WNNM) can enhance the sparsity of the singular vector by imposing different weights on each singular value [28]. Specifically, smaller singular values that carry less valid information are assigned larger weights and then are greatly reduced or even reduced to zero. On the contrary, larger values are assigned smaller weights, so that the main information is well preserved. The weighted nuclear norm Xb,∗ is adopted to improve the robustness and consistency of model (8), which is expressed as
(9)Xb,∗=∑i=1rbiσiX,where b=b1,b2,⋯ is the nonnegative weight vector and bi is the weight related to the ith singular value of the matrix X. At the kth iteration, bi is updated as follows:
(10)bik=1σiXk+ε,where ε>0 is a small constant added to avoid singularity and Xk represents the abundance matrix X at the kth iteration.In summary, the low-rank and spectral-spatial sparse unmixing model (LRSSU) is expressed as follows:(11)minX12Y−AXF2+λHspa⊙X2,1+τXb,∗s.t.X≥0.


### 3.2. Solution of the Optimization Problem

In this section, an iterative optimization scheme of internal and external dual circulation based on the alternating direction method of multipliers (ADMM) [29] is proposed. The external circulation updates the spatial weight matrix following equation (6), and the internal circulation solves the abundance matrix in the framework of ADMM. The flow of solving model (11) is described in detail as follows.The auxiliary matrix variablesU, V1, V2, V3, and V4 are introdued, and model (11) is equivalently converted to the following form:
(12)minU,V1,V2,V3,V412Y−V1F2+λHspa⊙V22,1+τV3b,∗+ιR+V4s.t.U=X,V1=AU,V2=U,V3=U,V4=U,where ιR+X=∑i=1nιR+xi is the indicator function; i.e., ιR+xi is zero if xi is nonnegative and +∞ otherwise. Equation (12) can be expressed compactly as
(13)minU,VgU,Vs.t.GU+BV=0,where gU,V=1/2Y−V1F2+λHspa⊙V22,1+τV3b,∗+ιR+V4, G=A,I,I,IT, B=diag−I, and V=V1,V2,V3,V4T.The augmented Lagrangian of equation (13) is written as
(14)LU,V,D=gU,V+μ2GU+BV−DF2,where μ is a nonnegative penalty parameter, D=D1,D2,D3,D4T is the Lagrange multipliers with regards to GU+BV=0.The optimization problem (13) can be solved by minimizing LU,V,D in the framework of ADMM. We update the variables U and V sequentially by solving the following subproblems:
(15)Uk+1=argminULU,Vk,Dk,Vk+1=argminVLUk+1,V,Dk.TheU-subproblem has the closed-form solution, which is obtained by taking the partial derivative of U. The variable U is calculated as follows:
(16)Uk+1=argminUμ2AU−V1k−D1kF2+μ2U−V2k−D2kF2+μ2U−V3k−D3kF2+μ2U−V4k−D4kF2=ATA+3I−1ATV1k+D1k+V2k+D2k+V3k+D3k+V4k+D4k.TheV-subproblem is split into four independent optimization problems with regard to V1,V2,V3, and V4. The solution of V1 is calculated as follows:
(17)V1k+1=argminV112Y−V1F2+μ2AUk+1−V1−D1kF2=11+μY+μAUk+1−D1k.The optimization problem related toV2 is
(18)V2k+1=argminV2λHspa⊙V22,1+μ2Uk+1−V2−D2kF2.Therth row of V2 is written as V2,r; then, the closed-form solution of problem (18) is
(19)V2,rk+1=vect‐softUrk+1−D2,rk,λμHspa,r,where vect‐softy,τ=ymaxy2−τ,0/maxy2−τ,0+τ denotes the soft threshold function for row vector. V3 is updated according to the following optimization problem:
(20)V3k+1=argminV3τV3b,∗+μ2Uk+1−V3−D3kF2=SVTb,τ/μUk+1−D3k,where SVTb,ω· denotes the singular value threshold function and the weight vector b is calculated according to equation (10). The singular value decomposition (SVD) of X is written as X=U^Diagσ1,⋯,σrV∧T, then SVTb,ωX=U^Diagmaxσ1−ωb1,0,⋯,maxσr−ωbr,0V∧T. In the same vein, we get the solution of V4 as follows:
(21)V4k+1=argminV4ιR+V4+μ2Uk+1−V4−D4kF2=maxUk+1−D4k,0.Finally, the Lagrange multipliers are updated as follows:(22)Dk+1=Dk−GUk+1+BVk+1.To sum up, the complete scheme for the proposed LRSSU model is presented in Algorithm1, which consists of internal and external dual circulations. The external circulation updates the spatial weight matrix in step 4, which is set up to 100 iterations. The internal circulation updates the variables and Lagrange multipliers in ADMM framework from step 6 to step 15 and sets 5 iterations empirically [10, 30, 31].The convergence of Algorithm1 is difficult to prove. Nevertheless, the residual of the optimization problem decreases with the increase of iterations. Figure 1 shows the residual GUt+BVtF as a function of the number of iterations for the complete algorithm. It can be observed that the residual arrived at a certain level close to zero after about 20 iterations of the external circulation. Therefore, the stopping criterion of the algorithm is set as the residual error reaches the error tolerance or the number of iterations reaches the upper limit, which guarantees the convergence of the algorithm.Algorithm 1: Pseudocode of the LRSSU algorithm.
1:Initialization:2:k,t=0, choose μ,λ,τ,ε>0, U0, V10, V20, V30, V40, D10, D20, D30, D403:Repeat:4:Hspa,ijt+1=1/fNpxpt+ε, xp=x1j,x2j,⋯,xmjT is the jth column of Ut−V2t5:Repeat:6:Uk+1⟵ATA+3I−1ATV1k+D1k+V2k+D2k+V3k+D3k+V4k+D4k7:V1k+1⟵1/1+μY+μAUk+1−D1k8:V2,rk+1⟵vect‐softUrk+1−D2,rk,λ/μHspa,r9:V3k+1⟵SVTb,τ/μUk+1−D3k10:V4k+1⟵maxUk+1−D4k,011:Update Lagrange multipliers:12:D1k+1⟵D1k−AUk+1+V1k+113:D2k+1⟵D2k−Uk+1+V2k+114:D3k+1⟵D3k−Uk+1+V3k+115:D4k+1⟵D4k−Uk+1+V4k+116:Update iteration:k⟵k+117:Ut+1⟵Uk18:D2t+1⟵D2k19:Update iteration:t⟵t+120:Until some stopping criterion is satisfied.Figure 1
The residualGUt+BVtF as a function of the number of iterations for the complete algorithm.


## 3.1. The Proposed LRSSU Model

Inspired by the idea of iterative weighting [27], the spatial weight based on neighborhood is introduced into ℓ2,1 norm sparse regularization term, which integrates the spatial and spectral information. Due to the advantage of the low-rank representation model in maintaining the low-dimensional structure of the image, a low-rank regularization term is incorporated to explore the global spatial correlation of the hyperspectral data. Combining the sparse and low-rank constraints, a new sparse unmixing model with low-rank constraint is proposed as follows:
(5)minX12Y−AXF2+λHspa⊙X2,1+τrankXs.t.X≥0,where Hspa∈ℝm×n is the spatial weight matrix, ⊙ is the Hadamard product, rankX is the rank of the abundance matrix X, and λ and τ are nonnegative regularization parameters. Hspa is updated iteratively as follows:
(6)Hspa,ijt+1=1fNpxpt+ε,where Hspa,ijt represents the element in the ith row and jth column of Hspa at the tth iteration. p represents the pixel corresponding to fractional abundance xij, xp=x1j,x2j,⋯,xmjT is the abundance vector for pixel p, Np represents the neighbor pixel set centered on p in the image, and fN·· denotes a function that explores the spatial correlation of pixels in a neighborhood system N·. In this paper, the 8 neighbor pixels of p constitute the neighborhood system N8p; then, the function f· is defined as
(7)fN8pxp=∑q∈N8pρqxq∑q∈N8pρq,where q is any pixel in the 8 neighbor pixel set of p, xq is the abundance vector corresponding to q, ρq=1/dEp,q denotes the neighborhood weight with regard to q, and dEp,q denotes the Euclidean distance between pixel q and its central pixel p.The rank ofX is closely relevant to the similarity of abundance vectors, and the reduced rank abundance matrix consists of fewer linearly independent vectors accordingly. When ℓ2,1 norm is used to promote the row sparsity of X, the low-rank regularizer helps to determine the endmembers from the spectral library more accurately. Nevertheless, it is difficult to calculate the rank of the matrix directly. The nuclear norm X∗ is the tightest convex relaxation of the rank of X. Therefore, we replace rankX with X∗ and rewrite model (5) as
(8)minX12Y−AXF2+λHspa⊙X2,1+τX∗s.t.X≥0,where X∗=∑iσiX and σiX is the ith singular value of X.The weighted nuclear norm minimization (WNNM) can enhance the sparsity of the singular vector by imposing different weights on each singular value [28]. Specifically, smaller singular values that carry less valid information are assigned larger weights and then are greatly reduced or even reduced to zero. On the contrary, larger values are assigned smaller weights, so that the main information is well preserved. The weighted nuclear norm Xb,∗ is adopted to improve the robustness and consistency of model (8), which is expressed as
(9)Xb,∗=∑i=1rbiσiX,where b=b1,b2,⋯ is the nonnegative weight vector and bi is the weight related to the ith singular value of the matrix X. At the kth iteration, bi is updated as follows:
(10)bik=1σiXk+ε,where ε>0 is a small constant added to avoid singularity and Xk represents the abundance matrix X at the kth iteration.In summary, the low-rank and spectral-spatial sparse unmixing model (LRSSU) is expressed as follows:(11)minX12Y−AXF2+λHspa⊙X2,1+τXb,∗s.t.X≥0.


## 3.2. Solution of the Optimization Problem

In this section, an iterative optimization scheme of internal and external dual circulation based on the alternating direction method of multipliers (ADMM) [29] is proposed. The external circulation updates the spatial weight matrix following equation (6), and the internal circulation solves the abundance matrix in the framework of ADMM. The flow of solving model (11) is described in detail as follows.The auxiliary matrix variablesU, V1, V2, V3, and V4 are introdued, and model (11) is equivalently converted to the following form:
(12)minU,V1,V2,V3,V412Y−V1F2+λHspa⊙V22,1+τV3b,∗+ιR+V4s.t.U=X,V1=AU,V2=U,V3=U,V4=U,where ιR+X=∑i=1nιR+xi is the indicator function; i.e., ιR+xi is zero if xi is nonnegative and +∞ otherwise. Equation (12) can be expressed compactly as
(13)minU,VgU,Vs.t.GU+BV=0,where gU,V=1/2Y−V1F2+λHspa⊙V22,1+τV3b,∗+ιR+V4, G=A,I,I,IT, B=diag−I, and V=V1,V2,V3,V4T.The augmented Lagrangian of equation (13) is written as
(14)LU,V,D=gU,V+μ2GU+BV−DF2,where μ is a nonnegative penalty parameter, D=D1,D2,D3,D4T is the Lagrange multipliers with regards to GU+BV=0.The optimization problem (13) can be solved by minimizing LU,V,D in the framework of ADMM. We update the variables U and V sequentially by solving the following subproblems:
(15)Uk+1=argminULU,Vk,Dk,Vk+1=argminVLUk+1,V,Dk.TheU-subproblem has the closed-form solution, which is obtained by taking the partial derivative of U. The variable U is calculated as follows:
(16)Uk+1=argminUμ2AU−V1k−D1kF2+μ2U−V2k−D2kF2+μ2U−V3k−D3kF2+μ2U−V4k−D4kF2=ATA+3I−1ATV1k+D1k+V2k+D2k+V3k+D3k+V4k+D4k.TheV-subproblem is split into four independent optimization problems with regard to V1,V2,V3, and V4. The solution of V1 is calculated as follows:
(17)V1k+1=argminV112Y−V1F2+μ2AUk+1−V1−D1kF2=11+μY+μAUk+1−D1k.The optimization problem related toV2 is
(18)V2k+1=argminV2λHspa⊙V22,1+μ2Uk+1−V2−D2kF2.Therth row of V2 is written as V2,r; then, the closed-form solution of problem (18) is
(19)V2,rk+1=vect‐softUrk+1−D2,rk,λμHspa,r,where vect‐softy,τ=ymaxy2−τ,0/maxy2−τ,0+τ denotes the soft threshold function for row vector. V3 is updated according to the following optimization problem:
(20)V3k+1=argminV3τV3b,∗+μ2Uk+1−V3−D3kF2=SVTb,τ/μUk+1−D3k,where SVTb,ω· denotes the singular value threshold function and the weight vector b is calculated according to equation (10). The singular value decomposition (SVD) of X is written as X=U^Diagσ1,⋯,σrV∧T, then SVTb,ωX=U^Diagmaxσ1−ωb1,0,⋯,maxσr−ωbr,0V∧T. In the same vein, we get the solution of V4 as follows:
(21)V4k+1=argminV4ιR+V4+μ2Uk+1−V4−D4kF2=maxUk+1−D4k,0.Finally, the Lagrange multipliers are updated as follows:(22)Dk+1=Dk−GUk+1+BVk+1.To sum up, the complete scheme for the proposed LRSSU model is presented in Algorithm1, which consists of internal and external dual circulations. The external circulation updates the spatial weight matrix in step 4, which is set up to 100 iterations. The internal circulation updates the variables and Lagrange multipliers in ADMM framework from step 6 to step 15 and sets 5 iterations empirically [10, 30, 31].The convergence of Algorithm1 is difficult to prove. Nevertheless, the residual of the optimization problem decreases with the increase of iterations. Figure 1 shows the residual GUt+BVtF as a function of the number of iterations for the complete algorithm. It can be observed that the residual arrived at a certain level close to zero after about 20 iterations of the external circulation. Therefore, the stopping criterion of the algorithm is set as the residual error reaches the error tolerance or the number of iterations reaches the upper limit, which guarantees the convergence of the algorithm.Algorithm 1: Pseudocode of the LRSSU algorithm.
1:Initialization:2:k,t=0, choose μ,λ,τ,ε>0, U0, V10, V20, V30, V40, D10, D20, D30, D403:Repeat:4:Hspa,ijt+1=1/fNpxpt+ε, xp=x1j,x2j,⋯,xmjT is the jth column of Ut−V2t5:Repeat:6:Uk+1⟵ATA+3I−1ATV1k+D1k+V2k+D2k+V3k+D3k+V4k+D4k7:V1k+1⟵1/1+μY+μAUk+1−D1k8:V2,rk+1⟵vect‐softUrk+1−D2,rk,λ/μHspa,r9:V3k+1⟵SVTb,τ/μUk+1−D3k10:V4k+1⟵maxUk+1−D4k,011:Update Lagrange multipliers:12:D1k+1⟵D1k−AUk+1+V1k+113:D2k+1⟵D2k−Uk+1+V2k+114:D3k+1⟵D3k−Uk+1+V3k+115:D4k+1⟵D4k−Uk+1+V4k+116:Update iteration:k⟵k+117:Ut+1⟵Uk18:D2t+1⟵D2k19:Update iteration:t⟵t+120:Until some stopping criterion is satisfied.Figure 1
The residualGUt+BVtF as a function of the number of iterations for the complete algorithm.


## 4. Experimental Analysis

In this section, the performance of the proposed LRSSU algorithm will be discussed in simulated and real hyperspectral data experiments. In the simulated data experiment, the signal-to-reconstruction error (SRE) is used to quantitatively evaluate the unmixing accuracy. The SRE (dB) is defined as follows:(23)SREdB=10·log10Ex22/Ex−x^22,where E· is the expectation function, x^ is the estimated abundance vector, x is the true abundance vector. In addition, the probability of success ps is another quantitative evaluation indicator, which appraises the stability of the results [5] . ps is defined as follows: ps≡Px∧−x2/x2≤threshold. The results will be considered successful, if the relative error power of results is less than an appropriate threshold, usually 3.16 (5 dB). The larger the SRE (dB) or ps, the better the performance of the algorithm. For fair comparison, the parameters of all algorithms in the experiment are tuned to the optimal in advance, and the maximum number of iterations is set to 500. The proposed LRSSU algorithm is compared with the classic SUnSAL, CLSUnSAL, and SUnSAL-TV algorithms, as well as the advanced DRSU algorithm and the JSPBLRU algorithm considering the low-rank constraint.


### 4.1. Experiment on Simulated Data

The spectral libraryA∈ℝ221×222 used in the simulated experiment is composed of 222 spectral curves randomly selected from the splib06 spectral library provided by the United States Geological Survey (USGS) (http://speclab.cr.usgs.gov/spectral.lib06.), with 221 spectral bands. The reflectance values are evenly distributed in the range of 0.4-2.5 μm. The abundance map containing 100×100 pixels is generated form the HyperMix standard dataset, where ANC and ASC are implemented. The data imitates the distribution features of surface objects in nature, such as surface water, coastline, vegetation, and mountains, so it is widely used to test the unmixing performance of different algorithms. The simulated data cube is synthesized by nine endmembers selected from the library A, and then, white Gaussian noise with different signal-to-noise ratios (SNR) is added. The true fractional abundances of the nine endmembers are displayed in Figure 2.Figure 2
True fractional abundances of the endmembers in the simulated data.Table1 lists the value of SRE (dB) and ps obtained by different algorithms on the simulated data at SNR=30 dB, 40 dB, and 50 dB. From Table 1, we observe that the proposed LRSSU algorithm achieves higher SRE (dB) and ps than other comparison algorithms at all noise levels. Compared with SUnSAL, CLSUnSAL, and SUnSAL-TV, the LRSSU algorithm shows great advantages in stability and accuracy, which proves its superior unmixing performance. In the case of high noise, the LRSSU and JSPBLRU algorithms maintain high ps, which indicates that minimizing the rank of the abundance matrix can effectively improve the robustness of the unmixing algorithm. The SRE (dB) obtained by the LRSSU algorithm is higher than that of the DRSU and JSpBLRU algorithms, which demonstrates that the local spatial weighting strategy plays an important role in improving the accuracy of unmixing.Table 1
SRE (B) andps results after applying different unmixing algorithms to the simulated data (the optimal parameter settings for which the reported values were achieved are shown in the parentheses).
AlgorithmSNR=30 dBSNR=40 dBSNR=50 dBSRE (dB)psSRE (dB)psSRE (dB)psSUnSAL6.39670.630311.55620.887118.92220.9991(λ=8e-3)(λ=2e-3)(λ=3e-4)CLSUnSAL6.61570.726914.78190.999926.68271(λ=3e-1)(λ=2e-2)(λ=7e-3)SUnSAL-TV9.03780.782215.44330.986125.36561(λ=4e-3; λTV=2e-3)(λ=6e-5; λTV=9e-4)(λ=5e-5; λTV=9e-5)DRSU14.32960.946526.0656134.49671(λ=2e-3)(λ=6e-4)(λ=1e-4)JSpBLRU15.05460.992824.8702134.00731(λ=2.3; τ=1.4)(λ=4e-1; τ=2e-1)(λ=2e-1; τ=4e-2)LRSSU19.83800.998927.9738136.47911(λ=4e-2; τ=8e-2)(λ=2e-2; τ=1e-2)(λ=4e-3; τ=2e-3)To illustrate the effect of the proposed LRSSU algorithm on endmember identification, Figure3 exhibits estimated abundances obtained by each algorithm on the simulated data with SNR=30 dB. For visual clarity, Figure 3 displays only 1000 pixels selected from the image. Each colored line in the figure represents a nonzero row vector in the abundance matrix, corresponding to an actual endmember. It can be seen from the figure that the abundance map obtained by the LRSSU algorithm is consistent with the ground-truth abundance. More specifically, the endmembers identified by the LRSSU algorithm from the spectral library are the same as the true ones, with the least outliers. The abundance maps estimated by the SUnSAL, CLSUnSAL, and SUnSAL-TV algorithms emerge a mass of false endmember lines and interference noises, which do not match the actual endmembers well. The DRSU algorithm confuses some actual endmembers with other highly similar entries in the spectral library, resulting in endmember mismatches, and it is susceptible to noise. The results obtained by the JSpBLRU algorithm are disturbed by abundance vectors with small values corresponding to false endmembers. It can be concluded that the proposed LRSSU algorithm further enhances the row sparsity of the abundance matrix and improves the accuracy of identifying endmembers from the spectral library.Figure 3
The ground-truth abundance and the estimated abundance maps obtained from the simulated data under noise ratio ofSNR=30 dB.Figure4 exhibits the abundance maps of the ninth endmember estimated by different algorithms from the simulated data at SNR=30 dB. The other endmembers exhibit similar behavior, so they are not shown. For a more intuitive comparison, Figure 4 also shows the difference maps between the estimated abundance obtained by each algorithm and the true abundance. The abundance map obtained by the proposed LRSSU algorithm is close to the true distribution, which reveals a significant advantage than other algorithms. The SUnSAL, CLSUnSAL, and SUnSAL-TV algorithms estimate the abundances inaccurately, and the unmixing results are obviously different from the true abundance map. Compared with the DRSU and JSpBLRU algorithms, the LRSSU algorithm retains more spatial detail information of the abundance map, showing better spatial smoothness and less noise. It is proved that the combination of spectral-spatial collaborative sparse and low-rank constraints enhances the sparsity of the solution, preserves the finer structure and texture information of the image, promotes the spatial consistency between pixels, and improves the accuracy of unmixing.Figure 4
Abundance maps obtained for the ninth endmember from the simulated data under noise ratio ofSNR=30 dB. Difference maps between the ground-truth abundance and the results obtained by (a) SUnSAL, (b) CLSUnSAL, (c) SUnSAL-TV, (d) DRSU, (e) JSpBLRU, and (f) LRSSU.
(a)(b)(c)(d)(e)(f)The setting of regularization parameters has an important impact on the unmixing performance of the algorithm. To analyze the sensitivity of parameters in the proposed LRSSU algorithm, Figure5 shows the SRE (dB) as functions of parameters λ and τ. We observe that when the noise level decreases, the value of the optimal parameter decreases accordingly. Moreover, the parameter λ affects the accuracy of solutions more than the parameter τ. Overall, when the regularization parameters λ and τ vary in a large range, SRE (dB) is still optimal or suboptimal. This is helpful to achieve the best performance of the algorithm and improve its practical application.Figure 5
SRE (B) as functions of parametersλ and τ for the simulated data at SNR=30 dB.


### 4.2. Experiment on Real Data

The Cuprite mining dataset collected by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS) (http://aviris.jpl.nasa.gov/html/aviris.freedata.html.) is used in the real data experiment to evaluate the unmixing performance of the proposed LRSSU algorithm. The part utilized in the experiment is a subset of the scene, with 100×100 pixels, including 224 bands in the range of 0.4-2.5 μm, with a spectral resolution of 10 nm. The experiment excludes the low SNR and low water absorption bands; that is, the 1-2, 105-115, 150-170, and 223-224 bands are removed, leaving 188 spectral bands. The spectral library utilized in the experiment is selected from the USGS library, which contains 240 spectral signatures representing different minerals. Low SNR and low water absorption bands are removed from the library in the same way, and 188 bands are left. Since the true abundance maps of the Cuprite dataset are not available, we will refer to the mineral classification map (as shown in Figure 6) generated by Tricorder software (http://speclab.cr.usgs.gov/PAPER/tetracorder.) to qualitatively analyze the performance of different unmixing algorithms [32]. The experiment meets the maximum number of iterations to ensure the convergence of each unmixing algorithm.Figure 6
USGS mineral map of Cuprite mining district in Nevada [5].The regularization parameters involved in SUnSAL, CLSUnSAL, and DRSU algorithms are empirically set toλ=0.001, λ=0.01, and λ=0.0001, and the parameters related to SUnSAL-TV, JSpBLRU, and LRSSU algorithms are empirically set to λ=0.001,λTV=0.001 and λ=0.05,τ=0.02, and λ=0.08,τ=0.02, respectively. Figure 7 takes three typical minerals: alunite, buddingtonite, and chalcedony, as examples to show the abundance maps obtained by unmixing the Cuprite dataset with the above algorithms. From Figure 7, we can observe that the abundance maps obtained by the SUnSAL and CLSUnSAL algorithms are noisy, especially in buddingtonite. The abundance maps obtained by the SUnSAL-TV algorithm are obviously oversmooth, such as alunite and buddingtonite. The abundance maps obtained by the DRSU algorithm have poor spatial consistency, such as chalcedony. In addition, the abundance maps estimated by the LRSSU algorithm are less disturbed by noise than the JSPBLRU algorithm that also adopts the low-rank constraint, such as buddingtonite, and exhibit better spatial consistency, such as chalcedony. The abundance maps of the proposed LRSSU algorithm are closest to the Tricorder reference maps among all comparison algorithms. The qualitative analysis for the real data experiment confirms that the LRSSU algorithm is effective in unmixing real hyperspectral data and gets better unmixing results.Figure 7
Fractional abundance maps estimated by different unmixing algorithms for the Cuprite subscene.


## 4.1. Experiment on Simulated Data

The spectral libraryA∈ℝ221×222 used in the simulated experiment is composed of 222 spectral curves randomly selected from the splib06 spectral library provided by the United States Geological Survey (USGS) (http://speclab.cr.usgs.gov/spectral.lib06.), with 221 spectral bands. The reflectance values are evenly distributed in the range of 0.4-2.5 μm. The abundance map containing 100×100 pixels is generated form the HyperMix standard dataset, where ANC and ASC are implemented. The data imitates the distribution features of surface objects in nature, such as surface water, coastline, vegetation, and mountains, so it is widely used to test the unmixing performance of different algorithms. The simulated data cube is synthesized by nine endmembers selected from the library A, and then, white Gaussian noise with different signal-to-noise ratios (SNR) is added. The true fractional abundances of the nine endmembers are displayed in Figure 2.Figure 2
True fractional abundances of the endmembers in the simulated data.Table1 lists the value of SRE (dB) and ps obtained by different algorithms on the simulated data at SNR=30 dB, 40 dB, and 50 dB. From Table 1, we observe that the proposed LRSSU algorithm achieves higher SRE (dB) and ps than other comparison algorithms at all noise levels. Compared with SUnSAL, CLSUnSAL, and SUnSAL-TV, the LRSSU algorithm shows great advantages in stability and accuracy, which proves its superior unmixing performance. In the case of high noise, the LRSSU and JSPBLRU algorithms maintain high ps, which indicates that minimizing the rank of the abundance matrix can effectively improve the robustness of the unmixing algorithm. The SRE (dB) obtained by the LRSSU algorithm is higher than that of the DRSU and JSpBLRU algorithms, which demonstrates that the local spatial weighting strategy plays an important role in improving the accuracy of unmixing.Table 1
SRE (B) andps results after applying different unmixing algorithms to the simulated data (the optimal parameter settings for which the reported values were achieved are shown in the parentheses).
AlgorithmSNR=30 dBSNR=40 dBSNR=50 dBSRE (dB)psSRE (dB)psSRE (dB)psSUnSAL6.39670.630311.55620.887118.92220.9991(λ=8e-3)(λ=2e-3)(λ=3e-4)CLSUnSAL6.61570.726914.78190.999926.68271(λ=3e-1)(λ=2e-2)(λ=7e-3)SUnSAL-TV9.03780.782215.44330.986125.36561(λ=4e-3; λTV=2e-3)(λ=6e-5; λTV=9e-4)(λ=5e-5; λTV=9e-5)DRSU14.32960.946526.0656134.49671(λ=2e-3)(λ=6e-4)(λ=1e-4)JSpBLRU15.05460.992824.8702134.00731(λ=2.3; τ=1.4)(λ=4e-1; τ=2e-1)(λ=2e-1; τ=4e-2)LRSSU19.83800.998927.9738136.47911(λ=4e-2; τ=8e-2)(λ=2e-2; τ=1e-2)(λ=4e-3; τ=2e-3)To illustrate the effect of the proposed LRSSU algorithm on endmember identification, Figure3 exhibits estimated abundances obtained by each algorithm on the simulated data with SNR=30 dB. For visual clarity, Figure 3 displays only 1000 pixels selected from the image. Each colored line in the figure represents a nonzero row vector in the abundance matrix, corresponding to an actual endmember. It can be seen from the figure that the abundance map obtained by the LRSSU algorithm is consistent with the ground-truth abundance. More specifically, the endmembers identified by the LRSSU algorithm from the spectral library are the same as the true ones, with the least outliers. The abundance maps estimated by the SUnSAL, CLSUnSAL, and SUnSAL-TV algorithms emerge a mass of false endmember lines and interference noises, which do not match the actual endmembers well. The DRSU algorithm confuses some actual endmembers with other highly similar entries in the spectral library, resulting in endmember mismatches, and it is susceptible to noise. The results obtained by the JSpBLRU algorithm are disturbed by abundance vectors with small values corresponding to false endmembers. It can be concluded that the proposed LRSSU algorithm further enhances the row sparsity of the abundance matrix and improves the accuracy of identifying endmembers from the spectral library.Figure 3
The ground-truth abundance and the estimated abundance maps obtained from the simulated data under noise ratio ofSNR=30 dB.Figure4 exhibits the abundance maps of the ninth endmember estimated by different algorithms from the simulated data at SNR=30 dB. The other endmembers exhibit similar behavior, so they are not shown. For a more intuitive comparison, Figure 4 also shows the difference maps between the estimated abundance obtained by each algorithm and the true abundance. The abundance map obtained by the proposed LRSSU algorithm is close to the true distribution, which reveals a significant advantage than other algorithms. The SUnSAL, CLSUnSAL, and SUnSAL-TV algorithms estimate the abundances inaccurately, and the unmixing results are obviously different from the true abundance map. Compared with the DRSU and JSpBLRU algorithms, the LRSSU algorithm retains more spatial detail information of the abundance map, showing better spatial smoothness and less noise. It is proved that the combination of spectral-spatial collaborative sparse and low-rank constraints enhances the sparsity of the solution, preserves the finer structure and texture information of the image, promotes the spatial consistency between pixels, and improves the accuracy of unmixing.Figure 4
Abundance maps obtained for the ninth endmember from the simulated data under noise ratio ofSNR=30 dB. Difference maps between the ground-truth abundance and the results obtained by (a) SUnSAL, (b) CLSUnSAL, (c) SUnSAL-TV, (d) DRSU, (e) JSpBLRU, and (f) LRSSU.
(a)(b)(c)(d)(e)(f)The setting of regularization parameters has an important impact on the unmixing performance of the algorithm. To analyze the sensitivity of parameters in the proposed LRSSU algorithm, Figure5 shows the SRE (dB) as functions of parameters λ and τ. We observe that when the noise level decreases, the value of the optimal parameter decreases accordingly. Moreover, the parameter λ affects the accuracy of solutions more than the parameter τ. Overall, when the regularization parameters λ and τ vary in a large range, SRE (dB) is still optimal or suboptimal. This is helpful to achieve the best performance of the algorithm and improve its practical application.Figure 5
SRE (B) as functions of parametersλ and τ for the simulated data at SNR=30 dB.


## 4.2. Experiment on Real Data

The Cuprite mining dataset collected by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS) (http://aviris.jpl.nasa.gov/html/aviris.freedata.html.) is used in the real data experiment to evaluate the unmixing performance of the proposed LRSSU algorithm. The part utilized in the experiment is a subset of the scene, with 100×100 pixels, including 224 bands in the range of 0.4-2.5 μm, with a spectral resolution of 10 nm. The experiment excludes the low SNR and low water absorption bands; that is, the 1-2, 105-115, 150-170, and 223-224 bands are removed, leaving 188 spectral bands. The spectral library utilized in the experiment is selected from the USGS library, which contains 240 spectral signatures representing different minerals. Low SNR and low water absorption bands are removed from the library in the same way, and 188 bands are left. Since the true abundance maps of the Cuprite dataset are not available, we will refer to the mineral classification map (as shown in Figure 6) generated by Tricorder software (http://speclab.cr.usgs.gov/PAPER/tetracorder.) to qualitatively analyze the performance of different unmixing algorithms [32]. The experiment meets the maximum number of iterations to ensure the convergence of each unmixing algorithm.Figure 6
USGS mineral map of Cuprite mining district in Nevada [5].The regularization parameters involved in SUnSAL, CLSUnSAL, and DRSU algorithms are empirically set toλ=0.001, λ=0.01, and λ=0.0001, and the parameters related to SUnSAL-TV, JSpBLRU, and LRSSU algorithms are empirically set to λ=0.001,λTV=0.001 and λ=0.05,τ=0.02, and λ=0.08,τ=0.02, respectively. Figure 7 takes three typical minerals: alunite, buddingtonite, and chalcedony, as examples to show the abundance maps obtained by unmixing the Cuprite dataset with the above algorithms. From Figure 7, we can observe that the abundance maps obtained by the SUnSAL and CLSUnSAL algorithms are noisy, especially in buddingtonite. The abundance maps obtained by the SUnSAL-TV algorithm are obviously oversmooth, such as alunite and buddingtonite. The abundance maps obtained by the DRSU algorithm have poor spatial consistency, such as chalcedony. In addition, the abundance maps estimated by the LRSSU algorithm are less disturbed by noise than the JSPBLRU algorithm that also adopts the low-rank constraint, such as buddingtonite, and exhibit better spatial consistency, such as chalcedony. The abundance maps of the proposed LRSSU algorithm are closest to the Tricorder reference maps among all comparison algorithms. The qualitative analysis for the real data experiment confirms that the LRSSU algorithm is effective in unmixing real hyperspectral data and gets better unmixing results.Figure 7
Fractional abundance maps estimated by different unmixing algorithms for the Cuprite subscene.


## 5. Conclusions and Future Work

In this paper, we propose a new sparse reduced-rank regression model for hyperspectral unmixing, which merges spatial and spectral information into the regularizer to enhance the interpretation of hyperspectral data. A local spatial weighting factor is introduced in the collaborative sparse unmixing model, which promotes the spatial smoothness of the image. Moreover, a low-rank constraint is imposed on the abundance estimation to enhance the row sparsity of the abundance matrix, which to some extent overcomes the influence of the high correlation of spectral signatures in the library. The proposed LRSSU model is solved by an ADMM-based the inner and outer loop scheme, the inner loop updates fractional abundances, and the outer loop calculates the spatial weights. The scheme accelerates the convergence of the algorithm. Simulated and real hyperspectral data experiments reveal that the proposed LRSSU algorithm enhances the identification of actual endmembers in the spectral library. Compared with SUnSAL, CLSUnSAL, SUnSAL-TV, DRSU, and JSPBLRU algorithms, the proposed LRSSU algorithm improves the accuracy of abundance estimation and shows the advantages of the spectral-spatial multiple constraint strategy. In the future, we will extend the constraint strategy to blind unmixing and further explore the spectral-spatial structure information of the hyperspectral data cube represented by a third-order tensor.



---

*Source: 9374908-2021-10-31.xml*