Keywords

1 Introduction

Hyperspectral Unmixing methods have been recently developed as the powerful techniques to characterize mixed pixels of the hyperspectral spectrum into a set of constituent spectral signatures called endmembers, and the corresponding set of fraction of these endmembers, called abundances [2, 3, 16]. In a linear spectral mixture analysis fashion, we can first identify a collection of pure constituent spectra and then represent the measured spectrum of each mixed pixel with a linear combination of endmembers weighted by their fractional abundances. Using spectral unmixing approach aims to extract endmembers from the hyperspectral image and to compute the corresponding abundances under certain constrains.

In spite of two stages-based unmixing methods that employ both the endmember extraction methods (e.g., N-FINDR [25] and vertex component analysis (VCA) [17]) and mixed pixels decomposition methods (i.e., \(\ell _2\)-norm, \(\ell _1\)-norm, and \(\ell _0\)-norm approximation e.g., [5–7]), one-step methods such as nonnegative matrix factorization (NMF) [14, 19] based approach is highly interested because of its noticeable supports. First of all, the nonnegativity constraints for both spectral signatures and their fractional abundances (called the abundance nonnegativity constraint (ANC)) due to the physical consideration are automatically included in the NMF-based methods. Furthermore, it can make decomposition matrices to be more intractable because of a part-based representation of the data, see [4, 21] and references therein. The classical NMF problem is a NP-hard optimization problem [24] beyond the very large feasible set of solution without any further constraints over that. Due to the nonconvexity of the corresponding cost function, the algorithm prone to noise corruption and computationally demanding [21]. Although there are various proposed methods based on the NMF approach for unmixing purpose such as \(\ell _{1/2}\)-norm sparsity constrained [21], substance dependence constrained [26] and manifold regularization into the sparsity constraint [15], they have their own drawbacks and the researches are still working to introduce better (sparser) constrained terms into the cost function in order to improve the current methods.

In this paper, we are initiated to achieve the fact that the fractional abundances of endmembers can be effected by imposing sparsity among the endmembers collaboratively for all hyperspectral pixels. In fact, in hyperspectral images, all the pixels share the same set of spectral signatures of materials lying into a lower dimensional subspace [12]. Although the motivation of simultaneous sparse technique [23] is used in [12] for the unmixing purpose, the spectral library was assumed to be known. Besides, the averaging is applied as an \(\ell _1\)-norm term. Inspired by this motivation, we introduce a new minimization problem that considers the collaborative \(\ell _{2,q}\)-norm term along with a sparse approximation term for the abundances when the spectral library is also unknown. To solve the acquired minimization problem, we apply the multiplicative updating rule used for the standard NMF problem [13]. We also show that the updating rules will be guarantee to reach a local minimum. The simulation results demonstrate the effectiveness of our proposed method and outperform the other state-of-the-art methods for both metrics of spectral angle distance (SAD) and abundance angle distance (AAD).

The rest of the paper is organized as follows. Problem formulation for the linear mixing model (LMM) and NMF method are briefly described in Sect. 2. Our proposed method for the spectral unmixing is presented in Sect. 3 and we evaluate it by the simulations in Sect. 4. Section 5 concludes the paper.

2 Problem Formulation

In the LMM, the received spectral signature of a pixel in any given spectral band is assumed to be a linear combination of all of the endmembers’ spectra present in the pixel at the respective spectral band. The measured reflectance values of the pixel can be provided by a L-dimensional vector \(\mathbf {y}\) and can be expressed as \(\mathbf {y = Ax + e}\), where \(\mathbf {y}\) is an \(L \times 1\) column vector, \(\mathbf {A} = [\mathbf {a}_1, \mathbf {a}_2, ..., \mathbf {a}_N] \in \mathbb {R}_+^{L\times N}\) is the mixing matrix contains N pure spectral signatures (endmembers) with \(L>>N\), \(\mathbf{x}\) is a \(N \times 1\) vector with the corresponding fractional abundances of the endmembers and \(\mathbf {e}\) is an \(L \times 1\) additive noise vector due to the errors affecting the measurements at each spectral band. The fractional abundances vector \(\mathbf {x}\) must satisfy two additional constraints because of the physical considerations as the abundance nonnegativity constraint (ANC), \(\mathbf {x}\ge 0\), and the abundance sum-to-one constraint (ASC), \(\mathbf {1}^T \mathbf {x} = 1\). Matrix representation for M pixels can be considered as

$$\begin{aligned} \mathbf {Y = AX + E}, \end{aligned}$$
(1)

where, \(\mathbf {Y}\in \mathbb {R}^{L\times M}\), \(\mathbf {X} = [\mathbf {x}_1, \mathbf {x}_2, ..., \mathbf {x}_M] \in \mathbb {R}_+^{N\times M}\) and \(\mathbf {E}\in \mathbb {R}^{L\times M}\).

2.1 NMF-Based Methods: An Overview

The goal of NMF is to find two nonnegative matrices \(\mathbf {A}\) and \(\mathbf {X}\) to approximate a matrix \(\mathbf {Y}\) with the size of \(L\times M\) as \(\mathbf {Y \approx AX}\), where \(N<\min (L,M)\). Then, the loss function for the classic NMF problem can be represented by \(\frac{1}{2}||\mathbf {Y - AX}||_2^2\). Finding the global optimal solution for the corresponding minimization problem is difficult due to the nonconvexity of the problem with respect to both \(\mathbf {A}\) and \(\mathbf {X}\). The regularization method is the natural approach to tackle this problem and to promote the constraints into the cost function. Hence, the \(\ell _2\)-norm regularizer (e.g., [20]) and the \(\ell _1\)-norm regularizer (e.g., [10, 11]) are the most common choices in which the former focuses the smooth solution rather than the sparse result and the latter one achieves the sparsity of the spectral signatures dictionary matrix and/or the fractional abundances matrix.

Thinking about the \(\ell _0\)-norm regularizer as the sparsest preference is always challenging since it is an NP-hard optimization problem and cannot be achieved in practice. Thus, some \(\ell _p (0<p<1)\) regularization methods specifically when \(p=\frac{1}{2}\) are the recent interests for unmixing purpose [21]. Accordingly, the cost function can be considered as follows:

$$\begin{aligned} c(\mathbf {A,X}) = \frac{1}{2}||\mathbf {Y - AX}||_2^2 + \alpha ||\mathbf {X}||_\frac{1}{2}, \end{aligned}$$
(2)

where \(||\mathbf {X}||_\frac{1}{2} = \sum _{i=1}^{N} \sum _{j=1}^{M} x_{ij}^{\frac{1}{2}}\) and \(x_{ij}\) is the corresponding abundance for the i-th endmember at the j-th pixel and \(\alpha >0\) is the Lagrangian parameter. More \(\ell _\frac{1}{2}\)-based methods are also studied for hyperspectral unmixing, e.g., [15, 22, 26].

The multiplicative updating rules proposed for the standard NMF is used to solve the corresponding minimization problem (2) as follows [21]:

$$\begin{aligned} \mathbf {A}\leftarrow & {} \mathbf {A.*YX}^T./\mathbf {AXX}^T \end{aligned}$$
(3)
$$\begin{aligned} \mathbf {X}\leftarrow & {} \mathbf {X.*A}^T\mathbf {Y}./(\mathbf {A}^T\mathbf {AX}+\frac{\alpha }{2}\mathbf {X}^{-\frac{1}{2}}) \end{aligned}$$
(4)

where \(.*\) and . /  are the element-wise multiplication and division, respectively.

3 Our Proposed Method

The idea of collaborative sparse technique was already used for the unmixing approach in [12] by assuming that the spectral library is known. However, this assumption is not reliable in many application as a blind source separation. On the other hand, using NMF-based methods are fairly attractive because of their obvious superiorities e.g., [15, 21, 26] as mentioned earlier.

Here, we first define the following \(\ell _{2,q}\)-norm based cost function along with the sparsity of the element-wisely fractional abundances term of \(\ell _{\frac{1}{2}}\)-norm

$$\begin{aligned} f(\mathbf {A,X}) = \frac{1}{2}||\mathbf {Y - AX}||_2^2 + \alpha ||\mathbf {X}||_\frac{1}{2} + \beta ||\mathbf {X}||_{2,q}^q, \end{aligned}$$
(5)

where

$$\begin{aligned} ||\mathbf {X}||_{2,q} = \Big (\sum _{i=1}^{N}||\mathbf {x}^i||_2^q \Big )^{\frac{1}{q}}, \end{aligned}$$
(6)

and \(||\mathbf {x}^i||_2^q = \big (\sum _{j=1}^{M} |x_{ij}|^2 \big )^{\frac{q}{2}}\) with the similar definition in [18] and \(\alpha >0\) and \(\beta >0\) are the Lagrangian regularizes. In fact, the term of \(||\mathbf {X}||_{2,q}^q\) in (5) is the q-th power of \(||\mathbf {X}||_{2,q}\) and results in the sum of \(||\mathbf {x}^i||_2^q\) over i’s. It is obvious that we have the exact averaging of the \(\ell _2\)-norm of vectors \(\{\mathbf {x}^i\}\) when \(q=1\), i.e., \(\ell _{2,1}\)-norm.

The updating rule for the spectral signatures \(\mathbf {A}\) keeps the same with (3). To find the updating rule for \(\mathbf {X}\), we first take the partial derivative of (5) with respect to \(\mathbf {X}\) which gives the following result:

$$\begin{aligned} \frac{\partial f(\mathbf {A,X})}{\partial \mathbf {A}} = \mathbf {A}^T\mathbf {AX} - \mathbf {A}^T\mathbf {Y} + \frac{\alpha }{2}\mathbf {X}^{-\frac{1}{2}} + \beta \mathbf {F}.* \mathbf {|X|}, \end{aligned}$$
(7)

where

$$\begin{aligned} \mathbf {F}= q \left( \begin{array}{cccc} (||\mathbf {x}^1||_2^2)^{\frac{q}{2}-1} &{} (||\mathbf {x}^1||_2^2)^{\frac{q}{2}-1} &{} \cdots &{} (||\mathbf {x}^1||_2^2)^{\frac{q}{2}-1} \\ (||\mathbf {x}^2||_2^2)^{\frac{q}{2}-1} &{} (||\mathbf {x}^2||_2^2)^{\frac{q}{2}-1} &{} \cdots &{} (||\mathbf {x}^2||_2^2)^{\frac{q}{2}-1} \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ (||\mathbf {x}^N||_2^2)^{\frac{q}{2}-1} &{} (||\mathbf {x}^N||_2^2)^{\frac{q}{2}-1} &{} \cdots &{} (||\mathbf {x}^N||_2^2)^{\frac{q}{2}-1} \\ \end{array} \right) .* \mathrm {sign}(\mathbf {X}) \end{aligned}$$
(8)

and \(\mathrm {sign}\) function operates element-wisely over \(\mathbf {X}\).

Following the Karush-Kuhn-Tucker (KKT) conditions and applying the transposition and division, the updating rule can be determined as follows:

$$\begin{aligned} \mathbf {X}\leftarrow & {} \mathbf {X.*A^TY./(A^TAX}+ \frac{\alpha }{2}\mathbf {X}^{-\frac{1}{2}} + \beta \mathbf {F}.*\mathbf {|X|}). \end{aligned}$$
(9)

The value of Lagrangian parameter \(\alpha \) depends on the degree of sparseness for the fractional abundances of endmembers and it can be estimated based on the proportion of \(\ell _1\)-norm of the observed signals to the corresponding \(\ell _2\)-norm as the similar way in [11]. In our simulations, we set a fraction of the value defined in [11] for our proposed method as follows:

$$\begin{aligned} \alpha = \frac{\eta }{\sqrt{L}} \sum _{k=1}^L\frac{\sqrt{M}-\frac{||\mathbf {y}_k||_1}{||\mathbf {y}_k||_2}}{\sqrt{M-1}}, \end{aligned}$$
(10)

where \(\mathbf {y}_k\) is the k-th band of the observed hyperspectral images and \(0<\eta \le 1\) is a known constant value.

To impose the ASC over the fractional abundances, we append \(\delta \mathbf {1}^T\) to \(\mathbf {Y}\) and \(\mathbf {A}\), respectively, during the iterative process as the same way used in [9] and many literatures afterwards, e.g., [15, 21, 22, 26]. Then, the matrices \(\mathbf {\bar{Y}}\) and \(\mathbf {\bar{A}}\) take the place of \(\mathbf {Y}\) and \(\mathbf {A}\) as follows:

$$\begin{aligned} \mathbf {\bar{Y}}= \begin{bmatrix} \mathbf {Y} \\ \delta \mathbf {1}_M^T \end{bmatrix} , \mathbf {\bar{A}}= \begin{bmatrix} \mathbf {A} \\ \delta \mathbf {1}_N^T \end{bmatrix} \end{aligned}$$
(11)

where \(\delta \) is a known constant value that controls the effect of ASC and \(\mathbf {1}_l\) is the \(l\times 1\) column vector with all elements equal to one.

Now, we show that our proposed method is converging.

Proposition 1

The loss function in (5) is nonincreasing under (3) and (9).

Proof

The convergence of the method under the updating rule for \(\mathbf {A}\) in (3) can be shown similar to [13]. For the updating rule of \(\mathbf {X}\) in (9), we can follow up the same procedure in [21] and shows the cost function in (5) is nonincreasing function under (9).

Our proposed \(\ell _{2,q}\)-NMF based unmixing method is summarized as the following algorithm.

figure a

4 Experimental Results

In this section, we evaluate our proposed method through different experiments. First, we use the USGS library [1] to generate the synthetic data as follows.

We select \(n_N\) random spectral signatures from [1] for all of the following experiments. It should be noted that these signatures must be linearly independent. Then, we consider \(M = r^2\times r^2\) pixels of entire image to produce linear mixtures where r is an integer known value. In fact, these pixels are divided into \(r\times r\) patches. Each patch is assigned randomly by an integer value, say \(n_E\) that is between 2 (to avoid to generate the pure pixels) and \(m_E\) shows the maximum number of endmember for pixels involved in the patch. Moreover, the fractional abundances for these spectral signatures are generated randomly with the Dirichlet distribution [8] based on the assigned number of materials for each patch, i.e., \(n_E\). In order to make sure that the number of spectral signatures constructs a mixed pixel has enough contribution to build such pixel, we replace the abundances of all pixels whose their fractional abundances are larger than a threshold, say \(n_T\), with the equal contributions, i.e., \(\frac{1}{n_E}\).

Afterwards, we passed the generated pixels through the Additive White Gaussian Noise (AWGN) and observed the outputs. For the evaluation purpose, we use the SAD metric to measure the similarity between the recovered spectral signatures and the ground-truth samples. Besides, we use the AAD metric to find the similarity of the estimated fractional abundances with their ground-truth values. They are defined as follows:

$$\begin{aligned} \mathrm {SAD}_i&= \arccos \Big (\frac{\hat{\mathbf {a}}_i^T\mathbf {a}_i}{||\hat{\mathbf {a}}_i||||\mathbf {a}_i||}\Big ),\end{aligned}$$
(12)
$$\begin{aligned} \mathrm {AAD}_i&= \arccos \Big (\frac{\hat{\mathbf {x}}_i^T\mathbf {x}_i}{||\hat{\mathbf {x}}_i||||\mathbf {x}_i||}\Big ), \end{aligned}$$
(13)

where \(\hat{\mathbf {a}}_i\) and \(\hat{\mathbf {x}}_i\) represent the estimated spectral signature of i-th material and the corresponding estimated fractional abundances.

Implementation Setting: In our simulations, we set \(q=0.01\) in (5). Also, we set \(\eta =0.5\) and \(\beta = 0.2\alpha \) in (10) due to the performance consideration and set \(\eta =1\) for \(\ell _{\frac{1}{2}}\)-NMF method as recommended in [21]. Selecting the larger value of \(\delta \) gives the closer the columns of \(\mathbf {X}\) to the full additivity constraint e.g., [21, 26] which leads to more time for unmixing process in simulations. We choose \(\delta =5\) for our experiments over all unmixing methods. Moreover, we set \(r=8\), \(n_N=12\), \(m_E=5\) and \(n_T=0.7\) to generate \(M=4096\) pixels from 12 spectral signatures chosen from [1] with \(L=224\) spectral bands. Thus, each pixel has at most 5 mixed materials where the maximum purity is set to 0.7. Finally, we set \(I_{max}=1000\) and \(\varepsilon =10^{-4}\) in Algorithm 1 as well as \(\ell _{\frac{1}{2}}\)-NMF method proposed in [21].

Fig. 1.
figure 1

Comparison of the estimated spectral signatures of four sample materials with the corresponding ground-truth values (SNR = 25 dB) (a) Almandine (b) Bronzite (c) Biotite (d) Chlorite

Full size image
Fig. 2.
figure 2

Simulation results of unmixing methods of (a) SAD (b) AAD as functions of SNRs.

Full size image

We did several and different experiments and we only report some of the simulation results due to the space limit. First, the unmixing results for a fixed value of signal-to-noise-ratio (SNR) is given. Hence, Fig. 1 compares the estimated spectral signatures by our proposed method and the other state-of-the-art methods for 4 sample materials out of 12 while the SNR is set to 25 dB. We can observe that our proposed method gives the closest results to the original spectral signatures of materials compared with the other two unmixing methods.

We also evaluate the performance of our proposed method for different values of the noise power. Figure 2 shows the values of SAD and AAD as the functions of SNRs for our proposed method as well as the other two state-of-the-art- methods. The value of inf represents noise-free observation. Our proposed method outperforms the other methods form lower SNRs to the noise-free environments.

5 Conclusion

In this paper, we proposed a method of unmixing hyperspectral images based on the collaborative property of the fractional abundances of endmembers through the NMF problem. We introduced a new cost function based on the NMF and \(\ell _{2,q}\)-norm to formulate our minimization problem. Then, we applied the multiplicative updating rules to solve the desire objective function. We showed that our proposed method converges and evaluated it by different experiments over the USGS spectral library. Our simulation results illustrated that the proposed method outperformed the other state-of-the-art methods in terms of SAD and AAD metrics. Although the experimental results obtained for the synthetic data analysis are promising, more experiments with real hyperspectral data are also interesting to fully support our contributions.