Weighted Group Sparsity-Constrained Tensor Factorization for Hyperspectral Unmixing

Abstract

Recently, unmixing methods based on nonnegative tensor factorization have played an important role in the decomposition of hyperspectral mixed pixels. According to the spatial prior knowledge, there are many regularizations designed to improve the performance of unmixing algorithms, such as the total variation (TV) regularization. However, these methods mostly ignore the similar characteristics among different spectral bands. To solve this problem, this paper proposes a group sparse regularization that uses the weighted constraint of the $L_{2,1}$ norm, which can not only explore the similar characteristics of the hyperspectral image in the spectral dimension, but also keep the data smooth characteristics in the spatial dimension. In summary, a non-negative tensor factorization framework based on weighted group sparsity constraint is proposed for hyperspectral images. In addition, an effective alternating direction method of multipliers (ADMM) algorithm is used to solve the algorithm proposed in this paper. Compared with the existing popular methods, experiments conducted on three real datasets fully demonstrate the effectiveness and advancement of the proposed method.

1. Introduction

Hyperspectral remote sensing images have a wide range of applications in the field of earth observation because of the rich spectral information [1,2,3]. However, the increase in spectral resolution usually results in a decrease in spatial resolution due to the influence of spectral splitting [4]. As a result, hyperspectral images generally have low spatial resolution, which leads to the inclusion of multiple substances in the image pixels. Therefore, mixed pixels are widely present in the hyperspectral images (HSI). The purpose of hyperspectral unmixing is to extract all pure materials (called endmembers) in the HSI and the percentage of endmembers in each pixel (called abundances) [5].

To solve the above-mentioned unmixing problem, numerous models have been proposed, the most popular of which is the linear mixed model (LMM). Under the LMM, each of the mixed pixels is regarded as the linear expression form of the endmembers, and the expression coefficients of endmembers constitute the abundance matrix [6] of the image. Among them, nonnegative matrix factorization (NMF) has been favored by researchers for its clear physical meaning and simple optimization process. The purpose of NMF is to find a base matrix of the original data and the expression coefficients under this base matrix [7]. In unmixing task, the optimized base matrix is the endmember matrix while the expression coefficient matrix is the abundance matrix. However, the solution of the simple NMF model is a non-convex problem, which leads to unstable unmixing results [8]. Therefore, to improve the performance of unmixing, many improved NMF-based unmixing methods have been proposed, for instance, sparse-based and spatial-based.

According to the prior of the distribution of ground objects, the endmembers in the mixed pixels in HSI are much smaller than the endmember spectral library [9]. Therefore, the abundance matrix is reasonably sparse. As a result, the sparsity constraint of abundance matrix becomes one of the most popular constraints. The $L_0$ regularization is the most ideal sparse regularizer [10]. Unfortunately, the NMF-based unmixing algorithm with $L_0$ sparsity constraint is an NP-hard problem [11]. As an improvement, the $L_p$ norm is widely used in the unmixing frameworks based on NMF to improve the stability and effectiveness of the algorithm. Guo et al. [12] uses the $L_1$ regularization to relax the $L_0$ regularization while [13] uses the $L_2$ regularization to constrain the abundance matrix. $L_{1/2}$ has also been proven to make the abundance matrix get a good sparse result [14,15]. Some methods also use the $L_{2,1}$ regularization to simultaneously explore the sparse structure in the abundance map [16], which has been proven to have good performance in the unmixing algorithm.

Besides, the spatial distribution of ground objects has the characteristics of local similarity and global similarity [17,18]. That is, there is a certain spatial similarity structure in HSI. Therefore, spatial similarity constraints are also really popular technics to improve the effectiveness of NMF-based methods. Lu et al. [19] propose a manifold regularization, which learns the global similarity structure of data to constrain the abundance matrix and force the abundance vectors of similar pixels to be similar. Manifold constraints have also proved by [20,21] to be a really effective way to explore the spatial structures. In [22], it is considered that the abundance of a pixel is the consistence as the average value of its neighboring pixels to maintain the smooth characteristics of HSI. In addition to these, the total variation regularization [23] and its variation [24] are also applied to the NMF-based unmixing algorithms, which can effectively retain the smoothing relationship of the image in the abundance matrix. In [25], the clustering algorithm is used to explore the global similarity relationship of pixels in the HSI, and it is combined with the NMF framework to improve the performance of unmixing.

The NMF-based methods mentioned above have significantly improved the performance of unmixing. Unfortunately, this kind of algorithms has a great flaw. The NMF-based framework reshapes the original 3-D HSI into a 2-D matrix before unmixing. This method of image reshaping severely destroys the 3-D cube structure of HSI, and it is difficult to completely explore the missing spatial information even if spatial regularization constraints are subsequently used [26]. To slow down the loss of information and make better use of all the spatial and spectral information in HSI, many researchers have begun to use nonnegative tensor factorization (NTF) instead of NMF framework to carry out unmixing tasks. The NTF framework directly treats HSI as a 3-D tensor, and regards the image as a whole during the optimization process, which can completely transfer all the spatial relationships of the original image to the abundance tensor.

Zhang et al. [27] innovatively applies tensor factorization (TF) to unmixing for the first time, which obtains performance comparable to NMF framework. The canonical polyadic decomposition-based (CPD) method is used to extract endmember matrix in [28]. However, this kind of methods decomposes the image into a finite number of rank-1 tensors, the number of which is generally much larger than the actual number of endmembers. Another TF method, tucker factorization, is also applied to the unmixing task in [29]. Regrettably, tucker factorization often fails to correspond the unmixing result with the actual distribution of objects. To make the result of TF method has a clear physical meaning, Qian et al. [26] propose a matrix-vector-based NTF (MV-NTF) method, which decomposes the HSI into multiple component tensors, then each component tensor is expressed by the outer product of a matrix and a vector. This method greatly promotes the application of NTF framework in unmixing task. In addition, many researchers add some effective regularizations to improve the effectiveness of the MV-NTF method. For instance, a spatial low-rank tensor regularization is proposed by [30,31] to maintain the sparsity of the abundance tensor. Sun et al. [32] uses a local spatial similarity constraint to explore the local smooth features of the data. A regularization for exploring the low-rank spatial structure of data is also introduced in [33]. The purpose is to preserve the HSI properties completely in the abundance tensor to improve the performance of unmixing. In addition, TV regularization has also been added to the NTF framework, such as [23,34].

Although the above methods have good performance in hyperspectral unmixing, the NTF-based models still have room for improvement in combining spatial and spectral information. According to [24], the difference image obtained after the image is differentiated along the horizontal and vertical directions is also sparse, which not only has zero elements, but also has zero connected regions in space. Therefore, a weighted group sparsity-constrained tensor factorization (WSCTF) method is proposed in this paper. Specifically, the HSI is differentiated along the horizontal and vertical directions at the same time to obtain two three-dimensional difference images firstly, which can explore the sparse structure of the HSI. Subsequently, the difference images are folded along the mode-1 and mode-2 to obtain two sparse matrices. Then, the $L_{2,1}$ norm is used to constrain the sparse structure of difference matrix, which can not only maintain the sparseness of the difference image, but also constrain the zero connected region. To enhance the sparsity, the proposed method WSCTF also provides a weight for the sparsity constraint in vertical and horizontal directions. In addition, the $L_{2,1}$ regularization is used to constrain the sparsity of the abundance tensor. In summary, the contributions of this paper can be concluded as:

• A sparse constrained regularization is proposed to explore the sparse structure of the HSI differential image in the horizontal and vertical directions. In addition, weight coefficients are used to enhance sparsity.
• The $L_{2,1}$ norm is utilized to constrain the sparse regularizer, and the $L_{2,1}$ norm is embedded in the NTF framework to maintain the sparsity of the abundance tensor.
• The proposed algorithm WSCTF uses the alternating direction method of multipliers to iteratively optimize. Experiments on three real data prove the effectiveness of our algorithm.

The paper is arranged as follows. Section 2 introduces the concepts related to tensor and the unmixing framework based on NTF. The proposed method WSCTF and its optimization are described in detail in Section 3. The synthetic and real experimental results are analyzed in Section 4. Section 5 is the conclusion of this paper.

2. Related Work

2.1. Notation

According to experiences, a matrix is a two-dimensional array, and a tensor is an extension of the matrix. In this paper, an italic letter is used to represent a scalar, e.g., y and I. A vector is expressed by a regular lowercase letter, e.g., $\mathrm{y}\in {\mathbb{R}}^{I\times 1}$. A matrix is given using a bold captial letter, e.g., $\mathbf{M}\in {\mathbb{R}}^{I_1\times I_2}$, and a tensor is defined as a calligraphic letter, e.g., $\mathcal{Y}\in {\mathbb{R}}^{I_1\times I_2\times \dots \times I_N}$.

Tensor Mode: Each dimension of a tensor is called a mode. For instance, the tensor $\mathcal{Y}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$ has 3 modes. Mode-n Unfolding: It can reshape a tensor to a matrix form. For example, given the tensor $\mathcal{Y}\in {\mathbb{R}}^{I\times J\times K}$, the mode-n unfolding of $\mathcal{Y}$ is represented as:

$$\begin{array}{cc}\hfill & {\left({\mathbf{Y}}_{<1>}\right)}_{(j-1)K+k,i}=y_{ijk},\hfill \\ \hfill & {\left({\mathbf{Y}}_{<2>}\right)}_{(k-1)I+i,j}=y_{ijk},\hfill \\ \hfill & {\left({\mathbf{Y}}_{<3>}\right)}_{(i-1)J+j,k}=y_{ijk}.\hfill \end{array}$$

(1)

in which ${\mathbf{Y}}_{<1>}$ is the mode-1 unfolding of tensor $\mathcal{Y}$. ${\mathbf{Y}}_{<2>}$ and ${\mathbf{Y}}_{<3>}$ are the mode-2 and mode-3 unfolding, respectively.

Frobenius Norm: The Frobenius Norm of tensor $\mathcal{Y}\in {\mathbb{R}}^{I\times J\times K}$ is given by:

$$\left|\right|\mathcal{Y}|{|}_F={\left(\sum _{i=1}^I\sum _{j=1}^J\sum _{k=1}^K|y_{ijk}{|}^2\right)}^{1/2}.$$

(2)

Mode-n Product: Given a tensor $\mathcal{Y}$, it can be expressed as $\mathcal{Y}={\mathcal{A}}_{\times n}\mathbf{M}$, in which ${\mathcal{A}}_{\times n}$ denotes the mode-n unfolding of tensor $\mathcal{A}$.

2.2. NTF Unmixing Method

According to [35], the linear spectral mixing model can be represented as a linear combination by the endmember matrix and abundance maps. Specifically, given hyperspectral data $\mathcal{Y}\in {\mathbb{R}}^{I\times J\times K}$, then it can be written as:

$$\mathcal{Y}=\mathcal{A}{\times }_3\mathbf{M}+\mathcal{N},$$

(3)

where $\mathcal{A}\in {\mathbb{R}}^{I\times J\times P}$ is the abundance tensor. $\mathcal{A}{\times }_3$ denotes the mode-3 unfolding of tensor $\mathcal{A}$. And $\mathbf{M}\in {\mathbb{R}}^{I\times J}$ and $\mathcal{N}\in {\mathbb{R}}^{I\times J\times K}$ are the endmember matrix and noise tensor existing in HSI, respectively. Besides, I, J and K are the dimensions of hyperspectral data, and P is the number of endmembers. As a result, the unmixing model based on tensor factorization can be given by

$$\begin{array}{c}\hfill J(\mathcal{A},\mathbf{M})=\frac{1}{2}{∥\mathcal{Y}-\mathcal{A}{\times }_3\mathbf{M}∥}_F^2,\end{array}$$

(4)

in which $\left|\right|·|{|}_F$ denotes the Frobenius norm. Since hyperspectral unmixing is a practical image processing process, there must be non-negative constraint here. In addition, the sum of the abundance coefficients inside each pixel should be one. These two properties create two constraints in the function (4). Specifically, the objective function (4) can be updated as:

$$\begin{array}{c}\hfill J(\mathcal{A},\mathbf{M})=\frac{1}{2}{∥\mathcal{Y}-{\mathcal{A}}_{\times 3}\mathbf{M}∥}_F^2,\\ \hfill s.t.\mathcal{A}\ge \mathbf{0},\mathcal{A}{\times }_1{\mathbf{1}}_P={\mathbf{1}}_{I\times J},\end{array}$$

(5)

in which $\mathcal{A}{\times }_1$ is mode-1 unfolding of the abundance tensor $\mathcal{A}$.

Unfortunately, the above function (5) is a non-convex problem, which needs to be optimized for both the abundance $\mathcal{A}$ and the endmember matrix $\mathbf{M}$. To improve the performance of the algorithm, Ref. [36] adds the $L_{2,1}$ sparsity constraint to the NTF unmixing model. Therefore, it can be given as:

$$\begin{array}{c}J(\mathcal{A},\mathbf{M})=\frac{1}{2}{∥\mathcal{Y}-{\mathcal{A}}_{\times 3}\mathbf{M}∥}_F^2+{\lambda }_S{∥\mathcal{A}∥}_{2,1}\\ s.t.\mathcal{A}\ge \mathbf{0},\mathcal{A}{\times }_1{\mathbf{1}}_P={\mathbf{1}}_{I\times J}.\end{array}$$

(6)

In addition, some spatial constraints are also added to the NTF unmixing model to explore the spatial information of HSI, especially TV regularization [34]. Specifically, the NTF unmixing model based on the TV regularization can be described as follows:

$$\begin{array}{c}J(\mathcal{A},\mathbf{M})=\frac{1}{2}{∥\mathcal{Y}-{\mathcal{A}}_{\times 3}\mathbf{M}∥}_F^2+{\lambda }_S{∥\mathcal{A}∥}_{2,1}+{\lambda }_{TV}{∥\mathcal{A}∥}_{TV}\\ s.t.\mathcal{A}\ge \mathbf{0},\mathcal{A}{\times }_1{\mathbf{1}}_P={\mathbf{1}}_{I\times J},\end{array}$$

(7)

where ${∥\mathcal{A}∥}_{TV}$ is the TV regularization, and ${\lambda }_{TV}$ denotes the weight parameter. Since the TV regularization has achieved excellent results in the spatial relationships of data, especially in maintaining the smoothness of the data, it has achieved comparable performance on the hyperspectral unmixing, which has been proved in [37,38]. However, the existing methods are still unable to fully explore some important characteristics of data. For instance, few work consider continuity in the direction of the spectral dimension. In addition, according to [39], the three-dimensional difference images obtained after also have group sparsity. Based on these ideas, this paper explores the three-dimensional spatial structure information of the image, including the spectral similarity of the difference images and the sparsity of the spatial group.

3. Proposed Method

In response to the questions mentioned above, a new regularization is designed to explore the spectral continuity and group sparsity of difference images of HSI. As a result, a weighted group sparsity constraint non-negative tensor factorization (WSCTF) method for hyperspectral unmixing is proposed. The proposed algorithm and its optimization method are described in detail in this section.

3.1. WSCTF Model

According to the existing experience, the TV regularization is widely used to explore the spatial smooth structure in HSI. Specifically, the TV regularizer can be expressed as follows:

$${\parallel \mathcal{A}\parallel }_{TV}={\parallel \mathcal{D}\mathcal{A}\parallel }_1={∥{\mathcal{D}}_1\mathcal{A}∥}_1+{∥{\mathcal{D}}_2\mathcal{A}∥}_1,$$

(8)

where $\mathcal{D}$ denotes the difference operator. ${\mathcal{D}}_1$ and ${\mathcal{D}}_2$ are vertical and horizontal difference, respectively. $\left|\right|·{\left|\right|}_1$ is the $l_1$ norm. In detail, given the spatial coordinate of a abundance tensor $(i,j,p)$, then the difference results in the vertical and horizontal directions can be described as:

$$\left\{\begin{array}{c}{\mathcal{D}}_1\mathcal{A}\left(i,j,p\right)=\mathcal{A}\left(i,j+1,p\right)-\mathcal{A}\left(i,j,p\right),\hfill \\ {\mathcal{D}}_2\mathcal{A}\left(i,j,p\right)=\mathcal{A}\left(i+1,j,p\right)-\mathcal{A}\left(i,j,p\right).\hfill \end{array}\right.$$

(9)

Based on the above analysis, it can be concluded that the difference images in the horizontal and vertical directions are also three-dimensional tensors. According to the prior, the difference image has continuity in the spectral dimension. In addition, the difference image has extremely high sparsity. In the difference image, not only many elements are zeros, but also some connected regions with zeros elements. Therefore, a new weighted sparsity regularization is designed to explore these properties. Specifically, the regularization can be given as:

$$\Re \left(\mathcal{A}\right)=\left|\right|W_1\odot \left({\mathcal{D}}_1{\mathcal{A}}_{\times 1}\right)|{|}_{2,1}+\left|\right|W_2\odot \left({\mathcal{D}}_2{\mathcal{A}}_{\times 2}\right)|{|}_{2,1},$$

(10)

in which $W_1$ and $W_2$ are aiming at enhancing the sparsity of difference images. ${\mathcal{A}}_{\times 1}$ and ${\mathcal{A}}_{\times 2}$ are the mode-1 and mode-2 unfolding of tensor $\mathcal{A}$. $\left|\right|·|{|}_{2,1}$ denotes the $l_{2,1}$ norm, which can be used to explore the group sparsity of difference image. Therefore, the objective function (7) can be improved by:

$$\begin{array}{c}J(\mathcal{A},\mathbf{M})=\frac{1}{2}{∥\mathcal{Y}-{\mathcal{A}}_{\times 3}\mathbf{M}∥}_F^2+{\lambda }_S{∥\mathcal{A}∥}_{2,1}+\lambda \Re \left(\mathcal{A}\right)\\ s.t.\mathcal{A}\ge \mathbf{0},\mathcal{A}{\times }_1{\mathbf{1}}_P={\mathbf{1}}_{I\times J},\end{array}$$

(11)

where $\lambda$ is regularization parameter. In addition, to be more prepared to explore the sparse structure of the abundance tensor $\mathcal{A}$, a weight constraint is also added to the sparse regularization. Finally, the objective function of WSCTF proposed in this paper can be expressed as:

$$\begin{array}{c}J(\mathcal{A},\mathbf{M})=\frac{1}{2}{∥\mathcal{Y}-{\mathcal{A}}_{\times 3}\mathbf{M}∥}_F^2+{\lambda }_S{∥W_S\odot \mathcal{A}∥}_{2,1}+\lambda \Re \left(\mathcal{A}\right)\\ s.t.\mathcal{A}\ge \mathbf{0},\mathcal{A}{\times }_1{\mathbf{1}}_P={\mathbf{1}}_{I\times J},\end{array}$$

(12)

in which $W_s$ is the weight to promote the sparsity of abundance tensor. In the proposed method, not only the spatial neighborhood smoothing information is considered, but also the sparsity characteristics of the difference image, especially the group sparsity, are described. The optimization method of objective function (12) is explained in detail in the next subsection.

3.2. Optimization

The Formula (12) is a non-convex problem, which is directly difficult to solve. According to the exist experience, the ADMM algorithm is used to update each variable, which can be solved by decomposing the original problem into multiple sub-problems.

(1) Update M: For the function (12), it first considers the update of the endmember matrix $\mathbf{M}$. At this time, the abundance tensor $\mathcal{A}$ is fixed, so the objective function can be expressed as:

$$J(\mathcal{A},\mathbf{M})=\frac{1}{2}{∥\mathcal{Y}-{\mathcal{A}}_{\times 3}\mathbf{M}∥}_F^2.$$

(13)

It can directly obtain the updated rule of endmember $\mathbf{M}$ as:

$$\mathbf{M}\leftarrow \mathbf{M}.\ast \left({{\mathcal{Y}}_{\times 3}{\mathcal{A}}_{\times 3}}^{\mathrm{T}}\right)./\left({\mathbf{M}{\mathcal{A}}_{\times 3}{\mathcal{A}}_{\times 3}}^{\mathrm{T}}\right),$$

(14)

in which ${\mathcal{Y}}_{\times 3}$ and ${\mathcal{A}}_{\times 3}$ are the mode-3 unfolding of tensors. $.\ast$ and $./$ are multiplication and division operators.

In addition, by introducing six intermediate variables ${\mathcal{G}}_1$, ${\mathcal{G}}_2$, ${\mathcal{G}}_3$, ${\mathcal{G}}_4$, ${\mathcal{G}}_5$ and ${\mathcal{G}}_6$, the Formula (12) can be rewritten as:

$$\begin{array}{c}\begin{array}{c}\hfill \underset{\mathcal{A},\mathbf{M}}{min}\frac{1}{2}{∥{\mathcal{G}}_1-Y∥}_F^2+{\lambda }_{\mathrm{s}}{∥W_S\odot {\mathcal{G}}_2∥}_{2,1}+\lambda {∥W_1\odot {\mathcal{G}}_3∥}_{2,1}+\lambda {∥W_2\odot {\mathcal{G}}_4∥}_{2,1}+l_{R^{+}}\left({\mathcal{G}}_5\right)+l_{\left\{1\right\}}\left({\mathcal{G}}_6\right),\end{array}\\ s.t.{\mathcal{G}}_1={{\mathcal{A}}_{\times }}_3\mathbf{M},\\ {\mathcal{G}}_2=\mathcal{A},\\ {\mathcal{G}}_3={\mathcal{D}}_1{\mathcal{A}}_{\times 1},\\ {\mathcal{G}}_4={\mathcal{D}}_2{\mathcal{A}}_{\times 2},\\ {\mathcal{G}}_5=\mathcal{A},\\ {\mathcal{G}}_6=\mathcal{A}.\end{array}$$

(15)

In the above Formula (15), ${\mathcal{G}}_5$ and ${\mathcal{G}}_6$ are used to meet the abundance nonnegative constraint (ANC) and abundance sum-to-one constraint (ASC). They are defined as:

$$l_{R^{+}}\left(\mathcal{X}\right)=\left\{\begin{array}{c}0,if\mathcal{X}\ge 0,\\ +\infty ,otherwise.\end{array}\right.l_{\left\{\mathbf{1}\right\}}\left(\mathcal{X}\right)=\left\{\begin{array}{c}0,if\mathrm{sum}\left({\mathcal{X}}_{\times 3}\right)=\mathbf{1},\hfill \\ +\infty ,otherwise.\hfill \end{array}\right.$$

(16)

in which $\mathbf{1}$ denotes a matrix with all ones. Hence, the Lagrangian function of objective function (15) is:

$$\begin{array}{c}\hfill \underset{\mathcal{A},\mathbf{M},{\mathcal{G}}_i}{min}\frac{1}{2}{∥{\mathcal{G}}_1-Y∥}_F^2+{\lambda }_{\mathrm{s}}{∥W_S\odot {\mathcal{G}}_2∥}_{2,1}+\lambda {∥W_1\odot {\mathcal{G}}_3∥}_{2,1}+\lambda {∥W_2\odot {\mathcal{G}}_4∥}_{2,1}+l_{R^{+}}\left({\mathcal{G}}_5\right)+l_{\left\{1\right\}}\left({\mathcal{G}}_6\right)\\ \hfill +\frac{\mu }{2}{∥{\mathcal{G}}_1-\mathcal{A}{\times }_3\mathbf{M}+{\mathcal{H}}_1∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_2-\mathcal{A}+{\mathcal{H}}_2∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_3-{\mathcal{D}}_1{\mathcal{A}}_{\times 1}+{\mathcal{H}}_3∥}_F^2\\ \hfill +\frac{\mu }{2}{∥{\mathcal{G}}_4-{\mathcal{D}}_2{\mathcal{A}}_{\times 2}+{\mathcal{H}}_4∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_5-\mathcal{A}+{\mathcal{H}}_5∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_6-\mathcal{A}+{\mathcal{H}}_6∥}_F^2.\end{array}$$

(17)

(2) Update $\mathcal{A}$: In problem (17), the subproblem of $\mathcal{A}$ can become:

$$\begin{array}{c}\hfill \underset{\mathcal{A}}{min}\frac{\mu }{2}{∥{\mathcal{G}}_1-\mathcal{A}{\times }_3\mathbf{M}+{\mathcal{H}}_1∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_2-\mathcal{A}+{\mathcal{H}}_2∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_3-{\mathcal{D}}_1{\mathcal{A}}_{\times 1}+{\mathcal{H}}_3∥}_F^2\\ \hfill +\frac{\mu }{2}{∥{\mathcal{G}}_4-{\mathcal{D}}_2{\mathcal{A}}_{\times 2}+{\mathcal{H}}_4∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_5-\mathcal{A}+{\mathcal{H}}_5∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_6-\mathcal{A}+{\mathcal{H}}_6∥}_F^2\end{array}.$$

(18)

In problem (18), it can be directly solved by:

$${\mathcal{A}}_{\times 3}{\mathbf{M}}^T\mathbf{M}+\sum _{k=1}^2{\mathcal{A}}_{\times k}\left({\mathcal{D}}_k^T{\mathcal{D}}_k\right)+3\mathcal{A}=\mathcal{K}$$

(19)

where

$$\begin{array}{c}\hfill \mathcal{K}={({\mathcal{G}}_1+{\mathcal{H}}_1)}_{\times 3}{\mathbf{M}}^T+{({\mathcal{G}}_3+{\mathcal{H}}_3)}_{\times 1}{\mathcal{D}}_1^T+{({\mathcal{G}}_4+{\mathcal{H}}_4)}_{\times 2}{\mathcal{D}}_2^T+({\mathcal{G}}_2+H_2+{\mathcal{G}}_5+{\mathcal{H}}_5+{\mathcal{G}}_6+{\mathcal{H}}_6)\end{array}$$

(20)

(19) is equal as

$${\mathbf{A}}_{\left(3\right)}^T{\mathbf{M}}^T\mathbf{M}+{\mathbf{CA}}_{\left(3\right)}^T+3{\mathbf{A}}_{\left(3\right)}^T={\mathbf{K}}_{\left(3\right)}^T$$

(21)

where $\mathbf{C}=[\mathbf{I}\otimes {\mathbf{D}}_1^T{\mathbf{D}}_1]+[{\mathbf{D}}_2^T{\mathbf{D}}_2\otimes \mathbf{I}]$. C is a matrix with a locking loop structure. It can be easily concluded that ${\mathbf{M}}^T\mathbf{M}$ is a symmetric matrix. According to [39], 2-D fast Fourier transformation and SVD is employed to solve $\mathbf{C}$ and ${\mathbf{M}}^T\mathbf{M}$, respectively. Then it can be obtained that $\mathbf{C}={\mathbf{F}}_{\mathbf{2}}^{\mathbf{T}}{\mathbf{\Psi }}_{\mathbf{2}}{\mathbf{F}}_{\mathbf{2}}$ and ${\mathbf{M}}^{\mathbf{T}}\mathbf{M}={\mathbf{U}}_{\mathbf{2}}\sum {\mathbf{U}}_{\mathbf{2}}^{\mathbf{T}}$, in which ${\mathbf{F}}_{\mathbf{2}}$ is a 2-D discrete Fourier transformation matrix. Then the update rule of $\mathcal{A}$ is

$$\mathcal{A}=Fol{d}_3\left({\left[{\mathbf{F}}_2^T((1\otimes {\mathbf{T}}_2)\odot \left({\mathbf{F}}_2{K}_{\left(3\right)}^T{\mathbf{U}}_2\right)){\mathbf{U}}_2^T\right]}^T\right)$$

(22)

in which ${\mathbf{T}}_2=(diag\left({\mathsf{\Psi }}_2\right),diag\left({\mathsf{\Psi }}_2\right),\dots ,diag\left({\mathsf{\Psi }}_2\right))+{(diag(\sum ),diag(\sum ),\dots ,diag(\sum ))}^T$.

(3) Update the intermediate variables ${\mathcal{G}}_i$:

$$\begin{array}{c}\hfill \underset{{\mathcal{G}}_i}{min}\frac{1}{2}{∥{\mathcal{G}}_1-Y∥}_F^2+l_{R^{+}}\left({\mathcal{G}}_5\right)+l_{\left\{1\right\}}\left({\mathcal{G}}_6\right)+\frac{\mu }{2}{∥{\mathcal{G}}_1-\mathcal{A}{\times }_3\mathbf{M}+{\mathcal{H}}_1∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_2-\mathcal{A}+{\mathcal{H}}_2∥}_F^2\\ \hfill +\frac{\mu }{2}{∥{\mathcal{G}}_3-{\mathcal{D}}_1{\mathcal{A}}_{\times 1}+{\mathcal{H}}_3∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_4-{\mathcal{D}}_2{\mathcal{A}}_{\times 2}+{\mathcal{H}}_4∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_5-\mathcal{A}+{\mathcal{H}}_5∥}_F^2+\frac{\mu }{2}{∥{\mathcal{G}}_6-\mathcal{A}+{\mathcal{H}}_6∥}_F^2.\end{array}$$

(23)

It can be easily concluded that the update formula of ${\mathcal{G}}_1$

$${\mathcal{G}}_1\leftarrow \frac{1}{1+\mu }\left[\mathcal{Y}+\mu \left(\mathcal{A}{\times }_3\mathbf{M}-{\mathcal{H}}_1\right)\right].$$

(24)

According to [40], the update rules related to ${\mathcal{G}}_2$ can be expressed as:

$${\mathcal{G}}_2=\left\{\begin{array}{c}\frac{\left|\right|{\mathcal{H}}_2(i,j,:\left)\right|{|}_2-W_S(i,j)\frac{{\lambda }_s}{\mu }}{\left|\right|{\mathcal{H}}_2(i,j,:)|{|}_2}{\mathcal{H}}_2(i,j,:),if{W}_S(i,j)\frac{{\lambda }_s}{\mu }<\left|\right|{\mathcal{H}}_2(i,j,:\left)\right|{|}_2,\\ 0,\mathrm{otherwise}.\end{array}\right.$$

(25)

Observing the problem (19), it can be seen the intermediate variables ${\mathcal{G}}_2$, ${\mathcal{G}}_3$, and ${\mathcal{G}}_4$ have similar forms. Hence, their solving rules are also similar, then the closed solutions of ${\mathcal{G}}_3$ and ${\mathcal{G}}_4$ can be directly obtained as:

$${\mathcal{G}}_3=\left\{\begin{array}{c}\frac{\left|\right|{\mathcal{C}}_1(i,j,:)|{|}_2-W_1(i,j)\frac{\lambda }{\mu }}{\left|{|\mathcal{C}}_1(i,j,:)\right|{|}_2}{\mathcal{H}}_3(i,j,:),if{W}_1(i,j)\frac{\lambda }{\mu }<\left|\right|{\mathcal{C}}_1(i,j,:)|{|}_2,\\ 0,\mathrm{otherwise}.\end{array}\right.$$

(26)

where ${\mathcal{C}}_1={\mathcal{D}}_1{\mathcal{A}}_{\times 1}+{\mathcal{H}}_3$.

$${\mathcal{G}}_4=\left\{\begin{array}{c}\frac{\left|\right|{\mathcal{C}}_2(i,j,:)|{|}_2-W_2(i,j)\frac{\lambda }{\mu }}{\left|\right|{\mathcal{C}}_2(i,j,:)|{|}_2}{\mathcal{H}}_4(i,j,:),if{W}_2(i,j)\frac{\lambda }{\mu }<\left|\right|{\mathcal{C}}_2(i,j,:)|{|}_2,\\ 0,\mathrm{otherwise}.\end{array}\right.$$

(27)

where ${\mathcal{C}}_2={\mathcal{D}}_2{\mathcal{A}}_{\times 2}+{\mathcal{H}}_4$.

It is worth noting that both $W_s$, $W_1$ and $W_2$ are weight matrixes, which are used to increase the degree of the sparsity. Therefore, their update rules can be described as:

$$W_s(i,j,k)=\frac{1}{|{\mathcal{G}}_2(i,j,k)|+\epsilon },$$

(28)

$$W_1(i,j)=\frac{1}{\left|\right|{\mathcal{C}}_1(i,j,:)|{|}_2+\epsilon },$$

(29)

$$W_2(i,j)=\frac{1}{\left|\right|{\mathcal{C}}_2(i,j,:)|{|}_2+\epsilon },$$

(30)

in which $\epsilon$ is a small positive number.

Besides, ${\mathcal{G}}_5$ and ${\mathcal{G}}_6$ are aiming at ensuring ANC and ASC of abundance. Therefore, the update rule of them can be written as:

$${\mathcal{G}}_5\leftarrow max\left(\mathcal{A}-{\mathcal{H}}_5,0\right),$$

(31)

$$\begin{array}{c}{\mathcal{G}}_6\leftarrow \left(\mathcal{A}-{\mathcal{H}}_6\right)+\mathrm{repmat}(\mathbf{N},I\left(1\right),1)),\\ \mathbf{N}=\frac{1}{P}\left[\mathbf{1}-\mathrm{sum}\left((\mathcal{A}-{\mathcal{H}}_6)\right)\right].\end{array}$$

(32)

(4) Update the Lagrange multipliers ${\mathcal{H}}_i$:

$$\begin{array}{c}{\mathcal{H}}_1={\mathcal{H}}_1-{\mathcal{A}}_{\times 3}\mathbf{M}+{\mathcal{G}}_1,\hfill \\ {\mathcal{H}}_2={\mathcal{H}}_2-\mathcal{A}+{\mathcal{G}}_2,\hfill \\ {\mathcal{H}}_3={\mathcal{H}}_3-{\mathcal{D}}_1{\mathcal{A}}_{\times 1}+{\mathcal{G}}_3,\hfill \\ {\mathcal{H}}_4={\mathcal{H}}_4-{\mathcal{D}}_2{\mathcal{A}}_{\times 2}+{\mathcal{G}}_4,\hfill \\ {\mathcal{H}}_5={\mathcal{H}}_5-\mathcal{A}+{\mathcal{G}}_5,\hfill \\ {\mathcal{H}}_6={\mathcal{H}}_6-\mathcal{A}+{\mathcal{G}}_6.\hfill \end{array}$$

(33)

In general, the framework of the algorithm WSCTF proposed in this paper is summarized in Algorithm 1.

Algorithm 1 The proposed method SCLT

Input:       $\mathcal{Y}\in {\mathbb{R}}^{I\times J\times K}$—the mixing hyperspectral data; P—the number of endmembers;       Parameter ${\lambda }_S$, and $\lambda$. 1. Update the endmember matrix $\mathbf{M}$ by (14); 2. Update the abundance tensor $\mathcal{A}$ via (22); 3. Update the variables ${\mathcal{G}}_1$ by (24), ${\mathcal{G}}_2$ by (25), ${\mathcal{G}}_3$ by (26), ${\mathcal{G}}_4$ by (27), ${\mathcal{G}}_5$ and ${\mathcal{G}}_6$ via (31) and (32); 4. Update the Lagrangian multipliers ${\mathcal{H}}_1$, ${\mathcal{H}}_2$, ${\mathcal{H}}_3$, ${\mathcal{H}}_4$, ${\mathcal{H}}_5$, ${\mathcal{H}}_6$ by (33). Output:       $\mathcal{A}$—the abundance tensor; $\mathbf{M}$—the endmember matrix.

4. Experiments

To verify the effectiveness of the WSCTF algorithm proposed in this paper, this section conducts related experiments on two sythetic data and real data sets of Cuprite, Jasper Ridge and Samon, and discusses the experimental results in detail. Since the proposed WSCTF is based on a non-negative tensor decomposition algorithm, the five comparison algorithms used here are: ULTRA-V [30], SULRSR-TV [24], MV-NTF [26], NMF-QMV [41], NL-TSUn [33]. It is worth noting that ULTRA-V, MV-NTF and NL-TSUn are all based on NTF, while SULRSR-TV and NMF-QMV are based on NMF. In addition, the SULRSR-TV algorithm also uses the TV regularization to constrain the local smooth structure of the abundance. The purpose of setting up these comparative experiments is to comprehensively verify the effectiveness of our algorithm.

Spectral angle distance (SAD) and root-mean-square error (RMSE) are the most widely used evaluation indicators in unmixing tasks. Among them, SAD is used to describe the spectral angular distance between the estimated endmember and the ground truth. Correspondingly, the distance of the estimated abundance and the ground truth is described by RMSE. They can be expressed in mathematical formulas as:

$$\begin{array}{c}\mathrm{RMSE}=\sqrt{{\left|\mathcal{A}-\widehat{\mathcal{A}}\right|}^2},\\ \mathrm{SAD}={cos}^{-1}\left(\frac{{\widehat{\mathbf{M}}}^T\mathbf{M}}{∥\widehat{\mathbf{M}}∥∥\mathbf{M}∥}\right)\end{array}$$

(34)

in which $\widehat{\mathcal{A}}$ and $\widehat{\mathbf{M}}$ are the ground truth of abundances and endmembers, respectively. Besides, to ensure the fairness of comparison, all experiments use the same experimental platform and use the vertex component analysis (VCA) and fully constrained least squares (FCLS) algorithm to initialize the endmembers and abundances. In addition, it is worth noting that the average and variance of the results of all experiments have been presented after 20 runs, which can illustrate the stability of algorithms.

4.1. Sythetic Data

Firstly, experiments are conducted on two synthetic data. All pure endmember spectra in the data are randomly sampled from the USGS spectral library [42].

(1) Dataset 1 (DS1): Set the number of endmembers $p=5$, the size of the HSI is $75\times 75\times 224$, where 224 is the number of spectral band. The specific data synthesis method is similar to the paper [34]. It is worth noting that each pixel here should comply with the LMM and meet the abundance ANC and ASC at the same time. Figure 1a shows the false color image of DS1.

(2) Dataset 2 (DS2): DS2 can form a control group with DS1. Here, the method described in [34] is used to synthesize DS2. Set the image size of DS2 to $90\times 90\times 224$, where 224 is the spectral band. Figure 1b shows the false color image of DS2. According to the above, two clean synthetic data can be obtained. Subsequently, 10 dB, 20 dB and 30 dB Gaussian noise are added to DS1 and DS2 respectively to form six data interfered by different noises. All experiments are all performed on these six data, which can prove the effectiveness of algorithms under different image sizes, different noise conditions, and different numbers of endmembers.

The SAD and RMSE results of the proposed WSCTF and five comparison algorithms on these six data are shown in

Table 1. RMSE results on the sythetic data.

Data	SNR	SULRSR-TV	NMF-QMV	MV-NTF	NL-TSUn	ULTRA-V	WSCTF
	10 dB	0.0995 ± 0.0049	0.0960	0.1906 ± 0.0041	0.1447 ± 0.0621	0.2339 ± 0.0031	0.0827 ± 0.0019
DS1	20 dB	0.0805 ± 0.0077	0.0723	0.1500 ± 0.0420	0.1024 ± 0.0352	0.1701 ± 0.0214	0.0682 ± 0.0241
	30 dB	0.0797 ± 0.0051	0.0624	0.1456 ± 0.0164	0.0828 ± 0.0614	0.2291 ± 0.0510	0.0621 ± 0.0146
	10 dB	0.1433 ± 0.0143	0.1200	0.1256 ± 0.0091	0.1326 ± 0.0312	0.1980 ± 0.0312	0.1126 ± 0.0321
DS2	20 dB	0.1366 ± 0.0106	0.1193	0.1580 ± 0.0217	0.1273 ± 0.0318	0.1876 ± 0.0021	0.1004 ± 0.0021
	30 dB	0.1344 ± 0.0054	0.0942	0.1005 ± 0.0032	0.0987 ± 0.0021	0.1690 ± 0.0410	0.0911 ± 0.0039

and

Table 2. SAD results on the Sythetic data.

Data	SNR	SULRSR-TV	NMF-QMV	MV-NTF	NL-TSUn	ULTRA-V	WSCTF
	10 dB	0.1115 ± 0.0139	0.3624	0.5127 ± 0.0164	0.1326 ± 0.0052	0.1427 ± 0.0216	0.0920 ± 0.0031
DS1	20 dB	0.0453 ± 0.0014	0.0867	0.4023 ± 0.0121	0.0952 ± 0.0143	0.1214 ± 0.0126	0.0434 ± 0.0108
	30 dB	0.0134 ± 0.0004	0.0920	0.3210 ± 0.0114	0.0623 ± 0.2432	0.0737 ± 0.0601	0.0126 ± 0.0245
	10 dB	0.1977 ± 0.0928	0.2325	0.5232 ± 0.0088	0.4329 ± 0.0621	0.3172 ± 0.0241	0.1920 ± 0.0314
DS2	20 dB	0.1688 ± 0.0794	0.0861	0.5161 ± 0.0128	0.4061 ± 0.0051	0.3202 ± 0.0152	0.0827 ± 0.0021
	30 dB	0.1423 ± 0.0366	0.0184	0.3237 ± 0.0672	0.2846 ± 0.0601	0.3110 ± 0.0031	0.0094 ± 0.0017

. According to the results of Table 1 and Table 2, we can analyze the results: (1) On the DS1 and DS2 data, as the noise continues to increase, the performance of all algorithms, including the proposed WSCTF, is declining. However, the performance of WSCTF method changes less under the influence of different noises and the performance is more stable than other methods, since it fully considers the spatial and spectral similarity. (2) MV-NTF performs poorly on both data, because MV-NTF is based on the NTF basic algorithm, and other algorithms are mostly improved methods on it. (3) Since DS2 is more complicated than DS1 and the number of endmembers is more, all algorithms perform slightly worse on DS2. The ULTRA-V algorithm is mainly proposed for the variability with endmembers, but the absence of variable endmembers in the settings of DS1 and DS2 leads to poor unmixing performance. (4) The lack of pure pixels in the DS2 is the main reason for the low performance of the SULRSR-TV method. SULRSR-TV assumed that there must be pure pixels in the image when unmixing. In summary, compared with other comparison algorithms, the proposed WSCTF has certain advantages under the conditions of different noises, different image sizes and different numbers of endmembers. The above analysis can also be verified in the visualization results of abundance. As shown in Figure 2, the abundance and groundtruth obtained by all algorithms are displayed on DS1 with 20 dB. The abundance map estimated by the WSCTF method is obviously better than other comparison algorithms. Due to the large number of endmembers in DS2 with 20 dB, only five abundances are shown in Figure 3 here.

4.2. Real Datasets

Secondly, experiments are conducted on three real datasets.

(1) Jasper Ridge Data: Jasper Ridge data is a data set widely used to test the performance of unmixing algorithms, covering the Jasper Ridge Natural Reserve in California, USA. The data contains 224 spectral bands from 300 nm to 2500 nm, and its spectral resolution is 9.46 nm. Figure 1c is a pseudo-color image of Jasper Ridge data. Here we use an image with a size of 100 × 100 pixels to verify the performance of all algorithms. Because some bands are severely interfered by moisture and noise during the acquisition process, 188 bands are finally used for testing after removing these bands. According to a prior, there are four end members in this area, including soil, water, road and trees.

Table 3. SAD results on the Jasper ridge data.

	SULRSR-TV	NMF-QMV	MV-NTF	NL-TSUn	ULTRA-V	WSCTF
Tree	0.2012 ± 0.0187	0.2571	0.2416 ± 0.0136	0.2081 ± 0.0306	0.0502 ± 0.0003	0.0495 ± 0.0476
Soil	0.2336 ± 0.0011	1.1890	0.2601 ± 0.0089	0.2334 ± 0.0001	0.1422 ± 0.0342	0.1071 ± 0.1660
Water	0.6194 ± 0.1761	0.1386	0.1586 ± 0.0956	0.4929 ± 0.2715	0.5641 ± 0.0013	0.3977 ± 0.0072
Road	0.2943 ± 0.0722	0.1840	0.4544 ± 0.0706	0.3665 ± 0.1179	0.0362 ± 0.0033	0.2461 ± 0.0303
Mean	0.3863 ± 0.0357	0.4421	0.3024 ± 0.0396	0.3723 ± 0.0482	0.3002 ± 0.0010	0.2406 ± 0.0434

shows the SAD results of the endmembers extracted by the proposed WSCTF compared with ground truth. The bold font indicates the best result. It can be concluded that the WSCTF method proposed in this paper has great advantages in the extraction of two endmembers, and the final SAD result is also optimal. In addition, the endmember spectra extracted by all algorithms are shown in Figure 4. Compared with the reference spectrum, it can be seen that the proposed algorithm has considerable advantages in endmember extraction compared with other methods. Figure 5 is the comparison result of all abundances, which can be concluded that the WSCTF algorithm is the closest to the reference. In addition, compared with the distribution of ground objects in the original image in Figure 1c, the accuracy of the proposed algorithm can also be inferred.

(2) Samon Data: According to [43], the samon data set is collected by the Florida Environment Research Institute through the Samon sensor. The data contains 156 bands and the spectral resolution is 3.13 nm. Here we are experimenting on an image with a size of 100 × 100 pixels. There are three types of ground features in this area, including water, soil and trees. Figure 1d is a pseudo-color image of Samon data. The SAD and its variance results extracted by the six methods on the Samon data are introduced in

Table 4. SAD results on the Samson data.

	SULRSR-TV	NMF-QMV	MV-NTF	NL-TSUn	ULTRA-V	WSCTF
Tree	0.0247 ± 0.0121	0.0337	0.0433 ± 0.0129	0.0288 ± 0.0161	0.0340 ± 0.0272	0.0201 ± 0.0024
Soil	0.0495 ± 0.0011	0.0732	0.0953 ± 0.0011	0.0496 ± 0.0001	0.0566 ± 0.0090	0.0685 ± 0.0102
Water	0.1299 ± 0.0004	1.4831	0.2810 ± 0.0050	0.1289 ± 0.0003	0.2401 ± 0.0221	0.0675 ± 0.0021
Mean	0.0818 ± 0.0019	0.8575	0.1733 ± 0.0036	0.0825 ± 0.0026	0.1438 ± 0.0154	0.0611 ± 0.0105

. Obviously, it can be concluded that all unmixing algorithms, including the WSCTF proposed in this paper, perform slightly worse on the endmember of water. Since the edge of water when it is mixed with the land is not easy to determine, the abundance map of the water body endmember is greatly disturbed. Similar conclusions can also be drawn from the visualization results in Figure 6 and Figure 7. However, it is gratifying that the method proposed in this paper has a great advantage in the separation of the endmember of water compared to other algorithms, which proves the superiority and stability of WSCTF proposed in this paper.

(3) Cuprite Data: This data set is also widely used in unmixing task. Its coverage area is a mining area, so there are many highly mixed pixels. Here, a subgraph with a size of 200 × 150 × 188 is intercepted to test the effectiveness of the algorithm, where 188 is the number of bands. Figure 1e is a pseudo-color image of Cuprite data. According to a prior, there are 12 pure mineral endmembers in this area. Unlike the Samon and Jasper Ridge data, the Cuprite data does not provide true values, so the ground truth of endmembers here are set using the method in [34].

Cuprite data mainly contains 12 minerals, namely Alunite, Kaolinite #1, Kaolinite #2, Sphene, Andradite, Nontronite, Muscovite, Pyrope, Chalcedony, Dumortierite, Montmorillonite, and Buddingtonite. Table 3 shows the SAD comparison results of the endmembers extracted by the six algorithms and the reference endmembers. As shown in

Table 5. SAD results on the Aviris Cuprite data.

	SULRSR-TV	NMF-QMV	MV-NTF	NL-TSUn	ULTRA-V	WSCTF
Alunite	0.0966 ± 0.0739	0.0860	0.0748 ± 0.0875	0.0783 ± 0.0243	0.1897 ± 0.0722	0.0735 ± 0.0134
Andradite	0.1353 ± 0.0157	0.1212	0.0661 ± 0.0207	0.0958 ± 0.0376	0.0971 ± 0.0323	0.0631 ± 0.0369
Buddingtonite	0.0699 ± 0.0174	0.0797	0.0718 ± 0.0193	0.0982 ± 0.0183	0.1097 ± 0.0193	0.0669 ± 0.0217
Chalcedony	0.1482 ± 0.0230	0.1428	0.1652 ± 0.0304	0.1253 ± 0.0444	0.1212 ± 0.0307	0.0829 ± 0.0317
Dumortierite	0.0855 ± 0.0152	0.1108	0.0984 ± 0.0208	0.1709 ± 0.0240	0.0936 ± 0.0158	0.0944 ± 0.0258
Kaolinite#1	0.0718 ± 0.0155	0.0795	0.0841 ± 0.0070	0.0738 ± 0.0095	0.1300 ± 0.0137	0.0859 ± 0.0104
Kaolinite#2	0.0678 ± 0.0187	0.1014	0.0786 ± 0.0098	0.0625 ± 0.0126	0.0391 ± 0.0160	0.0638 ± 0.0154
Montmorillonite	0.0504 ± 0.0132	0.1114	0.0535 ± 0.0065	0.0720 ± 0.0124	0.0508 ± 0.0134	0.0501 ± 0.0143
Muscovite	0.1066 ± 0.0450	0.0811	0.2518 ± 0.0498	0.0925 ± 0.0457	0.1257 ± 0.0389	0.1719 ± 0.0461
Nontronite	0.0730 ± 0.0193	0.0715	0.1086 ± 0.0111	0.0813 ± 0.0156	0.0835 ± 0.0246	0.0702 ± 0.0271
Pyrope	0.1876 ± 0.0521	0.1812	0.1431 ± 0.0598	0.1501 ± 0.0464	0.1737 ± 0.0535	0.1326 ± 0.0574
Sphene	0.0663 ± 0.0564	0.1879	0.0700 ± 0.0673	0.1642 ± 0.0137	0.2306 ± 0.0506	0.0744 ± 0.0829
Mean	0.1042 ± 0.0062	0.1190	0.1176 ± 0.0135	0.1054 ± 0.0066	0.1317 ± 0.0089	0.0925 ± 0.0134

, the algorithm proposed in this paper has obvious advantages in each endmember and average endmember results, and there are 7 optimal results. In addition, Figure 8 and Figure 9 are the visualized results of endmembers and abundance obtained after unmixing, respectively. In Figure 8, compared with the other five methods, it is obvious that the endmember spectrum obtained by the WSCTF algorithm is closer to the reference spectrum. The abundance maps are shown in Figure 9. Compared with [34], it can be seen that the proposed WSCTF can also obtain more accurate abundance results.

4.3. Parameter Analysis

The WSCTF algorithm proposed in this paper involves two sparse parameters $\lambda$ and ${\lambda }_s$. Here we take Jasper Ridge data as an example for parameter setting discussion. First, set $\lambda$ as 0.015, and ${\lambda }_s$ to {0, 1 × 10${}^{-3}$, 2 × 10${}^{-3}$, 3 × 10${}^{-3}$, 4 × 10${}^{-3}$, 5 × 10${}^{-3}$, 6 × 10${}^{-3}$, 7 × 10${}^{-3}$, 8 × 10${}^{-3}$, 9 × 10${}^{-3}$, 1 × 10${}^{-2}$}. The setting of 0 can verify the advantages of the regularization proposed by this algorithm. It can be seen from Figure 10a that when ${\lambda }_s$=0.015, WSCTF achieves the best results on Jasper Ridge data. Then set ${\lambda }_s$ as 0.015, $\lambda$ to {0, 1 × 10${}^{-3}$, 2 × 10${}^{-3}$, 3 × 10${}^{-3}$, 4 × 10${}^{-3}$, 5 × 10${}^{-3}$, 6 × 10${}^{-3}$, 7 × 10${}^{-3}$, 8 × 10${}^{-3}$, 9 × 10${}^{-3}$, 1 × 10${}^{-2}$}. The result is shown in Figure 10b, and when $\lambda$ = 0.02, WSCTF achieves the best results on Jasper Ridge data. Similarly, apply the same method to other data sets to ensure the optimal parameters. Finally, the results of the parameters set in this paper are shown in

Table 6. Parameter settings.

Data	${\mathit{\lambda }}_{\mathit{s}}$	$\mathit{\lambda }$
DS1	0.015	0.010
DS2	0.015	0.010
Jasper	0.02	0.015
Samon	0.015	0.010
Cuprite	0.010	0.010

. It is worth noting that the same parameters are set here for DS1 and DS2 under different noise conditions since the synthetic data is much simpler than the real data.

4.4. Complexity Analysis

Here, three real datasets are taken as examples to analyze the time complexity of the proposed WSCTF and other comparison algorithms. All experiments in this paper are performed on Intel (R) Xeon (R) CPU E5-2697 v3 at 2.60 GHz and 64G RAM platform, which can guarantee the fairness of comparison. The running time of all algorithms is presented in

Table 7. Running Time (Seconds) of All Algorithms on Three Real Datasets.

	SULRSR-TV	NMF-QMV	MV-NTF	NL-TSUn	ULTRA-V	WSCTF
Jasper Ridge Data	6	10	170	55	26	5
Samon Data	5	9	135	29	14	4
Cuprite Data	34	293	904	194	209	17

. It can be seen from the table that the running time of MV-NTF will explode as the image size and the number of endmembers increase. The running time of the WSCTF algorithm proposed in this paper is very short, because it does not involve complex tensor operations. In summary, WSCTF also has certain advantages in terms of time complexity compared to other algorithms.

5. Conclusions

In this paper, a regularization for exploring the group sparse structure in difference images of HSI is proposed, which not only maintains the structural features of the raw data space in the abundance tensor, but also keeps the continuity in the spectral direction. First, the raw image can be constrained by the TV regularization to obtain the difference images, which are 3-D tensors. Then, the difference images with group sparse structures are separately reshaped in two spatial dimensions. In particular, the $L_{2,1}$ norm is used to constrain the sparsity, which can mine the global sparse connected structure. In addition, weights are also used to improve the sparsity. The final experiments prove the superiority of the proposed WSCTF compared with other state-of-the-art methods.