# Double Regression-Based Sparse Unmixing for Hyperspectral Images

Authors: Shuaiyang Zhang; Wenshen Hua; Gang Li; Jie Liu; Fuyu Huang; Qianghui Wang
Journal: Journal of Sensors
 (2021)

---

## Abstract

Sparse unmixing has attracted widespread attention from researchers, and many effective unmixing algorithms have been proposed in recent years. However, most algorithms improve the unmixing accuracy at the cost of large calculations. Higher unmixing accuracy often leads to higher computational complexity. To solve this problem, we propose a novel double regression-based sparse unmixing model (DRSUM), which can obtain better unmixing results with lower computational complexity. DRSUM decomposes the complex objective function into two simple formulas and completes the unmixing process through two sparse regressions. The unmixing result of the first sparse regression is added as a constraint to the second. DRSUM is an open model, and we can add different constraints to improve the unmixing accuracy. In addition, we can perform appropriate preprocessing to further improve the unmixing results. Under this model, a specific algorithm called double regression-based sparse unmixing viaK-means (DRSUMK−means) is proposed. The improved K-means clustering algorithm is first used for preprocessing, and then we impose single sparsity and joint sparsity (using l2,0 norm to control the sparsity) constraints on the first and second sparse unmixing, respectively. To meet the sparsity requirement, we introduce the row-hard-threshold function to solve the l2,0 norm directly. Then, DRSUMK−means can be efficiently solved under alternating direction method of multipliers (ADMM) framework. Simulated and real data experiments have proven the effectiveness of DRSUMK−means.


---

## Body

## 1. Introduction

The hyperspectral remote sensing technology has shown great development and application potential in recent years [1–5]. However, due to spatial resolution and spatial complexity limitations, there are often many mixed pixels in the images, which seriously affect the processing and application of hyperspectral data. Therefore, spectral unmixing technology for mixed pixels has received widespread attention, which aims to identify the spectral features (endmembers) and estimate the corresponding abundance [6]. The linear unmixing model [7], which assumes that mixed pixels are linear combinations of endmembers, has been widely used for its simplicity. Under this model, unmixing algorithms are mainly categorized into three: statistic based, geometry based, and sparse regression based.Sparse unmixing [8] is a semisupervised unmixing algorithm, which assumes that the mixed pixels are linear combinations of only a few endmembers from the spectral library. It does not need to extract endmembers and estimate the number of endmembers. Numerous effective sparse unmixing algorithms have been proposed in recent years. The sparse unmixing algorithm via variable splitting and augmented Lagrangian (SUnSAL) [8] is a classic sparse unmixing algorithm. SUnSAL introduces new ideas into spectral unmixing. But the unmixing performance of SUnSAL is affected by the mutual correlation of spectral libraries, and the sparsity degree is insufficient. Collaborative SUnSAL (CLSUnSAL) [9] imposes joint sparsity constraint on hyperspectral images, which uses the l2,0 norm to control sparsity. However, to avoid the nondeterministic polynomial hard (NP-hard) problem, CLSUnSAL uses convex relaxation strategy to solve the l2,0 norm, which makes the sparsity insufficient in many scenarios. Reweighted sparse unmixing algorithm adopts the reweighting strategy to enhance the sparsity [10]. Collaborative sparse unmixing algorithm using l0 norm (CSUnL0) [11] solves the l2,0 norm directly and gets better sparsity than many convex relaxation algorithms.In addition to sparsity, how to make full use of spatial information in the hyperspectral images is also an important direction to improve unmixing performance. Sparse unmixing via variable splitting augmented Lagrangian and total variation (SUnSAL-TV) [12] introduces a total variation regularizer to utilize the spatial information. Spectral–spatial weighted sparse unmixing (S2WSU) [13] introduces spatial information by spatial weighting factor. Reference [14] combines total variation regularizer and joint-sparse-blocks to improve unmixing accuracy. Reference [15] proposed row-sparsity spectral unmixing via total variation (RSSUn-TV) and achieved promising results. Reference [16] uses the l1−l2 sparsity model and total variation regularizer to depict the sparse characteristic. The high spatial correlation of hyperspectral images is associated with the low column rank, and the low-rank constraint can make full use of the global structure of data. Therefore, the low-rank property can also be used to exploit spatial information. The alternating direction sparse and low-rank unmixing algorithm (ADSpLRU) [17] imposes the low-rank constraint on single sparsity model. Reference [18] incorporates the low-rank constraint into sparse unmixing to suppress noises and preserve more details. Reference [19] proposed nonconvex joint-sparsity and low-rank unmixing with dictionary pruning. Reference [20] combines the joint-sparse-blocks with low-rank constraint. Reference [21] proposed nonlocal low-rank prior sparse unmixing algorithm and achieved promising results.The above algorithms can all obtain promising results. However, higher unmixing accuracy often leads to higher computational complexity. To solve this problem, we propose a novel double regression-based sparse unmixing model (DRSUM). DRSUM decomposes the complex unmixing process into two simple sparse regressions, which can reduce the computational complexity. An important contribution of this paper is that the proposed DRSUM is an open model and framework. Therefore, we can add different constraints on the first and second sparse regressions to improve unmixing accuracy. In addition, we can perform appropriate preprocessing to further improve the unmixing results. Under this model, a specific algorithm called double regression-based sparse unmixing viaK-means (DRSUMK−means) is proposed. We first use the improved K-means clustering algorithm [22] for preprocessing, which can classify pixels into different categories. Pixels belonging to the same category have the same active endmember set. Therefore, we use the mean vector of each category instead of all pixels for the first sparse unmixing. In other words, the number of pixels participating in the first sparse unmixing is the number of categories, which can greatly reduce the computational complexity. Then, we impose single sparsity constraint to unmix these mean vectors. The unmixing results of mean vectors are considered as abundance estimate of all pixels in the corresponding categories. After the first unmixing, we take the result as a constraint and impose joint sparsity (using l2,0 norm to control the sparsity) constraint on the original hyperspectral image for the second unmixing. However, it is NP-hard to solve the l2,0 norm directly [23]. Most algorithms adopt the convex relaxation strategy to solve the l2,0 norm at the expense of unmixing accuracy and sparsity. To address this issue, we introduce the row-hard-threshold function to solve the l2,0 norm directly, which can improve the unmixing accuracy and sparsity. Simulated and real data experiments have proven that the proposed DRSUMK−means algorithm can obtain better unmixing results.


## 2. Sparse Unmixing

Sparse unmixing is aimed at finding an optimal linear combination of endmembers from the spectral library for mixed pixels. LetY∈RL×n represent a hyperspectral image, which contains n pixels and L spectral bands. The mixed model can be written as
(1)Y=AX+N,where A∈RL×m represents the known spectral library, which contains a total of m atoms (endmembers). X∈Rm×n denotes the abundance matrix, which represents the abundance of these m endmembers in the hyperspectral image. N∈RL×n is the error term. However, the number of endmembers in the real hyperspectral image is far less than that of the endmembers in the spectral library. Most endmembers of the spectral library do not appear in the hyperspectral image. Therefore, most of the elements in the abundance matrix X are zero. In other words, there are only a few elements in X that are nonzero [24]. Therefore, the abundance matrix X is sparse, and we can impose sparsity constraint to solve equation (1). In addition, the abundance nonnegative constraint (ANC) and the abundance sum-to-one constraint (ASC) should also be met. It has been proven in [8] that the ANC will automatically impose a generalized ASC. Therefore, we only impose ANC and relax the ASC. Then, the sparse unmixing problem can be written as:
(2)minX12Y−AXF2+λX1,0s.t.X≻¯0,where ⋅F represents the Frobenius norm, and λ denotes the regularization parameter. X1,0 represents the number of nonzero elements in X. ≻¯ denotes element-wise comparison. It is well known that solving equation (2) directly is NP-hard. The SUnSAL algorithm relaxes the X1,0 norm to X1,1 norm to solve it. SUnSAL is single-pixel sparse unmixing algorithm, and it is seriously affected by the mutual correlation of spectral libraries. CLSUnSAL imposes collaborative sparsity constraint to improve unmixing accuracy, and the unmixing problem is
(3)minX12Y−AXF2+λX2,1s.t.X≻¯0,where X2,1=∑i=1mXir2 and Xir represent the i-th row of X. It is worth noting that X2,1 is a convex relaxation of X2,0 to avoid NP-hard problem. However, this will reduce the sparsity constraint and unmixing accuracy. Although SUnSAL and CLSUnSAL are faster than many algorithms, their unmixing accuracy is insufficient. To obtain higher unmixing accuracy with low computational complexity, the double regression-based sparse unmixing model is proposed.


## 3. Double Regression-Based Sparse Unmixing Model

### 3.1. DRSUM

The core idea of DRSUM is to decompose the complex unmixing process into two simple sparse regressions. We can add different constraints on the first and second sparse unmixing to improve unmixing accuracy. In addition, we can perform appropriate preprocessing to further improve the unmixing results. The first sparse unmixing problem is as follows:(4)minX112Y−AX1F2+gX1s.t.X1≻¯0,where Y denotes the hyperspectral image, and X1 is the corresponding abundance matrix. gX1 denotes the constraint item. After the first sparse unmixing, the obtained abundance matrix X1 is added as a constraint to the second unmixing process. The second sparse unmixing problem is
(5)minX12Y−AXF2+α2X1−XF2+fXs.t.X≻¯0,,where α denotes the regularization parameter, and fX denotes the constraint for X. There are no specific restrictions on gX1 and fX. In different application scenarios, gX1 and fX can be different types of constraints to guarantee the unmixing accuracy, such as sparsity constraints, low-rank constraints, graph regularization constraints, total variation regularizer, nonlocal similarity, and image edge information. Therefore, DRSUM is flexible and has wide applicability.


### 3.2.DRSUMK−means

#### 3.2.1. Preprocessing

DRSUM is an open model and framework. Under this model, we can perform appropriate preprocessing and add different constraints to get different algorithms. In this paper, we first use the improvedK-means clustering algorithm to classify the hyperspectral image. The K-means algorithm is simple and fast. It has only one parameter to be determined: the number of categories K. It should be noted that the original K-means algorithm classifies pixels only based on spectral distance, which ignores spatial information in the hyperspectral image. To improve classification accuracy, we consider both spectral distance d1 and spatial distance d2 to classify pixels. Therefore, the distance calculation formula for pixels i and j becomes
(6)d1i,j=yi−yj22,(7)d2i,j=ai−aj2+bi−bj2,(8)d1′i,j=d1i,jmaxd1,d2′i,j=d2i,jmaxd2,(9)Di,j=d1′i,j2+ρd2′i,j2,where yi and yj denote the spectral vectors. ai,bi and aj,bj denote the position of pixels i and j, respectively. To effectively control the weight, d1i,j and d2i,j are normalized to obtain d1′i,j and d2′i,j. maxd1 and maxd2 denote the maximum value of spectral distance and spatial distance for any two pixels. ρ denotes the weight between spectral distance and spatial distance. To simplify the calculation, ρ is set to 1 in this paper. In each iteration, every pixel is assigned to the nearest cluster center. After each round of classification is completed, the cluster centers are updated by the average of the spectral vectors and position coordinates of all pixels in the corresponding categories. The classification process is completed when the cluster centers no longer change. More details about the K-means algorithm can be seen in reference [22].After preprocessing, the pixels are classified intoK categories. Pixels belonging to the same category have the same active endmember set. Therefore, we can use the mean vector of each category instead of all pixels for the first unmixing, which can greatly reduce the computational complexity. In other words, the number of pixels participating in the first sparse unmixing is the number of categories. In addition, the improved K-means clustering algorithm considers both spectral and spatial information for classification, which can improve the information utilization of hyperspectral images.


#### 3.2.2. First Sparse Unmixing

After preprocessing, we calculate the mean vector of each category and impose single sparsity constraint to unmix these mean vectors. Then, we can get the following unmixing problem:(10)minXM12YM−AXMF2+λ1XM1,1s.t.XM≻¯0,where YM=y1,y2,⋯,yK∈RL×K denotes mean vectors of the K categories. XM denotes the corresponding abundance matrix and XM1,1=∑i=1m∑j=1Kxi,j. λ1 denotes the regularization parameter. Before solving equation (10), we introduce the indicator function to denote the nonnegative constraint, and the unmixing problem becomes
(11)minXM12YM−AXMF2+λ1XM1,1+lR+XM,where lR+XM represents the indicator function, i.e., lR+XM=0 if XM∈R+ and lR+XM=+∞ otherwise. Then, equation (11) can be efficiently solved under alternating direction method of multipliers (ADMM) framework [25]. Set V1=AXM, V2=XM, V3=XM, and then equation (11) can be written as
(12)minXM,V1,V2,V312YM−V1F2+λ1V21,1+lR+V3s.t.V1=AXM,V2=XM,V3=XM.Set(13)gXM,V=12YM−V1F2+λ1V21,1+lR+V3,(14)V≡V1,V2,V3,(15)B=diag−I,(16)G=A,I,IT.Then, equation (12) becomes
(17)minXM,VgXM,Vs.t.GXM+BV=0.The following Lagrangian function can be constructed by equation (17)
(18)LXM,V,D≡gXM,V+μ2GXM+BV−DF2,where μ>0 denotes a positive constant, and D/μ denotes the Lagrangian multiplier of the constraint GXM+BV=0. Under the ADMM framework, each variable is independent when updated. First, expand equation (18) to get
(19)LXM,V1,V2,V3,D1,D2,D3=12YM−V1F2+λ1V21,1+lR+V3+μ2AXM−V1−D1F2+μ2XM−V2−D2F2+μ2XM−V3−D3F2.When updating each variable, only the items related to it need to be considered. The objective optimization functions of variablesXM, V1, V2, and V3 are
(20)XMi+1=argminXMμ2AXM−V1i−D1iF2+μ2XM−V2i−D2iF2+μ2XM−V3i−D3iF2,(21)V1i+1=argminV112YM−V1F2+μ2AXMi+1−V1−D1iF2,(22)V2i+1=argminV2λ1V21,1+μ2XMi+1−V2−D2iF2,(23)V3i+1=argminV3lR+V3+μ2XMi+1−V3−D3iF2.The corresponding solutions are(24)XMi+1=ATA+2I−1ATV1i+D1i+V2i+D2i+V3i+D3i,(25)V1i+1=11+μYM+μAXMi+1−D1i,(26)V2i+1=softXMi+1−D2i,λ1μ,(27)V3i+1=maxXMi+1−D3i,0,where soft denotes the soft-threshold function softy,τ=signymaxy−τ,0. Then, the final algorithm flow for the first sparse unmixing can be seen in Algorithm 1.Algorithm 1: ADMM pseudocode for the first sparse unmixing.
1. Input: mean vectors matrixYM, spectral library A, constant μ.2. Select parameters:λ1.3. Initialization:XM0,V0,D0 and set i=0.4. Repeat:XMi+1=ATA+2I−1ATV1i+D1i+V2i+D2i+V3i+D3iV1i+1=1/1+μYM+μAXMi+1−D1iV2i+1=softXMi+1−D2i,λ1/μV3i+1=maxXMi+1−D3i,0D1i+1=D1i−AXMi+1+V1i+1D2i+1=D2i−XMi+1+V2i+1D3i+1=D3i−XMi+1+V3i+1i=i+15. Until the Stopping Condition Is Met


#### 3.2.3. Second Sparse Unmixing

After the first sparse unmixing, we get the abundance estimateXM of these mean vectors. Then, the abundance matrix XM∈Rm×K for these mean vectors is converted to the abundance matrix X1∈Rm×n for the entire image. Specifically, the pixels in the same category have the same abundance estimate as the category’s mean vector. Then, X1 is added as a constraint for the second sparse unmixing, and we impose joint sparsity constraint on the original hyperspectral image. It is well known that solving the l2,0 norm directly will cause NP-hard problem. Most algorithms use the l2,1 norm instead of l2,0 norm to avoid NP-hard problem. However, the l2,1 norm just represents the sum of row-norms, which cannot meet the sparsity requirement sometimes. To address this issue, we introduce the row-hard-threshold function to solve the l2,0 norm directly. Then, we can get the following unmixing problem:
(28)minX12Y−AXF2+α2X1−XF2+λ2X2,0s.t.X≻¯0,where λ2 and α represent the regularization parameters, and the l2,0 norm is defined as follows
(29)X2,0=∑i=1m1Xir2>0,where Xir denotes the i-th row of X. We also use the ADMM method to solve equation (28). Set V1=AX, V2=X, and V3=X. Then equation (28) can be written as
(30)minX,V1,V2,V312Y−V1F2+α2X1−XF2+λ2V22,0+lR+V3s.t.V1=AX,V2=X,V3=X.We can construct the following Lagrangian function using the same method as equation (19)
(31)LX,V1,V2,V3,D1,D2,D3=12Y−V1F2+α2X1−XF2+λ2V22,0+lR+V3+μ2AX−V1−D1F2+μ2X−V2−D2F2+μ2X−V3−D3F2.The iterative update formulas of variablesX and V1 are as follows
(32)Xi+1=ATA+α+μI−1ATV1i+D1i+V2i+D2i+V3i+D3i+αX1,(33)V1i+1=11+μY+μAXi+1−D1i.Before solving variableV2, we first remove items independent of it and get
(34)V2i+1=argminV2λ2V22,0+μ2Xi+1−V2−D2iF2.Inspired by CSUnL0 algorithm [11], we introduce the row-hard-threshold function to solve equation (34). Define
(35)RHtXi,j=xi,j∑jxi,j2>t,0∑jxi,j2≤t,where RHtXi,j represents the j-th element in the i-th row of the row-hard-threshold function RHtX. xi,j denotes the j-th element in the i-th row of X and t denotes the threshold. It has been proved in reference [11] that H=RHtX is the closed-form solution for
(36)minHtH2,0+X−HF2.Equation (34) can be written as follows
(37)V2i+1=argminV2μ22λ2μV22,0+Xi+1−V2−D2iF2.Then, we can get the iterative update formula ofV2 from equation (36) and (37).
(38)V2i+1=RH2λ2/μXi+1−D2i.The iterative update formulas of other variables are similar to those in the first sparse unmixing. Then, the final algorithm flow for the second sparse unmixing can be seen in Algorithm2.Algorithm 2: ADMM pseudocode for the second sparse unmixing.
1. Input: hyperspectral imageY, spectral library A, abundance matrix X1, constant μ.2. Select parameters:λ2 and α.3. Initialization:X0,V0,D0 and set i=0.4. Repeat:Xi+1=ATA+α+μI−1ATV1i+D1i+V2i+D2i+V3i+D3i+αX1V1i+1=1/1+μY+μAXi+1−D1iV2i+1=RH2λ2/μXi+1−D2iV3i+1=maxXi+1−D3i,0D1i+1=D1i−AXi+1+V1i+1D2i+1=D2i−Xi+1+V2i+1D3i+1=D3i−Xi+1+V3i+1i=i+15:until the stopping condition is met.


## 3.1. DRSUM

The core idea of DRSUM is to decompose the complex unmixing process into two simple sparse regressions. We can add different constraints on the first and second sparse unmixing to improve unmixing accuracy. In addition, we can perform appropriate preprocessing to further improve the unmixing results. The first sparse unmixing problem is as follows:(4)minX112Y−AX1F2+gX1s.t.X1≻¯0,where Y denotes the hyperspectral image, and X1 is the corresponding abundance matrix. gX1 denotes the constraint item. After the first sparse unmixing, the obtained abundance matrix X1 is added as a constraint to the second unmixing process. The second sparse unmixing problem is
(5)minX12Y−AXF2+α2X1−XF2+fXs.t.X≻¯0,,where α denotes the regularization parameter, and fX denotes the constraint for X. There are no specific restrictions on gX1 and fX. In different application scenarios, gX1 and fX can be different types of constraints to guarantee the unmixing accuracy, such as sparsity constraints, low-rank constraints, graph regularization constraints, total variation regularizer, nonlocal similarity, and image edge information. Therefore, DRSUM is flexible and has wide applicability.


## 3.2.DRSUMK−means

### 3.2.1. Preprocessing

DRSUM is an open model and framework. Under this model, we can perform appropriate preprocessing and add different constraints to get different algorithms. In this paper, we first use the improvedK-means clustering algorithm to classify the hyperspectral image. The K-means algorithm is simple and fast. It has only one parameter to be determined: the number of categories K. It should be noted that the original K-means algorithm classifies pixels only based on spectral distance, which ignores spatial information in the hyperspectral image. To improve classification accuracy, we consider both spectral distance d1 and spatial distance d2 to classify pixels. Therefore, the distance calculation formula for pixels i and j becomes
(6)d1i,j=yi−yj22,(7)d2i,j=ai−aj2+bi−bj2,(8)d1′i,j=d1i,jmaxd1,d2′i,j=d2i,jmaxd2,(9)Di,j=d1′i,j2+ρd2′i,j2,where yi and yj denote the spectral vectors. ai,bi and aj,bj denote the position of pixels i and j, respectively. To effectively control the weight, d1i,j and d2i,j are normalized to obtain d1′i,j and d2′i,j. maxd1 and maxd2 denote the maximum value of spectral distance and spatial distance for any two pixels. ρ denotes the weight between spectral distance and spatial distance. To simplify the calculation, ρ is set to 1 in this paper. In each iteration, every pixel is assigned to the nearest cluster center. After each round of classification is completed, the cluster centers are updated by the average of the spectral vectors and position coordinates of all pixels in the corresponding categories. The classification process is completed when the cluster centers no longer change. More details about the K-means algorithm can be seen in reference [22].After preprocessing, the pixels are classified intoK categories. Pixels belonging to the same category have the same active endmember set. Therefore, we can use the mean vector of each category instead of all pixels for the first unmixing, which can greatly reduce the computational complexity. In other words, the number of pixels participating in the first sparse unmixing is the number of categories. In addition, the improved K-means clustering algorithm considers both spectral and spatial information for classification, which can improve the information utilization of hyperspectral images.


### 3.2.2. First Sparse Unmixing

After preprocessing, we calculate the mean vector of each category and impose single sparsity constraint to unmix these mean vectors. Then, we can get the following unmixing problem:(10)minXM12YM−AXMF2+λ1XM1,1s.t.XM≻¯0,where YM=y1,y2,⋯,yK∈RL×K denotes mean vectors of the K categories. XM denotes the corresponding abundance matrix and XM1,1=∑i=1m∑j=1Kxi,j. λ1 denotes the regularization parameter. Before solving equation (10), we introduce the indicator function to denote the nonnegative constraint, and the unmixing problem becomes
(11)minXM12YM−AXMF2+λ1XM1,1+lR+XM,where lR+XM represents the indicator function, i.e., lR+XM=0 if XM∈R+ and lR+XM=+∞ otherwise. Then, equation (11) can be efficiently solved under alternating direction method of multipliers (ADMM) framework [25]. Set V1=AXM, V2=XM, V3=XM, and then equation (11) can be written as
(12)minXM,V1,V2,V312YM−V1F2+λ1V21,1+lR+V3s.t.V1=AXM,V2=XM,V3=XM.Set(13)gXM,V=12YM−V1F2+λ1V21,1+lR+V3,(14)V≡V1,V2,V3,(15)B=diag−I,(16)G=A,I,IT.Then, equation (12) becomes
(17)minXM,VgXM,Vs.t.GXM+BV=0.The following Lagrangian function can be constructed by equation (17)
(18)LXM,V,D≡gXM,V+μ2GXM+BV−DF2,where μ>0 denotes a positive constant, and D/μ denotes the Lagrangian multiplier of the constraint GXM+BV=0. Under the ADMM framework, each variable is independent when updated. First, expand equation (18) to get
(19)LXM,V1,V2,V3,D1,D2,D3=12YM−V1F2+λ1V21,1+lR+V3+μ2AXM−V1−D1F2+μ2XM−V2−D2F2+μ2XM−V3−D3F2.When updating each variable, only the items related to it need to be considered. The objective optimization functions of variablesXM, V1, V2, and V3 are
(20)XMi+1=argminXMμ2AXM−V1i−D1iF2+μ2XM−V2i−D2iF2+μ2XM−V3i−D3iF2,(21)V1i+1=argminV112YM−V1F2+μ2AXMi+1−V1−D1iF2,(22)V2i+1=argminV2λ1V21,1+μ2XMi+1−V2−D2iF2,(23)V3i+1=argminV3lR+V3+μ2XMi+1−V3−D3iF2.The corresponding solutions are(24)XMi+1=ATA+2I−1ATV1i+D1i+V2i+D2i+V3i+D3i,(25)V1i+1=11+μYM+μAXMi+1−D1i,(26)V2i+1=softXMi+1−D2i,λ1μ,(27)V3i+1=maxXMi+1−D3i,0,where soft denotes the soft-threshold function softy,τ=signymaxy−τ,0. Then, the final algorithm flow for the first sparse unmixing can be seen in Algorithm 1.Algorithm 1: ADMM pseudocode for the first sparse unmixing.
1. Input: mean vectors matrixYM, spectral library A, constant μ.2. Select parameters:λ1.3. Initialization:XM0,V0,D0 and set i=0.4. Repeat:XMi+1=ATA+2I−1ATV1i+D1i+V2i+D2i+V3i+D3iV1i+1=1/1+μYM+μAXMi+1−D1iV2i+1=softXMi+1−D2i,λ1/μV3i+1=maxXMi+1−D3i,0D1i+1=D1i−AXMi+1+V1i+1D2i+1=D2i−XMi+1+V2i+1D3i+1=D3i−XMi+1+V3i+1i=i+15. Until the Stopping Condition Is Met


### 3.2.3. Second Sparse Unmixing

After the first sparse unmixing, we get the abundance estimateXM of these mean vectors. Then, the abundance matrix XM∈Rm×K for these mean vectors is converted to the abundance matrix X1∈Rm×n for the entire image. Specifically, the pixels in the same category have the same abundance estimate as the category’s mean vector. Then, X1 is added as a constraint for the second sparse unmixing, and we impose joint sparsity constraint on the original hyperspectral image. It is well known that solving the l2,0 norm directly will cause NP-hard problem. Most algorithms use the l2,1 norm instead of l2,0 norm to avoid NP-hard problem. However, the l2,1 norm just represents the sum of row-norms, which cannot meet the sparsity requirement sometimes. To address this issue, we introduce the row-hard-threshold function to solve the l2,0 norm directly. Then, we can get the following unmixing problem:
(28)minX12Y−AXF2+α2X1−XF2+λ2X2,0s.t.X≻¯0,where λ2 and α represent the regularization parameters, and the l2,0 norm is defined as follows
(29)X2,0=∑i=1m1Xir2>0,where Xir denotes the i-th row of X. We also use the ADMM method to solve equation (28). Set V1=AX, V2=X, and V3=X. Then equation (28) can be written as
(30)minX,V1,V2,V312Y−V1F2+α2X1−XF2+λ2V22,0+lR+V3s.t.V1=AX,V2=X,V3=X.We can construct the following Lagrangian function using the same method as equation (19)
(31)LX,V1,V2,V3,D1,D2,D3=12Y−V1F2+α2X1−XF2+λ2V22,0+lR+V3+μ2AX−V1−D1F2+μ2X−V2−D2F2+μ2X−V3−D3F2.The iterative update formulas of variablesX and V1 are as follows
(32)Xi+1=ATA+α+μI−1ATV1i+D1i+V2i+D2i+V3i+D3i+αX1,(33)V1i+1=11+μY+μAXi+1−D1i.Before solving variableV2, we first remove items independent of it and get
(34)V2i+1=argminV2λ2V22,0+μ2Xi+1−V2−D2iF2.Inspired by CSUnL0 algorithm [11], we introduce the row-hard-threshold function to solve equation (34). Define
(35)RHtXi,j=xi,j∑jxi,j2>t,0∑jxi,j2≤t,where RHtXi,j represents the j-th element in the i-th row of the row-hard-threshold function RHtX. xi,j denotes the j-th element in the i-th row of X and t denotes the threshold. It has been proved in reference [11] that H=RHtX is the closed-form solution for
(36)minHtH2,0+X−HF2.Equation (34) can be written as follows
(37)V2i+1=argminV2μ22λ2μV22,0+Xi+1−V2−D2iF2.Then, we can get the iterative update formula ofV2 from equation (36) and (37).
(38)V2i+1=RH2λ2/μXi+1−D2i.The iterative update formulas of other variables are similar to those in the first sparse unmixing. Then, the final algorithm flow for the second sparse unmixing can be seen in Algorithm2.Algorithm 2: ADMM pseudocode for the second sparse unmixing.
1. Input: hyperspectral imageY, spectral library A, abundance matrix X1, constant μ.2. Select parameters:λ2 and α.3. Initialization:X0,V0,D0 and set i=0.4. Repeat:Xi+1=ATA+α+μI−1ATV1i+D1i+V2i+D2i+V3i+D3i+αX1V1i+1=1/1+μY+μAXi+1−D1iV2i+1=RH2λ2/μXi+1−D2iV3i+1=maxXi+1−D3i,0D1i+1=D1i−AXi+1+V1i+1D2i+1=D2i−Xi+1+V2i+1D3i+1=D3i−Xi+1+V3i+1i=i+15:until the stopping condition is met.


## 3.2.1. Preprocessing

DRSUM is an open model and framework. Under this model, we can perform appropriate preprocessing and add different constraints to get different algorithms. In this paper, we first use the improvedK-means clustering algorithm to classify the hyperspectral image. The K-means algorithm is simple and fast. It has only one parameter to be determined: the number of categories K. It should be noted that the original K-means algorithm classifies pixels only based on spectral distance, which ignores spatial information in the hyperspectral image. To improve classification accuracy, we consider both spectral distance d1 and spatial distance d2 to classify pixels. Therefore, the distance calculation formula for pixels i and j becomes
(6)d1i,j=yi−yj22,(7)d2i,j=ai−aj2+bi−bj2,(8)d1′i,j=d1i,jmaxd1,d2′i,j=d2i,jmaxd2,(9)Di,j=d1′i,j2+ρd2′i,j2,where yi and yj denote the spectral vectors. ai,bi and aj,bj denote the position of pixels i and j, respectively. To effectively control the weight, d1i,j and d2i,j are normalized to obtain d1′i,j and d2′i,j. maxd1 and maxd2 denote the maximum value of spectral distance and spatial distance for any two pixels. ρ denotes the weight between spectral distance and spatial distance. To simplify the calculation, ρ is set to 1 in this paper. In each iteration, every pixel is assigned to the nearest cluster center. After each round of classification is completed, the cluster centers are updated by the average of the spectral vectors and position coordinates of all pixels in the corresponding categories. The classification process is completed when the cluster centers no longer change. More details about the K-means algorithm can be seen in reference [22].After preprocessing, the pixels are classified intoK categories. Pixels belonging to the same category have the same active endmember set. Therefore, we can use the mean vector of each category instead of all pixels for the first unmixing, which can greatly reduce the computational complexity. In other words, the number of pixels participating in the first sparse unmixing is the number of categories. In addition, the improved K-means clustering algorithm considers both spectral and spatial information for classification, which can improve the information utilization of hyperspectral images.


## 3.2.2. First Sparse Unmixing

After preprocessing, we calculate the mean vector of each category and impose single sparsity constraint to unmix these mean vectors. Then, we can get the following unmixing problem:(10)minXM12YM−AXMF2+λ1XM1,1s.t.XM≻¯0,where YM=y1,y2,⋯,yK∈RL×K denotes mean vectors of the K categories. XM denotes the corresponding abundance matrix and XM1,1=∑i=1m∑j=1Kxi,j. λ1 denotes the regularization parameter. Before solving equation (10), we introduce the indicator function to denote the nonnegative constraint, and the unmixing problem becomes
(11)minXM12YM−AXMF2+λ1XM1,1+lR+XM,where lR+XM represents the indicator function, i.e., lR+XM=0 if XM∈R+ and lR+XM=+∞ otherwise. Then, equation (11) can be efficiently solved under alternating direction method of multipliers (ADMM) framework [25]. Set V1=AXM, V2=XM, V3=XM, and then equation (11) can be written as
(12)minXM,V1,V2,V312YM−V1F2+λ1V21,1+lR+V3s.t.V1=AXM,V2=XM,V3=XM.Set(13)gXM,V=12YM−V1F2+λ1V21,1+lR+V3,(14)V≡V1,V2,V3,(15)B=diag−I,(16)G=A,I,IT.Then, equation (12) becomes
(17)minXM,VgXM,Vs.t.GXM+BV=0.The following Lagrangian function can be constructed by equation (17)
(18)LXM,V,D≡gXM,V+μ2GXM+BV−DF2,where μ>0 denotes a positive constant, and D/μ denotes the Lagrangian multiplier of the constraint GXM+BV=0. Under the ADMM framework, each variable is independent when updated. First, expand equation (18) to get
(19)LXM,V1,V2,V3,D1,D2,D3=12YM−V1F2+λ1V21,1+lR+V3+μ2AXM−V1−D1F2+μ2XM−V2−D2F2+μ2XM−V3−D3F2.When updating each variable, only the items related to it need to be considered. The objective optimization functions of variablesXM, V1, V2, and V3 are
(20)XMi+1=argminXMμ2AXM−V1i−D1iF2+μ2XM−V2i−D2iF2+μ2XM−V3i−D3iF2,(21)V1i+1=argminV112YM−V1F2+μ2AXMi+1−V1−D1iF2,(22)V2i+1=argminV2λ1V21,1+μ2XMi+1−V2−D2iF2,(23)V3i+1=argminV3lR+V3+μ2XMi+1−V3−D3iF2.The corresponding solutions are(24)XMi+1=ATA+2I−1ATV1i+D1i+V2i+D2i+V3i+D3i,(25)V1i+1=11+μYM+μAXMi+1−D1i,(26)V2i+1=softXMi+1−D2i,λ1μ,(27)V3i+1=maxXMi+1−D3i,0,where soft denotes the soft-threshold function softy,τ=signymaxy−τ,0. Then, the final algorithm flow for the first sparse unmixing can be seen in Algorithm 1.Algorithm 1: ADMM pseudocode for the first sparse unmixing.
1. Input: mean vectors matrixYM, spectral library A, constant μ.2. Select parameters:λ1.3. Initialization:XM0,V0,D0 and set i=0.4. Repeat:XMi+1=ATA+2I−1ATV1i+D1i+V2i+D2i+V3i+D3iV1i+1=1/1+μYM+μAXMi+1−D1iV2i+1=softXMi+1−D2i,λ1/μV3i+1=maxXMi+1−D3i,0D1i+1=D1i−AXMi+1+V1i+1D2i+1=D2i−XMi+1+V2i+1D3i+1=D3i−XMi+1+V3i+1i=i+15. Until the Stopping Condition Is Met


## 3.2.3. Second Sparse Unmixing

After the first sparse unmixing, we get the abundance estimateXM of these mean vectors. Then, the abundance matrix XM∈Rm×K for these mean vectors is converted to the abundance matrix X1∈Rm×n for the entire image. Specifically, the pixels in the same category have the same abundance estimate as the category’s mean vector. Then, X1 is added as a constraint for the second sparse unmixing, and we impose joint sparsity constraint on the original hyperspectral image. It is well known that solving the l2,0 norm directly will cause NP-hard problem. Most algorithms use the l2,1 norm instead of l2,0 norm to avoid NP-hard problem. However, the l2,1 norm just represents the sum of row-norms, which cannot meet the sparsity requirement sometimes. To address this issue, we introduce the row-hard-threshold function to solve the l2,0 norm directly. Then, we can get the following unmixing problem:
(28)minX12Y−AXF2+α2X1−XF2+λ2X2,0s.t.X≻¯0,where λ2 and α represent the regularization parameters, and the l2,0 norm is defined as follows
(29)X2,0=∑i=1m1Xir2>0,where Xir denotes the i-th row of X. We also use the ADMM method to solve equation (28). Set V1=AX, V2=X, and V3=X. Then equation (28) can be written as
(30)minX,V1,V2,V312Y−V1F2+α2X1−XF2+λ2V22,0+lR+V3s.t.V1=AX,V2=X,V3=X.We can construct the following Lagrangian function using the same method as equation (19)
(31)LX,V1,V2,V3,D1,D2,D3=12Y−V1F2+α2X1−XF2+λ2V22,0+lR+V3+μ2AX−V1−D1F2+μ2X−V2−D2F2+μ2X−V3−D3F2.The iterative update formulas of variablesX and V1 are as follows
(32)Xi+1=ATA+α+μI−1ATV1i+D1i+V2i+D2i+V3i+D3i+αX1,(33)V1i+1=11+μY+μAXi+1−D1i.Before solving variableV2, we first remove items independent of it and get
(34)V2i+1=argminV2λ2V22,0+μ2Xi+1−V2−D2iF2.Inspired by CSUnL0 algorithm [11], we introduce the row-hard-threshold function to solve equation (34). Define
(35)RHtXi,j=xi,j∑jxi,j2>t,0∑jxi,j2≤t,where RHtXi,j represents the j-th element in the i-th row of the row-hard-threshold function RHtX. xi,j denotes the j-th element in the i-th row of X and t denotes the threshold. It has been proved in reference [11] that H=RHtX is the closed-form solution for
(36)minHtH2,0+X−HF2.Equation (34) can be written as follows
(37)V2i+1=argminV2μ22λ2μV22,0+Xi+1−V2−D2iF2.Then, we can get the iterative update formula ofV2 from equation (36) and (37).
(38)V2i+1=RH2λ2/μXi+1−D2i.The iterative update formulas of other variables are similar to those in the first sparse unmixing. Then, the final algorithm flow for the second sparse unmixing can be seen in Algorithm2.Algorithm 2: ADMM pseudocode for the second sparse unmixing.
1. Input: hyperspectral imageY, spectral library A, abundance matrix X1, constant μ.2. Select parameters:λ2 and α.3. Initialization:X0,V0,D0 and set i=0.4. Repeat:Xi+1=ATA+α+μI−1ATV1i+D1i+V2i+D2i+V3i+D3i+αX1V1i+1=1/1+μY+μAXi+1−D1iV2i+1=RH2λ2/μXi+1−D2iV3i+1=maxXi+1−D3i,0D1i+1=D1i−AXi+1+V1i+1D2i+1=D2i−Xi+1+V2i+1D3i+1=D3i−Xi+1+V3i+1i=i+15:until the stopping condition is met.


## 4. Simulated Data Experiments

In this section, we use two sets of simulated hyperspectral data to validate the effectiveness of the proposedDRSUMK−means algorithm. The famous USGS spectral library [26], which contains 498 substances and 224 bands, was used in this paper. Its spectral range is from 0.4 to 2.5 μm. For simplicity, we only use a subset A1∈R224×240 of the USGS spectral library. The angle between any endmembers (Em) in A1 is greater than 4.44 degrees. The DRSUMK−means algorithm was compared with five state-of-the-art algorithms: SUnSAL [8], CLSUnSAL [9], SUnSAL-TV [12], CSUnL0 [11], and S2WSU-W1 [13]. The metrics signal-to-reconstruction error (SRE) and reconstruction overall root-mean-square error (rRMSE) [27] were used to evaluate the unmixing effect of these algorithms. Let xi denote the true abundance vector and x^i an estimate of xi. The metric SRE can be defined as follows
(39)SREdB=10log10Exi22Exi−x^i22,where E⋅ represents the expectation function. The larger the SRE (dB), the higher the unmixing accuracy. The metric rRMSE is defined by
(40)rRMSE=1LnY−Y^F2,where Y^=AX^ represents the reconstruction image, and X^ represents the estimated abundance matrix. n and L represent the number of pixels and spectral bands, respectively. The smaller the rRMSE, the higher the unmixing accuracy.


### 4.1. Generation of Simulated Data

#### 4.1.1. Simulated Data Cube 1 (DC1)

DC1 contains75×75 pixels. The abundance satisfied the ANC and ASC. Five endmembers were randomly selected from A1 for this experiment: Jarosite GDS101, Anorthite HS349.3B, Calcite WS272, Microcline HS82.3B, and Howlite GDS155. Figure 1 shows the true abundance map of each endmember. To be more realistic and detect the antinoise ability of each algorithm, DC1 was contaminated with white Gaussian noise for three signal-to-noise ratio (SNR) levels: 20, 30, and 40 dB.Figure 1
True abundance maps of the endmembers (Em) in DC1.


#### 4.1.2. Simulated Data Cube 2 (DC2)

A more complex dataset (DC2), which contains100×100 pixels, was used for this experiment. We also used the A1 spectral library and randomly selected 9 spectral signatures from it as the endmembers: Wollastonite HS348.3B, Heulandite GDS3, Mordenite GDS18, Almandine WS478, Nepheline HS19.3, Barite HS79.3B, Grossular WS485, Brookite HS443.2B, and Albite HS324.3B. Their corresponding abundance conforms to the Dirichlet distribution and is shown in Figure 2. We can see from Figure 2 that the endmember distribution in DC2 is more complicated than that in DC1. Both steep and smooth transitions exist in DC2, which is more in line with the real distribution of objects. In addition, DC2 was also contaminated by white Gaussian noise with SNRs of 20, 30, and 40 dB.Figure 2
True abundance maps of the endmembers (Em) in DC2.


### 4.2. Three Processes Analysis ofDRSUMK−means

In this section, we will analyze the three processes ofDRSUMK−means: preprocessing, the first sparse unmixing, and the second sparse unmixing. We first use the improved K-means clustering algorithm to classify pixels into different categories. The purpose of preprocessing is to improve the unmixing accuracy and reduce the computational complexity. The choice of K value is an important factor of preprocessing. If the value of K is too small, each category will contain more pixels. There may be large differences between pixels belonging to the same category. The standard deviation of each category will become larger. Then, pixels of the same category may have different active endmember set. In other words, it is inappropriate to use the mean vector of each category instead of all pixels for the first sparse unmixing. If the value of K is too large, the utilization of spatial information will decrease. Pixels with the same active endmember set may be classified into different categories. This will also affect the unmixing accuracy. Therefore, there should be a balance for the parameter K.After preprocessing, we calculate the mean vector of each category and impose single sparsity constraint on these mean vectors for the first sparse unmixing. If there are large differences between pixels belonging to the same category, the results of the first sparse unmixing will be quite different from those of the real situation. Generally speaking, the larger the number of categories, the smaller the difference between pixels in the same category. However, as the number of categories increases, the utilization of spatial information will decrease and the calculation time will increase. Therefore, the classification results have a great influence on the first sparse unmixing.After the first sparse unmixing, we get the abundance estimateXM of these mean vectors. The abundance matrix XM∈Rm×K for these mean vectors is converted to the abundance matrix X1∈Rm×n for the entire image. Then, X1 is added as a constraint for the second sparse unmixing, and we impose joint sparsity constraint on the original hyperspectral image. The abundance matrix X1 can be regarded as a reference for the second sparse unmixing. The closer X1 is to the true abundance, the higher the unmixing accuracy. Therefore, the first unmixing results can directly affect the second unmixing accuracy.


### 4.3. Parameters Selection

TheDRSUMK−means algorithm contains three processes. We can know from Section 4.2 that the previous process has a great influence on the latter one. However, these three processes still have a certain degree of independence. Therefore, their parameters can be selected independently. When selecting some parameters, other parameters maintain appropriate values. The parameters that lead to the highest SRE values are selected. According to the algorithm flow, the number of categories K is first selected from the interval [5, 200] with a step length of 5 for simulated data experiments. For real data experiment, the step length and interval are set to 10 and [200, 700], respectively. The parameters λ1, λ2, and α are temporarily set to 0.01, 0.01, and 10, respectively. After K is selected, the parameter λ1 is selected from the range {0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1}. Finally, the parameters λ2 and α are selected from the ranges {0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1} and {5, 10, 15, 20, 25, 30, 35, 40, 45, 50}, respectively. In addition, the constant μ is set to 0.01 for all experiments. Figure 3 presents the influence of parameters on the performance of DRSUMK−means for DC1 with SNR 30 dB. We can see that the parameters K, λ1, λ2, and α all have a great influence on the SRE value, which indicates that the preprocessing, the first, and second sparse unmixing can affect the unmixing effect greatly.Figure 3
SRE (dB) as a function of parameters inDRSUMK−means for DC1 with SNR 30 dB. (a) Parameter K, (b) parameter λ1, (c) parameters λ2 and α.
(a)(b)(c)


### 4.4. Comparison with Other Algorithms

In this section, we will compare the proposedDRSUMK−means algorithm with several state-of-the-art algorithms: SUnSAL [8], CLSUnSAL [9], SUnSAL-TV [12], CSUnL0 [11], and S2WSU-W1 [13]. The parameters of all algorithms were constantly adjusted to obtain their optimal performance. The SRE and rRMSE values, the optimal parameter values, and the required time for all algorithms are shown in Tables 1 and 2. We can see that the DRSUMK−means algorithm has the highest SRE values and lowest rRMSE values under all the SNRs. Both CSUnL0 and S2WSU-W1 algorithms can also achieve promising results. However, they are sensitive to noise and their unmixing accuracy is low when SNR is 20 dB. The SUnSAL-TV algorithm takes much more time than DRSUMK−means, and its unmixing accuracy is still worse than that of DRSUMK−means. Although the SUnSAL and CLSUnSAL algorithms are faster than DRSUMK−means, their unmixing accuracy is much worse than that of DRSUMK−means.Table 1
SRE (dB), rRMSE, time (s), and parameters of each algorithm for DC1.
SNRCriteriaSUnSALCLSUnSALSUnSAL-TVCSUnL0S2WSU-W1DRSUMK−means20 dBSRE (dB)4.56108.54359.42395.35147.862114.5899rRMSE0.00530.00470.00450.00520.00490.0039Time (s)3.954.1488.0728.2530.325.33Parametersλ=0.58λ=16λ=0.05λTV=0.05a0=0.04λ=0.2K=35, λ1=0.005α=40, λ2=0.530 dBSRE (dB)8.916811.126314.440813.583415.517317.2511rRMSE0.00460.00430.00400.00410.00370.0035Time (s)3.045.5887.6426.7129.175.81Parametersλ=0.09λ=2.8λ=0.007λTV=0.01a0=0.01λ=0.01K=90, λ1=0.005α=20, λ2=0.0540 dBSRE (dB)12.621416.308521.175923.044624.736626.0369rRMSE0.00420.00360.00310.00300.00280.0026Time (s)3.285.6285.9233.3131.456.19Parametersλ=0.05λ=0.1λ=0.001λTV=0.005a0=0.005λ=0.005K=120, λ1=0.001α=10, λ2=0.01Table 2
SRE (dB), rRMSE, time (s), and parameters of each algorithm for DC2.
SNRCriteriaSUnSALCLSUnSALSUnSAL-TVCSUnL0S2WSU-W1DRSUMK−means20 dBSRE (dB)4.26733.401110.97916.83229.164118.7322rRMSE0.00440.00460.00370.00420.00390.0019Time (s)5.8013.51114.6249.8452.3610.82Parametersλ=0.1λ=2.2λ=0.01λTV=0.03a0=0.05λ=0.05K=35, λ1=0.005α=40, λ2=0.530 dBSRE (dB)10.45297.350517.807817.452119.278521.3605rRMSE0.00380.00410.00220.00220.00170.0013Time (s)7.7416.04115.1751.0653.4912.36Parametersλ=0.01λ=0.3λ=0.005λTV=0.007a0=0.005λ=0.03K=100, λ1=0.005α=5, λ2=0.0140 dBSRE (dB)13.372412.723320.652221.215323.863825.1639rRMSE0.00310.00330.00150.00130.000110.0009Time (s)6.3311.28116.3446.5951.0812.68Parametersλ=0.005λ=0.05λ=0.001λTV=0.005a0=0.001λ=0.007K=150, λ1=0.001α=5, λ2=0.005Figures4 and 5 present the estimated abundance maps for DC1 and DC2 with an SNR of 30 dB, respectively. For space considerations, we only selected five endmembers from DC2 for comparison. It can be seen from Figure 4 that the SUnSAL, CLSUnSAL, CSUnL0, and S2WSU-W1 algorithms have more impurities in the background region than SUnSAL-TV and DRSUMK−means algorithms. The number of square region of SUnSAL-TV is less than that of DRSUMK−means, especially for Em2. We can see from Figure 5 that the estimated abundance maps of SUnSAL and CLSUnSAL algorithms are quite different from the real ones. The estimated abundance maps of SUnSAL-TV, CSUnL0, S2WSU-W1, and DRSUMK−means algorithms are similar. However, it can still be seen that DRSUMK−means obtains higher unmixing accuracy than other algorithms, especially for Em1 and Em2. Therefore, the simulated data experiments can validate the effectiveness of the proposed DRSUMK−means algorithm.Figure 4
Abundance maps for DC1 (from top to bottom: Em1 to Em5) withSNR=30dB.Figure 5
Abundance maps for DC2 (from top to bottom: Em1 to Em5) withSNR=30dB.


### 4.5. Computational Cost Analysis

To analyze the computational complexity of the proposedDRSUMK−means algorithm, all algorithms were run in MATLAB R2018a on a laptop equipped with an Intel Core 7 processor (2.3 GHz main frequency) and 16 GB of memory. The maximum number of iterations and the tolerance for all algorithms are set to 1000 and 10−4, respectively. Tables 1 and 2 show the running time of all algorithms for DC1 and DC2. We can see that SUnSAL-TV takes the most time among all algorithms. S2WSU-W1 and CSUnL0 algorithms take much more time than DRSUMK−means. Although DRSUMK−means takes a little more time than SUnSAL and CLSUnSAL algorithms, the unmixing effect of DRSUMK−means is much better than that of them. Therefore, the DRSUMK−means algorithm can obtain promising results with low computational complexity. It should be noted that the time required by SUnSAL-TV is several times longer than that of DRSUMK−means. Therefore, though DRSUMK−means has more parameters to be set, the time required by DRSUMK−means for selecting parameters is still less than that of SUnSAL-TV.


### 4.6. Convergence Analysis

After preprocessing, we impose single sparsity constraint on these mean vectors for the first sparse unmixing. From equation (10), we can know that the first sparse unmixing is similar to the SUnSAL algorithm. Reference [8] has proved that the SUnSAL algorithm is convergent. Therefore, we can conclude that the first sparse unmixing is also convergent. From equation (28), we can know that the result of the first sparse unmixing and the l2,0 norm is introduced as constraints for the second sparse unmixing. It is difficult to give the theoretical convergence analysis of equation (28) directly. Figure 6 presents the convergence curves of equation (28) for DC1 and DC2 with different SNRs. We can see from Figure 6 that the SRE value no longer changes when the number of iterations reaches a certain value. Therefore, the residuals and the maximum number of iterations can be used as iteration stop conditions of DRSUMK−means. In all experiments, the tolerance and the maximum number of iterations for DRSUMK−means are set to 10−4 and 1000, respectively. In addition, it is worth noting that though convergence is not theoretically guaranteed, and the DRSUMK−means algorithm exhibits a robust convergence behavior.Figure 6
Convergence curves ofDRSUMK−means for DC1 and DC2. (a) DC1 and (b) DC2.
(a)(b)


## 4.1. Generation of Simulated Data

### 4.1.1. Simulated Data Cube 1 (DC1)

DC1 contains75×75 pixels. The abundance satisfied the ANC and ASC. Five endmembers were randomly selected from A1 for this experiment: Jarosite GDS101, Anorthite HS349.3B, Calcite WS272, Microcline HS82.3B, and Howlite GDS155. Figure 1 shows the true abundance map of each endmember. To be more realistic and detect the antinoise ability of each algorithm, DC1 was contaminated with white Gaussian noise for three signal-to-noise ratio (SNR) levels: 20, 30, and 40 dB.Figure 1
True abundance maps of the endmembers (Em) in DC1.


### 4.1.2. Simulated Data Cube 2 (DC2)

A more complex dataset (DC2), which contains100×100 pixels, was used for this experiment. We also used the A1 spectral library and randomly selected 9 spectral signatures from it as the endmembers: Wollastonite HS348.3B, Heulandite GDS3, Mordenite GDS18, Almandine WS478, Nepheline HS19.3, Barite HS79.3B, Grossular WS485, Brookite HS443.2B, and Albite HS324.3B. Their corresponding abundance conforms to the Dirichlet distribution and is shown in Figure 2. We can see from Figure 2 that the endmember distribution in DC2 is more complicated than that in DC1. Both steep and smooth transitions exist in DC2, which is more in line with the real distribution of objects. In addition, DC2 was also contaminated by white Gaussian noise with SNRs of 20, 30, and 40 dB.Figure 2
True abundance maps of the endmembers (Em) in DC2.


## 4.1.1. Simulated Data Cube 1 (DC1)

DC1 contains75×75 pixels. The abundance satisfied the ANC and ASC. Five endmembers were randomly selected from A1 for this experiment: Jarosite GDS101, Anorthite HS349.3B, Calcite WS272, Microcline HS82.3B, and Howlite GDS155. Figure 1 shows the true abundance map of each endmember. To be more realistic and detect the antinoise ability of each algorithm, DC1 was contaminated with white Gaussian noise for three signal-to-noise ratio (SNR) levels: 20, 30, and 40 dB.Figure 1
True abundance maps of the endmembers (Em) in DC1.


## 4.1.2. Simulated Data Cube 2 (DC2)

A more complex dataset (DC2), which contains100×100 pixels, was used for this experiment. We also used the A1 spectral library and randomly selected 9 spectral signatures from it as the endmembers: Wollastonite HS348.3B, Heulandite GDS3, Mordenite GDS18, Almandine WS478, Nepheline HS19.3, Barite HS79.3B, Grossular WS485, Brookite HS443.2B, and Albite HS324.3B. Their corresponding abundance conforms to the Dirichlet distribution and is shown in Figure 2. We can see from Figure 2 that the endmember distribution in DC2 is more complicated than that in DC1. Both steep and smooth transitions exist in DC2, which is more in line with the real distribution of objects. In addition, DC2 was also contaminated by white Gaussian noise with SNRs of 20, 30, and 40 dB.Figure 2
True abundance maps of the endmembers (Em) in DC2.


## 4.2. Three Processes Analysis ofDRSUMK−means

In this section, we will analyze the three processes ofDRSUMK−means: preprocessing, the first sparse unmixing, and the second sparse unmixing. We first use the improved K-means clustering algorithm to classify pixels into different categories. The purpose of preprocessing is to improve the unmixing accuracy and reduce the computational complexity. The choice of K value is an important factor of preprocessing. If the value of K is too small, each category will contain more pixels. There may be large differences between pixels belonging to the same category. The standard deviation of each category will become larger. Then, pixels of the same category may have different active endmember set. In other words, it is inappropriate to use the mean vector of each category instead of all pixels for the first sparse unmixing. If the value of K is too large, the utilization of spatial information will decrease. Pixels with the same active endmember set may be classified into different categories. This will also affect the unmixing accuracy. Therefore, there should be a balance for the parameter K.After preprocessing, we calculate the mean vector of each category and impose single sparsity constraint on these mean vectors for the first sparse unmixing. If there are large differences between pixels belonging to the same category, the results of the first sparse unmixing will be quite different from those of the real situation. Generally speaking, the larger the number of categories, the smaller the difference between pixels in the same category. However, as the number of categories increases, the utilization of spatial information will decrease and the calculation time will increase. Therefore, the classification results have a great influence on the first sparse unmixing.After the first sparse unmixing, we get the abundance estimateXM of these mean vectors. The abundance matrix XM∈Rm×K for these mean vectors is converted to the abundance matrix X1∈Rm×n for the entire image. Then, X1 is added as a constraint for the second sparse unmixing, and we impose joint sparsity constraint on the original hyperspectral image. The abundance matrix X1 can be regarded as a reference for the second sparse unmixing. The closer X1 is to the true abundance, the higher the unmixing accuracy. Therefore, the first unmixing results can directly affect the second unmixing accuracy.


## 4.3. Parameters Selection

TheDRSUMK−means algorithm contains three processes. We can know from Section 4.2 that the previous process has a great influence on the latter one. However, these three processes still have a certain degree of independence. Therefore, their parameters can be selected independently. When selecting some parameters, other parameters maintain appropriate values. The parameters that lead to the highest SRE values are selected. According to the algorithm flow, the number of categories K is first selected from the interval [5, 200] with a step length of 5 for simulated data experiments. For real data experiment, the step length and interval are set to 10 and [200, 700], respectively. The parameters λ1, λ2, and α are temporarily set to 0.01, 0.01, and 10, respectively. After K is selected, the parameter λ1 is selected from the range {0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1}. Finally, the parameters λ2 and α are selected from the ranges {0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1} and {5, 10, 15, 20, 25, 30, 35, 40, 45, 50}, respectively. In addition, the constant μ is set to 0.01 for all experiments. Figure 3 presents the influence of parameters on the performance of DRSUMK−means for DC1 with SNR 30 dB. We can see that the parameters K, λ1, λ2, and α all have a great influence on the SRE value, which indicates that the preprocessing, the first, and second sparse unmixing can affect the unmixing effect greatly.Figure 3
SRE (dB) as a function of parameters inDRSUMK−means for DC1 with SNR 30 dB. (a) Parameter K, (b) parameter λ1, (c) parameters λ2 and α.
(a)(b)(c)


## 4.4. Comparison with Other Algorithms

In this section, we will compare the proposedDRSUMK−means algorithm with several state-of-the-art algorithms: SUnSAL [8], CLSUnSAL [9], SUnSAL-TV [12], CSUnL0 [11], and S2WSU-W1 [13]. The parameters of all algorithms were constantly adjusted to obtain their optimal performance. The SRE and rRMSE values, the optimal parameter values, and the required time for all algorithms are shown in Tables 1 and 2. We can see that the DRSUMK−means algorithm has the highest SRE values and lowest rRMSE values under all the SNRs. Both CSUnL0 and S2WSU-W1 algorithms can also achieve promising results. However, they are sensitive to noise and their unmixing accuracy is low when SNR is 20 dB. The SUnSAL-TV algorithm takes much more time than DRSUMK−means, and its unmixing accuracy is still worse than that of DRSUMK−means. Although the SUnSAL and CLSUnSAL algorithms are faster than DRSUMK−means, their unmixing accuracy is much worse than that of DRSUMK−means.Table 1
SRE (dB), rRMSE, time (s), and parameters of each algorithm for DC1.
SNRCriteriaSUnSALCLSUnSALSUnSAL-TVCSUnL0S2WSU-W1DRSUMK−means20 dBSRE (dB)4.56108.54359.42395.35147.862114.5899rRMSE0.00530.00470.00450.00520.00490.0039Time (s)3.954.1488.0728.2530.325.33Parametersλ=0.58λ=16λ=0.05λTV=0.05a0=0.04λ=0.2K=35, λ1=0.005α=40, λ2=0.530 dBSRE (dB)8.916811.126314.440813.583415.517317.2511rRMSE0.00460.00430.00400.00410.00370.0035Time (s)3.045.5887.6426.7129.175.81Parametersλ=0.09λ=2.8λ=0.007λTV=0.01a0=0.01λ=0.01K=90, λ1=0.005α=20, λ2=0.0540 dBSRE (dB)12.621416.308521.175923.044624.736626.0369rRMSE0.00420.00360.00310.00300.00280.0026Time (s)3.285.6285.9233.3131.456.19Parametersλ=0.05λ=0.1λ=0.001λTV=0.005a0=0.005λ=0.005K=120, λ1=0.001α=10, λ2=0.01Table 2
SRE (dB), rRMSE, time (s), and parameters of each algorithm for DC2.
SNRCriteriaSUnSALCLSUnSALSUnSAL-TVCSUnL0S2WSU-W1DRSUMK−means20 dBSRE (dB)4.26733.401110.97916.83229.164118.7322rRMSE0.00440.00460.00370.00420.00390.0019Time (s)5.8013.51114.6249.8452.3610.82Parametersλ=0.1λ=2.2λ=0.01λTV=0.03a0=0.05λ=0.05K=35, λ1=0.005α=40, λ2=0.530 dBSRE (dB)10.45297.350517.807817.452119.278521.3605rRMSE0.00380.00410.00220.00220.00170.0013Time (s)7.7416.04115.1751.0653.4912.36Parametersλ=0.01λ=0.3λ=0.005λTV=0.007a0=0.005λ=0.03K=100, λ1=0.005α=5, λ2=0.0140 dBSRE (dB)13.372412.723320.652221.215323.863825.1639rRMSE0.00310.00330.00150.00130.000110.0009Time (s)6.3311.28116.3446.5951.0812.68Parametersλ=0.005λ=0.05λ=0.001λTV=0.005a0=0.001λ=0.007K=150, λ1=0.001α=5, λ2=0.005Figures4 and 5 present the estimated abundance maps for DC1 and DC2 with an SNR of 30 dB, respectively. For space considerations, we only selected five endmembers from DC2 for comparison. It can be seen from Figure 4 that the SUnSAL, CLSUnSAL, CSUnL0, and S2WSU-W1 algorithms have more impurities in the background region than SUnSAL-TV and DRSUMK−means algorithms. The number of square region of SUnSAL-TV is less than that of DRSUMK−means, especially for Em2. We can see from Figure 5 that the estimated abundance maps of SUnSAL and CLSUnSAL algorithms are quite different from the real ones. The estimated abundance maps of SUnSAL-TV, CSUnL0, S2WSU-W1, and DRSUMK−means algorithms are similar. However, it can still be seen that DRSUMK−means obtains higher unmixing accuracy than other algorithms, especially for Em1 and Em2. Therefore, the simulated data experiments can validate the effectiveness of the proposed DRSUMK−means algorithm.Figure 4
Abundance maps for DC1 (from top to bottom: Em1 to Em5) withSNR=30dB.Figure 5
Abundance maps for DC2 (from top to bottom: Em1 to Em5) withSNR=30dB.


## 4.5. Computational Cost Analysis

To analyze the computational complexity of the proposedDRSUMK−means algorithm, all algorithms were run in MATLAB R2018a on a laptop equipped with an Intel Core 7 processor (2.3 GHz main frequency) and 16 GB of memory. The maximum number of iterations and the tolerance for all algorithms are set to 1000 and 10−4, respectively. Tables 1 and 2 show the running time of all algorithms for DC1 and DC2. We can see that SUnSAL-TV takes the most time among all algorithms. S2WSU-W1 and CSUnL0 algorithms take much more time than DRSUMK−means. Although DRSUMK−means takes a little more time than SUnSAL and CLSUnSAL algorithms, the unmixing effect of DRSUMK−means is much better than that of them. Therefore, the DRSUMK−means algorithm can obtain promising results with low computational complexity. It should be noted that the time required by SUnSAL-TV is several times longer than that of DRSUMK−means. Therefore, though DRSUMK−means has more parameters to be set, the time required by DRSUMK−means for selecting parameters is still less than that of SUnSAL-TV.


## 4.6. Convergence Analysis

After preprocessing, we impose single sparsity constraint on these mean vectors for the first sparse unmixing. From equation (10), we can know that the first sparse unmixing is similar to the SUnSAL algorithm. Reference [8] has proved that the SUnSAL algorithm is convergent. Therefore, we can conclude that the first sparse unmixing is also convergent. From equation (28), we can know that the result of the first sparse unmixing and the l2,0 norm is introduced as constraints for the second sparse unmixing. It is difficult to give the theoretical convergence analysis of equation (28) directly. Figure 6 presents the convergence curves of equation (28) for DC1 and DC2 with different SNRs. We can see from Figure 6 that the SRE value no longer changes when the number of iterations reaches a certain value. Therefore, the residuals and the maximum number of iterations can be used as iteration stop conditions of DRSUMK−means. In all experiments, the tolerance and the maximum number of iterations for DRSUMK−means are set to 10−4 and 1000, respectively. In addition, it is worth noting that though convergence is not theoretically guaranteed, and the DRSUMK−means algorithm exhibits a robust convergence behavior.Figure 6
Convergence curves ofDRSUMK−means for DC1 and DC2. (a) DC1 and (b) DC2.
(a)(b)


## 5. Real Data Experiment

The public Cuprite dataset (http://aviris.jpl.nasa.gov/html/aviris.freedata.html) was used for this experiment. It is representative and has been widely used for hyperspectral unmixing experiments. The Cuprite dataset was collected by the Airborne Visible Infrared Imaging Spectrometer in 1997. To simplify the calculation, we used a 250×191 subset of the Cuprite dataset, similarly as in [28–30]. It contains 224 spectral bands, but after removing the low-SNR and high-water-absorption bands (1–2, 105–115, 150–170, 223–224), 188 spectral bands were retained. The whole USGS spectral library A∈R224×498 was used. Figure 7 shows the distribution of different substances in Cuprite dataset. It was produced in 1995 by the USGS, and it is an important reference to qualitatively evaluate the unmixing performance. We selected three representative endmembers (buddingtonite, alunite, and chalcedony) from this dataset and compared the abundance maps estimated by each algorithm with those generated by USGS software Tricorder 3.32, which are presented in Figure 8. The parameters of SUnSAL, CLSUnSAL, SUnSAL-TV, CSUnL0, and S2WSU-W1 are set according to references [9, 11, 13]. Table 3 presents the parameter values, processing time, and rRMSE values of each algorithm for the Cuprite dataset.Figure 7
Distribution map of different substances for Cuprite dataset.Figure 8
Abundance maps estimated for (from top to bottom) alunite, buddingtonite, and chalcedony.Table 3
The rRMSE values, time (s), and parameters of each algorithm for Cuprite dataset.
CriteriaSUnSALCLSUnSALSUnSAL-TVCSUnL0S2WSU-W1DRSUMK−meansrRMSE8.7e−58.3e−55.9e−57.8e−53.6e−52.1e−5Time (s)217.3194.51441.4631.5824.7322.5Parametersλ=0.001λ=0.01λ=0.001λTV=0.001a0=6e−5λ=0.002K=390, λ1=0.001α=5, λ2=0.001We can see from Figure8 that the backgrounds of SUnSAL and SUnSAL-TV algorithms are not as pure and clear as those of other algorithms, especially for the mineral buddingtonite. There are also more oversmoothing and blurred edges phenomenon in the abundance maps of SUnSAL-TV. CLSUnSAL and CSUnL0 algorithms have worse abundance maps than other algorithms for the mineral chalcedony. The unmixing effects of S2WSU-W1 and DRSUMK−means algorithms are similar, but it can still be seen that DRSUMK−means has higher abundance for the mineral chalcedony and is closer to the true distribution. In addition, it can be seen from Table 3 that the DRSUMK−means algorithm has the lowest rRMSE value among all algorithms. That is to say, DRSUMK−means has the highest unmixing accuracy for the real dataset. Therefore, the real dataset can also prove the effectiveness of the proposed DRSUMK−means algorithm.


## 6. Conclusions

In this paper, a novel double regression-based sparse unmixing model (DRSUM) was proposed. DRSUM decomposes the complex unmixing process into two simple sparse regressions, which can reduce the computational complexity. DRSUM is an open model and framework. In different application scenarios, we can add different types of constraints into DRSUM to guarantee the unmixing accuracy. In addition, we can perform appropriate preprocessing to further improve the unmixing results. Under this model, a specific algorithm called double regression-based sparse unmixing viaK-means (DRSUMK−means) is proposed. By appropriate preprocessing, DRSUMK−means can comprehensively utilize spatial information, spectral information, and sparsity. Simulated and real data experiments have proven that DRSUMK−means can obtain better unmixing results with lower computational complexity.



---

*Source: 5575155-2021-09-06.xml*