Double Reweighted Sparse Regression and Graph Regularization for Hyperspectral Unmixing

: Hyperspectral unmixing, aiming to estimate the fractional abundances of pure spectral signatures in each mixed pixel, has attracted considerable attention in analyzing hyperspectral images. Plenty of sparse unmixing methods have been proposed in the literature that achieved promising performance. However, many of these methods overlook the latent geometrical structure of the hyperspectral data which limit their performance to some extent. To address this issue, a double reweighted sparse and graph regularized unmixing method is proposed in this paper. Speciﬁcally, a graph regularizer is employed to capture the correlation information between abundance vectors, which makes use of the property that similar pixels in a spectral neighborhood have higher probability to share similar abundances. In this way, the latent geometrical structure of the hyperspectral data can be transferred to the abundance space. In addition, a double weighted sparse regularizer is used to enhance the sparsity of endmembers and the fractional abundance maps, where one weight is introduced to promote the sparsity of endmembers as a hyperspectral image typically contains fewer endmembers compared to the overcomplete spectral library and the other weight is exploited to improve the sparsity of the abundance matrix. The weights of the double weighted sparse regularizer used for the next iteration are adaptively computed from the current abundance matrix. The experimental results on synthetic and real hyperspectral data demonstrate the superiority of our method compared with some state-of-the-art approaches.


Introduction
Hyperspectral imaging has the capacity to record the same scene over many contiguous spectral bands including the visible, nearinfrared and shortwave infrared spectral bands, and it has been widely used in practical applications such as terrain classification, mineral detection and exploration [1,2], military target discrimination, environmental monitoring and pharmaceutical counterfeiting [3]. Mixed pixels are commonly present in hyperspectral data due to the low spatial resolution of the sensors and the complexity of the terrain. The occurrence of mixed pixels severely degrades the application of hyperspectral data [4]. Thus, hyperspectral unmixing (HU) is an important process for hyperspectral data exploitation.
The task of HU is to identify a set of endmembers, and to estimate their corresponding abundances in the formation of each hyperspectral image (HSI) pixel [5][6][7]. HU algorithms mainly depend on the expected type of mixing, which can be classified as either a nonlinear or linear model [8]. Nonlinear mixing model assumes that part of the source radiation is affected by multiple scattering effects through several endmembers in the scene before being collected at the sensor. Therefore, the observed mixed pixel spectrum is generated from nonlinear function of the abundance associated with the endmembers. On the other hand, the linear mixture model (LMM) assumes that the mixing only occurs within the measuring instrument itself and interactions among distinct endmembers are negligible. Hence, the observed mixed pixel spectrum can be expressed as a linear combination of a collection of endmember signatures, which is weighted by their corresponding fractional abundances. Even though LMM is not always a precise model to characterize many real scenarios, it is generally recognized as an acceptable approximation of the light scattering and instrument mechanism. It has been widely used for HU due to its simplicity, effectiveness in different applications and computational tractability.
Based on the LMM, many methods have been proposed, such as geometry and statistics based approaches [9]. Geometry based methods usually require that pixel observations of hyperspectral data are within a simplex and their vertices correspond to the endmember [10]. This class of methods includes N-FINDR [11], vertex component analysis [12], minimum-volume simplex analysis [13], simplex growing algorithm [14], etc. In general, these methods require the estimation of the number of the endmembers or the presence of pure pixels. These assumptions may not hold in many hyperspectral data owing to the relatively low spatial resolution of imaging spectrometers and the presence of the mixture phenomena at different scales.
Statistical approaches utilize a statistical characteristic of the hyperspectral data, such as independent component analysis (ICA) [15], and nonnegative matrix factorization (NMF) [16]. NMF can learn a parts-based representation of the data and decompose the high-dimensional data into two nonnegative matrices. Consequently, one can get the endmember signatures and their corresponding abundances simultaneously. This method does not require the pure pixel assumption. Nevertheless, NMF could lead to many local minima owing to the nonconvexity of the objective function. To alleviate this disadvantage, many methods have been proposed by adding some extra constraints under the NMF framework. For example, spatial and spectral regularization is considered in [4,17], the sparsity-constrained of abundances is considered in [18], and the manifold structure of the hyperspectral data is considered in [19,20]. All these methods are unsupervised methods, i.e., both the endmember signatures and the abundances are unknown.
As more spectral libraries become available to the public, many semisupervised methods have been proposed recently. In [21], Bioucas-Dias et al. proposed a sparse unmixing method, named SUnSAL, which is based on variable splitting and augmented Lagrangian. This method assumes that the observed image signatures can be expressed by a linear combination of a small number of pure spectral signature. The high mutual coherence of the hyperspectral library limits the performance of SunSAL. To mitigate the impact of the mutual coherence, Iordache et al. proposed a collaborative sparse regression method for hyperspectral unmixing (CLSUnSAL) [22]. To further improve the performance of the CLSUnSAL, a HU method based on weighted sparse regression with weighted l d,1 (d = 1 or 2) mixed norm was proposed by Zheng et al. [23]. In [24], Tang et al. assumed that some materials in the spectral library are known to exist in the hyperspectral scene in many situations and they propose the SunSPI algorithm which employs the spectral a priori information. This method improves the performance of SunSAL and CLSUnSAL. However, these methods do not exploit the spatial-contextual information of the hyperspectral images. In [5], Iordache et al. proposed the SUnSAL-TV method, which utilizes the spatial information of the first-order pixel neighborhood system of the abundance map through total variation (TV) on sparse unmixing formulation. In [25], Zhao et al. studied total variation regularization in deblurring and sparse unmixing of hyperspectral images and their method gets better performance than SUnSAL-TV. Taking the non-local spatial information of whole abundance image into account, the non-local mean unmixing method which utilizes the similar patterns and structures in the abundance map is proposed in [26]. Allowing for the low rankness and sparsity property of abundance matrices, Giampouras et al. proposed a method which combines the weighted l 1 -norm and the weighted nuclear norm of the abundance matrix as a penalty term [27]. These methods ignore the underlying structure of the hyperspectral images. To take advantage of this property, some graph-based methods are proposed [28,29], which employ the graph topology and sparse group lasso regularization.
To fully utilize the underlying structure of hyperspectral images (HSIs) and sparsity property of abundance matrix, we propose a double reweighted sparse and graph regularized hyperspectral unmixing model. The motivation is based on two priors. First, the similar pixels in a spectral neighborhood are more likely to share the same abundances according to the fundamental spatial-spectral correlation property of HSI. In our method, this latent structure is modeled by a graph regularizer, which groups the highly related neighboring pixels on the graph. In this way, it transfers the structure of hyperspectral data to the unmixing results. Second, a hyperspectral image typically contains fewer endmembers comparing to the spectral library and the abundance maps are naturally sparse. Thus, we employ a double reweighted sparse regularizer to capture this prior. Particularly, we introduce one weight to promote the column sparsity of fractional abundance, while another weight is employed to enhance the sparsity of the abundance along the abundance vector corresponding to each endmember. Moreover, we design a symmetric alternating direction method of multipliers based algorithm to solve the proposed model.
There are some works related to us (e.g., [4,29]). In [4], Wang et al. introduced a double reweighted sparse unmixing and TV (DRSU-TV) algorithm for hyperspectral unmixing. The difference is that they employ the total variation regularization which utilizes local neighbourhood pixels on the abundances rather than graph regularization which can exploit the structures hidden in the hyperspectral data and allow connecting a pixel with as many other pixels in the image as long as they are similar. In [29], the authors incorporated a graph-based total variation framework within the unmixing problem, which takes the strength of graph TV regularization and l 2,1 -norm regularization. Unlike [29], in which the graph TV regularization is imposed on the reconstructed spectra, we exploit the graph Laplacian directly on abundances. In addition, we adopt a double reweighted l 1 -norm rather than l 2,1 -norm regularization employed in work [29].
The rest of this paper is organized as follows. In Section 2, we briefly review the linear mixture model and symmetric alternating direction method of multipliers (symmetric ADMM). In Section 3, we present our formulation, the newly developed algorithm, called double reweighted sparse and graph regularized hyperspectral unmixing (DRSGHU). Section 4 describes the experimental results on simulated hyperspectral data and real hyperspectral data. Section 5 discusses different choices of the graph Laplacian matrix L in graph regularization term and analyzes the influence of the parameters. Finally, some conclusions are provided in Section 6.

Preliminary
In this section, we briefly review the LMM and symmetric ADMM algorithm for subsequent discussion. The former is the basic model of the linear hyperspectral unmixing and the latter is one of the representative optimization methods for unmixing.

Linear Mixture Model (LMM)
LMM is a popular model used in hyperspectral image unmixing, which assumes that the spectral response of a pixel is a linear combination of spectral signatures (called endmembers) [19,20,[30][31][32]. Given an observed spectrum vector of a mixed pixel y ∈ R l×1 , it can be approximated by a nonnegative linear combination of m endmembers, i.e., where A ∈ R l×m represents a spectral library containing l spectral bands and m spectral signatures, x ∈ R m×1 is the abundance vector, and e ∈ R l×1 is the error term (i.e., the noise affecting the measurements). More generally, if the observed data set contains n mixed pixels organized in a matrix Y = [y 1 , y 2 , · · · , y n ] ∈ R l×n , the linear mixture model can be written as where X = [x 1 , x 2 , · · · , x n ] ∈ R m×n is the abundance fractional matrix (each column x j ∈ R m×1 corresponds with the abundance fractions of a mixed pixel y j ∈ R l×1 ), and E = [e 1 , e 2 , · · · , e n ] ∈ R l×n is the noise matrix. In general, two restrictions are imposed on the fractional abundances X in LMM due to physical constraints [5,21,24,26,33]. On the one hand, we often assume X ij 0 where X ij denotes the (i, j)-th entry of X, termed as the abundance non-negative constraint (ANC). On the other hand, it is common to assume 1 T x j = 1 where 1 T ∈ R 1×m is a line vector of 1's compatible with x j , termed as the abundance sum-to-one constraint (ASC). Under some circumstances, the sum-to-one constraint cannot be exactly satisfied because of the variability of the signature [24]. Hence, we do not consider the sum-to-one constraint.
Consider the following linearly separable convex minimization problem of the form: where θ 1 : R n 1 → R, θ 2 : R n 2 → R are convex functions, b ∈ R m , Ω 1 ⊂ R n 1 and Ω 2 ⊂ R n 2 are closed convex sets, A ∈ R m×n 1 and B ∈ R m×n 2 are given matrices. The augmented Lagrangian function of Equation (3) can be written as: where λ is the Lagrangian multipliers and β > 0 is a penalty parameter. Then, the symmetric ADMM scheme is:

Double Reweighted Sparse and Graph Regularized Unmixing Method
In this section, we present the details of the proposed DRSGHU method, including the model formulation and its optimization by symmetric ADMM.

Double Reweighted l 1 Sparse Prior
As each pixel can be expressed as the linear combinations of a few number of endmembers, the abundance vector of a pixel in HSI is usually assumed to be sparse. In other words, most of the entries of the abundance matrix are supposed to be zero. Many methods have been proposed to utilize the sparse property of abundances in recent years. A very straightforward method is to introduce l 0 quasi-norm regularizer to count the nonzero entries in the abundance matrix. However, this is a NP-hard problem. To circumvent this disadvantage, a common surrogate regularizer is l 1 -norm, which is a convex substitution for l 0 regularizer. To further improve the estimation performance over the l 1 -norm and promote the sparsity of the solution, some other regularizers, such as l p -norm regularizer, collaborative sparse regularizer [22], reweighted collaborative sparse regularizer [23] are proposed to solve unmixing problems. In this paper, we propose to use a weighted formulation of l 1 minimization to improve the estimation performance over the l 1 -norm, inspired by [4,44,45]. Considering the fact that a hyperspectral image usually includes fewer endmembers compared with the overcomplete spectral library and the abundance map is inherently sparse, we employ the double reweighted l 1 sparse constraint which defined as Then, we introduce the double reweighted sparse regularizer to unmixing problem, and the new objective function is where the operator denotes the Hadamard product of two variables, and W 1 and W 2 are weighted matrices. The first weight W 1 is introduced to promote the sparsity of endmembers in the spectral library as a hyperspectral image usually contains fewer endmembers compared with the overcomplete spectral library. Due to the unavailability of X, we use the abundance matrix of the current solution to compute the weight matrix W 1 for the next iteration where (W 1 ) k ij ∈ R m×n denotes the (i, j)-th element of W 1 estimated in the k-th iteration and X k (i, :) 1 denotes the l 1 -norm of X k (i, :) (the i-th row of X estimated in the k iteration). > 0 is a small positive constant to avoid singularity. Then, where X 1,0 is the l 0 quasi-norm of vector [ X(1, :) 1 , X(2, :) 1 , ·, X(m, :) 1 ], equivalent to the number of nonzero rows of X [23]. In this way, W 1 X 1,1 has the similar sparsity-promoting property as X 1,0 , which is used to measure the number of nonzero row of abundance matrix X. W 1 X 1,1 equally treats all of the entries in each abundance vector. In other words, the same material appears in every pixel. However, the fraction abundance matrix is intrinsically sparse. Hence, we introduce another weight W 2 to improve sparsity along the abundance vector corresponding to each pixel. Similar to W 1 , we compute W 2 iteratively, i.e., the iteration weighting coefficients of W 2 are based on the current estimation of X k : It suggests that the small weights of W 2 encourage nonzero entries in the estimated abundance matrix, and vice versa. It can be seen that when weighted matrices W 1 and W 2 are all-ones matrices, formulation (P1) become to sparse unmixing model proposed in [21].

Graph Regularizer
The sparse constraint only studies sparsity characteristic of the abundance matrix. However, it ignores the latent geometrical structure of the HSI data. As we know, similar pixels in a spatial or spectral neighborhood have higher probability to share similar abundances [46][47][48]. Therefore, fully making use of this property would enhance the accuracy of the unmixing. We consider the pixel similarity relationship among the neighboring pixels in the spectral space, that is, the abundance vector data set maintains the manifold structure of the spectral feature data set, and incorporates this information into the unmixing problem by the graph Laplacian regularizer. Given a hyperspectral data Y ∈ R l×n with n pixels, we view each pixel y i as a node, then we find its p nearest neighbors for each data point y i and put a edge between y i and its neighbors to construct a weighted graph W on the n data points. There are many different methods to define the weighted graph matrix. Here, we exploit heat kernel weighting, i.e., if nodes i and j are connected, set where σ is the kernel's bandwidth W [48]. Otherwise, if two nodes y i and y j are not connected, we set the weight W ij = 0. Obviously, when y i and y j are very close, the weight W ij tends to be 1, i.e., a larger W ij value indicates that pixels y i and y j are similar.
Through above process, we can get the weight matrix W of size n × n that contains the latent structures in hyperspectral data. We aim to transfer the structure of Y to the abundance space. As the corresponding abundances of any two points y i and y j are expected to maintain the same structure, we adopt the following graph Laplacian regularizer: where Tr(·) represents the trace of a matrix and D is a diagonal matrix whose entries are column (or row, since W is symmetric) sums of W, D ii = ∑ j W ij . L = D − W is called graph Laplacian [46]. By minimizing R 2 (X), we expect that if two data points y i and y j in the HSI data space are close, then x i and x j in the abundance space also should be close.

The Double Reweighted Sparse and Graph Regularized Unmixing Model
To simultaneously utilize structure information of HSI data and the sparse property of abundance matrix, we propose the double reweighted sparse and graph regularized unmixing model which can be expressed as follows: where α 1 > 0 and α 2 > 0 are regularization parameters. The first term is the reconstruction error, the second term is introduced to capture the correlation information between the abundance vectors that is hidden in the HSI data, and the third term is exploited to promote the sparseness of the abundance matrix.
The optimization problem (P2) can be reformulated as Then, Equation (13) can be further expressed as where , and B and C are, respectively, given by By introducing the Lagrangian multiplier Λ = (Λ 1 ; Λ 2 ) ∈ R 2m×n , the augmented Lagrangian function of Equation (13) can be written as where β > 0 is a penalty parameter. Note that for any fixed value of X, the variables V 1 and V 2 are decoupled. Therefore ,we can solve V 1 and V 2 on their corresponding subproblems separately.
In the framework of symmetric ADMM, we use the following iteration scheme to approach the solution of Equation (15): For the X-subproblem, we get that The solution of X k+1 is explicitly given by: For the V 1 -subproblem, we solve the following equation: the solution of which is Denoting W = W 2 W 1 , then the V 2 -subproblem can be written as As the V 2 subproblem is component-wise separable, the solution of V 2 is , "•" represents the point-wise product and sign(x) is the signum function.
In summary, the proposed double reweighted sparse graph regularized hyperspectral unmixing algorithm (DRSGHU) is described in Algorithm 1.

Experimental Results
In this section, we test DRSGHU for hyperspectral unmixing using two simulated hyperspectral data and one real data. Both simulated and real hyperspectral data are generated to evaluate the effectiveness of the proposed algorithms. In the synthetic data experiments, we compare our method with some representative methods: (1) SUnSAL [21]; (2) CLSUnSAL [22]; (3) DRSU-TV [4]; (4) GLUP-Lap [28]; (5) Graph-regularized unmixing method (without sparse constraint on abundances), abbreviated as GraphHU; and (6) Sparse graph-regularized unmixing method (set W 1 = W 2 as all ones matrices in problem (P2)), abbreviated as SGHU. For quantitative comparison, the signal-to-reconstruction error (SRE) is introduced to measure the quality of the estimated fractional abundances of endmembers. SRE is defined as: , measured in dB: SRE(dB) ≡ 10 log 10 (SRE), where E denotes the expected value, see [22]. Here, X true is the true fractional abundance of endmembers andX is the estimated abundance of endmembers. Generally speaking, when the SRE becomes bigger, the estimation approximates the truth better. For initialization, we set V 0 1 = 0, V 0 2 = 0, Λ 0 1 = 0, Λ 0 2 = 0, where 0 is a zero matrix with size m × n.
The parameters r and s are set to be r = 0.9, and s = 0.9, the same as [51]. The stopping criteria is set as or the maximum iteration number is reached, which is set to be 500 in our experiments.
To achieve the best performance, we test the proposed method for different combinations of α 1 , α 2 and β according to different noise levels.
All experiments are performed under windows 10 and Matlab R2016b running on a desktop with an Intel(R) Core(TM) i7-7700 CPU at 3.60 GHz and 16.0 GB of memory.

Experiments on Simulated Data
In the synthetic data experiments, we evaluate the unmixing performance of DRSGHU in situations of different signal-to-noise ratios (SNR ≡ 10 log( AX 2 F / N 2 F )), different spectral libraries, and different simulated data. The simulated data are contaminated by Gaussian noise with different noise levels of the SNR: 20, 30 and 40 dB. We consider the following two spectral libraries in our simulated data experiments: • A 1 ∈ R 100×120 is generated from a library of 262 spectral library signatures generally found on satellites from the National Aeronautics and Space Administration Johnson Space Center (NASA JSC) Spacecraft Materials Spectral Database [52], with 100 spectral bands. • A 2 ∈ R 224×240 is generated from a random selection of 240 materials from the U.S. Geological Survey (USGS) digital spectral library (splib06) [24]. The reflectance values are measured for 224 spectral bands uniformly distributed in the interval 0.4-2.5 µm.
In our simulation experiments, we test two different data sets. The details are presented in the following examples.

Example 1.
In this example, we test the performance on simulated data cube 1 (DC1). This simulated data cube contains 100 × 100 pixels generated from nine randomly selected signature from A 1 , and each pixel has 100 bands. Then, a Gaussian white noise with different SNR (20 dB, 30 dB or 40 dB) is added to DC1. The true fractional abundances of nine piecewise smooth with sharp transitions endmembers are shown in Figure 1. Table 1 lists the SRE results for all simulated data sets. From this table, we observe that DRSGHU gets higher SRE than GLUP-Lap, GraphHU, SGHU, SUnSAL for all noise levels. It indicates that both the double reweighted l 1 sparse prior and graph regularizer mutually enhance the abundance estimation accuracy. As the level of noise decreases, the rate at which DRSGHU improves comparing to SUnSAL increases. This is because that the observed HSI data contains less noise, the graph Laplacian matrix becomes more reliable. In most cases, DRSGHU achieves the best performance and DRSU-TV achieves second best results except when the level of noise is 20 dB. One possible reason is that the graph Laplacian L is obtained from the very noisy hyperspectral image rather than the clean image, making the graph Laplacian less reliable. Figure 2 shows the estimated abundance maps by DRSGHU for a noise level of 30 dB. It can be seen that the unmixing results by DRSGHU are very close to the true fractional abundances presented in Figure 1. We further present the relative errors between the estimated fractional abundances and the true fractional abundances of all the endmembers in Table 2 for detailed performance evaluation of our method. In Table 2, the relative errors for the each endmember obtained by DRSGHU is smallest in most cases, and the average relative error of all the endmembers is also smallest. It indicates that the estimation by DRSGHU approximates the truth best.
Furthermore, we present the estimated fractional abundances of all endmembers in the library A 1 for the first 150 pixels under a noise level of 30dB in Figure 3 to appreciate the sparsity of the fractional abundances. We observe that DRSUTV and DRSGHU produce less outliers than other methods. As both methods employ the double reweighted sparse prior, they can get relative sparser solution. This suggests that the double weighted matrices promote the sparsity of abundance maps. We also see that the unmixing results by SGHU are much better than SUnSAL. This is because the graph regularizer is imposed onto the sparsity constraint, where the latent structure of HSI data can be utilized. It indicates that incorporating graph regularizer in the unmixing problem improves the abundance estimation accuracy.

Example 2.
In this example, we test our method on simulated data cube 2 (DC2). This simulated data cube contains 75 × 75 pixels generated from five randomly selected signature from A 2 . DC2 has 100 bands per pixel [5]. The true abundances of five endmembers are shown in Figure 4. To illustrate the effectiveness of the proposed method, we present the estimated abundance of endmember 4 (EM4) by different methods in Figure 5. We can observe that, when the noise level is very high (SNR = 20 dB), all methods performs not very good and the unmixing results still contain some noise. With the decreasing of noise level, the unmixing results look much more satisfactory. Although the unmixing results by DRSU-TV are much better than SUnSAL and CLSUnSAL, some unexpected transitions occurred and some squares are not correctly estimated, especially when the noise level is high as shown in Figure 5. This is possibly because that DRSU-TV links a pixel with only four neighbors and tends to produce piecewise smooth results. On the contrary, the unmixing results by DRSGHU is visually much more acceptable and perfectly recover the abundances. For quantitative comparison, we present the relative errors between the estimated fractional abundance EM4 and the true one at different noise levels in Table 3. It can be seen that DRSGHU achieves smallest relative errors of EM4.

Experiments on Real Data
In this subsection, we use the well-known Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) Cuprite data (http://aviris.jpl.nasa.gov/html/aviris.freedata.html) set in our real data experiments. We use the portion corresponding to a 350 × 350-pixel subset with 188 spectral bands in this experiment (see [5] for more details). A mineral map produced in 1995 by USGS, in which the Tricorder 3.3 software product was used to map different minerals present in the Cuprite mining district (http://speclab.cr.usgs.gov/spectral.lib06). Essential calibration was performed to mitigate the mismatches between the real image spectra and the spectra available in the library before unmixing the AVIRIS Cuprite hyperspectral data [5].The regularization parameters α 1 and α 2 are set to be 10, 0.1, and the penalty parameter β is set to 0.01 for the proposed method.
The estimated fractional abundance maps for Alunite, Muscovite, and Chalcedony by DRSU-TV, GraphHU, SGHU and DRSGHU are shown in Figure 6. In this figure, we also present the spatial distribution of these materials extracted from the Tricorder software product providing a reference for visual assessment. All compared algorithms produce consistent unmixing results. GraphHU, SGHU and DRSGHU reveal more spatial details for Alunite and Muscovite, and produce fewer outliers for Chalcedony comparing to DRSU-TV.

The Graph Laplacian Matrix L
From our numerical experiments, we notice that it is very important to choose appropriate graph Laplacian matrix L to achieve good performance, as L reflects the similarity relationship among the spectral pixels. There are many methods to construct the graph Laplacian matrix L, most of which used the observed data directly to define the graph Laplacian, such as the ones mentioned in Section 2. However, when the observed data are highly corrupted by the noise, the graph Laplacian matrix directly constructed by observed data may not be reliable. Therefore, it is necessary to preprocess the observed data Y to get better estimation of L. For example, some denoising methods can be applied to observed data to get a good graph Laplacian matrix. Table 4 lists unmixing results by DRSGHU after applying the gaussian filter or bilateral filter [53] to the observed data when the SNR = 20 dB. The "reference" column represents the SREs of unmixing results using clean data (not contaminated by the noise) to construct the graph Laplacian matrix. It clearly demonstrates the advantage of pre-denoising before constructing L. Some other denoising methods could also be applied to get better performance, such as non-local total variation, BM3D, and sparse representation [2,26].

Parameter Analysis
In this section, we analyze the influence of the parameters α 1 , α 2 and β for DRSGHU. For simplicity, we take the simulation data DC1 under different noise levels as an example to illustrate the influence of parameters. We separate these three parameters into two groups, one is the regularization parameters α 1 and α 2 , and another group is penalty parameter β. To analyze the influence of the regularization parameters, we fix the penalty parameter β. Figure 7 presents SRE as a function of the regularization parameters of α 1 and α 2 . A general trend observed is that the optimal regularization parameters tend to be bigger when the noise level is higher. In addition, one can see that the surface is not notably steep, which suggests that the algorithm is quite robust around the optimal parameters. To analyze the influence of the penalty parameter, the regularization parameters α 1 and α 2 are fixed. Figure 8 presents SRE as a function of the penalty parameter β. It can be observed that the optimal penalty parameter becomes bigger with higher noise level. Moreover, we can see that SREs are relatively stable around the optimal penalty parameter under different noise level, also suggesting the robustness of our algorithm.

Convergence Analysis
The theoretical convergence analysis for the reweighted norm minimization problem is difficult to establish [44]. Nevertheless, we can observe that the relative error between restored abundance X and the truth abundance X true becomes smaller with the increasing of iteration through the experiments. We take simulated data DC1 as an example to illustrate the convergence of our method. In this test, we set the maximum iteration number as 1000. In Figure 9, one can see that the relative error changes little after 200 iterations, which indicates the stability of our method.

Conclusions
In this paper, we propose a double reweighted sparse and graph regularization based formulation for hyperspectral unmixing. Then, we solve the proposed optimization problem by symmetric alternating direction method with multipliers. The proposed method takes both the structure information of HSI data and sparse property of abundance matrix into consideration. The graph regularization term is introduced to capture the structure information of HSI data, and the double reweighted l 1 -norm is used to describe the sparsity property of endmembers and fractional abundance maps. The experiments indicate that graph regularizer and weighted sparse regularizer both improve the unmixing performance. It is worth pointing out that the graph Laplacian matrix obtained from the observed data is important to the performance of unmixing result. When the noise level is very high, we could first pre-denoise the hyperspectral data before computing the graph Laplacian matrix to get better performance.
Some issues are worth being investigated in the future. We only consider the Gaussian noise in this work. In real applications, HSIs are usually degraded by different types of noise, such as hybrid Gaussian and stripe noise. We will study how to simultaneously remove these types of noise and unmixing. In addition, we utilize the graph regularizer to capture the structure information of hyperspectral images. Some other regularizers, such as nonlocal TV and low rank tensor, may also be well-suited for revealing spatial and spectral information of hyperspectral images. It also needs to be studied in further.