Robust Dual Spatial Weighted Sparse Unmixing for Remotely Sensed Hyperspectral Imagery

Abstract

Sparse unmixing plays a crucial role in the field of hyperspectral image unmixing technology, leveraging the availability of pre-existing endmember spectral libraries. In recent years, there has been a growing trend in incorporating spatial information from hyperspectral images into sparse unmixing models. There is a strong spatial correlation between pixels in hyperspectral images (that is, the spatial information is very rich), and many sparse unmixing algorithms take advantage of this to improve the sparse unmixing effect. Since hyperspectral images are susceptible to noise, the feature separability of ground objects is reduced, which makes most sparse unmixing methods and models face the risk of degradation or even failure. To address this challenge, a novel robust dual spatial weighted sparse unmixing algorithm (RDSWSU) has been proposed for hyperspectral image unmixing. This algorithm effectively utilizes the spatial information present in the hyperspectral images to mitigate the impact of noise during the unmixing process. For the proposed RDSWSU algorithm, which is based on ${\mathsf{\ell }}_1$ sparse unmixing framework, a pre-calculated superpixel spatial weighting factor is used to smooth the noise, so as to maintain the original spatial structure of hyperspectral images. The RDSWSU algorithm, which builds upon the ${\mathsf{\ell }}_1$ sparse unmixing framework, employs a pre-calculated spatial weighting factor at the superpixel level. This factor aids in noise smoothing and helps preserve the inherent spatial structure of hyperspectral images throughout the unmixing process. Additionally, another spatial weighting factor is utilized in the RDSWSU algorithm to capture the local smoothness of abundance maps at the sub-region level. This factor helps enhance the representation of piecewise smooth variations within different regions of the hyperspectral image. Specifically, the combination of these two spatial weighting factors in the RDSWSU algorithm results in an enhanced sparsity of the abundance matrix. The RDSWSU algorithm, which is a sparse unmixing model, offers an effective solution using the alternating direction method of multiplier (ADMM) with reduced requirements for tuning the regularization parameter. The proposed RDSWSU method outperforms other advanced sparse unmixing algorithms in terms of unmixing performance, as demonstrated by the experimental results on synthetic and real hyperspectral datasets.

1. Introduction

Hyperspectral imaging spectrometers have the ability to capture images of specific regions in the electromagnetic spectrum, enabling the detection of objects that are unidentifiable in wide-band imaging remote sensing [1]. The use of hyperspectral remote sensing images enables the acquisition of data containing abundant information about the Earth surface, thereby providing the potential for achieving quantitative and refined ground observation [2]. Hyperspectral remote sensing imagery has been widely utilized in diverse fields, encompassing marine exploration [3], environmental monitoring [4], forestry monitoring [5], geological exploration [6], modern military affairs [7], and precision agriculture [4]. Owing to the limited spatial resolution of imaging sensors and the intricate variability of ground features, hyperspectral remote sensing images often exhibit mixed pixels. This phenomenon poses challenges in accurately interpreting the information related to ground objects. This seriously restricts the development and application of quantitative analysis in hyperspectral remote sensing [8]. To cope with the problem of mixed pixels, spectral unmixing [9] techniques have been introduced to effectively identify the spectral characteristics of different ground features in each pixel, called endmembers, and estimate their corresponding percentages, called abundance. This is a key method for analyzing hyperspectral images.

The linear mixture model (LMM) becomes an active research line in hyperspectral imaging analysis owing to its simplicity and mathematical tractability [10]. LMM assumes that the spectral signature of a mixed pixel can be represented as a linear combination of the spectra of multiple endmembers. Building upon LMM, various unsupervised spectral unmixing methods have been developed, including geometry-based approaches [11,12,13] and statistics-based approaches [14,15]. However, most of these methods require extracting endmembers directly from the given scene, which can lead to obtaining virtual endmembers that have no physical meaning from mixed hyperspectral images [16]. Sparse unmixing methods have emerged as a solution to tackle this issue [17,18,19]. These methods adopt a semi-supervised approach, aiming to identify the optimal subset of entries from a known endmember spectral library that effectively represents the mixed pixels in a hyperspectral image. By leveraging sparsity constraints, they seek to achieve a sparse representation of the mixed pixels using a reduced set of endmembers [20,21].

In the context of sparse unmixing, the abundance vector associated with a mixed pixel is characterized by containing far fewer elements than the total number of spectral signatures present in a given spectral library. As a result, the abundance vector exhibits sparsity as a prominent feature [22]. To enhance the sparsity of solutions, several typical methods have been proposed. For example, the renowned SUnSAL [17] algorithm utilizes the ${\mathsf{\ell }}_1$ regularizer to promote sparsity in the fractional abundances, yielding promising outcomes. Moreover, the collaborative SUnSAL algorithm, known as CLSUnSAL [23], incorporates the ${\mathsf{\ell }}_{2,1}$ regularizer to leverage collaborative row sparsity across all pixels. In order to further punish the abundance coefficient, some studies introduce the weight factor into the sparse unmixing models [24,25,26]. The DRSU [27] algorithm employs dual weight factors to efficiently boost the sparsity of the fractional abundance. While the aforementioned sparse unmixing algorithms demonstrate promising potential, they primarily emphasize the utilization of spectral information from hyperspectral images while overlooking the incorporation of spatial information, which impacts the accuracy of unmixing.

Spatial-contextual information is a valuable resource that can be harnessed to enhance the effectiveness of unmixing techniques in practical applications [28]. The incorporation of spatial information into sparse unmixing models has become a prominent approach, leading to the development of numerous sparse unmixing algorithms that leverage spatial information. For example, the SUnSAL-TV [29] algorithm utilizes anisotropic total variation regularization to encourage the smoothing of spatial neighborhoods in images. The NLSU [30] algorithm utilizes a non-local spatial regularization technique to identify and capture resemblances in patterns and structures present in the abundance map. Furthermore, a S${}^2$WSU method is introduced, incorporating spectral and spatial weight factors to enhance the sparsity of solutions within the ${\mathsf{\ell }}_1$ sparse regularization framework. The proposed DRSU-TV [31] method is developed, which utilizes a combination of double weighted factors and TV regularization to improve the sparsity of solutions while simultaneously achieving smoothness in the abundance map. A novel approach, based on superpixels and graph regularization, has been proposed for sparse unmixing [32]. This method leverages spectral weighting factors to encourage sparsity in the solutions, while utilizing a graph regularization term to enhance the spatial correlation of the resulting images. These methods can enhance the accuracy of abundance estimation to a certain degree. Nevertheless, in scenarios where hyperspectral images are heavily corrupted by high levels of noise, these methods often struggle and achieving satisfactory accuracy in unmixing becomes challenging. For this problem, a self-supervised robust deep unmixing model (SSRDMF) is presented [33], which decreases the influence of Gaussian noise and sparse noise through the deep autoencoder, thereby enhancing unmixing capabilities. The spectral-spatial robust NMF unmixing model (SSRNMF) is introduced [34], leveraging double regularization norms to attain noise robustness, resulting in favorable outcomes. The SSNPNMF unmixing model [35] is presented, building a graph based on spatial-spectral information to modify the reconstruction weight values of pixels and neighboring pixels, mitigating noise effects.

In recent years, superpixel segmentation technology has gained widespread adoption in the field of hyperspectral image processing due to its effectiveness in enhancing data interpretation performance [36], such as hyperspectral image classification [37,38], spectral unmixing [39], denoising [40], fusion [41] and other fields. For example, a superpixel-based weighting method was proposed in [42] to exploit the spatial correlation of images. In [43], a sparse unmixing model incorporating superpixel reweighted low-rank constraint and TV regularization is introduced, aiming to promote the smoothness of the image.

In order to effectively deal with spectral unmixing task under low signal-to-noise ratio (SNR) [44], in view of the advantages of superpixel segmentation strategy, a novel robust dual spatial weighted sparse unmixing (RDSWSU) algorithm is proposed in this paper. The main objective of the proposed RDSWSU algorithm is to acquire more precise spatial context information in order to mitigate the adverse impact of noise on the unmixing process. Specifically, for the proposed RDSWSU algorithm, which is based on the ${\mathsf{\ell }}_1$ sparse regularization framework, on the one hand, the pre-calculated spatial weighting factor based on superpixel is used to smooth the noise. On the other hand, the utilization of a spatial weighting factor enhances the spatial coherence of the abundance map.

The efficient solution to the optimization problem of the proposed RDSWSU algorithm can be achieved by leveraging the ADMM technique [45]. To sum up, we can summarize the contribution of this paper. For the newly proposed RDSWSU method, based on the ${\mathsf{\ell }}_1$ sparse unmixing framework, it introduced two key enhancements: the utilization of pre-calculated superpixel spatial weighting factor and spatial neighborhood weighting factor. These improvements aimed to reduce the impact of noise on unmixing and capture the localized smoothness characteristics of the abundance map, respectively. The RDSWSU unmixing model, as proposed, demonstrates a notable enhancement in its ability to handle high levels of noise in challenging scenes. Furthermore, it effectively preserves the precise spatial information of images, thereby leading to improved accuracy in the unmixing process. Moreover, our proposed model has a simple formulation with only one regularization parameter, which makes it easy to tune.

The remaining sections of the paper are outlined as follows. Section 2 is dedicated to the introduction of the linear spectral unmixing model. In Section 3, a detailed description of our proposed robust dual spatial weighted sparse unmixing model is provided. Section 4 specifically presents the experimental results using synthetic hyperspectral data, whereas Section 5 focuses on showcasing experimental results using actual hyperspectral data. Finally, Section 6 concludes the paper by providing a summary, remarks, and suggestions for future research topics.

2. Linear Spectral Unmixing

Suppose we have an observed hyperspectral image denoted by $\mathbf{Y}\in {\mathbb{R}}^{L\times n}$, where n represents the total number of pixel vectors and L represents the number of spectral bands. Let $\mathbf{A}\in {\mathbb{R}}^{L\times m}$ be a comprehensive spectral library comprising m spectral signatures. The aim of sparse unmixing is to identify the most suitable subset of endmembers from the spectral library $\mathbf{A}$ in order to represent mixed pixels in a hyperspectral image. To sum up, the LMM can be represented by the following equation:

$$\mathbf{Y}=\mathbf{A}\mathbf{X}+\mathbf{N}\mathrm{s}.\mathrm{t}.\text{:}\mathbf{X}\ge 0,$$

(1)

where $\mathbf{X}\in {\mathbb{R}}^{m\times n}$ denotes the fractional abundance matrix, and $\mathbf{N}\in {\mathbb{R}}^{L\times n}$ denotes the noise or the model error. The abundance nonnegativity constraint (ANC) states that $\mathbf{X}$ must be greater than or equal to zero (i.e., $\mathbf{X}\ge 0$), while the abundance sum-to-one constraint (ASC) requires that the sum of the elements in each column of $\mathbf{X}$ equals one (i.e., ${\mathbf{1}}^T\mathbf{x}=\mathbf{1}$). In principle, the abundance matrix $\mathbf{X}$ should satisfy these two physical constraints. However, it should be noted that, for some reasons regarding ASC in [17], the ASC constraint is not explicitly enforced here, only ANC constraint is enforced. In practical applications, the number of endmembers constituting mixed pixels is usually much less than that of the large spectral library $\mathbf{A}$, which means that the abundance matrix $\mathbf{X}$ contains many zero values [46].

In practical scenarios, the number of endmembers present in mixed pixels is often significantly smaller compared to the size of the extensive spectral library $\mathbf{A}$. As a consequence, the abundance matrix $\mathbf{X}$ tends to have numerous zero values (i.e., the abundance matrix $\mathbf{X}$ exhibits sparsity.) [46]. Taking these factors into account, the problem of sparse unmixing can be formulated as follows.

$$\underset{\mathbf{X}}{min}\frac{1}{2}\left|\right|\mathbf{A}\mathbf{X}-{\mathbf{Y}\left|\right|}_F^2+{\lambda \left|\right|\mathbf{X}\left|\right|}_0\mathrm{s}.\mathrm{t}.\text{:}\mathbf{X}\ge 0,$$

(2)

where ${\parallel ·\parallel }_F$ represents the Frobenius norm, and $\lambda \ge 0$ serves as a regularization parameter in the equation. The minimization of the ${\mathsf{\ell }}_0$ norm in objective function (2) is a computationally difficult problem, known to be NP-hard [47]. However, an alternative approach is possible by leveraging the Restricted Isometric Property (RIP) [48,49], which allows us to relax the problem and approximate it using the ${\mathsf{\ell }}_1$ norm instead. The formal expression for this approximation is as follows.

$$\underset{\mathbf{X}}{min}\frac{1}{2}\left|\right|\mathbf{A}\mathbf{X}-{\mathbf{Y}\left|\right|}_F^2+{\lambda \left|\right|\mathbf{X}\left|\right|}_{1,1}\mathrm{s}.\mathrm{t}.\text{:}\mathbf{X}\ge 0,$$

(3)

where ${\left|\right|\mathbf{X}\left|\right|}_{1,1}={\sum }_{j=1}^n\left|\right|{\mathbf{X}}_j{\left|\right|}_1$, the jth column of the fractional abundance matrix $\mathbf{X}$ is represented by ${\mathbf{X}}_j$. Then the optimization problem in (3) can be solved efficiently by the SUnSAL algorithm under the ADMM.

Earlier approaches to sparse unmixing typically emphasized the examination of spectral information while neglecting the integration of spatial information within the image. To achieve greater consistency in abundance estimation, a Total Variation (TV) regularizer is utilized to encourage the piecewise smoothness of fractional abundance. The optimization problem is expressed in an alternative formulation.

$$\begin{array}{c}\hfill \underset{\mathbf{X}}{min}\frac{1}{2}\left|\right|\mathbf{A}\mathbf{X}-{\mathbf{Y}\left|\right|}_F^2+{\lambda \left|\right|\mathbf{X}\left|\right|}_{1,1}+{\lambda }_{TV}\mathbf{T}\mathbf{V}\left(\mathbf{X}\right)\mathrm{s}.\mathrm{t}.\text{:}\mathbf{X}\ge 0.\end{array}$$

(4)

where ${\lambda }_{TV}⩾0$ denotes the regularization parameter. $\mathbf{T}\mathbf{V}\left(\mathbf{X}\right)\equiv {\sum }_{\left\{p,q\right\}\in \varphi }\left|\right|{\mathbf{x}}_p-{\mathbf{x}}_q{\left|\right|}_1$ represents an extended version of nonisotropic TV [29], where ${\mathbf{x}}_p$ refers to the vector of abundances for pixel p, and ${\mathbf{x}}_q$ denotes the vector of abundances for pixel q. $\varphi$ represents the collection of neighboring pixels in the abundance map, including both horizontal and vertical directions. TV regularization can smooth the abundance map in some regions, but it usually leads to oversmoothing and boundary blurring in the abundance maps. Moreover, the estimation results of SUnSAL-TV rely on the selection of two regularization parameters.

3. The Robust Dual Spatial Weighted Sparse Unmixing Model

3.1. Proposed RDSWSU Method

Drawing on the achievements of the sparse unmixing technique utilizing reweighted ${\mathsf{\ell }}_1$-norm, it can effectively impose a penalty on the abundance coefficients. This paper presents a novel approach called the robust dual spatial weighted sparse unmixing algorithm (RDSWSU) to address the challenges of noise reduction in unmixing and improve the spatial smoothness of abundance mapping. The flowchart of the proposed RDSWSU method is shown in Figure 1.

To enhance the robustness of the unmixing algorithm against noise, the superpixel spatial weight ${\mathbf{H}}_1$ is constructed based on a coarse hyperspectral image, serving as a guidance for abundance estimation. It should be noted that the coarse hyperspectral image is reconstructed by a superpixel segmentation algorithm, which considered the similarity of pixels in terms of both spatial location and spectral characteristics. This reconstruction process resembles a mean filter and effectively suppresses noise in the image. Subsequently, a coarse-scale weight guide matrix is established predictively using the reconstructed hyperspectral image. ${\mathbf{H}}_1$ is then calculated row by row from the coarse-scale weight guide matrix, making it insensitive to noise. Overall, the RDSWSU algorithm does not involve any denoising process. Instead, it mitigates the impact of noise on unmixing by performing superpixel segmentation on the original hyperspectral image and constructing superpixel spatial weights based on the segmented image. The optimization objective of the proposed RDSWSU algorithm can be defined as follows.

$$\begin{array}{cc}\hfill \underset{\mathbf{X}}{min}& \frac{1}{2}\left|\right|\mathbf{A}\mathbf{X}-{\mathbf{Y}\left|\right|}_F^2+\lambda {\parallel \left({\mathbf{H}}_1{\mathbf{H}}_2\right)\odot \mathbf{X}\parallel }_{1,1}\mathrm{s}.\mathrm{t}.\text{:}\mathbf{X}\ge 0.\hfill \end{array}$$

(5)

where the operator ⊙ represents the elementwise multiplication of two variables. It should be noted that ${\mathbf{H}}_1\in {\mathbb{R}}^{m\times m}$ is the pre-calculated superpixel spatial weighting matrix, and ${\mathbf{H}}_2\in {\mathbb{R}}^{m\times n}$ is the spatial neighborhood weighting matrix.

As stated earlier, superpixel segmentation technology [43,50] refers to a method utilized for extracting spatial information from hyperspectral images. This technique allows for the analysis and characterization of the spatial structure within the image.

The superpixel spatial weighting matrix ${\mathbf{H}}_1$ is designed to preserve the original spatial contextual information of hyperspectral images, particularly in scenarios with high levels of noise. Its purpose is to ensure that the spatial details in the image are accurately represented, even under challenging conditions with significant noise interference. To obtain the superpixel spatial weighting matrix ${\mathbf{H}}_1$, the initial hyperspectral image $\mathbf{Y}\in {\mathbb{R}}^{L\times n}$ is segmented into g superpixel blocks using the SLIC algorithm [36,51]. Importantly, the SLIC algorithm is fast, simple to implement, and produces superpixels of good quality. It balances segmentation performance against computational complexity, serving as a useful and efficient image oversegmentation approach. This segmentation process divides the image into compact and homogeneous regions, each referred to as a superpixel block. By utilizing the SLIC algorithm, the hyperspectral image is effectively partitioned into g coherent and visually meaningful superpixels. Let ${\mathbf{Y}}_t\in {\mathbb{R}}^{L\times n_t}(t=1,\dots ,g,{\sum }_{t=1}^g{n}_t=n)$ represent the t-th superpixel block, which consists of $n_t$ pixels. Within this block, we can refer to any pixel (column) as ${\mathbf{y}}_k\in {\mathbb{R}}^{L\times 1}(k=1,\dots ,n_t)$. This notation allows us to refer to and manipulate individual pixels within specific superpixel blocks for further analysis or computations. Then, based on the segmentation result, the pixels within each superpixel block are averaged to obtain the coarse image $\tilde{\mathbf{Y}}\in {\mathbb{R}}^{L\times n}$. The statement is reformulated as follows.

$$\begin{array}{c}\hfill {\tilde{\mathbf{y}}}_k=\frac{1}{n_t}\sum _{k=1}^{n_t}{\mathbf{y}}_k\end{array}$$

(6)

In the given context, when referring to ${\tilde{\mathbf{y}}}_k$, it represents a column from the coarse image ${\tilde{\mathbf{Y}}}_t$ that corresponds to a superpixel block in ${\mathbf{Y}}_t$. It is important to note that all elements within a single superpixel block possess identical values.

The SUnSAL method [17] is utilized to derive the coarse fractional abundance matrix $\tilde{\mathbf{X}}$. The corresponding optimization model is expressed as follows.

$$\underset{\tilde{\mathbf{X}}}{min}\frac{1}{2}\left|\right|\mathbf{A}\tilde{\mathbf{X}}-\tilde{\mathbf{Y}}{\left|\right|}_F^2+\lambda \left|\right|\tilde{\mathbf{X}}{\left|\right|}_{1,1},\mathrm{s}.\mathrm{t}.\text{:}\tilde{\mathbf{X}}\ge 0.$$

(7)

where the low-resolution fractional abundance matrix $\tilde{\mathbf{X}}$ can keep the spatial information among the pixels. Finally, the superpixel-based weighting factor ${\mathbf{H}}_1$ can be calculated as follows:

$$\begin{array}{c}\hfill {\mathbf{H}}_{1,ij}=\frac{1}{\left|\right|\tilde{\mathbf{X}}(i,:){\left|\right|}_2+\epsilon },i=1,\dots ,m,j=1,\dots ,n\end{array}$$

where $\epsilon >0$ represents a small constant. The ith row vector of the abundance matrix $\tilde{\mathbf{X}}$ is represented as $\tilde{\mathbf{X}}(i,:)$, where $(i=1,\dots ,m)$. ${\mathbf{H}}_1$ is influenced by the coarse abundance matrix $\tilde{\mathbf{X}}$ and can increase the sparsity of the row in the fractional abundance matrix.

Moreover, all the entries in each homogeneous region are averaged by the coarse hyperspectral image, so the impact of noise can be mitigated effectively. Therefore, the superpixel-based weighting matrix ${\mathbf{H}}_1$ is not sensitive to noise.

To enhance the smoothness of abundance maps, the spatial weighting matrix ${\mathbf{H}}_2$ takes into account the correlation between neighboring regions, following a similar approach as described in [52]. Let $h_{2,ij}^{t+1}$ represent the entry in the $t+1$ iteration of the weight matrix ${\mathbf{H}}_2$ at the ith row and jth column. The spatial weighting factor ${\mathbf{H}}_2$ can be iteratively updated as follows:

$$h_{2,ij}^{(t+1)}=\frac{1}{f_{\mathcal{N}\left(x_{ij}\right)}\left(x_{ij}^{\left(t\right)}\right)+\epsilon }$$

(8)

where $\mathcal{N}\left(x_{ij}\right)$ is the adjacent set for elements $x_{ij}$, $f(·)$ is a function that calculates the fractional abundance of pixels using a neighborhood system, with the goal of utilizing spatial correlations in images.

In this work, a new value of a pixel can be calculated by using the values of its 8-neighboring pixels and their importance (i.e., Euclidean distance). The function $f(·)$ can be defined as follows:

$$f_{\mathcal{N}\left(x_{ij}\right)}\left(x_{ij}\right)=\frac{{\sum }_{x_a\in {\mathcal{N}}_8\left(x_{ij}\right)}{\Lambda }_a{x}_a}{{\sum }_{x_a\in {\mathcal{N}}_8\left(x_{ij}\right)}{\Lambda }_a}$$

(9)

where ${\mathcal{N}}_8\left(x_{ij}\right)$ represents adjacent elements ($x_a$) of the central pixel $x_{ij}$. $\Lambda =\frac{1}{im(x_a,x_{ij})}$ represents the neighboring importance. Where the function $im(·)$ represents the importance of the elements $x_a$ and $x_{ij}$. Considering the European distance, ${\Lambda }_a$ is written as follows: ${\Lambda }_a=1/{\left({(u-m)}^2+{(v-n)}^2\right)}^{1/2}$, where $(u,v)$ and $(m,n)$ represent the spatial coordinates of $x_a$ and $x_{ij}$, respectively.

3.2. Optimization Model with the ADMM

In order to address the optimization difficulties of model (5), we adopt the ADMM technique as a solution method. The central concept of the ADMM method involves introducing a set of additional variables and subsequently transforming it into a sequence of more manageable sub-problems derived from the original problem. The optimization problem expressed as Equation (5) can be reformulated in the following manner.

$$\begin{array}{cc}\hfill \underset{\mathbf{X},{\mathbf{V}}_1,{\mathbf{V}}_2,{\mathbf{V}}_3}{min}& \frac{1}{2}\parallel {\mathbf{V}}_1{-\mathbf{Y}\parallel }_F^2+\lambda {\parallel \left({\mathbf{H}}_1{\mathbf{H}}_2\right)\odot {\mathbf{V}}_2\parallel }_{1,1}+{\iota }_{R+}\left({\mathbf{V}}_3\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.\text{:}& {\mathbf{V}}_1=\mathbf{AX},{\mathbf{V}}_2=\mathbf{X},{\mathbf{V}}_3=\mathbf{X}\hfill \end{array}$$

(10)

the equation states that ${\iota }_{R+}\left(\mathbf{X}\right)={\sum }_{i=1}^n{\iota }_{R+}\left({\mathbf{x}}_i\right)$ is an indicator function, and ${\iota }_{R+}\left({\mathbf{x}}_i\right)=0$ if ${\mathbf{x}}_i$ belongs to the nonnegative orthant and ${\iota }_{R+}\left({\mathbf{x}}_i\right)=+\infty$ otherwise.

Define $\mathbf{V}\equiv {({\mathbf{V}}_{\mathbf{1}},{\mathbf{V}}_{\mathbf{2}},{\mathbf{V}}_{\mathbf{3}})}^{\mathbf{T}}$, $\mathbf{G}={\left[\mathbf{A},\mathbf{I},\mathbf{I}\right]}^{\mathbf{T}}$ and $\mathbf{B}=\mathrm{diag}(-\mathbf{I})$, where $\mathbf{I}$ is the identity matrix, and

$$\begin{array}{cc}\hfill g\left(\mathbf{V}\right)\equiv & \frac{1}{2}\parallel {\mathbf{V}}_1{-\mathbf{Y}\parallel }_F^2+\lambda {\parallel \left({\mathbf{H}}_1{\mathbf{H}}_2\right)\odot {\mathbf{V}}_2\parallel }_{1,1}+{\iota }_{R+}\left({\mathbf{V}}_3\right)\hfill \end{array}$$

(11)

We can express Equation (10) in a condensed form as follows.

$$\begin{array}{c}\hfill \underset{\mathbf{X},\mathbf{V}}{min}g\left(\mathbf{V}\right)\mathrm{subject}\mathrm{to}\mathbf{GX}+\mathbf{BV}=\mathbf{0},\end{array}$$

(12)

We can denote the augmented Lagrangian function as given below.

$$\mathcal{L}\left(\mathbf{X},\mathbf{V},\mathbf{D}\right)\equiv g\left(\mathbf{V}\right)+\frac{\mu }{2}{\parallel \mathbf{GX}+\mathbf{BV}-\mathbf{D}\parallel }_F^2$$

(13)

where $\mu >0$ represents a positive constant, and $\mathbf{D}={({\mathbf{D}}_{\mathbf{1}},{\mathbf{D}}_{\mathbf{2}},{\mathbf{D}}_{\mathbf{3}})}^{\mathbf{T}}$ is the Lagrange multiplier related to the constraint $\mathbf{GX}+\mathbf{BV}=\mathbf{0}$. Subsequently, we can minimize the augmented Lagrangian function $\mathcal{L}(\mathbf{X},\mathbf{V},\mathbf{D})$ to obtain the values for $\mathbf{X}$ and $\mathbf{V}$. Simultaneously, we update $\mathbf{D}$ according to the following procedure.

$$\left\{\begin{array}{cc}\hfill {\mathbf{X}}^{(k+1)}& =\underset{\mathbf{X}}{argmin}\mathcal{L}\left(\mathbf{X},{\mathbf{V}}^{\left(\mathbf{k}\right)},{\mathbf{D}}^{\left(\mathbf{k}\right)}\right)\hfill \\ \hfill {\mathbf{V}}^{(k+1)}& =\underset{\mathbf{V}}{argmin}\mathcal{L}\left({\mathbf{X}}^{(\mathbf{k}+\mathbf{1})},\mathbf{V},{\mathbf{D}}^{\left(\mathbf{k}\right)}\right)\hfill \\ \hfill {\mathbf{D}}^{(k+1)}& ={\mathbf{D}}^{\left(k\right)}-\mathbf{G}{\mathbf{X}}^{(k+1)}-\mathbf{B}{\mathbf{V}}^{(k+1)}\hfill \end{array}\right.$$

Afterwards, we can utilize the aforementioned framework to calculate the $\mathbf{X}$ and $\mathbf{V}$-problems. First, the equivalent solution of the $\mathbf{X}$ subproblem is as follows.

$$\begin{array}{cc}\hfill \underset{\mathbf{X}}{min}& \frac{\mu }{2}\parallel \mathbf{AX}-{\mathbf{V}}_1^{\left(k\right)}-{\mathbf{D}}_1^{\left(k\right)}{\parallel }_F^2+\frac{\mu }{2}\parallel \mathbf{X}-{\mathbf{V}}_2^{\left(k\right)}-{\mathbf{D}}_2^{\left(k\right)}{\parallel }_F^2+\frac{\mu }{2}{\parallel \mathbf{X}-{\mathbf{V}}_3^{\left(k\right)}-{\mathbf{D}}_3^{\left(k\right)}\parallel }_F^2\hfill \end{array}$$

(14)

The solution for $\mathbf{X}$ is denoted as

$$\begin{array}{c}\hfill {\mathbf{X}}^{(k+1)}={({\mathbf{A}}^T\mathbf{A}+2\mathbf{I})}^{-1}({\mathbf{A}}^T{\xi }_1+{\xi }_2+{\xi }_3)\end{array}$$

(15)

where ${\mathbf{A}}^T$ denotes the transpose of matrix $\mathbf{A}$, and ${\xi }_1={\mathbf{V}}_1^{\left(k\right)}+{\mathbf{D}}_1^{\left(k\right)}$, ${\xi }_2={\mathbf{V}}_2^{\left(k\right)}+{\mathbf{D}}_2^{\left(k\right)}$, ${\xi }_3={\mathbf{V}}_3^{\left(k\right)}+{\mathbf{D}}_3^{\left(k\right)}$.

The subproblem of $\mathbf{V}$ decouples into three subparts about ${\mathbf{V}}_1,{\mathbf{V}}_2$, and ${\mathbf{V}}_3$. ${\mathbf{V}}_1$ is calculated as follows:

$$\begin{array}{c}\hfill \underset{{\mathbf{V}}_{\mathbf{1}}}{min}\frac{1}{2}\parallel {\mathbf{V}}_1-\mathbf{Y}{\parallel }_F^2+\frac{\mu }{2}{\parallel \mathbf{A}{\mathbf{X}}^{(k+1)}-{\mathbf{V}}_1-{\mathbf{D}}_1^{\left(k\right)}\parallel }_F^2\end{array}$$

(16)

It clearly follows that

$$\begin{array}{c}\hfill {\mathbf{V}}_1^{(k+1)}=\frac{1}{1+\mu }\left(\mathbf{Y}+\mu ({\mathbf{AX}}^{(k+1)}-{\mathbf{D}}_1^{\left(k\right)})\right)\end{array}$$

(17)

${\mathbf{V}}_2$ is obtained by solving the problem using the following method.

$$\begin{array}{cc}\hfill \underset{{\mathbf{V}}_2}{min}& \lambda \parallel \left({\mathbf{H}}_1{\mathbf{H}}_2\right)\odot {\mathbf{V}}_2{\parallel }_{1,1}+\frac{\mu }{2}{\parallel {\mathbf{X}}^{(k+1)}-{\mathbf{V}}_2-{\mathbf{D}}_2^{\left(k\right)}\parallel }_F^2\hfill \end{array}$$

(18)

The solution is

$$\begin{array}{c}\hfill {\mathbf{V}}_2^{(k+1)}=\mathrm{soft}({\mathbf{X}}^{(k+1)}-{\mathbf{D}}_2^{\left(k\right)},\frac{\lambda }{\mu }{\mathbf{H}}_1{\mathbf{H}}_2)\end{array}$$

(19)

the function $\mathrm{soft}(·,\tau )$ represents a soft-threshold function, which is defined as $\mathrm{soft}(y,\tau )=\mathrm{sign}\left(y\right)\max \left\{\right|y|-\tau ,0\}$.

For the ${\mathbf{V}}_3$, it can be obtained by solving as

$$\begin{array}{c}\hfill \underset{{\mathbf{V}}_3}{min}{\iota }_{R+}\left({\mathbf{V}}_3\right)+\frac{\mu }{2}{\parallel {\mathbf{X}}^{(k+1)}-{\mathbf{V}}_3-{\mathbf{D}}_3^{\left(k\right)}\parallel }_F^2\end{array}$$

(20)

The solution is

$$\begin{array}{c}\hfill {\mathbf{V}}_3^{(k+1)}=\max ({\mathbf{X}}^{(k+1)}-{\mathbf{D}}_3^{\left(k\right)},0)\end{array}$$

(21)

Lastly, the Lagrange multipliers are updated in the following manner.

$$\left\{\begin{array}{cc}\hfill {\mathbf{D}}_1^{(k+1)}& ={\mathbf{D}}_1^{\left(k\right)}-\mathbf{A}{\mathbf{X}}^{(k+1)}+{\mathbf{V}}_1^{(k+1)}\hfill \\ \hfill {\mathbf{D}}_2^{(k+1)}& ={\mathbf{D}}_2^{\left(k\right)}-{\mathbf{X}}^{(k+1)}+{\mathbf{V}}_2^{(k+1)}\hfill \\ \hfill {\mathbf{D}}_3^{(k+1)}& ={\mathbf{D}}_3^{\left(k\right)}-{\mathbf{X}}^{(k+1)}+{\mathbf{V}}_3^{(k+1)}\hfill \end{array}\right.$$

To summarize, Algorithm 1 outlines the use of the ADMM method to solve the optimization problem in the RDSWSU model. The RDSWSU algorithm, suggested in the study, is composed of an inner loop and an outer loop. In the inner loop, the Lagrange multipliers in the ADMM method are updated. On the other hand, the weight coefficients are updated in the outer loop. The maximum number of iterations in the inner loop is set to 5, while the number of iterations in the outer loop is set to 120. The iterative optimization of the model is similar to the well-known SUnSAL method in complexity due to the simplicity of the model.

Proving formal convergence for Algorithm 1 is difficult, however empirically the residual of the original problem is seen to decrease with more iterations. As illustrated in Figure 2, the residual $\parallel {\mathbf{GX}}^{\left(\mathbf{t}\right)}+{\mathbf{BV}}^{\left(\mathbf{t}\right)}{\parallel }_F\le 1\times {10}^{-5}$ as a function of outer loop counts reaches a stable near-zero level after only 40 iterations. The dual loop structure exhibits good empirical convergence and speed overall. The full algorithm halts when hitting the maximum iterations or once the residual satisfies a tolerance threshold.

Algorithm 1: The pseudocode representation for RDSWSU

1: Input: Y, A, and parameters of the SLIC algorithm 2: Perform: 3: Use the SLIC algorithm to segment image Y. 4: Compute Equation (6) to reconstruct the $\tilde{\mathbf{Y}}$ 5: $\tilde{\mathbf{X}}$ obtained by solving the optimization problem (7) 6: ${\mathbf{H}}_{1,ij}\leftarrow \frac{1}{\left|\right|\tilde{\mathbf{X}}(i,:){\left|\right|}_2+\epsilon }$ 7: Initialization: 8: set $k,t=0$, choose $\lambda ,\mu ,\epsilon >0,{\mathbf{X}}^{\left(0\right)},{\mathbf{V}}_1^{\left(0\right)},{\mathbf{V}}_2^{\left(0\right)},{\mathbf{V}}_3^{\left(0\right)},{\mathbf{D}}_1^{\left(0\right)},{\mathbf{D}}_2^{\left(0\right)},{\mathbf{D}}_3^{\left(0\right)}$ 9: Repeat: 10: ${\mathbf{H}}_2^{(t+1)}\leftarrow \left[\begin{array}{ccc}{h}_{2,11}^{\left(t\right)}& \dots & h_{2,1n}^{\left(t\right)}\\ ⋮& h_{2,ij}^{\left(t\right)}& ⋮\\ h_{2,m1}^{\left(t\right)}& \dots & h_{2,mn}^{\left(t\right)}\end{array}\right]$, 11: where $h_{2,ij}^{(t+1)}=\frac{1}{f_{\mathcal{N}\left(x_{ij}\right)}(x_{ij}^{\left(t\right)}-{d_{2,ij}}^{\left(t\right)})+\epsilon }$, $x_{ij}$ and $d_{2,ij}$ are respectively the elements in the ith row and jth column of $\mathbf{X}$ and ${\mathbf{D}}_2$. 12: Repeat: 13: ${\mathbf{X}}^{(k+1)}\leftarrow {({\mathbf{A}}^{\mathbf{T}}\mathbf{A}+2\mathbf{I})}^{-1}\left({\mathbf{A}}^{\mathbf{T}}({\mathbf{V}}_1^{\left(k\right)}+{\mathbf{D}}_1^{\left(k\right)})+{\mathbf{V}}_2^{\left(k\right)}+{\mathbf{D}}_2^{\left(k\right)}+{\mathbf{V}}_3^{\left(k\right)}+{\mathbf{D}}_3^{\left(k\right)}\right)$ 14: ${\mathbf{V}}_1^{(k+1)}\leftarrow \frac{1}{1+\mu }\left(\mathbf{Y}+\mu (\mathbf{A}{\mathbf{X}}^{(k+1)}-{\mathbf{D}}_1^{\left(k\right)})\right)$ 15: ${\mathbf{V}}_2^{(k+1)}\leftarrow \mathrm{soft}({\mathbf{X}}^{(k+1)}-{\mathbf{D}}_2^{\left(k\right)},\frac{\lambda }{\mu }{\mathbf{H}}_1{\mathbf{H}}_2^{\left(t\right)})$ 16: ${\mathbf{V}}_3^{(k+1)}\leftarrow \max ({\mathbf{X}}^{(k+1)}-{\mathbf{D}}_3^{\left(k\right)},0)$ 17: Revise the Lagrange multipliers: 18: ${\mathbf{D}}_1^{(k+1)}\leftarrow {\mathbf{D}}_1^{\left(k\right)}-\mathbf{A}{\mathbf{X}}^{(k+1)}+{\mathbf{V}}_1^{(k+1)}$ 19: ${\mathbf{D}}_2^{(k+1)}\leftarrow {\mathbf{D}}_2^{\left(k\right)}-{\mathbf{X}}^{(k+1)}+{\mathbf{V}}_2^{(k+1)}$ 20: ${\mathbf{D}}_3^{(k+1)}\leftarrow {\mathbf{D}}_3^{\left(k\right)}-{\mathbf{X}}^{(k+1)}+{\mathbf{V}}_3^{(k+1)}$ 21: Refresh the iteration: $k\leftarrow k+1$ 22: ${\mathbf{X}}^{(t+1)}\leftarrow {\mathbf{X}}^{\left(k\right)}$ 23: ${\mathbf{D}}_2^{(t+1)}\leftarrow {\mathbf{D}}_2^{\left(k\right)}$ 24: Refresh the iteration: $t\leftarrow t+1$ 25: until certain stopping criteria are satisfied.

4. Experiments with Synthetic Data

In these experiments, we will showcase the performance of the proposed RDSWSU method by conducting two simulated hyperspectral datasets. In our experiments, we will conduct a comparison between the proposed RDSWSU method and five other advanced sparse unmixing methods: SUnSAL [17], SUnSAL-TV [29], DRSU [27], DRSU-TV [31] and S${}^2$WSU [52]. It is worth mentioning that the last three algorithms utilize weighted sparse regression techniques.

Furthermore, two commonly utilized quality metrics are employed to assess the efficacy of sparse unmixing. The signal reconstruction error ratio (SRE), measured in decibels (dB), is a commonly used metric in quantitative analysis. It serves as an indicator to evaluate the accuracy of unmixing and gauge the performance of various unmixing algorithms. The calculation formula of SRE (dB) is:

$$\mathrm{SRE}\left(\mathrm{dB}\right)=10·{log}_{10}{\displaystyle \frac{E\left(\right|\left|\mathbf{x}\right|{|}_2^2)}{E\left(\right||\mathbf{x}-\widehat{\mathbf{x}}|{|}_2^2)}}$$

(22)

where $\widehat{\mathbf{x}}$ denotes the estimated fractional abundance, $\mathbf{x}$ denotes the true fractional abundance, and $E(·)$ indicates the expected function. The greater the value of SRE (dB), the higher the quality of unmixing achieved by this method.

In addition, we also adopt another measure, that is, the probability of success $p_s$. It represents an estimation of the likelihood that the relative error in power is below a certain threshold. The $p_s$ can be obtained by formally calculating as follows: $p_s\equiv P(\parallel \widehat{\mathbf{x}}{-\mathbf{x}\parallel }^2/{\parallel \mathbf{x}\parallel }^2\le threshold)$. In previous studies [17], it has been shown that the abundance estimation is deemed successful when $\parallel \widehat{\mathbf{x}}{-\mathbf{x}\parallel }^2/{\parallel \mathbf{x}\parallel }^2\le 3.16$ (5 dB). The unmixing performance improves as the $p_s$ value increases.

4.1. Simulated Data Sets

In this paper, we utilized two spectral libraries in the synthetic data analysis, both of which were extracted from the mineral dictionary of the US Geological Survey (USGS) spectral library (Available online at http://speclab.cr.usgs.gov/spectral.lib06 (accessed on 21 June 2023)). The spectral library ${\mathbf{A}}_1\in {\mathbb{R}}^{224\times 240}$ comprises $m=240$ spectral signatures with reflectance values across $L=224$ spectral bands, evenly distributed within the wavelength range of 0.4–2.5 $\mathsf{\mu }$m. Similar to ${\mathbf{A}}_1$, the spectral library ${\mathbf{A}}_{\mathbf{2}}\in {\mathbb{R}}^{221\times 222}$ consists of $m=222$ spectral signatures with $L=221$ spectral bands. According to LMM, two distinct simulated datasets were generated from spectral libraries ${\mathbf{A}}_{\mathbf{1}}$ and ${\mathbf{A}}_{\mathbf{2}}$. The fractional abundances were enforced with ANC and ASC constraints. Subsequently, a detailed description of the two simulated datasets will be provided.

• Simulated Data Cube 1 (DC1): DC1 comprises $100\times 100$ pixels, with each pixel containing 224 spectral bands. DC1 was generated by randomly selecting nine spectral signatures from spectral library ${\mathbf{A}}_{\mathbf{1}}$. Figure 3 illustrates the actual abundance maps associated with the nine endmembers. The simulated hyperspectral data generation process involved the introduction of independent and identically distributed (i.i.d.) Gaussian noise at five different SNR levels: 10 dB, 20 dB, 30 dB, 40 dB, and 50 dB.
• Simulated Data Cube 2 (DC2): DC2 has $100\times 100$ pixels, with each pixel containing 221 spectral bands. Simulated data DC2 was generated by randomly selecting nine spectral signatures from spectral library ${\mathbf{A}}_{\mathbf{2}}$, as detailed in [53,54]. Figure 4 displays the genuine abundance maps corresponding to the nine endmembers. In the experiments, similar to DC1, the simulated data for DC2 was also subjected to independent and i.i.d. Gaussian noise at five different SNR levels: 10 dB, 20 dB, 30 dB, 40 dB, and 50 dB.

4.2. Results and Discussion

In this section, we use two simulated datasets to test the proposed unmixing performance based on a dual spatial weighting strategy. For the parameter settings in the experiments, each method was tested with different combinations of parameter settings to ensure the optimal results.

Table 1. The SRE (dB) values achieved by various unmixing algorithms on DC1 and DC2 are as follows, with the optimal parameters indicated in parentheses.

Date Cube	SNR	SUnSAL	SUnSAL-TV	DRSU	DRSU-TV	S${}^2$WSU	RDSWSU
	10	1.6582	3.9092	2.1998	4.3212	3.4088	5.6387
		($\lambda =2$)	($\lambda =4·{10}^{-1}$;	($\lambda =6·{10}^{-1}$)	($\lambda =1·{10}^{-1}$;	($\lambda =2·{10}^{-1}$)	($\lambda =1·{10}^{-1}$)
			${\lambda }_{TV}=6·{10}^{-1}$)		${\lambda }_{TV}=8·{10}^{-3}$)
	20	4.1966	5.3175	4.5620	8.4389	7.8752	14.2101
		($\lambda =2·{10}^{-1}$)	($\lambda =1·{10}^{-1}$;	($\lambda =1·{10}^{-1}$)	($\lambda =2·{10}^{-2}$;	($\lambda =9·{10}^{-3}$)	($\lambda =3·{10}^{-2}$)
			${\lambda }_{TV}=2·{10}^{-2}$)		${\lambda }_{TV}=3·{10}^{-2}$)
DC1	30	8.5108	11.9497	15.2985	18.9425	21.8718	22.6783
		($\lambda =2·{10}^{-2}$)	($\lambda =1·{10}^{-2}$;	($\lambda =6·{10}^{-4}$)	($\lambda =2·{10}^{-3}$;	($\lambda =4·{10}^{-3}$)	($\lambda =6·{10}^{-3}$)
			${\lambda }_{TV}=4·{10}^{-3}$)		${\lambda }_{TV}=2·{10}^{-3}$)
	40	15.2388	17.8119	29.5464	31.0163	32.1470	32.3118
		($\lambda =5·{10}^{-3}$)	($\lambda =5·{10}^{-3}$;	($\lambda =1·{10}^{-3}$)	($\lambda =3·{10}^{-3}$;	($\lambda =2·{10}^{-3}$)	($\lambda =1·{10}^{-3}$)
			${\lambda }_{TV}=1·{10}^{-3}$)		${\lambda }_{TV}=1·{10}^{-3}$)
	50	23.1066	26.0519	41.2077	40.8223	41.5089	41.9386
		($\lambda =1·{10}^{-3}$)	($\lambda =2·{10}^{-3}$;	($\lambda =3·{10}^{-3}$)	($\lambda =6·{10}^{-4}$;	($\lambda =4·{10}^{-4}$)	($\lambda =2·{10}^{-4}$)
			${\lambda }_{TV}=2·{10}^{-4}$)		${\lambda }_{TV}=1·{10}^{-4}$)
	10	0.6993	2.4822	0.9698	2.6009	1.9869	6.1135
		($\lambda =1$)	($\lambda =1·{10}^{-1}$;	($\lambda =6·{10}^{-1}$)	($\lambda =6·{10}^{-2}$;	($\lambda =2·{10}^{-1}$)	($\lambda =1·{10}^{-1}$)
			${\lambda }_{TV}=4·{10}^{-1}$)		${\lambda }_{TV}=1·{10}^{-1}$)
	20	2.6569	4.9553	2.7730	7.6542	4.5408	11.9055
		($\lambda =2·{10}^{-1}$)	($\lambda =2·{10}^{-2}$;	($\lambda =7·{10}^{-2}$)	($\lambda =1·{10}^{-2}$;	($\lambda =9·{10}^{-3}$)	($\lambda =8·{10}^{-3}$)
			${\lambda }_{TV}=2·{10}^{-2}$)		${\lambda }_{TV}=4·{10}^{-2}$)
DC2	30	5.9799	9.5503	15.6040	18.8693	19.8190	20.1171
		($\lambda =8·{10}^{-3}$)	($\lambda =4·{10}^{-3}$;	($\lambda =2·{10}^{-3}$)	($\lambda =2·{10}^{-3}$;	($\lambda =5·{10}^{-3}$)	($\lambda =2·{10}^{-3}$)
			${\lambda }_{TV}=2·{10}^{-3}$)		${\lambda }_{TV}=2·{10}^{-3}$)
	40	11.8340	15.2535	26.9881	28.2618	28.6033	28.6337
		($\lambda =2·{10}^{-3}$)	($\lambda =6·{10}^{-5}$;	($\lambda =6·{10}^{-4}$)	($\lambda =9·{10}^{-4}$;	($\lambda =7·{10}^{-4}$)	($\lambda =5·{10}^{-4}$)
			${\lambda }_{TV}=9·{10}^{-4}$)		${\lambda }_{TV}=4·{10}^{-4}$)
	50	18.7115	24.6150	35.3070	36.5517	36.7807	37.0450
		($\lambda =3·{10}^{-4}$)	($\lambda =5·{10}^{-5}$;	($\lambda =1·{10}^{-4}$)	($\lambda =2·{10}^{-4}$;	($\lambda =5·{10}^{-5}$)	($\lambda =2·{10}^{-5}$)
			${\lambda }_{TV}=9·{10}^{-5}$)		${\lambda }_{TV}=9·{10}^{-5}$)

and

Table 2. The $p_s$ values achieved by various unmixing algorithms on DC1 and DC2 are as follows, with the optimal parameters indicated in parentheses.

Date Cube	SNR	SUnSAL	SUnSAL-TV	DRSU	DRSU-TV	S${}^2$WSU	RDSWSU
	10	0.2933	0.5206	0.3777	0.5899	0.5058	0.6810
		($\lambda =2$)	($\lambda =4·{10}^{-1}$;	($\lambda =6·{10}^{-1}$)	($\lambda =1·{10}^{-1}$;	($\lambda =2·{10}^{-1}$)	($\lambda =1·{10}^{-1}$)
			${\lambda }_{TV}=6·{10}^{-1}$)		${\lambda }_{TV}=8·{10}^{-3}$)
	20	0.5594	0.6406	0.6185	0.8387	0.8005	0.9716
		($\lambda =2·{10}^{-1}$)	($\lambda =1·{10}^{-1}$;	($\lambda =1·{10}^{-1}$)	($\lambda =2·{10}^{-2}$;	($\lambda =9·{10}^{-3}$)	($\lambda =3·{10}^{-2}$)
			${\lambda }_{TV}=2·{10}^{-2}$)		${\lambda }_{TV}=3·{10}^{-2}$)
DC1	30	0.7992	0.9492	0.9782	0.9994	1	1
		($\lambda =2·{10}^{-2}$)	($\lambda =1·{10}^{-2}$;	($\lambda =6·{10}^{-4}$)	($\lambda =2·{10}^{-3}$;	($\lambda =4·{10}^{-3}$)	($\lambda =6·{10}^{-3}$)
			${\lambda }_{TV}=4·{10}^{-3}$)		${\lambda }_{TV}=2·{10}^{-3}$)
	40	0.9902	0.9996	1.0000	1.0000	1.0000	1
		($\lambda =5·{10}^{-3}$)	($\lambda =5·{10}^{-3}$;	($\lambda =1·{10}^{-3}$)	($\lambda =3·{10}^{-3}$;	($\lambda =2·{10}^{-3}$)	($\lambda =1·{10}^{-3}$)
			${\lambda }_{TV}=1·{10}^{-3}$)		${\lambda }_{TV}=1·{10}^{-3}$)
	50	1	1	1	1	1	1
		($\lambda =1·{10}^{-3}$)	($\lambda =2·{10}^{-3}$;	($\lambda =3·{10}^{-3}$)	($\lambda =6·{10}^{-4}$;	($\lambda =4·{10}^{-4}$)	($\lambda =2·{10}^{-4}$)
			${\lambda }_{TV}=2·{10}^{-4}$)		${\lambda }_{TV}=1·{10}^{-4}$)
	10	0.1812	0.3029	0.1903	0.4015	0.3451	0.6368
		($\lambda =1$)	($\lambda =1·{10}^{-1}$;	($\lambda =6·{10}^{-1}$)	($\lambda =6·{10}^{-2}$;	($\lambda =2·{10}^{-1}$)	($\lambda =1·{10}^{-1}$)
			${\lambda }_{TV}=4·{10}^{-1}$)		${\lambda }_{TV}=1·{10}^{-1}$)
	20	0.3766	0.5031	0.4359	0.7186	0.6013	0.9076
		($\lambda =2·{10}^{-1}$)	($\lambda =2·{10}^{-2}$;	($\lambda =7·{10}^{-2}$)	($\lambda =1·{10}^{-2}$;	($\lambda =9·{10}^{-3}$)	($\lambda =8·{10}^{-3}$)
			${\lambda }_{TV}=2·{10}^{-2}$)		${\lambda }_{TV}=4·{10}^{-2}$)
DC2	30	0.5900	0.8023	0.9795	0.9976	1	1
		($\lambda =8·{10}^{-3}$)	($\lambda =4·{10}^{-3}$;	($\lambda =2·{10}^{-3}$)	($\lambda =2·{10}^{-3}$;	($\lambda =5·{10}^{-3}$)	($\lambda =2·{10}^{-3}$)
			${\lambda }_{TV}=2·{10}^{-3}$)		${\lambda }_{TV}=2·{10}^{-3}$)
	40	0.9193	0.9884	1.0000	1.0000	1.0000	1
		($\lambda =2·{10}^{-3}$)	($\lambda =6·{10}^{-5}$;	($\lambda =6·{10}^{-4}$)	($\lambda =9·{10}^{-4}$;	($\lambda =7·{10}^{-4}$)	($\lambda =5·{10}^{-4}$)
			${\lambda }_{TV}=9·{10}^{-4}$)		${\lambda }_{TV}=4·{10}^{-4}$)
	50	0.9999	1	1	1	1	1
		($\lambda =3·{10}^{-4}$)	($\lambda =5·{10}^{-5}$;	($\lambda =1·{10}^{-4}$)	($\lambda =2·{10}^{-4}$;	($\lambda =5·{10}^{-5}$)	($\lambda =2·{10}^{-5}$)
			${\lambda }_{TV}=9·{10}^{-5}$)		${\lambda }_{TV}=9·{10}^{-5}$)

exhibit the SRE (dB) and $p_s$ values obtained by each tested unmixing method on the two simulated datasets at different SNR levels. The tables show the optimal values obtained within the specified parameter range.

Table 1 and Table 2 demonstrate a trend where the performance of the listed algorithms for unmixing worsens as the noise level increases. Compared with other methods, the RDSWSU method achieves the highest SRE (dB) values across all noise levels. It is noteworthy that the proposed algorithm exhibits a significant advantage in low signal-to-noise ratio situations, and can achieve more accurate separation performance. For example, in DC1 and DC2, when the SNR is 10 dB (i.e., in high noise conditions), our method has increased the SRE (dB) by more than 2 dB compared to other methods, indicating that the RDSWSU method has stronger noise resistance. Furthermore, compared with several other algorithms, the algorithm presented in this paper also achieves superior results in terms of $p_s$ value, especially in low SNR conditions, which validates the high robustness of the dual-spatial weighting strategy. In summary, based on the dual-spatial weighting strategy, it has provided a good prospect for improving the unmixing accuracy in two different scenes.

To facilitate a direct comparison between different unmixing algorithms, we select DC1 and DC2 as examples at an SNR of 20 dB. This selection aims to demonstrate the improved sparsity of the solution achieved by the proposed RDSWSU algorithm. Figure 5 and Figure 6 demonstrate the corresponding fractional abundances estimated by each unmixing algorithm and the ground-truth abundances. It should be noted that each line represents a nonzero vector of the abundance matrix (with 100 randomly selected pixels shown here for display). It is evident that both the SUnSAL and SUnSAL-TV algorithms produce inaccurate estimations of the number of recovered endmembers. In contrast, the RDSWSU algorithm exhibits a remarkable ability to accurately estimate the number of endmembers, surpassing the performance of the DRSU, DRSU-TV, and S${}^2$WSU algorithms in this regard. Therefore, compared with SUnSAL-TV, DRSU, DRSU-TV, and S${}^2$WSU, the proposed RDSWSU algorithm can accurately identify the optimal subset of endmembers from the spectral library. Compared with other unmixing methods, the RDSWSU algorithm improved the ability to identify endmembers from a spectral library while also enhancing the sparsity of the abundance matrix.

Figure 7 and Figure 8 present the abundance maps of the first endmember achieved using various unmixing algorithms in DC1 and DC2 when the SNR is set to 20 dB. For a more intuitive comparison, the difference maps, which are generated by subtracting the estimated abundances from the ground truth abundances, are also shown in the figures. From Figure 7 and Figure 8, it can be observed that the abundance maps derived from SUnSAL and SUnSAL-TV lack accuracy. The abundance maps generated by the SUnSAL algorithm exhibit noticeable noise, indicating a relatively high level of noise present in the results. The SUnSAL-TV algorithm, on the one hand, enhances the quality of abundance estimation to some extent by incorporating spatial information. On the other hand, the abundance maps produced by the SUnSAL-TV algorithm appear over-smoothness phenomenon and blurred in their visual representation. The unmixing outcomes achieved by RDSWSU exhibit a remarkable resemblance to the true abundance maps, distinguishing it from other algorithms. Additionally, its abundance maps demonstrate a greater amount of spatial detail information across various homogeneous regions. In Figure 7 and Figure 8, it is evident that the abundance map obtained by the RDSWSU algorithm is smoother at the boundary compared to the results of the SUnSAL, DRSU, and S${}^2$WSU algorithms. However, it does not exhibit over-smoothing like the SUnSAL-TV and DRSU-TV algorithms. The RDSWSU algorithm integrates the dual spatial weighting factor with the estimated abundance through element-wise multiplication, and minimizes the ${\mathsf{\ell }}_1$ norm of their product. Thereinto, the superpixel spatial weighting factor is established at the scale of superpixels, where each superpixel corresponds to a homogeneous region in the image. During the unmixing process, the pixels within the same superpixel block are assigned identical weights, while the weights for pixels in different superpixel blocks differ. In a sense, the superpixel spatial weighting factor helps to preserve the spatial heterogeneity of the image. Furthermore, even within a superpixel, the final weight of each pixel can vary because there is another spatial neighborhood weighting factor that affects the final weight. The spatial neighborhood weighting factor is determined by the abundance vectors of its neighboring pixels within an 8-neighborhood. Each pixel is assigned a unique weight and dynamically changes with iterations. Unlike the TV regularization that directly constrains the differences in the abundance vectors of each pixel in the horizontal and vertical directions, the spatial neighborhood weighting factor induces the pixel with lower neighboring abundance to have an even lower abundance, rather than forcing its abundance to be equal to that of the neighboring pixels. As a result, the dual spatial weighting factor effectively utilizes the spatial correlation between pixels while also maintaining the heterogeneity of different homogeneous regions and promoting sparsity of abundance, thus avoiding excessive smoothing. It is proved that the weight factors based on superpixels and spatial neighborhood can significantly improve the effect of sparse unmixing.

Table 3. Computation times of different algorithms on processing DC1 with SNR = 20 dB.

Algorithm	SUnSAL	SUnSAL-TV	DRSU	DRSU-TV	S${}^2$WSU	RDSWSU
Time(s)	23.29	231.13	72.45	250.00	79.53	80.84

shows the computational time taken by different algorithms to process DC1 at SNR = 20 dB. All algorithms were coded using MATLAB (R2018b) and tested in a standardized computing environment. The testing was conducted on a desktop computer equipped with an Intel Core i7 dual-core processor running at 3.6 GHz and 16 GB of RAM.

It is noteworthy that the RDSWSU method did not include the time for superpixel segmentation and spatial weighting computation in its computational time, as both of these steps were preprocessed and the generated weight values were introduced as constants. Consequently, there is no necessity to recompute or make adjustments to the superpixel segmentation and spatial weighting while performing the unmixing procedure. By referring to Table 3, it can be noted that the SUnSAL algorithm exhibits the shortest processing time. On the other hand, the remaining algorithms are comparatively slower as they are derived from the SUnSAL algorithm. The SUnSAL-TV and DRSU-TV methods, which employ TV regularization, exhibit longer processing times due to the high computational cost associated with solving the TV spatial regularization term. This results in a longer duration compared to other algorithms. The processing times of the three algorithms without TV regularization are nearly identical, indicating that the proposed RDSWSU algorithm achieves superior unmixing performance in comparison to other algorithms without introducing additional computational complexity.

Furthermore, we examined the effects of varying the two ADMM parameters—regularization $\lambda$ and Lagrangian multiplier $\mu$- on unmixing performance. Using the DC1 dataset at 30 dB SNR as an example, Figure 9a shows our proposed RDSWSU algorithm can readily achieve suboptimal tuning by first fixing $\mu$ randomly, then selecting $\lambda$ from a set. As Figure 9b demonstrates, $\mu$ selection is robust over a wide range, producing minimal variability in the results.

5. Experiments with Real Hyperspectral Data

We employ two real hyperspectral data sets in this section—the Cuprite data and Houston data—to assess the performance of the various algorithms.

5.1. Experiment on Cuprite Data

The performance of the proposed method was evaluated using the well-known Cuprite dataset from the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS). The Cuprite site holds significant recognition in the field of mineralogy. The scene has been extensively utilized as a benchmark for validating the effectiveness of unmixing algorithms due to several exposed minerals. The experimental part focuses on a subset of the scene that measures $350\times 350$ pixels. This subset includes 224 spectral bands with wavelengths spanning from 0.4 to 2.5 $\mathsf{\mu }$m, and each band has a spectral resolution of 10 nm. Before conducting the analysis, certain bands were excluded due to their low water absorption and inadequate SNR. Specifically, bands 1–2, 105–115, 150–170, and 223–224 were removed from consideration. As a result, there were 188 spectral bands remaining for the analysis after the exclusion of the aforementioned bands. More specifically, the spectral library utilized in this experiment is equivalent to the spectral library ${\mathbf{A}}_{\mathbf{1}}$ that was employed in our simulation experiments. In addition, bands with a low SNR and those exhibiting low water absorption were also excluded from ${\mathbf{A}}_{\mathbf{1}}$. Furthermore, Tricorder software was utilized to generate classification maps depicting the distribution of these minerals. Unfortunately, the classification map was producted in 1995, and the publicly accessible Cuprite data was collected in 1997, a direct comparison between the classification map and the AVIRIS Cuprite data is not feasible. Therefore, the classification accuracy is not provided. Nonetheless, the classification map proves to be a valuable tool for qualitatively assessing the fractional abundance maps generated by the unmixing algorithms. Figure 10 illustrates a mineral map drawn by the USGS utilizing Tricorder 3.3 software [55] to represent the spatial distribution of various minerals within the Cuprite mining area.

Figure 11 presents a visual comparison of the classification maps produced by the USGS Tricorder algorithm in real hyperspectral scenes, along with three distinct mineral abundance maps (Alunite, Buddingtonite, and Chalcedony) generated by the SUnSAL, SUnSAL-TV, DRSU, DRSU-TV, S${}^2$WSU, and RDSWSU algorithms. In this experiment, the regularization parameters were set as $\lambda =0.001$, $\lambda =0.0001$, $\lambda =0.002$, and $\lambda =0.003$ for SUnSAL, DRSU, S${}^2$WSU, and RDSWSU, respectively. For SUnSAL-TV and DRSU-TV, the regularization parameters were set to $\lambda =0.001,{\lambda }_{TV}=0.001$ and $\lambda =0.002,{\lambda }_{TV}=0.0001$, respectively. The reasonable unmixing results achieved by all methods in Figure 11 demonstrate the effectiveness of sparse unmixing algorithms in interpreting the considered real hyperspectral datasets. As shown in Figure 11, the abundance maps (such as Buddingtonite) produced by SUnSAL and SUnSAL-TV algorithms displayed many noisy points. The results obtained using the SUnSAL-TV algorithm display an over-smoothness phenomenon, particularly noticeable in the estimation of Alunite and Buddingtonite. Additionally, the abundance map generated by DRSU algorithm has a mediocre performance in terms of spatial consistency (e.g., Chalcedony). We can observe that compared to the DRSU-TV and S${}^2$WSU algorithms, the abundance map generated by RDSWSU algorithm (such as Buddingtonite) is more similar to the reference image and shows smaller noise, demonstrating good spatial detail information. Considering qualitative criteria, it can be inferred that the proposed RDSWSU method outperforms other methods in terms of interpretability for the given real hyperspectral scene. To quantitatively assess the unmixing results in the real hyperspectral data experiment, we introduce a new evaluation metric, known as sparsity. This metric quantifies the ratio of elements in the estimated abundance that exceed 0.005 to the total number of elements. A lower sparsity value indicates a more sparse unmixing outcome, aligning with the expectations of sparse regression-based unmixing algorithms. The sparsity achieved by SUnSAL, SUnSAL-TV, DRSU, DRSU-TV, S${}^2$WSU and RDSWSU is 0.0682, 0.0743, 0.0428, 0.0422, 0.0421, and 0.0408 respectively. From these minor differences, we can infer that the proposed method represents the data with fewer components, improving sparsity. This qualitative analysis indicates the potential of the newly proposed RDSWSU algorithm to boost unmixing capabilities in real-world applications.

5.2. Experiment on Houston Data

In addition, we have introduced another commonly used benchmark dataset, the Houston dataset [56], for real hyperspectral data experiments to further demonstrate the superiority of the proposed RDSWSU algorithm. The Houston dataset was acquired by the ITRES CASI-1500 sensor over the University of Houston campus in June 2012 (Available online at http://hyperspectral.ee.uh.edu (accessed on 21 June 2023)). It comprises $170\times 170$ pixels, with 144 spectral bands spanning 364 nm to 1046 nm. Four endmembers are present in the data: Running track, Parking lot1, Parking lot2 and Grass healthy. Figure 12 shows the Houston data along with the real abundance maps and spectral signatures of the four endmembers. Notably, the endmember spectral library ${\mathbf{A}}_3\in {\mathbb{R}}^{144\times 4}$ used here only contains the actual spectral signatures of the four endmembers.

Figure 13 presents the abundance maps estimated by the six unmixing methods for the four endmembers. The obtained SRE (dB) values are 20.5499, 21.5405, 20.5757, 21.5473, 21.5537, and 21.5805 for SUnSAL, SUnSAL-TV, DRSU, DRSU-TV, S${}^2$WSU and RDSWSU, respectively. The quantitative results demonstrate our proposed method achieves comparable or slightly better performance versus other approaches. Notably, the differences are marginal, largely owing to the endmember spectral library containing only the four actual spectral signatures, which limits the full potential of weighted sparse unmixing algorithms. Still, the stable results from our method signify its merits.

6. Conclusions

This paper introduces a novel algorithm called robust dual spatial weighted sparse unmixing (RDSWSU) to enhance the outcomes of sparse unmixing. The RDSWSU algorithm presented in this study leverages superpixel-based spatial weighting and spatial neighborhood weighting factors for effectively utilizing the spatial information embedded in hyperspectral images. This can achieve accurate and stable unmixing results in high-noise environments. Additionally, this method utilizes an inner and outer loop strategy aimed at accelerating the convergence speed of the iterative optimization process of the model. The experimental results obtained from simulated and real hyperspectral datasets have showcased the remarkable capabilities of the proposed RDSWSU method. It exhibits substantial potential in enhancing unmixing performance, particularly in high-noise scenarios. While the experimental results obtained in this study are highly promising, it is important to acknowledge the need for conducting experiments on additional datasets in order to comprehensively assess the effectiveness of the proposed algorithm. In future endeavors, novel techniques such as low-rank prior [57] and deep learning [58] will be employed to effectively leverage the intrinsic feature information in images. This approach aims to enhance data interpretability and extract more valuable insights from the data.