Hyperspectral image unmixing algorithm based on endmember-constrained nonnegative matrix factorization

Abstract

The objective function of classical nonnegative matrix factorization (NMF) is non-convexity, which affects the obtaining of optimal solutions. In this paper, we proposed a NMF algorithm, and this algorithm was based on the constraint of endmember spectral correlation minimization and endmember spectral difference maximization. The size of endmember spectral overall-correlation was measured by the correlation function, and correlation function was defined as the sum of the absolute values of every two correlation coefficient between the spectra. In the difference constraint of the endmember spectra, the mutation of matrix trace was slowed down by introducing the natural logarithm function. Combining the image decomposition error with the influences of endmember spectra, in the objective function the projection gradient was used to achieve NMF. The effectiveness of algorithm was verified by the simulated hyperspectral images and real hyperspectral images.

Key words： hyperspectral image; unmixing; nonnegative matrix factorization (NMF); correlation; logarithm function

Introduction

Owing to the limitation of spatial resolution and the complexity of real situation, some of the pixels in the hyperspectral images are composed of a variety of substances (endmember), and these pixels are called mixed pixels [ 1– 4]. To improve the accuracy in describing the real situation, the mixed pixels were decomposed, and the ratio (abundance) of a kind of ground object type (endmember) in the pixels were calculated. The topic of unmixing is important in the quantitative research of hyperspectral image. The main algorithms of endmember extraction include pixel purity index (PPI), N-FINDR, vertex component analysis (VCA) which are based on convex geometry, automated morphological endmember extraction (AMEE) which is based on mathematics morhology, and independent component correlation algorithm (ICA), by which endmember extraction and abundance inversion may be performed simultaneously.

Nonnegative matrix factorization (NMF) algorithm has powerful information processing and problem solving capabilities, and the nonnegative constraints are required to solve many practical problems. NMF algorithm is widely used in the engineering field [ 5]. The objective function of NMF algorithm has obvious non-convexity, possessing local minimum values. To make the NMF algorithm to be applicable in different fields, the corresponding constraint conditions are required to be more complex according to the physical characters of different applications. At present, the conditions of smoothness constraint (SC), minimum volume constraint (MVC) and sparsity constraint, especially the sparsity constraint that based on clustering regularization, have already been proposed [ 6– 9]. In this paper, in order to improve the unmixing precision, combining with NMF, an unmixing algorithm, which may be called endmember constraint nonnegative matrix factorization (EC-NMF), was presented. In this algorithm, the correlations and differences between endmembers were used as constraint conditions. The validity of algorithm was verified by the simulated hyperspectral images and real hyperspectral images.

Linear spectral mixing model

The unmixing models can be classified into two categories, i.e., linear spectral mixture model (LSMM) and nonlinear spectral mixture model (NLSMM). LSMM has been widely studied and applied. This model has several advantages, such as, simple principle and calculation, high efficiency and clear physical meaning. In the linear model, the mixed spectra are the linear combination of endmember spectra and its ratio. Assumed that X is pixel spectral vector, S=[s_{\mathrm{1}},s_{\mathrm{2}},\cdots ,s_N] is endmember spectra matrix, N is stochastic noise, A={[a_{\mathrm{1}},a_{\mathrm{2}},\cdots ,a_N]}^T is N-dimensional vector, and each of the components of A corresponds to an endmember abundance. LSMM can be denoted as

(1) X=SA+N,

(2) a_i\ge \mathrm{0},

(3) {\displaystyle \sum _{i=\mathrm{1}}^N{a}_i}=\mathrm{1.}

Nonnegative matrix factorization based on endmember constraints

Nonnegative matrix factorization

In the NMF algorithm, based on the objective function of minimizing Euclidean distance, the optimal solutions of S and A were obtained by taking into account of X,

(4) f(S,A)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\Vert X-SA\Vert }^{\mathrm{2}}.

The iterative formula can be given as follows,

(5) S\leftarrow S-{\beta }_{\mathrm{1}}(SA-X)A^T,

(6) A\leftarrow A-{\beta }_{\mathrm{2}}{S}^T(SA-X),

where {\beta }_{\mathrm{1}} and {\beta }_{\mathrm{2}} are weights.

Spectral unmixing algorithm is based on NMF. In this algorithm, it is not necessary to determine whether there is pure pixel, the corresponding abundance may be obtained when the endmember is extracted.

Endmember spectral correlation constraint

Owing to abundance sum-to-one constraint (ASC) in Eq. (3), the independence between endmembers is undermined [ 10]. The normalized correlation coefficient function of two spectral vectors may be expressed as

(7) \rho （{\widehat{S}}_i,{\widehat{S}}_j）={\displaystyle \sum _{k=\mathrm{1}}^L{\widehat{S}}_{ki}}{\widehat{S}}_{kj}.

{\widehat{S}}_i and {\widehat{S}}_j are normalized S_i and S_j ,respectively. It can be seen in Eq. (7), when the two vectors have the same variation trend, the correlation coefficient function is positive; on the contrary, the correlation coefficient function is negative. Whereas when the two vectors are not related, the correlation coefficient function is zero.

To measure the overall correlation of N spectra, the spectral correlation function may be defined as

(8) \mu (\widehat{S})={\displaystyle \sum _{i=\mathrm{1}}^N{\displaystyle \sum _{j=i+\mathrm{1}}^N\left|\rho （{\widehat{S}}_i,{\widehat{S}}_j）\right|}},

where |\cdot | is the absolute value function. Equation (8) denotes the sum of absolute values of correlation coefficient function of every two spectra. When the correlation between every two spectra decreases, the value of \left|\rho （{\widehat{S}}_i,{\widehat{S}}_j）\right| decreases, tending to zero. Consequently, the value of \mu (\widehat{S}) also decreases, and vice versa.

As the constraint condition, the spectral correlation function is adopted.

(9) J_{\mathrm{1}}(\widehat{S})={\displaystyle \sum _{i=\mathrm{1}}^N{\displaystyle \sum _{j=i+\mathrm{1}}^N\left|\rho （{\widehat{S}}_i,{\widehat{S}}_j）\right|}}.

The constraint should be minimized, its derivation can be written to be

(10) {\nabla }_{\widehat{s}}{J}_{\mathrm{1}}(\widehat{S})={\displaystyle \sum _{i=\mathrm{1}}^N{\displaystyle \sum _{j=i+\mathrm{1}}^N\sin \mathrm{g}(\rho （{\widehat{S}}_i,{\widehat{S}}_j）){\widehat{S}}_j^{\mathrm{T}}{\widehat{S}}_i}},

where sign(\cdot ) is sign function.

Difference constraint of endmember spectra

The endmember spectra difference constraint is needed to be introduced, or the ideal unmixing accuracy cannot be obtained by taking only Eq. (9) as constraint condition of NMF. Endmember spectra should have greater independence, and the differences between them should be the largest [ 11]. The differences of endmember spectra can be expressed to be the reciprocal of the trace of endmember spectra autocorrelation matrix.

(11) J_{\mathrm{2}}(S)={(\mathrm{Tr}(S^{T}S))}^{-\mathrm{1}},

where \mathrm{Tr}(\cdot ) is the matrix trace. The lower the value of J_{\mathrm{2}}(S) , the greater the differences between endmember spectra. Since S^{T}S is a symmetric positive definite matrix, all the eigenvalues of S^{T}S should not be negative, i.e., {\lambda }_i\ge \mathrm{0} , {\displaystyle \sum {\lambda }_i>\mathrm{0}} , \mathrm{Tr}(S^{T}S)={\displaystyle \sum {\lambda }_i} , as a result, J_{\mathrm{2}}(S)>\mathrm{0} . In view of the large impact of the matrix on the value of Eq. (11), to stabilize the algorithm, the optimization function item should be used, which is given to be

(12) J_{\mathrm{2}}(S)=\ln ({(\mathrm{Tr}(S^{T}S))}^{-\mathrm{1}})=-\ln (\mathrm{Tr}(S^{T}S)).

In Eq. (12), natural logarithm function is introduced to slow the mutation of trace. The stability of the algorithm is improved without changing the monotony of trace.

The constraint should be minimized, its derivation can be written to be

(13) {\nabla }_S{J}_{\mathrm{2}}(S)=-{\displaystyle \frac{S}{\mathrm{Tr}(S^{T}S)}}.

Endmember constraint nonnegative matrix factorization algorithm description

The objective function of EC-NMF is given to be

(14) J(S,A)=f(S,A)-{\alpha }_{\mathrm{1}}{J}_{\mathrm{1}}(\widehat{S})-{\alpha }_{\mathrm{2}}{J}_{\mathrm{2}}(S),

where {\alpha }_{\mathrm{1}} and {\alpha }_{\mathrm{2}} are the weights of J_{\mathrm{1}}(\widehat{S}) and J_{\mathrm{2}}(S) , respectively. The relations between f(S,A) and J_{\mathrm{1}}(\widehat{S}) or J_{\mathrm{2}}(S) may be adjusted via {\alpha }_{\mathrm{1}} and {\alpha }_{\mathrm{2}} .

The number of endmembers may be evaluated with virtual dimensionality (VD) methods [ 12].

Based on NMF, Eq. (2) was met. To meet Eq. (3), the method presented in Ref. [ 13] was adopted, which is given to be

(15) \bar{X}=\left[\begin{array}{l}X\\ \delta {\mathrm{1}}_N^T\end{array}\right],\bar{S}=\left[\begin{array}{l}S\\ \delta {\mathrm{1}}_N^T\end{array}\right],

where \mathrm{1} is a vector which all the components are 1, i.e., {\mathrm{1}}_N^T={[\mathrm{1},\mathrm{1},\cdots ,\mathrm{1}]}^T\in U^N . \delta is weight. When iterated, X is replaced by \bar{X} , and S is replaced by \bar{S} .

The projection gradient iterative formula of EC-NMF should be

{\overline{S}}^{k+1}\leftarrow P_{\Omega }(x)({\overline{S}}^k-{\beta }_{1k}(({\overline{S}}^k{\mathrm{A}}^k-{\overline{\mathrm{X}}}^k){\left({\mathrm{A}}^k\right)}^{\mathrm{T}}-{\alpha }_{1k}{\nabla }_{\overline{s}}{J}_1\left({\overline{S}}^k\right)-{\alpha }_{2k}{\nabla }_{\overline{s}}{J}_2\left({\overline{S}}^k\right))),

(16) A^{k+\mathrm{1}}\leftarrow P_{\Omega }(x)(A^k-{\beta }_{\mathrm{2}k}{\left({\bar{S}}^k\right)}^T({\bar{S}}^k{A}^k-{\bar{X}}^k\left)\right).

The projection function can be selected as [ 14]

(17) P_{\Omega }(x)=\max (\mathrm{0},x).

Simulation experiment

Performance indexes

The unmixing effects are measured by the two common indicators, i.e., spectral angle distance (SAD) and root mean square error (RMSE). They are often used to calculate the approximation degree of the estimated values of spectra and abundance unmixing. For the ith endmember, SAD is defined as

(18) SA{D}_i={\cos }^{-\mathrm{1}}\left({\displaystyle \frac{S_i^T{S^{\prime }}_i}{\Vert S_i^T\Vert \Vert {S^{\prime }}_i\Vert }}\right),

where S is the true endmember spectrum, and S^{\prime } is the estimated value of S .

RMSE is defined as

(19) RMS{E}_i={\left[{\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \sum _{j=\mathrm{1}}^N{(A_{ij}-{A^{\prime }}_{ij})}^{\mathrm{2}}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}},

where A_i is the real abundance of endmember, and {A^{\prime }}_i is estimated value of A_i .

Simulation data experiment

Five linearly independent endmember spectra (alunite, buddingtonite, calcite, kaolinite, muscovite) were selected from United States Geological Survey (USGS) mineral spectral library. The five spectra were mixed according to the Dirichlet distribution, the sum of endmember abundances were normalized, and different white noise were added to generate the simulation experimental data.

Anti-noise performance experiment: as shown in Fig. 1, the unmixing performances of the EC-NMF algorithm was compared with minimum volume constraint nonnegative matrix factorization (MVCNMF), smoothness constraint nonnegative matrix factorization (SCNMF) and constrained sparse nonnegative matrix factorization (CSNMF) in different noise conditions [6,7,9]. The experiments were conducted when the signal to noise ratios (SNR) were set to be ∞ (without noise), 30, 20 and 10 dB, respectively. \overline{SAD} and \overline{RMSE} denoted the mathematical expectation.

Fig.1 Algorithm performance comparisons in different noise intensities. (a) \overline{SAD} ; (b) \overline{RMSE}

Full size|PPT slide

Unmixing performance experiments for different pixel numbers: as shown in Fig. 2, when the numbers of pixels were different, the unmixing performance of EC-NMF algorithm was compared with those of MVCNMF, SCNMF and CSNMF. Pixel numbers were set to be 1600, 3600, 6400, and 10000; the SNR was set to 30 dB.

Fig.2 Algorithm performance comparisons for different pixel numbers. (a) \overline{SAD} ; (b) \overline{RMSE}

Full size|PPT slide

As can be seen from experimental results, the unmixing performance of EC-NMF algorithm is better than those of MVCNMF, SCNMF and CSNMF.

Real data experiment

Hyperspectral images of cuprite region in Nevada were collected by Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) in American JPL. Figure 3 is the 172th band image. When the low SNR and water vapor absorption band were removed, the remaining 49 valid bands were obtained. The real object distributions in this region were given by Ref. [ 15].

Fig.3 The 172th band image of Cuprite region

Full size|PPT slide

The unmixing results of EC-NMF algorithm were presented in Fig. 4, and the \overline{SAD} results of four algorithms were listed in Table 1. It can be seen that the EC-NMF performance was better than those of the other three algorithms.

Fig.4 Unmixing results in Cuprite region obtained in EC-NMF algorithm. (a) Alunite; (b) buddingtonite; (c) chalcedony; (d) jarosite; (e) kaolinite#1; (f) kaolinite#2; (g) kaolinite#3; (h) montmorillonite; (i) muscovite; (j) nontronite; (k) sphene

Full size|PPT slide

Tab.1 SAD comparisons for different algorithms in Cuprite region

endmember	EC-NMF	CSNMF	SCNMF	MVCNMF
alunite	0.0731	0.1053	0.0836	0.0772
buddingtonite	0.0966	0.0842	0.1197	0.1493
chalcedony	0.1425	0.1622	0.2921	0.1641
jarosite	0.1990	0.2013	—	—
kaolinite#1	0.2297	0.2342	0.3717	0.2562
kaolinite#2	0.3285	0.3310	0.3595	—
kaolinite#3	0.1337	0.1507	0.2546	—
montmorillonite	0.1273	0.0935	—	0.1353
muscovite	0.0819	0.1365	0.0874	0.0945
nontronite	0.0875	0.0881	—	—
sphene	0.0868	0.0954	0.3173	0.3528

Conclusions

In view of the characteristics in endmember independence, an unmixing algorithm based on endmember-constrained NMF was proposed. The experimental analysis of the simulated data and real data indicated that the proposed algorithm can improve the unmixing accuracy. The convergence rate of projected gradient algorithm was slow, and a more efficient optimized algorithm is needed to improve the unmixing efficiency. To further improve the unmixing accuracy of mixed pixels, other forms of constraint conditions should be adopted.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 61275010) and the Scientific Research Fund of Heilongjiang Provincial Education Department (No. 12541734).

References

1 Thouvenin P A, Dobigeon N, Tourneret J Y. Hyperspectral unmixing with spectral variability using a perturbed linear mixing model. IEEE Transactions on Signal Processing, 2016, 64(2): 525–538

2 Heylen R, Scheunders P. A multilinear mixing model for nonlinear spectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(1): 240–251

3 Zheng C Y, Li H, Wang Q, Chen C L P. Reweighted sparse regression for hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(1):479–488

4 Altmann Y, Pereyra M, Bioucas-Dias J. Collaborative sparse regression using spatially correlated supports—application to hyperspectral unmixing. IEEE Transactions on Image Processing, 2015, 24(12): 5800–5811

5 Guillamet D, Vitrià J, Schiele B. Introducing a weighted non-negative matrix factorization for image classification. Pattern Recognition Letters, 2003, 24(14): 2447–2454

6 Pauca V P, Piper J, Plemmons R J. Nonnegative matrix factorization for spectral data analysis. Linear Algebra and Its Applications, 2006, 416(1): 29–47

7 Miao L, Qi H. Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(3): 765–777

8 Hoyer P O. Non-negative matrix factorization with sparseness constraints. Journal of Machine Learning Research, 2004, 5(1): 1457–1469

9 Lu X, Wu H, Yuan Y. Double constrained NMF for hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(5): 2746–2758

10 Luo W F, Zhong L, Zhang B, Gao L R. Independent component analysis for spectral unmixing in hyperspectral remote sensing image. Spectroscopy and Spectral Analysis, 2010, 30(6): 1628–1633 (in Chinese)

PMID

11 Wu B, Zhao Y, Zhou X. Unmixing mixture pixels of hyperspectral imagery using endmember constrained nonnegative matrix factorization. Computer Engineering, 2008, 34(22): 229–231

12 Chang C, Du Q. Estimation of number of spectrally distinct signal sources in hyperspectral imagery. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(3): 608–619

13 Heinz D C, Chang C. Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(3): 529–545

14 Gillis N, Glineur F. Using underapproximations for sparse nonnegative matrix factorization. Pattern Recognition, 2010, 43(4): 1676–1687

15 Clark R N, Swayze G A. Evolution in imaging spectroscopy analysis and sensor signal-to-noise: an examination of how far we have come. In: Proceedings of The 6th Annual JPL Airborne Earth Science Workshop, 1996