Hyperspectral Unmixing with Gaussian Mixture Model and Low-Rank Representation

Gaussian mixture model (GMM) has been one of the most representative models for hyperspectral unmixing while considering endmember variability. However, the GMM unmixing models only have proper smoothness and sparsity prior constraints on the abundances and thus do not take into account the possible local spatial correlation. When the pixels that lie on the boundaries of different materials or the inhomogeneous region, the abundances of the neighboring pixels do not have those prior constraints. Thus, we propose a novel GMM unmixing method based on superpixel segmentation (SS) and low-rank representation (LRR), which is called GMM-SS-LRR. we adopt the SS in the first principal component of HSI to get the homogeneous regions. Moreover, the HSI to be unmixed is partitioned into regions where the statistical property of the abundance coefficients have the underlying low-rank property. Then, to further exploit the spatial data structure, under the Bayesian framework, we use GMM to formulate the unmixing problem, and put the low-rank property into the objective function as a prior knowledge, using generalized expectation maximization to solve the objection function. Experiments on synthetic datasets and real HSIs demonstrated that the proposed GMM-SS-LRR is efficient compared with other current popular methods.


Introduction
In the last few decades, Hyperspectral image (HSI) has received considerable attention in the field of earth observation and geoinformation science.With the wealth of spatial and spectral information, HSI has been successfully applied in many applications, such as spectral unmixing, environment monitoring, matching and object classification [1][2][3][4][5][6][7][8].However, due to the low spatial resolution of current HSI sensors and the mixed effects of the ground surface, these factors will seriously affect the accurate interpretation of the image content.In this case, spectral unmixing (SU) problem has become a major issue for the deep development of HSI.
The information of hyperspectral images can be simplified by the linear mixing model (LMM), which assumes that the physical region corresponding to a pixel contains several pure materials.Hence, the observed spectra y n ∈ R B , n = 1, ..., N (B is the number of wavelengths and N is the number of pixels) is a (non-negative) linear combination of the pure material (called endmember) spectra m j ∈ R B , j = 1, ..., M (M is the number of endmembers) [9], i.e., where α nj is the proportion (called abundance) for the jth endmember at the nth pixel (with the positivity and sum-to-one constraint) and n n ∈ R B is additive noise.Here, the endmember set m j : j = 1, ..., M is fixed for all the pixels.Many endmember detection algorithms are under the pixel purity assumption to extract the pure endmembers, such as the pixel purity index [10], successive projection algorithm [11], and vertex component analysis (VCA) [12].Other algorithms assume that all the pixels lie in a convex hull in a high-dimensional subspace: N-Finder [13] and iterative constrained endmembers (ICE) [14].
In our proposed abundance estimation method, the pure endmembers are assumed to be known or they can be identified from the target image by one of these endmember detection techniques.However, in practice, the LMM may not be an appropriate model in many real scenarios.Even for a pure pixel that only contains one material, over the whole image its spectrum may not be consistent.This is due to several factors such as intrinsic variability, atmospheric conditions and topography.Equation ( 1) can be generalized to a more abstract form, y n = S(m j , α nj : j = 1, ..., M), which is called nonlinear mixing models (NLMMs).For example, in [15], the generalized bilinear model (GBM) generalizes the LMM by introducing bilinear terms to handle the vegetation case with taking into account the multipath effects.In that case, various representative abundance estimation algorithms have been proposed based on different nonlinear model assumptions, such as least squares [16], kernel-based least squares [17,18], and Bayesian methods [15].A panoply of nonlinear models can be found in the review article [19].The reason we note those models is that the endmember set is still assumed to be fixed.
Although NLMM models have abounded recently, it is still not appropriate to account for all the scenarios.LMM has also been widely used due to its simplicity and physical meaning.To model real scenarios more accurately, researchers have taken another route by generalizing Equation (1) to where m nj ∈ R B : j = 1, ..., M, n = 1, ..., N. We note that the endmember set in that case is not fixed and the endmember spectra could be different for each pixel y n .This is called endmember variability.
Given y n , inferring m nj , α nj is a much more difficult problem than inferring m j , α nj in Equation (1).
In the review paper [20], the methods taking into account the endmember variability can be expressed in two categories: (1) endmembers represented as a discrete set; or (2) endmembers represented using a continuous distribution.In the first category, the endmember spectra are modeled as a series of spectra clusters, and expressed as discrete sets in mathematical form.The multiple endmember spectral mixture analysis (MESMA) is one of the widely used methods of this category [21], which tries every endmember combination to search the minimum mean square error in the discrete set as the final calculation result.There are many variations to the original MESMA, such as the multiple-endmember linear spectral unmixing model (MELSUM) [22], Auto-Monte Carlo Unmixing (AutoMCU) [23,24] and the Bayesian spectral mixture analysis (BSMA) [25].Besides those variants, there are also some other set based methods, e.g., endmember bundles [26] and band weighting or transformation approaches [27,28].However, those methods mentioned above have a common disadvantage that, when the spectral library is large, their complexity will increase exponentially, resulting in extreme computational inefficiency difficulties.
The second category usually takes statistical method to model the endmember distribution.Specifically, it is assumed that the endmembers of each pixel are sampled from a probability distribution, and hence embrace large libraries while being numerically tractable.In [29], Eches et al. proposed the normal compositional model (NCM).This model is the early representative work in this direction, which assumes the endmembers for each pixel are sampled from unimodal Gaussian distribution (primarily due to mathematical simplicity).Due to the complexity of the model's prior knowledge and hyperparameters, the maximum a posteriori (MAP) constructed is often a non-convex function.Then, NCM uses different optimization approaches, such as expectation maximization [30], sampling method [31], and particle swarm optimization [32], to determine the hyperparameters.However, the NCM allows for endmember samples to range outside the interval of [0, 1] and in the real scenarios the endmember distribution could be skewed.Hence, Du et al. proposed a Beta compositional model (BCM) to model the endmember variability in the HSI unmixing [33].However, the true distribution may not be well approximated by either a Gaussian or Beta distribution.In that case, Zhou et al. proposed the Gaussian mixture model (GMM) to solve the LMM by considering endmember variability [34].GMM method assumes that the endmember m nj follows the GMM distribution and noise n n follows the Gaussian distribution: p(n n ) := N (n n |0, D), where D is the noise covariance matrix, and with proper abundances constraint under the Bayesian framework to lead the conditional density function to a standard MAP problem.
The performance of those methods in this category is often dependent on the initial value of parameters, but does not rely on the large-scale spectral database, which is the research hotspot for the endmember variability problem.However, neither method mentioned above takes into account the possible spatial correlation between local pixels.Researchers in [2,35,36] proved that the spatial correlations between the observed pixels can enhance the performance of spectral unmixing.On the other hand, this strategy has great advantages of easily generalizing the Bayesian algorithms, which is introduced in [37,38].In that vein, in an attempt to achieve better abundance estimation results, Zhou's model (GMM) [34] uses the proper smoothness and sparsity prior constraints on the abundances, while Iordache et al. [39] tried to find the minimum difference of the abundance between the adjacent pixels by the total variation (TV) constraint.However, those prior constraint assumptions require a rather strict assumption that the abundance of the local pixel should be piecewise smooth, which means that the mixing pixel and its associated abundance should be similar for the adjacent pixels.When the pixels lie on the boundaries of different materials or the inhomogeneous region, the abundances in the neighboring pixels do not have those prior constraints.Driven by those considerations, to further exploit the spatial data structure, we take superpixel segmentation (SS) in the first principal component of HSI to get the homogeneous regions.Then, the hyperspectral image to be unmixed is partitioned into several homogeneous regions (or classes) in which the abundance vector has the same first and second order statistical moments (means and covariances).The shape and size can be accustomed on the basis of different spatial structures, and each superpixel can be seen as a non-overlapping region.
Furthermore, in those homogeneous regions, there is a high degree of correlation among the spectral signatures of the neighboring pixels.The high spatial similarity usually implies the low-rank property in the abundance map.Qu et al. [40] proved that the abundance matrix has low-rank property in homogeneous region and constructed an unmixing model by using low-rank property to verify the effectiveness in homogeneous region unmixing.Giampouras et al. [36] proposed a novel low-rank representation (LRR) for HSI unmixing, which jointly obtains a representation for all the data under a global lowest rank constraint.In that vein, we take LRR to further exploit the spatial information between the local pixels in the homogeneous region, using the Bayesian framework to construct the objective function together seeking the lowest rank.
Thus, in this paper, in an attempt to achieve better abundance estimation performance and fully consider the possible spatial correlation between local pixels, we propose a novel GMM unmixing method based on SS and LRR, which is called GMM-SS-LRR.Firstly, we adopt the SS algorithm to cut the HSI into different regions.In these regions, pixels are highly spatial correlated, the mixing pixel and its associated abundance should be similar for the adjacent pixels and the abundance map has low-rank property.Secondly, considering the endmember variability phenomenon, under the Bayesian framework, we use GMM to formulate the unmixing problem, and put the low-rank property into the objective function as a prior knowledge.Finally, we use generalized expectation maximization (GEM) to solve the objective function [41].
In summary, the main contribution in this paper are two folds: (1) We take SS in the first principal component of HSI to get the homogeneous regions, which can help the GMM method make pixels and their associated abundances similar for the adjacent pixels, since directly using the GMM method to unmixing will make the estimated distribution of the clusters cannot fit the ground truth well and using SS can enhance the performance of the endmember estimation.
(2) We take LRR for the abundance estimation problem, which can further exploit the spatial information between the local pixels.Because after taking SS, the regions will be homogeneous, and the high spatial correlation of the data implies the low-rank property of the abundance matrix, using the LRR property can better capture the spatial data structure by seeking the lowest rank representation, and also get better abundance estimation.
The rest of this paper is structured as follows.In Section 2, we briefly introduce the related GMM method.In Section 3, we describe the proposed model GMM-SS-LRR and the GEM for solving and optimization are introduced.In Section 4, we evaluate the performances of the proposed GMM-SS-LRR and the compared state-of-the-art algorithms on the synthetic datasets and real HSIs.We conclude this paper in Section 5.

Related Models
Since the proposed model is closely related with the GMM method, we briefly introduce the GMM in this section.
GMM method [34] is an LMM by considering the endmember variability, use a mixture of Gaussians to approximate any distribution of the endmember and can be classified in the second category mentioned above (endmembers represented using a continuous distribution).GMM model assumes the endmember p(m nj ) follows a mixture of Gaussians, and has the density function as follows: subject to π jk ≥ 0, ∑ K j k=1 π jk = 1, with K j being the number of components, π jk (µ jk ∈ R B ) or Σ jk ∈ R B×B being the weight (mean or covariance matrix) of its kth Gaussian component, Θ := π jk , µ jk , Σ jk : j = 1, ..., M, k = 1, ..., K j , m nj : j = 1, ..., M are independent.The noise n n also follows the Gaussian distribution, which has the density function p(n n ) := N (n n |0, D), where D is the noise covariance matrix and then obtains the distribution of y n by the equation y n = ∑ M j=1 m nj α nj + n n , which turns out to be another mixture of Gaussians, we can obtain the density function of y n as follows: where are defined by: More specifically, in the prior of the abundances of p(A), Zhou's model [34] assumes the abundances A have the proper smoothness and sparsity prior constraints.The density function of the abundances A can be generalized as follows: where L is a graph Laplacian matrix constructed from w nm , n, m = 1, . . ., N with w nm = e ||y n −y m || 2 /2Bη 2 for neighboring pixels and 0 otherwise.Tr(•) is the trace of the matrix, and with β 1 controlling smoothness and β 2 controlling sparsity of the abundance maps.

GMM Unmixing with Superpixel Segmentation and Low-Rank Representation
In this section, we describe the specific steps in implementing the superpixel segmentation at first, then introduce the GMM unmixing based on the low-rank representation, finally develop an GEM for solving the proposed unmixing method GMM-SS-LRR.

Formulation of the Proposed GMM-SS-LRR
When considering the endmember variability, the current popular unmixing methods which model the endmember following the probability distribution have been mentioned above: normal compositional model (NCM) [29], beta compositional model (BCM) [33] and Gaussian mixture model (GMM) [34].Those methods all ignore the local spatial correlation of HSI.However, the spatial correlation of HSI is very important in the analysis of HSI, and researchers in [2,35,36] proved that the spatial correlations between the observed pixels can enhance the performance of spectral unmixing.In [34], the abundance A has been modeled as assumed having proper smoothness and sparsity prior constraints.From Equation ( 6) we can get that the smoothness is modeled by the well known symmetric positive semidefinite graph Laplacian matrix.This constraint is based on the assumptions that both the mixing material and its associated abundance should be similar for the adjacent pixels.When the pixels lie on the boundaries of different materials or the pixels are concentrated on different components, the prior will not have the smoothness property.The pixels in that region are not pure; this also cannot conduct a sparse prior constraint.In that vein, to better incorporate the abundance prior constraints, we take SS to cut the HSI into several homogeneous regions.The SS methods are originally designed for the visible images, since the whole original HSI cube usually has hundreds of bands, and the SS cannot be directly used for the HSI cube; hence, we take PCA to get the first principle component of HSI which is used as the base image when using the SS.Since entropy rate (ER) possesses the advantage of the computational efficiency [42], and the superpixels incline to have common features and similar sizes [43].Therefore, in this paper, we adapt the ER to generate superpixels.
After adopting the PCA and SS, the HSI is segmented into several homogeneous regions, and, as shown in Figure 1, each superpixel can be seen as a non-overlapping region, the shape is adaptive and size can be accustomed on the basis of different spatial structures.Thus, each superpixel in these regions, pixels are highly spatial correlated, and the abundance matrix has the low-rank property.According to Horn and Johnson [44], we further propose to exploit its rank property by the following theorem: For the real HSI, the extracted endmembers E are generally distinct from each other and the number of bands L is usually larger than the total number of the endmember R, which makes the rank(E) = R, and rank(Z) = k ≤ min(R, N); according to Theorem 1, we can get rank(Z) = rank(A), since adopting the PCA and SS, the original HSI has been cut into different regions, and the columns of Z are highly correlated, which means that the matrix Z is low rank, thus the abundance matrix A in each homogeneous regions is also low rank.To use this property, together with the proper smoothness and sparse prior constraints, and after conducting principal component analysis (PCA) and SS on the HSI, the priors for abundance A based on the proposed GMM-SS-LRR can be mathematically written as follows: However, Equation ( 8) is a highly nonconvex optimization problem or even NP-hard, which is very difficult to work out.According to Liu et al. [45], we use the nuclear norm to substitute the rank function, and the nuclear norm is widely used as a surrogate to the rank function.Thus, the prior constraints on A can be reformulated as: where ||A|| * denotes the nuclear norm of the matrix A defined as Then, based on GMM method analysis in the Section 2, given all the abundances A := [α 1 , ..., α N ] ∈ R N×M (α N := [α 1 , ..., α nM ] T ) and GMM parameters Θ, we can model the conditional distribution of the pixels Y := [y 1 , ..., y N ] ∈ R N×B , together we also assume the noise D follows the Gaussian distribution, and assuming the conditional distributions of each y n are both independent, then, using the result in Equation ( 4), given A, Θ, D becomes: On the basis of the conditional density function, the priors and Bayes' theorem, we can obtain the posterior given by p(A, Θ|Y, D) ∝ p(Y|A, Θ, D)p(A)p(Θ), (12) where p(Θ) is assumed to be a uniform distribution.Since we have obtained the distribution of p(Y|A, Θ, D) (Equation ( 11)), p(A) (Equation ( 9)), and maximizing p(A, Θ|Y, D) is equivalent to minimizing − log p(A, Θ|Y, D), we can obtain the objective function as follows: where ||A|| * and µ nk , Σ nk are defined in Equation ( 5).

Optimization of the Proposed GMM-SS-LRR
There are numerical methods in estimating the parameters of GMM method, such as expectation maximization (EM) [46,47] and projection-based clustering [48].However, in our method, each pixel y n is generated to be a different GMM determined by the coefficients α n , which is a more challenging problem to estimate the parameters.Since EM is more flexible and can be seen as a special case of majoriziation-minimization algorithms [49], we take this method to approach.Because the parameters A, Θ both need to be updated in each M step, we adopt the method called generalized expectation maximization [41], where both parameters A, Θ are updated sequentially as long as the complete data log-likelihood increases.
Then, following the EM step, we calculate the posterior probability of the latent variable γ nk given the observed data and old parameters in the E step: In the M step, we maximize the expected value of the complete data log-likelihood.Putting the priors in the Bayesian formulation, the final objective function ε M we need to minimize becomes the following: where the ε prior are defined in Equation (13).The weight of the gaussian mixture π k can be updated as In the next step, we need to update the µ jk , Σ jk and A. Using Equation ( 5), we can obtain the derivatives of the objective function ε M in Equation (15) with respect to where δ lk j = 1 when l = k j and 0 otherwise, λ nk ∈ R B×1 and Ψ nk ∈ R B×B are given by and the K = L − β 2 β 1 I N (suppose β 1 = 0), U, V are the singular value decomposition of the abundance A.
where U ∈ R n×r and V ∈ R m×r are matrices with orthogonal columns and the singular values It is convenient to represent the derivatives in matrix forms.Considering the multiple summations in Equations ( 17)-( 19), we can write them as where If the optimization problem has no constraint condition, we can let ∂ε M ∂µ jl = 0, ∂ε M ∂Σ jl = 0 and ∂ε M ∂A = 0 to get the minimum of ε M .However, the Σ jk has the positive definite constraint and α nj has the non-negativity (ANC) and sum-to-one (ASC) constraint, which make minimizing ε M become very difficult.Thus, we follow the method [24].In each M step, using Equations ( 23)-( 25) and decreasing the objective function by project gradient decent.In estimating the number of components K j from the data, we use the Kullback-Leibler (KL) divergence to minimize difference with true density function and the estimated density function, and can be mathematically written as follows: where f m j (y|Θ j ) is the estimated density function with Θ j := π jk , µ j k, Σ j k : k = 1, ..., K j , g m j (y) is the true density function, N j denotes number of the jth endmember and y j n denotes nth element for the jth endmember.Then, we take the cross-validation-based information criterion CVIC [50,51] to correct for the bias.To sum up, the detailed procedure for solving Equation ( 15) is listed in Algorithm 1.

Algorithm 1 Solving Equation (15) with EM.
Input: Collect mixed pixel matrix Y, endmember E, the parameter of smoothness and sparse constraint β 1 , β 2 , and the parameter of low-rank property β 3 ; Output: The estimated abundance matrix A; 1: Implement PCA and set α nk ← (R k R T k + I M ) −1 R k y n as initialization; 2: Using the KL divergence to get the number of components K j ; 3: Using CVIC to correct for the bias; 4: while not converged do 5: E step: 7: Update γ nk and A; 8: end while 9: Return A.
Given the complexity of EM scheme in solving GMM-SS-LRR, the most time-consuming step is estimating the abundances A in each iteration.Thus, the EM complexity of GMM-SS-LRR has spatial complexity O(NB 2 ) and time complexity O(NB 3 ).The detailed procedure of the proposed GMM-SS-LRR is listed in Algorithm 2.

Algorithm 2 Proposed GMM-SS-LRR.
Input: Collect mixed pixel matrix Y, endmember E, the parameter of smoothness and sparse constraint β 1 , β 2 , and the parameter of low-rank property β 3 and the number of the superpixels S; Output: The estimated abundance matrix A; |Recover A with Algorithm 1; 5: end for 6: Obtain A from each homogeneous region; 7: Return A.

Experimental Result
We evaluated the performances of the proposed GMM-SS-LRR and compared it with the state-of-the-art algorithms on synthetic datasets and real HSIs.To demonstrate the efficiency of the proposed GMM-SS-LRR, we mainly compared it with three other algorithms, namely the normal compositional model (NCM), the beta compositional model (BCM) (BCM takes the spectral version with quadratic programming) and GMM, because these three methods are both based on LMM and model the endmember using continuous distribution.For the proposed method GMM-SS-LRR and GMM, the original image data were projected to a subspace with 10 dimensions to speed up the computation for abundance estimation.Since NCM is a supervised algorithm, we took the ground truth pixels as input, and modeled them by an unimodal Gaussian distribution and obtained the abundance maps by maximizing the log-likelihood.BCM was also implemented as a supervised algorithm; ground truth pure pixels were again taken as input and the results were the abundance maps.In the following experiments, we implemented the algorithm in MATLAB and experiments were performed on a laptop with 2.6-GHz Intel Core CPU, 16-GB memory.All parameters of β 1 (smoothness prior constraint) and β 2 (sparse prior constraint) were set following the GMM method [34] (in the synthetic dataset, the parameters were set to β 1 = 0.1, β 2 = 0.1; and in the real HSIs, the parameters were set to β 1 = 5, β 2 = 5 ), and, for the proposed GMM-SS-LRR method, the low-rank property in the synthetic dataset was set to β 3 = 0.1 and in the real HSIs was set to β 3 = 1.The number of superpixels for both synthetic dataset and real HSIs was set to S = 5.
For comparison of abundances, we calculated the root mean squared error (RMSE) for abundance error, which is defined as follows: where α GT nj are the ground truth abundances and α est nj are the estimated values.Since only some pure pixels were identified as ground truth in the real HSI dataset, we calculated error given the pure pixel index set I.

Synthetic Datasets
For the synthetic data experiments, the mean endmember spectra we used were randomly selected from the ASTER spectral library [24] with slight constant changes (their original spectra are shown in Figure 2b), which determined a spectral range from 0.4 µm to 14 µm.The size of the synthetic was 60 × 60 pixels and constructed from five endmembers (limestone, basalt, concrete, conifer and asphalt), whose spectral signatures are highly differentiable.The covariance matrices Were constructed by α 2 jk I B + b 2 jk µ jk µ T jk where µ jk is a unit vector controlling the major variation direction.We made the first material as background, and the procedure of generating the abundance maps followed Zhou et al. [52]: for each other material (not as background), 600 Gaussian blobs were randomly placed in the corner, and the shape, width and location were sampled from Gaussian distributions.The Gaussian noise was added to generated pixels; the noise sigma was set with σ Y = 0.001.Figure 2a shows the resulting color images by extracting the bands corresponding to wavelengths 488 nm, 556 nm, and 693 nm.The endmember spectra we used to generate the synthetic data are shown in Figure 3, where we can see the centers and the major variations of the Gaussians.
Figure 2c,d shows the first principle component of the synthetic dataset and the superpixel map we used, respectively.The synthetic image was correctly partitioned into several regions, whose shape and size were accustomed on the basis of different spatial structures, meaning there was a high degree of correlation among the spectral signatures of the neighboring pixels, hence we could take low-rank representation to further exploit the spatial information and improve the abundance estimation.The abundance maps of our method and other comparison algorithms for the synthetic dataset are shown in Figure 4; since the materials were randomly placed in the corner, the abundance maps of the last four materials (basalt, concrete, conifer and asphalt) look relatively clean and less cluttered.Comparing them with the ground truth, we can see that all algorithms correctly estimated the abundance maps scatter.However, for BCM and NCM, although the abundance maps performed relatively clean as well, the shape and size were quite different compared with the ground truth.This means that many pixels with original abundance of 1 were predicted to have abundance of 0, which caused a relatively large value of RMSE.As for the GMM method, although the abundance maps appeared to have clutter point, the shape and size were similar to the ground truth.Hence, the estimating accuracy was higher than NCM and BCM.The GMM-SS-LRR algorithm with SS and low-rank property was much closer to the ground truth map.The quantitative analysis of those four algorithms is shown in Table 1.
The histograms of the synthetic pure pixels for the five materials are shown in Figure 5. Since the BCM algorithm was not modeled as Gaussian distribution, the probability of each distribution was only compared among GMM-SS-LRR, GMM and NCM algorithm.The histogram was the statistical probability value of the pure materials for the five materials (when projected to one-dimensional space determined by performing PCA).The probability of each distribution was calculated by multiplying the value of the density function at each bin location with the bin size.Figure 5 shows that the graph of NCM is similar to GMM and GMM-SS-LRR when the distribution of pure pixels was generated by an unimodal Gaussian.Nevertheless, for basalt and concrete, GMM and GMM-SS-LRR provided a more accurate estimation because NCM was modeled by a single Gaussian hypothesis.When comparing GMM and GMM-SS-LRR, the distribution was similar for limestone, concrete, conifer and asphalt.However, for basalt, GMM-SS-LRR fitted the histograms better than GMM, as the first step in GMM approach was to separate the library into several groups, each of which representing a material and clustered into several centers.The size of each cluster affected the probability of picking its center to a large extent.When directly using the GMM method for unmixing, sometimes the estimated distribution of the cluster did not fit the ground truth well.Taking SS in the first principle component of HSI helped the GMM better separate the clusters.The quantitative analysis of these three algorithms is shown in Table 2.We calculated the probability value in each histograms between the ground truth and estimate value using Equation (29).

Mississippi Gulfport Datasets
The dataset was collected at the University of Southern Mississippi's Gulfpark Campus, and is a 289 × 284 image with 72 bands corresponding to wavelengths from 0.368 µm to 1.043 µm.The spatial resolution is 1 m/pixel.The scene contains several man-made and natural materials including sidewalks, roads, various types of building roofs, concretes, shrubs, trees and grasses.To better compare the performance between our method and GMM, we chose the same ROI area followed in [34].The ROI was a 58 by 65 pixel area containing five materials.The original RGB image and the selected ROI are shown in Figure 6a,c.
Figure 7 shows the abundance maps comparison in the Gulfport dataset.Comparing them with the ground truth (the first row of Figure 7), we can see that BCM failed to estimate the pure pixels of tree, even though ground truth pure pixels were used for training.For example, the fourth and fifth abundance maps of BCM showed that the pixels of tree were mixed with grass and asphalt.For BCM and NCM, we did not use PCA to get the main information while using the original HSI dataset as input.However, they still performed poorly on this dataset.The result of GMM-SS-LRR not only showed sparse abundances for the region but also interpreted the boundary as a combination of neighboring materials.Although the results of GMM method appeared to be good in general, the abundances in a pure material region were inconsistent.The result of these abundance maps also verified that the GMM-SS-LRR algorithm with SS and low-rank property could better capture the spatial data structure and enhance the performance in the abundance estimation.The errors of abundances for these algorithms are shown in Table 3, which demonstrates that GMM-SS-LRR performed best overall.
Figure 8 shows the GMM-SS-LRR result in the wavelength-reflectance space for the Gulfport.Figure 9 shows histograms of pure pixels for the five materials (when projected to one-dimensional space determined by performing PCA on the pure pixels of each material) of the Gulfport dataset.Although there were no multiple peaks in any of the histograms, NCM algorithm still did not fit the histograms of the shadow.In contrast, our method and GMM algorithm gave a much better fit for this pure pixel distribution.For the performance of all five materials, our method matched the ground truth best, and this was also verified in the quantitative analysis presented in Table 4.

Salinas-A Datasets
The dataset was collected by the AVIRIS sensor over Salinas Valley, California, and is characterized by high spatial resolution (3.7-m pixels).It is a 512 by 217 image with 224 bands.The original image includes vegetables, bare soils, and vineyard fields.Considering that the whole dataset contains many different objects, we only performed experiments on the exemplar ROI.It is a small sub scene of Salinas image, denoted Salinas-A, which is commonly used.It comprises 86 by 83 pixels and includes six classes: broccoli, corn, lettuce 4wk, lettuce 5wk, lettuce 6wk and lettuce 7wk.The RGB image, ground truth and the superpixel map are shown in Figure 10.
Figure 11 shows the abundance maps comparison in the Salinas-A dataset.Comparing them with the ground truth (the first row of Figure 11), we can see that BCM and NCM both failed to estimate the pure pixels of corn, and, as for the lettuce 4wk, the abundance maps of GMM, NCM and BCM were mixed with the endmember corn.GMM-SS-LRR matched the ground truth best followed by GMM.The errors of abundances for these algorithms are shown in Table 5, which also implies that GMM-SS-LRR performed best overall.Figure 12 shows GMM-SS-LRR result in the wavelength-reflectance space.Figure 13 shows histograms of pure pixels for the six materials and the estimated distributions of GMM-SS-LRR, GMM and NCM.When the distribution of pure pixels was a single peak, the NCM algorithm still did not closely approximate the ground truth.In the lettuce 5wk, the GMM algorithm failed to estimate the distribution since the pure pixels had multiple peaks.This also showed that our method helped the GMM better separate the clusters and improve the performance in the estimated distribution.The quantitative analysis is shown in Table 6.

Effects of the Size of the Superpixels
For the proposed method GMM-SS-LRR, the number of superpixels for both synthetic dataset and real HSIs was set to S = 5.To compare the effects of different sizes of superpixels on the unmixing accuracy, we experimented with the synthetic dataset by changing the size of the parameters S. From the first row of Figure 14b, we can see that, when S = 1, the GMM-SS-LRR abundance maps were similar to those of GMM (the second row of Figure 14a), which means that, when we did not take the SS, the pixels were not in the homogeneous region, thus only taking low-rank representation for prior constraint could not improve the performance of the abundance estimation.As the number S increased, the performance of the abundance estimation was improved.Compared with the ground truth (the first row of Figure 14a), we can see that, when the number of superpixels S was set to 5 and 7, the precision of unmixing became the best.That means only segmenting the pixel into a highly spatially related region, and then taking low-rank representation could greatly improve the abundance estimation.However, when the number of superpixels S was set to 10 and 15, the accuracy of the unmixing declined.This was because, when the superpixel was too small, there were not enough data for training, and those superpixels had no statistical significance.The quantitative analysis of different size of superpixels is shown in Table 7.

Figure 1 .
Figure 1.Algorithm flow of the proposed method.

Figure 2 .
Figure 2. (a) The synthetic color images by extracting the bands corresponding to wavelengths 488 nm, 556 nm, and 693 nm; (b) the original spectra from the ASTER library; (c) the first principle component of the synthetic dataset; and (d) the superpixel map used for the synthetic dataset.

Figure 3 .
Figure 3.Estimated GMM-SS-LRR in the wavelength-reflectance space for the synthetic dataset.The background gray image represents the histogram created by placing the pure pixel spectra into the reflectance bins at each wavelength.The different colors represent different components, where the solid curve is the center µ jk , the dashed curves are µ jk ± 2σ jk v jk (σ jk is the square root of the large eigenvalue of Σ jk while v jk is the corresponding eigenvector), and the legend shows the prior probabilities.

Figure 4 .
Figure 4. Abundance maps for the synthetic dataset.From top to bottom: Ground truth, GMM-SS-LRR, GMM, NCM, and BCM.The corresponding endmembers from left to right are: limestone, basalt, concrete, conifer and asphalt.

Figure 5 .
Figure 5. Histograms of pure pixels for the five materials (when projected to one-dimensional space determined by performing PCA on the pure pixels of each material) of the synthetic dataset.The probability of each distribution was calculated by multiplying the value of the density function at each bin location with the bin size.

Figure 6 .
Figure 6.(a) The original RGB image of the Mississippi Gulfport dataset; (b) the original ground truth map; (c) the selected ROI 58 by 65 pixels; (d) ground truth materials of asphalt, grass, shadow, tree and grey roof in the ROI; (e) the first principle component of the Gulfport dataset; (f) the superpixel map used for the Gulfport dataset; and (g) the mean spectra of the five materials.

Figure 7 .
Figure 7. Abundance maps for the Gulfport dataset.From top to bottom: Ground truth, GMM-SS-LRR, GMM, NCM, and BCM.The corresponding endmembers from left to right are: asphalt, shadow, roof, grass and tree.

Figure 8 .
Figure 8.Estimated GMM-SS-LRR in the wavelength-reflectance space for the Gulfport dataset.The background gray image and the curves have the same meaning as in Figure 3.

Figure 9 .
Figure 9. Histograms of pure pixels for the Gulfport dataset and the estimated distribution from GMM-SS-LRR, GMM and NCM when projected to one dimension.

Figure 10 .
Figure 10.(a) Original RGB image of the Salinas-A dataset; (b) ground truth materials of broccoli, corn, lettuce 4wk, lettuce 5wk, lettuce 6wk and lettuce 7wk; (c) the first principle component of the Salinas-A dataset; (d) the superpixel map used for the Salinas-A dataset: and (e) the mean spectra of the six materials.

Figure 12 .
Figure 12.Estimated GMM-SS-LRR in the wavelength-reflectance space for the Salinas-A dataset.The background gray image and the curves have the same meaning as in Figure 3.

Figure 13 .
Figure 13.Histograms of pure pixels for the Salinas-A dataset and the estimated distribution from GMM-SS-LRR, GMM and NCM when projected to one dimension.

Table 1 .
Abundance errors for synthetic dataset.

Table 2 .
The RMSE calculated between the probability value in each histogram and the estimated value at each bin location for the synthetic dataset.

Table 3 .
Abundance errors for Gulfport dataset.

Table 4 .
The RMSE calculated between the probability value in each histogram and the estimated value at each bin location for the Gulfport dataset.

Table 6 .
The RMSE calculated between the probability value in each histogram and the estimated value at each bin location for the Sainas-A dataset.