Distributed Unmixing of Hyperspectral Data With Sparsity Constraint

Spectral unmixing (SU) is a data processing problem in hyperspectral remote sensing. The significant challenge in the SU problem is how to identify endmembers and their weights, accurately. For estimation of signature and fractional abundance matrices in a blind problem, nonnegative matrix factorization (NMF) and its developments are used widely in the SU problem. One of the constraints which was added to NMF is sparsity constraint that was regularized by L 1/2 norm. In this paper, a new algorithm based on distributed optimization has been used for spectral unmixing. In the proposed algorithm, a network including single-node clusters has been employed. Each pixel in hyperspectral images considered as a node in this network. The distributed unmixing with sparsity constraint has been optimized with diffusion LMS strategy, and then the update equations for fractional abundance and signature matrices are obtained. Simulation results based on defined performance metrics, illustrate advantage of the proposed algorithm in spectral unmixing of hyperspectral data compared with other methods. The results show that the AAD and SAD of the proposed approach are improved respectively about 6 and 27 percent toward distributed unmixing in SNR=25dB.


INTRODUCTION
Hyperspectral remote sensors capture the electromagnetic energy emitted from materials and collect hyperspectral images as a data cube with two-dimensional spatial information over many contiguous bands of high spectral resolution. One of the challenges in hyperspectral imaging is mixed pixels, which are pixels containing more than one kind of materials. So, spectrum of a pixel is a mixture of spectrum of some materials in the scene, named endmembers. Each endmember in a pixel is weighted by its fractional abundance (Miao and Qi, 2007). Decomposition of the mixed pixels is known as spectral unmixing (SU) problem (Mei et al., 2011). Most of the spectral unmixing methods are based on linear mixing model (LMM), in which it is assumed that the recorded spectrum of a particular pixel is linearly mixed by endmembers, which exist in the pixel. If the number of endmembers that are present in the scene and its signatures, are unknown, the SU problem becomes a blind source separation (BSS) problem (Qian et al., 2011).
Several SU methods were proposed in different models. Pixel purity index (PPI) (Chang and Plaza, 2006), N-FINDR (Winter, 1999), simplex volume maximization (Chan et al., 2011), convex cone analysis (CCA) (Ifarraguerri and Chang, 1999), successive projections algorithm (SPA) (Araújo et al., 2001), principal component analysis (PCA) (Smith et al., 1985), vertex component analysis (VCA) (Nascimento and Dias, 2005), (Lopez et al., 2012), are some of convex geometric methods. They are based on the pure pixel assumption, that the simplex volume is considered as a criterion for detection of endmembers. Some of the methods such as Independent Component Analysis (ICA) (Bayliss et al., 1998) use statistical models to solve the SU problem.
Another method is nonnegative matrix factorization (NMF) (Paatero and Tapper, 1994), (Lee and Seung, 1999), which decomposes the data into two nonnegative matrices. Recently, this basic method was developed with adding constraints, such as the minimum volume constrained NMF (MVC-NMF) method (Miao and Qi, 2007), graph regularized NMF (GNMF) (Rajabi and Ghassemian, 2013), NMF with local smoothness constraint (NMF-LSC) (Yang et al., 2015), multilayer NMF (MLNMF) (Rajabi and Ghassemian, 2015), structured discriminative NMF (SDNMF) (Li et al., 2016), and regionbased structure preserving NMF (R-NMF) (Tong et al., 2017). One of the constraints that has been used to improve performance of NMF methods is sparsity constraint that can be applied to NMF cost function using Lq regularizers. The sparsity constraint means that most of the pixels composed of only a few of the endmembers in the scene (Iordache et al., 2010), and the fractional abundances of other endmembers are equal to zero. So, the abundance matrix has many zero elements, and it has a large degree of sparsity. Regularization methods have been used to provide updating equations for signatures and abundances. Using L1/2 regularization into NMF, which leads to an algorithm named L1/2-NMF, has been proposed in (Qian et al., 2011), that enforces the sparsity of fractional abundances.
As another approach, the distributed strategy has been used for utilization of neighborhood information. There are some distributed strategies such as consensus strategies (Tsitsiklis and Athans, 1984), incremental strategies (Bertsekas, 1997) and diffusion strategies . In this article, a diffusion strategy is used because it has high stability over adaptive networks . Diffusion LMS strategy has been proposed in (Cattivelli and Sayed, 2010).
To solve a distributed problem, a network is considered. There are three types of networks: 1) a single-task network, that nodes estimate a common unknown and optimum vector, 2) a multitask network, which each node estimate its own optimum vector and 3) a clustered multitask network includes clusters that each of them has to estimate a common optimum vector (Chen et al., 2014). Unmixing problem is a multitask problem where each pixel considered to be as a node. Here, we first consider the general case, where there is a clustered multitask network and each cluster has an optimum vector (fractional abundance vector) that should be estimated. Then we will reduce that to a multitask case.
In this paper, the distributed unmixing of hyperspectral images is considered in which neighborhood information and sparsity constraint are used. This paper is organized as follows. Section 2 presents the proposed method and optimization procedure. Section 3 includes introduction of datasets. Section 4 provides simulation results and the last section gives conclusions.

DISTRIBUTED UNMIXING OF HYPERSPECTRAL DATA WITH SPARSITY CONSTRAINT
In this section, a new method that utilizes sparsity constraint and neighborhood information is proposed. First, we will express linear mixing model in subsection 2.1, then we will formulate the distributed cost functions with sparsity constraint in 2.2, and finally, we will use them to solve SU problem in 2.3.

Linear Mixing Model (LMM)
To solve the SU problem, we focus on a simple but efficient model, named LMM. In this model, there exists a linear relation between the endmembers that weighted by their fractional abundances, in the scene. Mathematically, this model is described as: where y is an observed data vector, A is the signature matrix, s is the fractional abundance vector and ε is assumed as a noise vector.
In the SU problem, fractional abundance vectors have two constraints in each pixel, abundance sum to one constraint (ASC) and abundance nonnegativity constraint (ANC) (Ma et al., 2014), which are as follows, for p endmembers in a scene.
Where sk(n) is the fractional abundance of the n-th endmember in the k-th pixel of the image. Some methods have been proposed without applying ASC constraint. For instance, this constraint has not been used in (Zhang et al., 2012) for complicated ground scene with nonlinear interferences. However, according to (Heinz, 2001), that studied utilization or removing of this constraint, in this article ASC constraint is adopted to gain better results.

Distributed Cost Functions and Optimization
As explained earlier, three types of networks containing single task, multitask and clustered multitask networks are supposed. First, N nodes are considered in a clustered multitask network as pixels in a hyperspectral image, and a p×1 optimum vector sC(k) is estimated, and ak is the signature vector at node k. A cost function Jk(sC(k)), when node k is in cluster C(k), defined as follows from the mean square error criterion (Chen et al., 2014): Then, the following equation is written, using the iterative steepest-descent solution (Sayed, 2003): where µ>0 is a step-size parameter, then by computing complex gradient and substituting it into (5), the following iterative equation is obtained: This equation requires knowledge of the autocorrelations, and they are replaced by instantaneous approximations, in the LMS algorithm as follows: Then, the recursive equation is changed to: Equation (9) is not distributed, because it requires knowledge of {yk , ak,} from the entire network (Cattivelli and Sayed, 2010).
In a distributed network, relationships between nodes are used for improving the accuracy. In this article, we utilize the squared Euclidean distance (Chen et al., 2014): Then, the L1/2 regularizer for sparsity constraint is used (Qian et al., 2011): Combining (4), (10) and (11), the following cost function is obtained: where it is the cost function for abundances of Q clusters, and the symbol \ is the set difference, η>0 is a regularization parameter (Chen et al., 2014) that controls the effect of neighborhood term, λ is a scalar that weights the sparsity function (Qian et al., 2011), Nk shows nodes that are in the neighborhood of node k, and the nonnegative coefficients ρkƖ are normalized spectral similarity which is obtained from correlation of data vectors (Chen et al., 2014). The coefficients are computed as introduced in ( where Nk -includes neighbors of node k except itself, and θ is computed as (Chen et al., 2014): If sk o is considered as the minimizer of the cost function, this equation is denoted as: Rayleigh-Ritz characterization eigenvalues (Sayed, 2013), is a strategy that let us to simplify the second term of equation (16) to: where bƖk are some nonnegative scalar (Chen et al., 2014). Then the cost function changes as follows for one cluster: As explained earlier, the SU is a multitask problem that each cluster only has one node. Thus, above equation with adoption of LMS strategy, is simplified to: Therefore, the optimum vectors are computed with the recursive equation, using some initial values.

Spectral Unmixing Updating Equations
According to the NMF algorithm, the conventional least squares error should be minimized with respect to the signatures and abundances matrices, subject to the non-negativity constraint (Lee and Seung, 2001). So, we have: where A and S are signatures and fractional abundances matrices, respectively, and Y denotes Hyperspectral data matrix. Then, similar to the distributed unmixing with sparsity constraint, these terms are added to the NMF cost function as follows: This cost function has been minimized with respect to A, using multiplicative update rules (Lee and Seung, 2001), and then recursive equation for signatures matrix is obtained as: Then, according to equation (22), and using d||x||1/dx=sign(x), recursive equation for abundance vectors is obtained as follows: In this paper, the ASC and ANC constraints are adopted for abundance vectors, with using the operator that explained in (Chen and Ye, 2011). This operator projects vectors onto a simplex which size of its sides are equal to one.

DATASETS
The proposed algorithm is tested on synthetic and real data. This section introduces a real dataset that recorded with hyperspectral sensors and a synthetic data set that are generated using spectral libraries.

Synthetic Images
To generate synthetic data, some spectral signatures are chosen from a digital spectral library (USGS) (Clark et al., 2007), that include 224 spectral bands, with wavelengths from 0.38µm to 2.5µm. Size of intended images is 64×64, and one endmember has been contributed in spectral signature of each pixel, randomly. Pixels of this image are pure, so to have an image containing mixed pixels, a low pass filter is considered. It averages from abundances of endmembers in its window, so that the LMM would be confirmed. Size of the window controls degree of mixing (Miao and Qi, 2007). With smaller dimension of the window and more endmembers in the image, degree of sparsity is increased.

Real Data
The real dataset that the proposed method was applied on it, is hyperspectral data captured by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) over Cuprite, Nevada. This dataset has been used since the 1980s. AVIRIS spectrometer has 224 channels and covers wavelengths from 0.4µm to 2.5µm. Its spectral and spatial resolutions is about 10nm and 20m, respectively (Green et al., 1998). 188 bands of these 224 bands are used in the experiments. The other bands (covering bands 1, 2, 104-113, 148-167, and 221-224) have been removed which are related to water-vapor absorption or low SNR bands. Figure 1 illustrates a sample band (band #3) of this dataset.

EXPERIMENTS AND RESULTS
In this section, for quantity comparison between the proposed and the other methods, the performance metrics such as spectral angle distance (SAD) and abundance angle distance (AAD) (Miao and Qi, 2007) Figure 2. The SAD performance metric of 5 methods versus SNR, using random initialization and applied on synthetic data. SAD of the proposed algorithm is star-dashed line.   where â is the estimation of spectral signature vectors and ŝ is the estimation of abundance fraction vectors.
The proposed algorithm and some other algorithms such as VCA-FCLS, NMF, L1/2-NMF and distributed unmixing is applied on synthetic data, using a 3×3 low pass filter, with 4 endmembers, and without pure pixel. According to (Chen et al., 2014), this algorithm gain the best results with (µ,η) = (0.01,0.1). After generating data, noise is added to it with 5 different values of SNR, and then performance metrics are computed by averaging 20 Monte-Carlo runs. The simulation results are shown in figures 2, 3 and 4. Also Table 1 describe the average of running time of NMF, Lq-NMF, distributed unmixing and proposed method using MATLAB R2015b with Intel Core i5 CPU at 2.40 GHz and 4 GB memory. This table shows that one of the main advantages of sparse representation is its efficiency and improvement in running time. Afterwards, the proposed algorithm is applied on real data and simulation results are shown in figures 5 and 6. According to figures 2 and 3, the AAD and SAD measures decrease as the SNR values increase. The methods that eventuate lower level of SAD or AAD, have better performance. So the star-dashed lines illustrate that the proposed algorithm produces the best results in comparison with other algorithms.

CONCLUSION
Hyperspectral remote sensing is a distinguished research topic in data processing. The purpose of spectral unmixing is decomposition of pixels in the scene into their constituent materials. The proposed distributed unmixing method with sparsity constraint was developed that estimates signatures and their abundances, and improves fractional abundances estimation toward NMF method. Simulation results on real and synthetic dataset illustrated better performance of the proposed algorithm compared with NMF, L1/2-NMF, VCA and distributed unmixing. Furthermore, this method converges faster in comparison with the distributed method.