An Endmember Initialization Scheme for Nonnegative Matrix Factorization and Its Application in Hyperspectral Unmixing

Nonnegative matrix factorization (NMF) is a blind source separation (BSS) method often used in hyperspectral unmixing. However, it tends to converge to a local optimum. To overcome this limitation, we present a simple, but effective endmember initialization scheme for NMF, which is realized by improving initial values through the application of the automatic target generation process (ATGP) algorithm. The initial spectra and abundances of target endmembers are first obtained using the ATGP algorithm and nonnegative least squares (NNLS) method, respectively. The preliminary results are then optimized through iterative application of NMF. To validate the applicability and effectiveness of the proposed method, we analyzed the improvement of NMF by the ATGP algorithm, using the synthetic hyperspectral data and real hyperspectral images. The results from the proposed method are compared with those of the vertex component analysis (VCA)-NMF algorithm, which uses the VCA algorithm to perform initialization for NMF, the minimum volume constrained NMF (MVC-NMF) algorithm, the traditional two-step VCA-fully-constrained least squares (FCLS) algorithm, which uses the VCA to extract the endmember matrix, and the FCLS algorithm to estimate the abundance matrix. The comparison results prove that proper endmember initialization can help the NMF algorithm yield better estimation results. Through the optimization of target endmembers’ initial values, the proposed ATGP-NMF algorithm can consistently produce good results at a lower computational cost, especially in the case of a real hyperspectral image for which pure pixels do not exist and there is little prior knowledge. With its high applicability and effectiveness, the ATGP-NMF algorithm has a great potential to solve hyperspectral unmixing problems.


Introduction
Hyperspectral imagery can provide abundant spatial and spectral information for observed ground objects.It has attracted increasingly more attention and has been widely used.However, due to the low spatial resolution of hyperspectral sensors, mixed pixels prevail in hyperspectral images.The mixed pixel problem severely influences the precision of object recognition and classification; therefore, it has become an obstacle to quantification analysis of hyperspectral images.Hence, spectral unmixing has always been an important subject in hyperspectral data analysis [1].Traditional spectral unmixing models are usually developed based on endmember spectra and involve two steps, i.e., endmember extraction and spectral unmixing [2][3][4][5].These methods depend heavily on the precision of endmember extraction and require certain prior knowledge [6].Their effectiveness is greatly limited when there is insufficient prior information.Therefore, it is greatly important and necessary to develop data-driven blind unmixing methods for hyperspectral images.
Over the past decade, the development of a new technique in the signal processing field, blind source separation (BSS) [7], has offered novel methods for spectral unmixing.Different from traditional two-stage methods, the BSS method can simultaneously obtain the endmembers and their corresponding abundances.The current blind unmixing methods mainly include independent component analysis (ICA) [8], nonnegative matrix factorization (NMF) [9], complexity-based BSS methods [10] and sparse component analysis (SCA) [11,12].
Different from other methods, NMF can obtain nonnegative results with physical significance.Therefore, NMF and its extensions have been widely used for spectral unmixing.However, these NMF algorithms often provide only "local best" solutions due to the non-convexity of objective functions [13].Most NMF-based algorithms deal with this problem by adding different constraints, such as the smoothness constraint, the minimum volume constraint, the sparseness constraint and the regularized constraint.Representative algorithms of the smoothness constraint have been reported in [14,15].Miao et al. developed the minimum volume constrained nonnegative matrix factorization (MVC-NMF) [16], and Li et al. introduced the sparseness constraint of abundance matrix into MVC-NMF [17].Various sparse regression-based unmixing methods have been introduced by incorporating the sparseness constraints into spectral unmixing [18][19][20][21].Representative research on the sparseness constraints include the coupled nonnegative matrix factorization (CNMF) [22], the substance dependence constrained sparse nonnegative matrix factorization (SDSNMF) [23], the total variation regularized reweighted sparse NMF (TV-RSNMF) [24] and the simultaneously sparse and low-rank constrained NMF (DSPLR-NMF) [25].Recent research has also considered regularized constraints, such as the sparsity-constrained NMF with L 1/2 regularization [26], the sparsity-regularized robust NMF (RNMF) [27], the Arctan-NMF with smooth and sparse regularization [28], the robust collaborative nonnegative matrix factorization (R-CoNMF) [29] and the graph regularized nonnegative matrix factorization (GNMF) [13].These methods are mainly proposed to get rid of the local minima of NMF by adding constraints to the endmembers or the abundances.However, most statistical-based blind unmixing methods can still be trapped in local minima if the initial values are not properly provided.
Researchers have also tried to solve the local optimization problem by providing reasonable initialization and have found that the proper initialization of NMF and its extensions can help them to avoid the "local best" problem and refine their solutions [30,31].Several methods have been proposed to improve the initialization, such as singular value decomposition (SVD) [32], clustering-based initialization approaches [33] and vertex component analysis (VCA) [3,34].Although great progress has been made in improving the local optimization problem of the NMF algorithm, these studies have mainly focused on signal processing and related areas.Further studies still need to be done to verify their applicability in solving the mixed pixel problem of a hyperspectral image.Compared to traditional supervised endmember extraction methods, such as the pixel purity index (PPI), ATGP requires less prior information to extract pure pixel vectors as the endmember spectrum from the target image [35].Moreover, most of the ATGP-generated target pixels turn out to be endmembers [36].In other words, ATGP can provide more accurate initial points to identify endmembers.When it was used to generate initial values for another BSS method, ICA, ICA showed a better performance [37].Therefore, we endeavored to use the ATGP algorithm to extract better initial endmembers for the NMF algorithm.
In this paper, we propose a new method that utilizes the ATGP algorithm to perform endmember initialization for NMF, with the purpose of solving the NMF's local optimum problem.First, the ATGP algorithm is used to find initial target endmembers from a hyperspectral image, which are used to initialize the endmember matrices of NMF.The initial abundance matrices of NMF are extracted by using NNLS.Second, the endmember spectrum matrices and abundance matrices are iteratively updated by the multiplicative update algorithm until the optimization goal is achieved.Finally, the effectiveness and efficiency of the proposed method (ATGP-NMF) are verified by comparing its experimental results in different situations with the results of other methods.
The remainder of this paper is organized as follows.Section 2 briefly describes the linear mixture model (LMM), NMF and ATGP and introduces the proposed endmember initialization scheme for NMF.Section 3 presents experimental results and evaluates the performances of the proposed method and the other three algorithms on synthetic hyperspectral data and real hyperspectral images.A summary of this study and the conclusions are presented in Section 4.

Linear Mixture Model
LMM assumes that the observed pixel is a linear combination of all the endmembers weighted by their corresponding abundance fractions [38,39].In addition, LMM has been widely used for hyperspectral unmixing due to its physical effectiveness and mathematical simplicity.The LMM can be expressed as where L×N represents the observation matrix with L bands and N pixels, and the endmember matrix and the abundance matrix, respectively.P denotes the number of endmembers and E ∈ R L×N is the additive noise.In LMM, the abundance vectors are subjected to the abundance nonnegative constraint (ANC) and abundance sum-to-one constraint (ASC).

Nonnegative Matrix Factorization
NMF is a general matrix decomposition method, which decomposes the data matrix X as a product of two nonnegative matrices A and S. Note that L denotes the number of bands, N denotes the number of pixels, P denotes the number of endmembers, and P < min (L, N).In essence, NMF is an optimization problem.An objective function based on Euclidean distance is presently the most widely used method; in this case the goal of the NMF model is to find the minimum reconstruction error by two nonnegative matrices A and S: where a ij ≥ 0, and s ij ≥ 0, 1 2 • 2 F represents the Frobenius norm.For the objective function 1  2 X − AS 2 F , many learning algorithms, such as the multiplicative update algorithm, alternating least squares algorithm and gradient descent algorithm, are presently often used to solve the optimization problem of NMF.Among these algorithms, the multiplicative update algorithm can guarantee the non-negativity of both the calculation process and the final results [40].It also adopts an auto-adjusted step size in the iteration, rather than the traditional manual setting, which reduces the impacts of parameter selection.In this paper, we used the general iterative expression formula of the multiplicative update algorithm to solve the objective function.
For each matrix of A and S, the objective function 1 2 X − AS 2 F is a convex function, but if matrices A and S are considered as a whole, the objective function is no longer convex.This makes it difficult to converge to a global optimum.Tao et al. found that the NMF algorithm can determine the optimum and physically significant solution more quickly and easily if a better starting point is provided [41].

Automatic Target Generation Process
The selection of an initial target endmember has a great impact on spectral unmixing [42].Based on unsupervised and unconstrained orthogonal subspace projection (OSP) theory, the ATGP algorithm searches for the candidate endmembers with the maximal orthogonal projection in a sequence of orthogonal projection subspaces [43].ATGP is an unsupervised target detection technique.It implements an orthogonal subspace projector to find a set of potential targets with the maximal orthogonal projections (OPs) from the image.The orthogonal subspace projector of ATGP can be defined as: where the subspace U is composed of the targets previously generated by the ATGP algorithm.
It should be noted that before the ATGP algorithm is applied, the number of spectral endmembers should be determined, especially in cases for which there is not sufficient prior information.The Harsanyi-Farrand-Chang (HFC) method based on the Neyman-Pearson detection theory has recently been widely used to determine the number of spectral endmembers [44].Therefore, in this study, the HFC algorithm is used to determine the virtual dimensionality (VD), which is representative of the number of endmembers.
Detailed steps of the ATGP algorithm are described as follows: 1.
The observed pixel of the hyperspectral image is set to a vector x, and the VD determined by the HFC algorithm is q. 2.
Based on the convex geometry theory, we choose a target pixel with the maximum vector length, which corresponds to the brightest pixel in the image, as the initial target pixel vector denoted by

3.
Based on the orthogonal subspace projection theory, the ATGP algorithm begins with the initial target pixel vector t 0 by applying an orthogonal subspace projector P ⊥ t 0 to all pixel vectors in the image and finds the first target pixel vector denoted by t 1 with the maximum orthogonal projection in the orthogonal complement space t 1 = arg max 4. Suppose that we find the i th target t i generated at the i th step by t i = arg max where is the target pixel vector matrix generated at the (i − 1) st step.

5.
The orthogonal subspace of U i is given by indicates that the projector P ⊥ U i maps the observed pixel x into the range ⊥ U i , the orthogonal complement of U i ; and I is a unit matrix.

6.
The ATGP stops when i equals to the number of the endmembers q.The target pixel vector matrix [t 1 , t 2 , • • • , t q which contains q target pixel vectors apart from the initial target pixel vector t 0 , is then regarded as the endmember spectra.
Another widely used endmember finding algorithm, i.e., the VCA algorithm, extracts endmembers by the maximal OPs from a sequence of successive orthogonal projection subspaces, which is quite similar to that of the ATGP [45,46].Actually, the VCA algorithm can be considered as a variant of ATGP.However, the VCA algorithm does not have stable performance when random initial values are used [47].Scholars have proven that the performance of VCA can be improved by using ATGP-generated targets as its initial conditions [48].Compared to VCA, the ATGP algorithm is not constrained by initial conditions or dimension reduction transformations.Therefore, it serves as a more effective and accurate endmember extraction method when little or no prior information is available.

Proposed Endmember Initialization Scheme for NMF
The NMF blind unmixing method improved by the ATGP-generated initial target endmembers (ATGP-NMF) can be expressed as follows: 1.
Initialization: target endmembers extracted by the ATGP algorithm are used as the initial endmember matrix A of NMF.Since the NNLS algorithm has been used for hyperspectral unmixing, especially for the initialization of NMF [33,49], the corresponding abundances of target endmembers derived by NNLS are used as the initial abundance matrix S of NMF.

2.
Objective function and optimization algorithm: The objective functions are constructed based on the Euclidean distance, due to its simplicity, intuition and wide application.To solve the objective functions, the multiplicative update algorithm is selected as the optimization criterion.The update formula of endmember matrix A and abundance matrix S can be expressed as Equations ( 4) and ( 5): where S pn is the p th endmember abundance value of the n th mixed pixel, and A l p is the p th endmember spectrum value of the l th band.Here, a minimal positive value λ is used to ensure that the denominator is nonzero [14].

3.
Update: the rows of the weight coefficient matrix S are updated using Equation (4), while the columns of the basic matrix A are updated using Equation (5). 4.
Normalization: according to the ASC in the spectral mixing model, and following the methods suggested in [2,9,14], the weight coefficient matrix S is normalized after every iteration using Equation ( 6) to guarantee that the sum of each column vector in the weight coefficient matrix S is always equal to one.
Final endmember matrix A and endmember abundance matrix S: matrices A and S are repeatedly updated until the maximum number of iterations is reached or until the stopping condition is satisfied, i.e., the value of the objective function is equal to zero or the prescribed error threshold ε .
where X(k) is the reconstruction of X after the k th iteration.
The endmembers extracted by the blind unmixing method do not contain information about their categories; therefore, it is necessary to compare the endmember spectra extracted from hyperspectral data to the real endmember spectra and further identify their exact categories.In this study, we search for the best match of each endmember's spectral curve in the reference spectra set until the average correlation coefficient of the spectra set reaches the maximum value.
The flowchart of the proposed method is shown in Figure 1.

Objective function and optimization algorithm:
The objective functions are constructed based on the Euclidean distance, due to its simplicity, intuition and wide application.To solve the objective functions, the multiplicative update algorithm is selected as the optimization criterion.
The update formula of endmember matrix A and abundance matrix S can be expressed as Equations ( 4) and ( 5): A lp ←A lp (XS T ) lp (ASS T ) lp +λ (5) where S pn is the p th endmember abundance value of the n th mixed pixel, and A lp is the p th endmember spectrum value of the l th band.Here, a minimal positive value λ is used to ensure that the denominator is nonzero [14].
3. Update: the rows of the weight coefficient matrix S are updated using Equation ( 4), while the columns of the basic matrix A are updated using Equation ( 5). 4. Normalization: according to the ASC in the spectral mixing model, and following the methods suggested in [2,9,14], the weight coefficient matrix S is normalized after every iteration using Equation ( 6) to guarantee that the sum of each column vector in the weight coefficient matrix S is always equal to one.
5. Final endmember matrix A and endmember abundance matrix S : matrices A and S are repeatedly updated until the maximum number of iterations is reached or until the stopping condition is satisfied, i.e., the value of the objective function is equal to zero or the prescribed error threshold ε .
where X � (k) is the reconstruction of X after the k th iteration.
The endmembers extracted by the blind unmixing method do not contain information about their categories; therefore, it is necessary to compare the endmember spectra extracted from hyperspectral data to the real endmember spectra and further identify their exact categories.In this study, we search for the best match of each endmember's spectral curve in the reference spectra set until the average correlation coefficient of the spectra set reaches the maximum value.
The flowchart of the proposed method is shown in Figure 1.

Experimental Results
In this study, we carried out spectral unmixing experiments using the synthetic hyperspectral data and two real hyperspectral images to analyze the applicability and validity of the ATGP-NMF algorithm.The synthetic hyperspectral data were used to simulate a series of situations with different SNR levels, different purity levels and different numbers of endmembers.The real hyperspectral images representsituations in which pure pixels do not exist and there are strong spectral correlations between endmembers.The performance of the proposed method was compared with those of three representative algorithms, namely the VCA-NMF algorithm, the MVC-NMF algorithm and the VCA-FCLS algorithm.The VCA-NMF algorithm is included to test whether the ATGP approach is a better choice than the VCA to perform the initialization of NMF.The MVC-NMF algorithm, which adds the minimum volume constraint into NMF, is used to compare two different solutions for the NMF's local optimum problem, i.e., to impose certain constraints and provide better initialization.In addition, it is a representative of constraint-NMF algorithms and has been widely used in comparative studies of NMF-based algorithms [23,50,51].The VCA-FCLS algorithm, which uses VCA to identify endmembers and FCLS to calculate abundances, represents the traditional two-stage spectral unmixing method.

Performance Metrics
The performances of the abovementioned algorithms are evaluated from two perspectives, the similarity to real spectra and the accuracy of the abundance estimation.The similarity of the extracted endmember's spectra to the real spectra is evaluated via the indices of spectral angle distance (SAD) [52] and spectral information divergence (SID) [3].The accuracy of the abundance estimation is evaluated by the root mean square error (RMSE) [53].The equations of these indices are given as follows: • Spectral angle distance (SAD): where A is the actual spectral vector, B is the estimated spectral vector, and n is the number of bands.The real endmember spectra can be spectra in the spectral library, observed spectra from the laboratory or field, or spectra extracted from pure pixels in remote sensing imagery.
• Spectral information divergence (SID): where A is the actual spectral vector, and B is the estimated spectral vector.
• Root mean square error (RMSE): where S ij is the actual abundance of the j th endmember abundance of the i th mixed pixel, Ŝij is the estimated value of S ij and n is the number of bands.The RMSE metric can partly reflect how the inversion results reconstruct the original information.

Spectral Unmixing Using Synthetic Data
The hyperspectral synthetic data were generated using five spectral signatures randomly selected from the United States Geological Survey (USGS) digital spectral library [54], i.e., those for brucite, chabazite, olivine, spessartine and witherite (Figure 2).These spectral signatures are linearly independent and can be regarded as endmembers [16].All signatures are available online (http://speclab.cr.usgs.gov/spectral-lib.html).The spectral data contain 420 bands covering wavelengths from 395.1-2560 nm.To create the linear mixtures, the abundance of each endmember was randomly determined using the Dirichlet distribution under the ANC and ASC.Moreover, to simulate possible errors and sensor noise, zero mean white Gaussian noise with a 30-dB signal-to-noise ratio (SNR), as defined in [55], was added to the mixture dataset.was randomly determined using the Dirichlet distribution under the ANC and ASC.Moreover, to simulate possible errors and sensor noise, zero mean white Gaussian noise with a 30-dB signal-to-noise ratio (SNR), as defined in [55], was added to the mixture dataset.The abovementioned four algorithms, including the ATGP-NMF, VCA-NMF, MVC-NMF and VCA-FCLS algorithms, were first applied to the synthetic mixture data (SNR = 30 dB, purity level ρ = 0.8) to extract the endmember spectra (Figure 3) and calculate the corresponding abundances.The similarity of the extracted spectra to the real spectra was evaluated by SAD and SID (Figure 4).The accuracy of the abundance estimation was evaluated by RMSE (Table 1).The abovementioned four algorithms, including the ATGP-NMF, VCA-NMF, MVC-NMF and VCA-FCLS algorithms, were first applied to the synthetic mixture data (SNR = 30 dB, purity level = 0.8) to extract the endmember spectra (Figure 3) and calculate the corresponding abundances.The similarity of the extracted spectra to the real spectra was evaluated by SAD and SID (Figure 4).The accuracy of the abundance estimation was evaluated by RMSE (Table 1).
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 7 of 19 was randomly determined using the Dirichlet distribution under the ANC and ASC.Moreover, to simulate possible errors and sensor noise, zero mean white Gaussian noise with a 30-dB signal-to-noise ratio (SNR), as defined in [55], was added to the mixture dataset.The abovementioned four algorithms, including the ATGP-NMF, VCA-NMF, MVC-NMF and VCA-FCLS algorithms, were first applied to the synthetic mixture data (SNR = 30 dB, purity level ρ = 0.8) to extract the endmember spectra (Figure 3) and calculate the corresponding abundances.The similarity of the extracted spectra to the real spectra was evaluated by SAD and SID (Figure 4).The accuracy of the abundance estimation was evaluated by RMSE (Table 1).As shown in Figures 3 and 4, the spectra extracted by the proposed ATGP-NMF algorithm are quite similar to the real spectra, with smooth curves and little noise, followed by those extracted by the MVC-NMF, VCA-NMF and VCA-FCLS algorithms.In this experiment, although a pure pixel does not exist (purity level ρ = 0.8), the purity level of the synthetic mixture data is higher.Therefore, the endmember spectra extracted by VCA-FCLS also resemble the real endmember spectra, with small SAD and SID values.To further evaluate the robustness of these four algorithms, we compared their performances under different conditions, in terms of the SNR, purity level ρ and number of endmembers p, by using synthetic data.The comparison results are shown in Figures 5-7.As shown in Figures 3 and 4, the spectra extracted by the proposed ATGP-NMF algorithm are quite similar to the real spectra, with smooth curves and little noise, followed by those extracted by the MVC-NMF, VCA-NMF and VCA-FCLS algorithms.In this experiment, although a pure pixel does not exist (purity level = 0.8), the purity level of the synthetic mixture data is higher.Therefore, the endmember spectra extracted by VCA-FCLS also resemble the real endmember spectra, with small SAD and SID values.As shown in Figures 3 and 4, the spectra extracted by the proposed ATGP-NMF algorithm are quite similar to the real spectra, with smooth curves and little noise, followed by those extracted by the MVC-NMF, VCA-NMF and VCA-FCLS algorithms.In this experiment, although a pure pixel does not exist (purity level ρ = 0.8), the purity level of the synthetic mixture data is higher.Therefore, the endmember spectra extracted by VCA-FCLS also resemble the real endmember spectra, with small SAD and SID values.To further evaluate the robustness of these four algorithms, we compared their performances under different conditions, in terms of the SNR, purity level ρ and number of endmembers p, by using synthetic data.The comparison results are shown in Figures 5-7.To further evaluate the robustness of these four algorithms, we compared their performances under different conditions, in terms of the SNR, purity level and number of endmembers p, by using synthetic data.The comparison results are shown in Figures 5-7.
ISPRS Int.J. Geo-Inf.2018, 7, 195 9 of 19 (1) Robustness to noise interference: In this experiment, tests of the algorithms' sensitivities to noise were carried out by changing the SNR level from 50 to 10 dB.The synthetic data consist of five endmembers, i.e., the spectral signatures selected from the USGS spectral library, at the purity level of 0.8, which means pure pixels do not exist.The SAD, SID and RMSE values of the four algorithms are shown in Figure 5a-c, respectively.The ATGP-NMF algorithm performed well, yielding lower SAD and SID values than the other algorithms at different SNR levels.With the decrease in SNR, the RMSE values of ATGP-NMF decreased and stabilized; for MVC-NMF, the RMSE values kept increasing.
(2) Robustness to mixing degree: Following the method used in [56], we generated mixed spectra at different purity levels.In this experiment, we changed the purity level from 1 to 0.6 with an SNR of 30 dB.The comparisons of the SAD, SID and RMSE values derived from the four algorithms are shown in Figure 6a-c.The SAD, SID and RMSE values of these four algorithms are inversely proportional to the purity level.At the purity level of one, i.e., when pure pixels exist, the algorithms of four methods all good solutions, with low SAD, SID and RMSE values.With regard to the purity level, the RMSE values of ATGP-NMF were lower overall than those of the other algorithms.
(3) Generalization of the number of endmembers: In this experiment, we generated mixed spectra with different numbers of endmembers to test these algorithms.The endmember numbers were set to 3, 5 and 10, while the SNR was set to 30 dB and the purity level was set to 0.8.When the number of endmembers was set to three, the mixed spectra were generated using the spectral signatures of brucite, chabazite and olivine.When the number of endmembers was set to 10, another five spectral signatures were selected from the USGS digital spectral library, in addition to the previously-selected five signatures.Figure 7a-c shows that the performances of VCA-NMF and ATGP-NMF are relatively better and more robust than those of the other algorithms for the tested number of endmembers.
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 9 of 19 (1) Robustness to noise interference: In this experiment, tests of the algorithms' sensitivities to noise were carried out by changing the SNR level from 50 to 10 dB.The synthetic data consist of five endmembers, i.e., the spectral signatures selected from the USGS spectral library, at the purity level of 0.8, which means pure pixels do not exist.The SAD, SID and RMSE values of the four algorithms are shown in Figure 5a-c, respectively.The ATGP-NMF algorithm performed well, yielding lower SAD and SID values than the other algorithms at different SNR levels.With the decrease in SNR, the RMSE values of ATGP-NMF decreased and stabilized; for MVC-NMF, the RMSE values kept increasing.
(2) Robustness to mixing degree: Following the method used in [56], we generated mixed spectra at different purity levels.In this experiment, we changed the purity level ρ from 1 to 0.6 with an SNR of 30 dB.The comparisons of the SAD, SID and RMSE values derived from the four algorithms are shown in Figure 6a-c.The SAD, SID and RMSE values of these four algorithms are inversely proportional to the purity level.At the purity level of one, i.e., when pure pixels exist, the algorithms of four methods all produced good solutions, with low SAD, SID and RMSE values.With regard to the purity level, the RMSE values of ATGP-NMF were lower overall than those of the other algorithms.
(3) Generalization of the number of endmembers: In this experiment, we generated mixed spectra with different numbers of endmembers to test these algorithms.The endmember numbers were set to 3, 5 and 10, while the SNR was set to 30 dB and the purity level was set to 0.8.When the number of endmembers was set to three, the mixed spectra were generated using the spectral signatures of brucite, chabazite and olivine.When the number of endmembers was set to 10, another five spectral signatures were selected from the USGS digital spectral library, in addition to the previouslyselected five signatures.Figure 7a-c shows that the performances of VCA-NMF and ATGP-NMF are relatively better and more robust than those of the other algorithms for the tested number of endmembers.

Spectral Unmixing Using Berlin HyMap Image
To evaluate the performance of the proposed ATGP-NMF algorithm on the real image, we conducted experiments on a subset of the Berlin-Urban-Gradient dataset, which mainly comprises a HyMap reflectance image, a library of material spectra and the reference land cover information [57].More information on this dataset can be found in [58][59][60].For this work, a subset with 80 × 80 pixels of the Berlin HyMap image (Figure 8a), acquired on 20 August 2009 at 9 m spatial resolution, was selected as the study data.The image 111 bands after removing noisy bands with spectral sampling distances between 13 and 17 nm, and spectral coverage from 440 to 2500 nm.Based on the reference land cover information, we considered five urban land cover types in the study area, i.e., roof, pavement, grass, tree and soil.The reference abundances of each land cover type (Figure 9), defined as the corresponding proportions at each pixel, were calculated based on the reference land cover map (Figure 8b).We selected the pixels with an area percentage of one and calculated their average spectral values as the reference spectrum of the land cover type.Because there were no pure pixels of the grass category, four pixels with an area percentage greater than 0.7 were selected instead.
The comparisons of the VCA-FCLS, VCA-NMF, MVC-NMF and ATGP-NMF algorithms on the Berlin HyMap image are shown in Figure 10 and Table 2. Figure 10a,b shows that the proposed ATGP-NMF algorithm provided five endmembers with relatively low SAD and SID values.Figure 11 shows the estimated abundance maps by ATGP-NMF.Comparing to the reference abundance maps (Figure 9), it is demonstrated that the ATGP-NMF algorithm can produce satisfactory abundance separation.For all four algorithms, the estimated spectra of grass are in bad accordance with the reference endmembers.In the study area, grasses are scattered mainly around the buildings and trees, which are susceptible to the shadows of the buildings and trees and the spectral similarity between the grasses and trees.

Spectral Unmixing Using Berlin HyMap Image
To evaluate the performance of the proposed ATGP-NMF algorithm on the real image, we conducted experiments on a subset of the Berlin-Urban-Gradient dataset, which mainly comprises a HyMap reflectance image, a library of material spectra and the reference land cover information [57].More information on this dataset can be found in [58][59][60].For this work, a subset with 80 × 80 pixels of the Berlin HyMap image (Figure 8a), acquired on 20 August 2009 at 9 m spatial resolution, was selected as the study data.The image contains 111 bands after removing noisy bands with spectral sampling distances between 13 and 17 nm, and spectral coverage from 440 to 2500 nm.Based on the reference land cover information, we considered five urban land cover types in the study area, i.e., roof, pavement, grass, tree and soil.The reference abundances of each land cover type (Figure 9), defined as the corresponding proportions at each pixel, were calculated based on the reference land cover map (Figure 8b).We selected the pixels with an area percentage of one and calculated their average spectral values as the reference spectrum of the land cover type.Because there were no pure pixels of the grass category, four pixels with an area percentage greater than 0.7 were selected instead.
The comparisons of the VCA-FCLS, VCA-NMF, MVC-NMF and ATGP-NMF algorithms on the Berlin HyMap image are shown in Figure 10 and Table 2. Figure 10a,b shows that the proposed ATGP-NMF algorithm provided five endmembers with relatively low SAD and SID values.Figure 11 shows the estimated abundance maps by ATGP-NMF.Comparing to the reference abundance maps (Figure 9), it is demonstrated that the ATGP-NMF algorithm can produce satisfactory abundance separation.For all four algorithms, the estimated spectra of grass are in bad accordance with the reference endmembers.In the study area, grasses are scattered mainly around the buildings and trees, which are susceptible to the shadows of the buildings and trees and the spectral similarity between the grasses and trees.

Spectral Unmixing Using Berlin HyMap Image
To evaluate the performance of the proposed ATGP-NMF algorithm on the real image, we conducted experiments on a subset of the Berlin-Urban-Gradient dataset, which mainly comprises a HyMap reflectance image, a library of material spectra and the reference land cover information [57].More information on this dataset can be found in [58][59][60].For this work, a subset with 80 × 80 pixels of the Berlin HyMap image (Figure 8a), acquired on 20 August 2009 at 9 m spatial resolution, was selected as the study data.The image contains 111 bands after removing noisy bands with spectral sampling distances between 13 and 17 nm, and spectral coverage from 440 to 2500 nm.Based on the reference land cover information, we considered five urban land cover types in the study area, i.e., roof, pavement, grass, tree and soil.The reference abundances of each land cover type (Figure 9), defined as the corresponding proportions at each pixel, were calculated based on the reference land cover map (Figure 8b).We selected the pixels with an area percentage of one and calculated their average spectral values as the reference spectrum of the land cover type.Because there were no pure pixels of the grass category, four pixels with an area percentage greater than 0.7 were selected instead.
The comparisons of the VCA-FCLS, VCA-NMF, MVC-NMF and ATGP-NMF algorithms on the Berlin HyMap image are shown in Figure 10 and Table 2. Figure 10a,b shows that the proposed ATGP-NMF algorithm provided five endmembers with relatively low SAD and SID values.Figure 11 shows the estimated abundance maps by ATGP-NMF.Comparing to the reference abundance maps (Figure 9), it is demonstrated that the ATGP-NMF algorithm can produce satisfactory abundance separation.For all four algorithms, the estimated spectra of grass are in bad accordance with the reference endmembers.In the study area, grasses are scattered mainly around the buildings and trees, which are susceptible to the shadows of the buildings and trees and the spectral similarity between the grasses and trees.

Spectral Unmixing Using a Hyperion Image
We used a Hyperion image to further verify the applicability of the proposed ATGP-NMF algorithm in complicated situations.The Hyperion image (Figure 12a), acquired on 19 January 2011, covers parts of Guangzhou and Foshan, both of which are cities in Guangdong Province, China.The test image of this study is a subset of 101 × 101 pixels with 30-m spatial resolution.In addition, the image contains 242 bands acquired in the spectral range of 356-2577 nm.To improve the spectral unmixing performance, a series of pretreatments was performed on this image, including radiometric calibration and atmospheric correction.In addition, we removed the low SNR bands and the water vapor absorption bands.A total of 158 bands (including Bands 8-57, 79-120, 133-165 and 188-220) were used in this experiment.
For the Hyperion image, the number of endmembers was estimated using HFC mentioned in Section 2.3.Considering the actual situation of the study area, we finally determined that there are eight endmembers in the study area.In this study, the reference abundances were extracted from a high-resolution remote sensing image acquired from Google Earth (Figure 12b).The classification map of the eight land cover types of the study area (Figure 12c) was then acquired through supervised classification of the high-resolution remote sensing image using the k-nearest neighbor [61].Finally, the area percentage of each land cover type in one 30 × 30 m pixel was calculated and considered as its endmember abundance (Figure 13).In addition, due to the variability of the endmember spectra, we selected reference endmember spectra from the Hyperion image.We selected 30 samples for each land cover type and used the average spectra of these 30 samples as the reference spectra for each land cover type.
The comparison of the SAD and SID values of the four algorithms is shown in Figure 14, while the errors of the abundance estimation are given in Table 3.As shown in Figure 14a,b, the proposed ATGP-NMF algorithm provided eight endmembers with relatively low SAD and SID values.Compared to those of the other three algorithms, the spectra extracted by ATGP-NMF are closer to those of the actual land cover types.The ATGP-NMF algorithm can extract endmembers with high correlations, such as buildings with white roofs and blue roofs.Figure 15 shows the abundance maps estimated by ATGP-NMF.The distribution of the estimated endmembers is close to the distribution of their corresponding land cover types.Table 3 shows that the ATGP-NMF algorithm yielded the most accurate abundance estimation results with the lowest RMSE for the hyperspectral image.This is primarily because the ATGP-NMF algorithm can effectively and efficiently find the most representative endmember spectra from hyperspectral images.Among the four different algorithms, the traditional VCA-FCLS algorithm had the poorest performance, which again proves the advantages of blind unmixing methods.

Spectral Unmixing Using a Hyperion Image
We used a Hyperion image to further verify the applicability of the proposed ATGP-NMF algorithm in complicated situations.The Hyperion image (Figure 12a), acquired on 19 January 2011, covers parts of Guangzhou and Foshan, both of which are cities in Guangdong Province, China.The test image of this study is a subset of 101 × 101 pixels with 30-m spatial resolution.In addition, the image contains 242 bands acquired in the spectral range of 356-2577 nm.To improve the spectral unmixing performance, a series of pretreatments was performed on this image, including radiometric calibration and atmospheric correction.In addition, we removed the low SNR bands and the water vapor absorption bands.A total of 158 bands (including Bands 8-57, 79-120, 133-165 and 188-220) were used in this experiment.
For the Hyperion image, the number of endmembers was estimated using HFC mentioned in Section 2.3.Considering the actual situation of the study area, we finally determined that there are eight endmembers in the study area.In this study, the reference abundances were extracted from a high-resolution remote sensing image acquired from Google Earth (Figure 12b).The classification map of the eight land cover types of the study area (Figure 12c) was then acquired through supervised classification of the high-resolution remote sensing image using the k-nearest neighbor [61].Finally, the area percentage of each land cover type in one 30 × 30 m pixel was calculated and considered as its endmember abundance (Figure 13).In addition, due to the variability of the endmember spectra, we selected reference endmember spectra from the Hyperion image.We selected 30 samples for each land cover type and used the average spectra of these 30 samples as the reference spectra for each land cover type.
The comparison of the SAD and SID values of the four algorithms is shown in Figure 14, while the errors of the abundance estimation are given in Table 3.As shown in Figure 14a,b, the proposed ATGP-NMF algorithm provided eight endmembers with relatively low SAD and SID values.Compared to those of the other three algorithms, the spectra extracted by ATGP-NMF are closer to those of the actual land cover types.The ATGP-NMF algorithm can extract endmembers with high correlations, such as buildings with white roofs and blue roofs.Figure 15 shows the abundance maps estimated by ATGP-NMF.The distribution of the estimated endmembers is close to the distribution of their corresponding land cover types.Table 3 shows that the ATGP-NMF algorithm yielded the most accurate abundance estimation results with the lowest RMSE for the hyperspectral image.This is primarily because the ATGP-NMF algorithm can effectively and efficiently find the most representative endmember spectra from hyperspectral images.Among the four different algorithms, the traditional VCA-FCLS algorithm had the poorest performance, which again proves the advantages of blind unmixing methods.

Computational Efficiency
To evaluate the computational efficiency of the particular algorithm, the computing time was often used as an important measure in many related studies [35].The computational efficiency of an algorithm can be defined as a numerical function T(n): time versus the input size n.In this study, since the input data are the same, we used the average computing time (in seconds) of each algorithm's application on synthetic data, the Berlin HyMap image, and Hyperion image to evaluate their computational efficiency.The four algorithms were executed in Mathworks MATLAB R2014b on a computer with Intel(R) Xeon(R) E3-1246 CPU (3.50 GHz) and 12 GB RAM.For the three NMF derivative algorithms, i.e., VCA-NMF, MVC-NMF and ATGP-NMF, the endmember spectral and abundance estimation results were determined at the maximum iteration number of 300.
Table 4 shows that the computing time of the proposed ATGP-NMF algorithm is less than those of the traditional two-stage VCA-FCLS algorithm and the MVC-NMF algorithm.The MVC-NMF algorithm required more processing time than the other three algorithms, particularly for processing the real image.With the endmember initialization, the proposed ATGP-NMF and the existing VCA-NMF required less computation time.Table 4 also shows that ATGP-generated initial endmembers can significantly reduce the computational time of the spectral unmixing process.This is mainly because ATGP only needs to process previously-generated target pixels, not the whole dataset, to extract endmembers [42].

Computational Efficiency
To evaluate the computational efficiency of the particular algorithm, the computing time was often used as an important measure in many related studies [35].The computational efficiency of an algorithm can be defined as a numerical function T(n): time versus the input size n.In this study, since the input data are the same, we used the average computing time (in seconds) of each algorithm's application on synthetic data, the Berlin HyMap image, and Hyperion image to evaluate their computational efficiency.The four algorithms were executed in Mathworks MATLAB R2014b on a computer with Intel(R) Xeon(R) E3-1246 CPU (3.50 GHz) and 12 GB RAM.For the three NMF derivative algorithms, i.e., VCA-NMF, MVC-NMF and ATGP-NMF, the endmember spectral and abundance estimation results were determined at the maximum iteration number of 300.
Table 4 shows that the computing time of the proposed ATGP-NMF algorithm is less than those of the traditional two-stage VCA-FCLS algorithm and the MVC-NMF algorithm.The MVC-NMF algorithm required more processing time than the other three algorithms, particularly for processing the real image.With the endmember initialization, the proposed ATGP-NMF and the existing VCA-NMF required less computation time.Table 4 also shows that ATGP-generated initial endmembers can significantly reduce the computational time of the spectral unmixing process.This is mainly because ATGP only needs to process previously-generated target pixels, not the whole dataset, to extract endmembers [42].

Conclusions
The traditional NMF algorithm suffers from the "local best" problem, which greatly limits its applications.In this paper, via the use of ATGP-generated initial target endmembers, the ATGP-NMF algorithm is proposed to improve the performance of NMF in hyperspectral unmixing.The ATGP algorithm is introduced to provide more accurate initial points to identify endmembers, which can then be used by the NMF to obtain the global solutions for hyperspectral unmixing.Furthermore, the proposed ATGP-NMF algorithm requires little prior knowledge and can be applied in complex situations, even if there are no pure pixels.
Experiments were carried out to evaluate the performances of ATGP-NMF and three other algorithms, namely, VCA-NMF, MVC-NMF, and VCA-FCLS, using synthetic hyperspectral data and real hyperspectral images.Based on analysis of the experimental results, the following conclusions can be made: (1) The traditional two-stage spectral unmixing method, VCA-FCLS, can only achieve satisfactory results when pure pixels exist.When there are little prior knowledge and no pure pixels, the studied BSS methods, namely, VCA-NMF and ATGP-NMF, and MVC-NMF without the pure pixel assumption can achieve better unmixing results than VCA-FCLS.(2) Of all the algorithms, the ATGP-NMF algorithm demonstrated the best applicability and consistently produced good results in different situations.The extracted endmember spectra were closer to the real spectra, and the abundance estimation was more accurate.(3) Proper endmember initialization can help the NMF algorithm yield better estimation results with less computation time.

Figure 1 .
Figure 1.Flowchart of the proposed blind unmixing method improved by ATGP-generated initial target endmembers.

Figure 1 .
Figure 1.Flowchart of the proposed blind unmixing method improved by ATGP-generated initial target endmembers.

Figure 2 .
Figure 2. Reflectance spectra of the five materials selected from the USGS spectral library.

Figure 2 .
Figure 2. Reflectance spectra of the five materials selected from the USGS spectral library.

Figure 2 .
Figure 2. Reflectance spectra of the five materials selected from the USGS spectral library.

Figure 8 .Figure 8 .
Figure 8.(a) A subset of the Berlin HyMap image (R = 833 nm, G = 1652 nm, B = 632 nm); (b) the reference land cover map of the study area.

Figure 12 .
Figure 12.(a) Hyperion image of the study area (R = 936 nm, G = 2143 nm, B = 2345 nm); (b) the highresolution remote sensing image covering the study area; (c) the land use classification map of the study area.

Figure 12 .Figure 13 .Figure 14 .
Figure 12.(a) Hyperion image of the study area (R = 936 nm, G = 2143 nm, B = 2345 nm); (b) the high-resolution remote sensing image covering the study area; (c) the land use classification map of the study area.

Table 1 .
Spectral angle distance (SAD), spectral information divergence (SID) and RMSE values of the four algorithms on synthetic data.

Table 1 .
Spectral angle distance (SAD), spectral information divergence (SID) and RMSE values of the four algorithms on synthetic data.

Table 1 .
Spectral angle distance (SAD), spectral information divergence (SID) and RMSE values of the four algorithms on synthetic data.

Table 2 .
SAD, SID and RMSE values of different algorithms with the Berlin HyMap image.

Table 2 .
SAD, SID and RMSE values of different algorithms with the Berlin HyMap image.

Table 2 .
SAD, SID and RMSE values of different algorithms with the Berlin HyMap image.

Table 2 .
SAD, SID and RMSE values of different algorithms with the Berlin HyMap image.

Table 3 .
SAD, SID and RMSE values of different algorithms with the Hyperion image.

Table 3 .
SAD, SID and RMSE values of different algorithms with the Hyperion image.

Table 4 .
Computing time (seconds) of different methods on synthetic data, the Berlin HyMap image, and Hyperion Image.