- Email: [email protected]
A Gaussian elimination based fast endmember extraction algorithm for hyperspectral imagery
ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Contents lists available at SciVerse ScienceDirect
ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs
A Gaussian elimination based fast endmember extraction algorithm for hyperspectral imagery Xiurui Geng a,⇑, Zhengqing Xiao b, Luyan Ji a, Yongchao Zhao a, Fuxiang Wang c a
Key Laboratory of Technology in Geo-spatial Information Processing and Application System, Institute of Electronics, Chinese academy of Sciences, Beijing 100080, China College of Resources Science & Technology, Beijing Normal University, Beijing 100086, China c School of Electronics and Information Engineering, Beihang University, 37 XueYuan Road, Beijing 100191, China b
a r t i c l e
i n f o
Article history: Received 27 February 2012 Received in revised form 28 February 2013 Accepted 28 February 2013
Keywords: Hyperspectral data Endmember Gaussian elimination Simplex
a b s t r a c t A fast endmember-extraction algorithm based on Gaussian Elimination Method (GEM) is proposed in this paper under the fact that a pixel is an endmember if it has the maximum value in any spectral band of a hyperspectral image when based on linear mixing model. Applying Gaussian elimination is much like performing a lower triangular matrix to transform the hyperspectral image. As more endmembers have been extracted, fewer bands are needed to be involved in the Gaussian elimination process, thus greatly reducing the computing time. The experimental results with both simulated and real hyperspectral images indicate that the method proposed here is much faster than the vertex component analysis (VCA) method, and can provide a similar performance with VCA. Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved.
1. Introduction The mixed pixel problem has been a critical issue in the application of hyperspectral remotely sensed data. Using traditional classification approaches, it is difficult to quantitatively interpret mixed pixels. Recently, research based on mixed pixel analysis has made great progress (Somers et al., 2011). Most often a linear mixture model (LMM) is assumed due to its simplicity (Adams et al., 1989, 1986; Asner and Heidebrecht, 2003; Combe et al., 2008; Full et al., 1982; Gill and Phinn, 2009; Gong et al., 1994; Roberts et al., 1998; Smith et al., 1990). The LMM assumes that any pixel in an image is a linear combination of endmembers, and their abundance fractions satisfy the non-negative and sum-to-one constraints. When all endmembers are available, any pixel in the image can be solved quantitatively by linear unmixing. Therefore, the challenge is how to determine all the endmembers from an image based on LMM. If the LMM holds, all the pixels of an image are found within a simplex in feature space, and the vertexes of the simplex are endmembers corresponding to the pure pixels of surface objects. Based on this geometrical feature, a series of endmember extraction algorithms have been developed, such as Pixel Purity Index (PPI)(Boardman, 1993), NFINDR(Winter, 1999), Orthogonal Bases Algorithm (OBA)(Tao et al., 2009), IEA(Neville et al., 1999), Simplex Growing Algorithm (SGA) (Chein et al., 2006), Successive
⇑ Corresponding author. E-mail address: [email protected] (X. Geng).
Projection Algorithm (SPA) (Zhang et al., 2008), and Piecewise Convex Endmember Detection (PCE) (Zare and Gader, 2010). Nascimento and Dias presented a fast endmember extraction method, named vertex component analysis (VCA) (Nascimento and Dias, 2003, 2005). Compared to other methods, VCA is much faster in endmember extraction (Martin and Plaza, 2011; Shaohui et al., 2010). The VCA algorithm is also based on the fact that each endmember is located at a vertex of the simplex. By introducing a projection matrix, all the endmembers can be extracted sequentially. Xia et al. (2011) presented a method for maximizing the simplex volume, called Simplex Volume Analysis based on Triangular Factorization (SVATF). SVATF uses multiplicative operations to replace the volume calculation of N-FINDR, which reduces the computational complexity significantly. Unlike VCA, SVATF does not have to perform dimensionality reduction (DR). Liu and Zhang (2012) proposed a new endmember extraction algorithm, called maximum volume based on Householder transformation (MVHT). The MVHT algorithm improves the SGA algorithm by including the Householder transformation when searching the simplex vertex, which results in a tremendous reduction of computational complexity. Compared to VCA, the MVHT algorithm provides consistent results, but its computational complexity is a little higher (Liu and Zhang, 2012). In this paper, we introduce the Gaussian Elimination Method (GEM) for endmember extraction. Like most simplex-based endmember searching methods, our method also assumes the existence of pure pixels for all the endmembers in the image. It is based on the fact that under LMM, the pixel owning a maximum value in any band of the hyperspectral image is an endmember. In
0924-2716/$ - see front matter Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.isprsjprs.2013.02.020
212
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
order to extract all of the endmembers, GEM conducts a Gaussian elimination on the data set after each new endmember has been extracted. Compared to VCA, GEM uses one less band after each endmember is extracted, and also avoids the operation of matrix inversion. As a result, GEM has a lower computational complexity. Besides, the proposed GEM can provide consistent consults.
Component Analysis (ICA). Here the most commonly used DR technique, PCA is selected for GEM. The detailed procedures are given in Algorithm 1. Symbol [X]:,i:j stands for the ith to jth columns of matrix X; XT stands for the transposition of X. Step 1
Perform the PCA DR on the original data. M is the total number of vectors in the image. The projection matrix U is obtained by PCA from RRT/ M, where R = [r1, r2, . . . , rM]. The initialization of the endmember matrix E, the projection vector f, and the complete Gaussian elimination matrix G. In order to determine the ith endmember, the data matrix Z is firstly projected onto a vector f in Step 5. The vector f is equal to f = (uTG)T (Step 11), where G = GGi (Steps 9–10) and u is the band selection vector (Step 8). Thus y = uTGZ. When extracting the ith endmember, G is the product of i elementary transform matrixes, i.e. G = Gi G2 G1. GZ is to perform i1 times Gaussian elimination on the data matrix Z, so as to remove the influence of the i1 extracted endmembers in Z. And the function of u is to select the ith row of eliminated data matrix GZ. Then, in the Step 6, search the maximum position in |y| for the ith endmember. It is easy to validate that matrix G is a lower triangular matrix. When i = 1, f = [1,0, . . ., 0], Steps 3, 5 and 6 actually mean to search the maximum absolute value at the first band of the image after PC transform.
2. Gaussian elimination method 2.1. Linear mixture model Steps 2–4 If the assumption of LMM holds, the spectral vector r for a given pixel in an image can be expressed by a linear combination of several endmember vectors ei, i = 1, . . ., p with the constrained conditions (2):
r¼xþn¼
p X ci ei þ n ¼ Ec þ n
Steps 5–11
ð1Þ
i¼1 p X ci ¼ 1; ci P 0
ð2Þ
i¼1
where p is the number of endmembers in the image, and ci is a scalar representing the fractional abundance of endmember vector ei in the pixel r. E = [e1,e2, . . ., ep] is a L p mixing matrix (L is the number of bands of the original data). 2.2. Gaussian elimination method GEM is based on the fact that under the LMM, a pixel is considered as an endmember if it has a maximum value in any band of the hyperspectral image (referring to Appendix A). Thus we can extract the endmembers by simply searching the maximum position in each band. However, we may not find all of the endmembers in this way. Because some endmembers in the original image may have extreme values in more than one band or even across all bands, while some may have no maximum value in any band at all. To solve this problem we introduce the Gaussian elimination method as follows. The main idea of GEM is to perform the Gaussian elimination on the data set after extracting a new endmember. In the first step, we find the first endmember by searching the maximum position in the first band. Then we use Gaussian elimination on the hyperspectral image remove the first endmembers spectral signature from all other spectra. Next, the algorithm searches the position with maximum absolute value in the second band of the eliminated hyperspectral image, and then performs Gaussian elimination again on the eliminated hyperspectral image to remove the second endmember’s spectral values on all other spectra. During this Gaussian elimination, the first band of the eliminated hyperspectral image actually does not need to participate in the computation. And as more endmembers have been extracted, fewer bands are needed to be involved in the Gaussian elimination process, which will greatly reduce the computing time. This process iterates until all the endmembers have been found. Each step of Gaussian elimination is equivalent to conducting the linear transformation against the original hyperspectral image. The simplex will remain as a simplex after the linear transformation. Moreover, the Gaussian elimination in each iteration ensures that one endmember will be extracted only once. Therefore, all the endmembers can be obtained by GEM. Usually, the number of endmembers present in a given scene is much smaller than the number of bands. Therefore DR is necessary for many endmember extraction methods for both computational time saving and signal-to-noise ratio (SNR) improvements. Many algorithms can be used for DR, such as Principle Component Analysis (PCA), Minimum Noise Fraction (MNF), and Independent
Algorithm 1. Gaussian Elimination Method (GEM)
INPUT: p and R = [r1r2 . . . rM] (1) Z := UTR {U is obtained by PCA} (2) E := 0; {E is the p p mixing matrix}; (3) f :¼ ½1; 0; ; 0T ; {f is the projection vector} |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} p
(4) (5) (6)
G := Ipp;{G is the complete Gaussian elimination matrix} for i := 1to p do y := fTZ; k :¼ arg max jyj j; {find the projection extreme} j¼1;2;;M
(7) (8)
[E]:,i = [Z]:,k;
(9)
Let ei ¼ ½E:;i , g ¼ Gei , 2 1 0 .. .. . 6 6 0 .. . . 6. .. 6 .. 0 1 . 6 .. . .. Gi ¼ 6 g iþ1 6 .. . g . 6 i 6. .. .. 6. 4. . . 0 gp 0 0 0 g
u :¼ ½0; ; 0; 1 ; 0 ; ; 0 T 1
i iþ1 iþ2
i
p
.. . .. . .. . .. . 0
3 0 .. 7 .7 .. 7 .7 7 .. 7; fGi is the ith .7 7 7 7 05 1
elementary Gaussian elimination matrix} (10) Update G := GGi; {G is (i1)th complete Gaussian elimination matrix} (11) f := (uTG)T; endfor
213
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Compared to VCA, GEM is guaranteed to produce consistent results. The projection vector of VCA for each endmember is randomly selected in the orthogonal space spanned by endmembers already extracted, so the result of VCA is unrepeatable (Liu and Zhang, 2012). In contrast, the initial projection vector of GEM is fixed (referring to Step 3). And the rest projection vectors are equal to f=(uTG)T, where u is a fixed vector in the Algorithm 1. For a fixed data set, Gi in each iteration is fixed, so G is also fixed. Therefore, for a fixed data set, the result of GEM is constant for different runs, and even the extraction order is fixed. It should be pointed out that the endmember corresponding to one surface object type will be extracted only one time in the process of Gaussian elimination. There will be new endmember generated during each iteration. Therefore, all of the endmembers can be obtained automatically via GEM. Of course, like PPI, NFINDR, VCA and MVHT, the endmembers obtained by GEM are more mathematically than physically meaningful. Furthermore, GEM also makes the assumption that there are pure pixels for all the surface objects in the image. Since LMM is only an approximation for real imagery, the extracted endmembers are also the approximation of the real endmembers. 2.3. The time consumption of GEM Here, the computational complexity of VCA, MVHT and GEM is measured by the number of floating point operations (ignoring the lower order items). Ignoring the computation time of PCA, the time consumption of GEM during each iteration mainly includes two parts: (1) the projection of the data Z on the projection vector f (referring to the Step 5); (2) the generation of the projection vector f (referring to the Steps 10–11). For the first part, since the projection vector of GEM (noted as fGEM for clarity) has only i non-zero values in the ith iteration, the computational complexity of GEM is 2iM. However, since all of the elements of the projection vector of VCA (noted as fVCA) are non-zero, the computation of fVCA has a complexity of 2pM flops. Therefore, for all the p endmembers, the total computational complexity of the first part of GEM is p2M. And the total computational complexity of VCA is 2p2M, which is two times of that of GEM. For the second part, the calculation of fVCA needs to perform not only the matrix multiplication but also the pseudo-inversion operation, while the calculation of fGEM only needs to perform the multiplication of a lower triangular and a sparse matrix. After a little algebra, we can get the computational P 2 2 3 complexity of this part is p4 for GEM and 2 pk¼1 ðkp þ 2k p þ k Þ for VCA. Combining these two parts, we can calculate the total computational complexity tabulated in Table 1. It can be shown that computational savings using GEM are quite evident. The detailed computational analysis of MVHT can be found in Liu and Zhang (2012). 3. Experiment 3.1. Experiment with simulated data In the first experiment we compared the performance of the three methods, namely, VCA, GEM, and MVHT using simulated hyperspectral data. All of the algorithms were implemented using Matlab 6.5b; the system used was based on Intel (R) Core (TM) 2 Duo CPU T8300 @ 2.40 GHz, Hard Disk of 160 GB and Memory of
2 GB. To evaluate the performance of these methods, the metric of rmsSAE and rmsSID (Nascimento and Dias, 2005) was used to measure the similarity of the real endmembers and endmembers extracted from VCA, GEM and MVHT respectively. We used four reference spectra, Alunite, Calcite, Desert Varnish and Kaolinite from the U.S. Geological Survey (USGS) digital spectral library (Clark et al., 1993) included in ENVI (the Environment for Visualizing Images) software, as shown in Fig. 1, to generate the mixing data with 1000 pixels. The corresponding abundance fractions were generated by Dirichlet distribution (Nascimento and Dias, 2005). A simulated white Gaussian noise with SNR defined by SNR = 10 log10(E(xTx)/E(nTn)) was added to the synthetic data mentioned above, where E(.) denotes the expectation operator. Since the results of the three algorithms will all be influenced by the data set used, a total of 50 hyperspectral data sets were simulated by adding noise under different SNR levels, and the average of the 50 results was used for comparison. Before endmember extraction, the dimensionality of the simulated data was reduced to 4 by PCA. Fig. 2 shows the performance of the three methods as function of SNR. For each SNR, both rmsSAE and rmsSID were calculated from an average of 50 runs. It can be seen that GEM provides a similar result with VCA, and the performance of both GEM and VCA is much better than that of MVHT. Meanwhile, the presence of noise degrades the performance of all algorithms, whether in term of rmsSAE or rmsSID. It is noticeable that though the values of all bands of the Desert Varnish are minimum (which can be regarded as dark object), GEM can extract it under all SNR levels. The next experiment was used to evaluate run times for GEM, VCA and MVHT. It was assumed that the number of bands and endmembers were equal in this experiment. Two sets of simulated data were generated. For the first one, p vectors with p dimensions (p ranged from 3 to 52) were generated randomly as endmembers with M = 1000; and for the second one, M (M ranged from 4000 to 20,000) pixels based on the LMM were generated with p = 3. In order to make sure that every endmember exists in the data set, the top p pixels are replaced by the p endmembers. GEM, VCA and MVHT can all extract endmembers correctly from the data set. Fig. 3 shows the computing time of GEM, VCA and MVHT as function of the number of endmembers (Fig. 3a) and pixels (Fig. 3b). We can see that compared to VCA and MVHT, GEM has obvious advantages in terms of run time. This was especially true when the number of endmembers is larger. For example, when the number of endmembers was greater than 20, the run time for GEM was significantly lower than it was for VCA and MVHT. Fig. 3b presents the relationship between the computational time and the number of pixels for the three methods when the number of the endmembers was fixed at 3. It is shown that the computational time of MVHT was much more than run times for VCA and GEM, while the computational speed of GEM is slightly better than that of VCA. The last experiment was used to evaluate the effect of DR on the performance of GEM. The simulated data in the first experiment was used here. Though in Algorithm 1, DR was performed before endmember extraction, GEM can run on data without DR. Theoretically, GEM can extract 4 endmembers for any 4 linearly independent bands of the original data. In this case, GEM may produce different results when starting on different bands. Fig. 4 shows the performance of GEM running on data with DR (using the first 4 PC bands) and without DR (using the first 4 bands of original
Table 1 Computational complexity of VCA, MVHT, and GEM. Algorithm
VCA
Complexity (flops)
P 2 2 3 2p2 M þ 2 pk¼1 ðkp þ 2k p þ k Þ
MVHT P 2ð p1 4p 4k þ 3ÞM k¼2
GEM p2 M þ p4
214
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
bands. From Fig. 4 we can see that for all SNR levels, the performance of GEM running on PC bands is always better than that on original data. It should be pointed out that GEM can extract all the p endmembers from any p linearly independent bands if the LMM holds and the data are noiseless. However, real hyperspectral scenes usually cannot fully satisfy the LMM and are easily corrupted by sensor noise. Therefore, to improve endmember extraction accuracy, GEM needs to perform DR first, which can improve the data SNR and at the same time remove the redundancies between bands. 3.2. Experiment with real data
Fig. 1. The four endmembers (Alunite, Calcite, Desert Varnish and Kaolinite) from USGS library used to generate the simulated data.
data). Since the first p PC bands retain much more information than any p bands of the original data, the results of GEM running on PC bands are usually better than that of GEM running on original
From the above experiments, we found that the computational complexity of VCA is much closer to that of GEM. As a result, in this section, a more detailed comparative analysis between GEM and VCA is presented with measured hyperspectral data, collected by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS) sensor over Cuprite, Nevada. AVIRIS is a high quality low noise imaging spectrometer, which acquires data in 224 contiguous spectral bands from 370 to 2510 nm (Green et al., 1998). The selected subscene is shown in Fig. 5 with a size of 190 250 pixels. Due to water vapor absorption or low SNR, bands 1–3, 107–113, 153–
Fig. 2. The performances of GEM, VCA and MVHT as function of SNR. (a) rmsSAE and (b) rmsSID.
Fig. 3. Computing time for VCA, MVHT and GEM. (a) M = 1000 and (b) p = 3.
215
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Fig. 4. The effect of PCA DR to the performance of GEM. (a) rmsSAE and (b) rmsSID.
Table 2 Spectral angle (in degrees) between USGS reference spectra and extracted endmembers by GEM and VCA in 3 different runs (VCA1 for the 1st run, VCA2 for the 2nd run, and VCA3 for the 3rd run).
Fig. 5. The band 8 (468.71 nm) of the real hyperspectral subscene.
169, 221–224 were removed. The resampled spectra from the USGS Digital Spectral Library included in ENVI software were selected as ground references for the comparison. According to the ground truth presented in reference (Swayze et al., 1992) and endmember extraction results in references (Liu and Zhang, 2012; Nascimento and Dias, 2005, 2007), there are less than 18 materials in this subscene. Thereby, the estimated number of endmembers was set to p = 18, and the dimensionality of the image was reduced to 18 by PCA. As GEM may find more than one endmember for one mineral, they are numbered successively, such as Alunite #1 and Alunite #2. For each mineral, the library spectrum, which has the minimum spectral angle degree (SAD) with the first endmember of that mineral, was selected as the reference for all the endmembers belonging to that mineral. Taking Alunite for example, there are 6 library spectra of Alunite in USGS spectral library, the one having the minimum SAD with Alunite #1 was chosen as the reference spectrum for both Alunite #1 and 2. Since VCA may produce different results for different runs, we ran VCA randomly for 3 times, noted as VCA1-3 in Table 2. Table 2 compares the SAD between extracted endmembers and USGS spectra quantitatively. We found that only 10 endmembers could be extracted by both the GEM and the 3 VCA runs. Compared to GEM, there were three endmembers missed by VCA1-3. It should be noted that GEM also missed some of the remaining VCA endmembers. For a fair comparison, the average SAD in Table 2 was calculated for the 10 common endmembers only. Within these 10 endmembers, 2 endmembers corresponded to the same minerals, so only 8 endmembers’ spectral curves are shown in Fig. 6.
Substance
GEM
VCA1
VCA2
VCA3
Alunite #1 Alunite #2 Andradite Buddingtonite Chalcedony #1 Chalcedony #2 Dumortierite #1 Dumortierite #2 Dumortierite #3 Kaolinite #1 Kaolinite #2 Montmorillonite Muscovite #1 Muscovite #2 Nontronite #1 Nontronite #2 Nontronite #3 Pyrope Average
4.3262 16.0481 5.7873 3.7080 5.4079 5.7812 9.4846 4.5170 9.3771 4.7599 9.8903 7.6854 3.9345 7.9540 3.5958 6.2838 4.8177 3.0012 6.8701
7.3982 18.5156 – 4.4668 – 6.1915 9.4846 – 8.6622 4.6720 – 5.7613 3.9345 – 3.5958 4.5381 – – 7.2682
6.2043 18.1011 – 5.0350 – 4.8733 9.4846 – 8.6622 4.9172 11.7588 6.8411 3.9345 – 3.6867 – – 5.4229 7.1740
15.0414 7.5699 6.8183 4.4668 – 5.9371 9.4846 – 8.6622 5.6708 11.1889 6.5376 3.9345 8.4721 5.5266 – – – 7.2832
From Table 2 we can see that GEM has a better performance for 3 endmembers (Alunite #1, Buddingtonite and Nontronite #1), an equal performance for 2 endmembers (Dumortierite #1 and Muscovite #1) and a worse performance for 2 endmembers (Montmorillonite and Dumortierite #3). For the remaining 3 endmembers (Alunite #2, Chalcedony #2 and Kaolinite #1), GEM performed much better than the worst VCA. Overall GEM produced a smaller average angle than VCA. Fig. 7 shows the abundance images of the 18 endmembers by the fully constrained least squares method in (Chang and Heinz, 2000). The abundance maps for the rest of the endmembers are in good agreement with the ground truth presented in reference(Swayze, et al., 1992) and the results presented in (Liu and Zhang, 2012; Nascimento and Dias, 2005, 2007). Fig. 8 shows the computing time of GEM and VCA versus p ranging from 10 to 18 using the real hyperspectral data. Similar to Fig. 3a, the advantage of GEM in computing time is more and more obvious as p increases.
4. Discussion and conclusions Based on the pure pixel assumption, we propose an endmember extraction method using Gaussian elimination method for
216
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Fig. 6. Comparison of endmember signatures extracted by (a) GEM, (b–d) VCA1-3, and (e) the corresponding laboratory spectra.
hyperspectral imagery, named Gaussian elimination method (GEM). For simulated data with p endmembers, GEM can extract
all the p endmembers from any p linearly independent bands based on LMM. However, real hyperspectral images usually do not satisfy
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
217
Fig. 7. Abundance maps of different endmembers from GEM. (a) Muscovite #1, (b) Alunite #1, (c) Dumortierite #1, (d) Kaolinite #1, (e) Montmorillonite, (f) Kaolinite #2, (g) Buddingtonite, (h) Nontronite #1, (i) Alunite #2, (j) Andradite, (k) Pyrope, (l) Nontronite #2, (m) Chalcedony #1, (n) Dumortierite #2, (o) Nontronite #3, (p) Muscovite #2, (q) Dumortierite #3 and (r) Chalcedony #2.
LMM and are easily corrupted by noise, so we need to perform DR before endmember extraction for better results. Besides, GEM assumes the existence of pure pixels for all the endmembers in an image. However, GEM might perform slightly worse than those endmember generation methods (Chan et al., 2009; Miao and Qi, 2007) if the pure-pixel assumption is not satisfied. In addition, it has to be pointed out that the performance of GEM was evaluated
on a relative small area of hyperspectral image. Both the performance and processing speed of GEM for large images will be studied in the future. Compared to the orthogonal matrix of VCA, the introducing of the Gaussian elimination matrix, which is a lower triangular matrix, can reduce the computational complexity of GEM a lot (referring to Table 1, Figs. 3 and 8). Another difference between
218
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Fig. 8. Computing time (s) for VCA and GEM versus the number of endmembers, p when using real hyperspectral image.
GEM and VCA is that GEM uses the fixed projection vector for each endmembers, i.e. GEM always searches endmembers from the first PC band to the pth PC bands. Thereby, GEM can produce repeatable results. Moreover, the performance of GEM is comparable to that of VCA (see Fig. 2 and Table 2). Acknowledgements This work was subsidized by the National High Technology Research and Development Program (863 Program) of China (Nos. 2009AA122101, 2009AA12ZX1486531). Appendix A The proof of the GEM is given as follows. Theorem. Suppose Ri = [ri1, ri2, . . ., riM] represents all the M pixels in ith (1 6 i 6 L) band of a hyperspectral image, and there exists only one maximum absolute value in Ri and denoted its index as J, that is absðriJ Þ ¼ max ðabsðrij ÞÞ, then rJ ¼ ½r 1J ; r 2J ; ; r LJ T is one of the 16j6M
endmembers of the image. Proof. If rJ is not an endmember, it must be a linear combination P of the endmembers of the image, i.e. rJ ¼ pk¼1 ck ek , where Pp k¼1 ck ¼ 1, 0 6 ck 6 1, fek ; 1 6 k 6 pg are all the endmembers of the image. P In the ith band of image, we have riJ ¼ pk¼1 ck eik , which can be Pp transformed to 1 ¼ k¼1 ck erikiJ by dividing both sides with riJ . Since riJ is the only maximum absolute value, there is at least one e ik k0 ð1 6 k0 6 pÞ which meets the condition that abs riJ0 < 1. Then Pp P P we have 1 ¼ absð k¼1 ck erikiJ Þ 6 pk¼1 absðck ÞabsðerikiJ Þ < pk¼1 ck ¼ 1 which results in the contradiction. Therefore rJ ¼ ½r1J ; r 2J ; ; rLJ T must be one of the endmembers. h References Adams, J.B., Smith, M.O., Johnson, P.E., 1986. Spectral mixture modeling: a new analysis of rock and soil types at the Viking Lander I site. Journal of Geophysical Research 91, 8098–8812. Adams, J.B., Smith, M.O., Gillespie, A.R., 1989. Simple models for complex natural surfaces: a strategy for the hyperspectral era of remote sensing. In: Paper presented at the IGARSS’89, Vancouver, BC, Canada.
Asner, G.P., Heidebrecht, K.B., 2003. Imaging spectroscopy for desertification studies: comparing AVIRIS and EO-1 Hyperion in Argentina drylands. IEEE Transactions on Geoscience and Remote Sensing 41 (6), 1283–1296. Boardman, J.W., 1993. Automated spectral unmixing of AVIRIS data using convex geometry concepts. In: Paper presented at the 4th JPL Airborne Geoscience Workshop. Chan, T.-H., Chi, C.-Y., Huang, Y.-M., Ma, W.-K., 2009. A convex analysis-based minimum-volume enclosing simplex algorithm for hyperspectral ummixing. IEEE Transaction on Signal Processing 57, 4418–4432. Chang, C.I., Heinz, D.C., 2000. Constrained subpixel target detection for remotely sensed imagery. IEEE Transactions on Geoscience and Remote Sensing 38 (3), 1144–1159. Chein, I.C., Chao-Cheng, W., Wei-min, L., Yen-Chieh, O., 2006. A new growing method for simplex-based endmember extraction algorithm. IEEE Transactions on Geoscience and Remote Sensing 44 (10), 2804–2819. Clark, R.N., Swayze, G., Gallagher, A., King, T.V., Calvin, W.M., 1993. The U.S. geological survey digital spectral library: Version 1: 0.2 to 3.0 um, U.S. Geological Survey, Washington, DC, Open File Rep. 93–592. Combe, J.P., Le Mouélic, S., Sotin, C., Gendrin, A., Mustard, J.F., Le Deit, L., et al., 2008. Analysis of OMEGA/Mars Express data hyperspectral data using a MultipleEndmember Linear Spectral Unmixing Model (MELSUM): Methodology and first results. Planetary and Space Science 56 (7), 951–975. http://dx.doi.org/10.1016/ j.pss.2007.12.007. Full, W.E., R, E., Bezdek, J.C., 1982. Fuzzy Qmodel- a new model approach for linear mixing. Mathematical Geology 14 (3), 259–270. Gill, T.K., Phinn, S.R., 2009. Improvements to ASTER-derived fractional estimates of bare ground in a Savanna Rangeland. IEEE Transactions on Geoscience and Remote Sensing 47 (2), 662–670. Gong, P., Miller, J.R., Spanner, M., 1994. Forest canopy closure from classification and spectral unmixing of scene components-multisensor evaluation of an open canopy. IEEE Transactions on Geoscience and Remote Sensing 32 (5), 1067– 1080. Green, R.O., Eastwood, M.L., Sarture, C.M., Chrien, T.G., Aronsson, M., Chippendale, B.J., et al., 1998. Imaging spectroscopy and the airborne visible/infrared imaging spectrometer (AVIRIS). Remote Sensing of Environment 65, 227–248. Liu, J., Zhang, J., 2012. A new maximum simplex volume method based on householder transformation for endmember extraction. IEEE Transactions on Geoscience and Remote Sensing 50 (1), 104–118. Martin, G., Plaza, A., 2011. Region-based spatial preprocessing for endmember extraction and spectral unmixing. IEEE Geoscience and Remote Sensing Letters 99, 745–749. Miao, L., Qi, H., 2007. Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization. IEEE Transactions on Geoscience and Remote Sensing 45 (3), 765–777. Nascimento, J.M.P., Dias, J.M.B., 2003. Vertex component analysis: a fast algorithm to extract endmembers spectra from hyperspectral data. In: Paper presented at the Pattern Recognition and Image Analysis: First Iberian Conference. Nascimento, J.M.P., Dias, J.M.B., 2005. Vertex component analysis: a fast algorithm to unmix hyperspectral data. IEEE Transactions on Geoscience and Remote Sensing 43 (4), 898–910. Nascimento, J.M.P., Dias, J.M.B., 2007. Hyperspectral Signal Subspace Estimation. In: Paper presented at the Proceedings of IGARSS. Neville, R.A., Staenz, K., Szeredi, T., Lefebvre, J., Hauff, P., 1999. Automatic endmember extraction from hyperspectral data for mineral exploration. In: Paper presented at the Proceedings of 21st Canadian Symposim Remote Sensing. Roberts, D.A., Gardmer, M., Church, R., Ustin, S., Scheer, G., Green, R.O., 1998. Mapping chaparral in the Santa Monica mountains using multiple endmember spectral mixture models. Remote Sensing of Environment 65, 267–279. Shaohui, M., Mingyi, H., Zhiyong, W., Dagan, F., 2010. Spatial purity based endmember extraction for spectral mixture analysis. IEEE Transactions on Geoscience and Remote Sensing 48 (9), 3434–3445. Smith, M.O., Ustin, S.L., Adams, J.B., Gillespie, A.R., 1990. Vegetation in deserts: I. A regional measure of abundance from multispectral images. Remote Sensing of Environment 31, 1–52. Somers, B., Asner, G.P., Tits, L., Coppin, P., 2011. Endmember variability in spectral mixture analysis: a review. Remote Sensing of Environment 115, 1603–1616. Swayze, G., Clark, R., Sutley, S., Gallagher, A., 1992. Ground-truthing AVIRIS mineral mapping at cuprite, Navada. In: Paper presented at the Proceedings of Summaries Third ANNu. JPL Airborne Geoscience Workshop. Tao, X., Wang, B., Zhang, L., 2009. Orthogonal bases approach for the decomposition of mixed pixels in hyperspectral imagery. IEEE Geoscience and Remote Sensing Letters 6 (2), 219–223. Winter, M.E., 1999. Fast autonomous spectral endmember determination in hyperspectral data. In: Paper Presented at the Proceedings of SPIE. Xia, W., Wang, B., Zhang, L., Lu, Q., 2011. Simplex volume analysis based on triangular factorization: A framework for hyperspectral for hyperspectral unmixing. In: Paper presented at the Geoscience and Remote Sensing Symposium (IGARSS), IEEE International. Zare, A., Gader, P., 2010. PCE: piecewise convex endmember detection. IEEE Transactions on Geoscience and Remote Sensing 48 (6), 2620–2632. Zhang, J.K., Rivard, B., Rogge, D.M., 2008. The successive projection algorithm (SPA), an algorithm with a spatial constraint for the automatic search of endmembers in hyperspectral data. Sensors 8, 1321–1342.
Contents lists available at SciVerse ScienceDirect
ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs
A Gaussian elimination based fast endmember extraction algorithm for hyperspectral imagery Xiurui Geng a,⇑, Zhengqing Xiao b, Luyan Ji a, Yongchao Zhao a, Fuxiang Wang c a
Key Laboratory of Technology in Geo-spatial Information Processing and Application System, Institute of Electronics, Chinese academy of Sciences, Beijing 100080, China College of Resources Science & Technology, Beijing Normal University, Beijing 100086, China c School of Electronics and Information Engineering, Beihang University, 37 XueYuan Road, Beijing 100191, China b
a r t i c l e
i n f o
Article history: Received 27 February 2012 Received in revised form 28 February 2013 Accepted 28 February 2013
Keywords: Hyperspectral data Endmember Gaussian elimination Simplex
a b s t r a c t A fast endmember-extraction algorithm based on Gaussian Elimination Method (GEM) is proposed in this paper under the fact that a pixel is an endmember if it has the maximum value in any spectral band of a hyperspectral image when based on linear mixing model. Applying Gaussian elimination is much like performing a lower triangular matrix to transform the hyperspectral image. As more endmembers have been extracted, fewer bands are needed to be involved in the Gaussian elimination process, thus greatly reducing the computing time. The experimental results with both simulated and real hyperspectral images indicate that the method proposed here is much faster than the vertex component analysis (VCA) method, and can provide a similar performance with VCA. Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved.
1. Introduction The mixed pixel problem has been a critical issue in the application of hyperspectral remotely sensed data. Using traditional classification approaches, it is difficult to quantitatively interpret mixed pixels. Recently, research based on mixed pixel analysis has made great progress (Somers et al., 2011). Most often a linear mixture model (LMM) is assumed due to its simplicity (Adams et al., 1989, 1986; Asner and Heidebrecht, 2003; Combe et al., 2008; Full et al., 1982; Gill and Phinn, 2009; Gong et al., 1994; Roberts et al., 1998; Smith et al., 1990). The LMM assumes that any pixel in an image is a linear combination of endmembers, and their abundance fractions satisfy the non-negative and sum-to-one constraints. When all endmembers are available, any pixel in the image can be solved quantitatively by linear unmixing. Therefore, the challenge is how to determine all the endmembers from an image based on LMM. If the LMM holds, all the pixels of an image are found within a simplex in feature space, and the vertexes of the simplex are endmembers corresponding to the pure pixels of surface objects. Based on this geometrical feature, a series of endmember extraction algorithms have been developed, such as Pixel Purity Index (PPI)(Boardman, 1993), NFINDR(Winter, 1999), Orthogonal Bases Algorithm (OBA)(Tao et al., 2009), IEA(Neville et al., 1999), Simplex Growing Algorithm (SGA) (Chein et al., 2006), Successive
⇑ Corresponding author. E-mail address: [email protected] (X. Geng).
Projection Algorithm (SPA) (Zhang et al., 2008), and Piecewise Convex Endmember Detection (PCE) (Zare and Gader, 2010). Nascimento and Dias presented a fast endmember extraction method, named vertex component analysis (VCA) (Nascimento and Dias, 2003, 2005). Compared to other methods, VCA is much faster in endmember extraction (Martin and Plaza, 2011; Shaohui et al., 2010). The VCA algorithm is also based on the fact that each endmember is located at a vertex of the simplex. By introducing a projection matrix, all the endmembers can be extracted sequentially. Xia et al. (2011) presented a method for maximizing the simplex volume, called Simplex Volume Analysis based on Triangular Factorization (SVATF). SVATF uses multiplicative operations to replace the volume calculation of N-FINDR, which reduces the computational complexity significantly. Unlike VCA, SVATF does not have to perform dimensionality reduction (DR). Liu and Zhang (2012) proposed a new endmember extraction algorithm, called maximum volume based on Householder transformation (MVHT). The MVHT algorithm improves the SGA algorithm by including the Householder transformation when searching the simplex vertex, which results in a tremendous reduction of computational complexity. Compared to VCA, the MVHT algorithm provides consistent results, but its computational complexity is a little higher (Liu and Zhang, 2012). In this paper, we introduce the Gaussian Elimination Method (GEM) for endmember extraction. Like most simplex-based endmember searching methods, our method also assumes the existence of pure pixels for all the endmembers in the image. It is based on the fact that under LMM, the pixel owning a maximum value in any band of the hyperspectral image is an endmember. In
0924-2716/$ - see front matter Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.isprsjprs.2013.02.020
212
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
order to extract all of the endmembers, GEM conducts a Gaussian elimination on the data set after each new endmember has been extracted. Compared to VCA, GEM uses one less band after each endmember is extracted, and also avoids the operation of matrix inversion. As a result, GEM has a lower computational complexity. Besides, the proposed GEM can provide consistent consults.
Component Analysis (ICA). Here the most commonly used DR technique, PCA is selected for GEM. The detailed procedures are given in Algorithm 1. Symbol [X]:,i:j stands for the ith to jth columns of matrix X; XT stands for the transposition of X. Step 1
Perform the PCA DR on the original data. M is the total number of vectors in the image. The projection matrix U is obtained by PCA from RRT/ M, where R = [r1, r2, . . . , rM]. The initialization of the endmember matrix E, the projection vector f, and the complete Gaussian elimination matrix G. In order to determine the ith endmember, the data matrix Z is firstly projected onto a vector f in Step 5. The vector f is equal to f = (uTG)T (Step 11), where G = GGi (Steps 9–10) and u is the band selection vector (Step 8). Thus y = uTGZ. When extracting the ith endmember, G is the product of i elementary transform matrixes, i.e. G = Gi G2 G1. GZ is to perform i1 times Gaussian elimination on the data matrix Z, so as to remove the influence of the i1 extracted endmembers in Z. And the function of u is to select the ith row of eliminated data matrix GZ. Then, in the Step 6, search the maximum position in |y| for the ith endmember. It is easy to validate that matrix G is a lower triangular matrix. When i = 1, f = [1,0, . . ., 0], Steps 3, 5 and 6 actually mean to search the maximum absolute value at the first band of the image after PC transform.
2. Gaussian elimination method 2.1. Linear mixture model Steps 2–4 If the assumption of LMM holds, the spectral vector r for a given pixel in an image can be expressed by a linear combination of several endmember vectors ei, i = 1, . . ., p with the constrained conditions (2):
r¼xþn¼
p X ci ei þ n ¼ Ec þ n
Steps 5–11
ð1Þ
i¼1 p X ci ¼ 1; ci P 0
ð2Þ
i¼1
where p is the number of endmembers in the image, and ci is a scalar representing the fractional abundance of endmember vector ei in the pixel r. E = [e1,e2, . . ., ep] is a L p mixing matrix (L is the number of bands of the original data). 2.2. Gaussian elimination method GEM is based on the fact that under the LMM, a pixel is considered as an endmember if it has a maximum value in any band of the hyperspectral image (referring to Appendix A). Thus we can extract the endmembers by simply searching the maximum position in each band. However, we may not find all of the endmembers in this way. Because some endmembers in the original image may have extreme values in more than one band or even across all bands, while some may have no maximum value in any band at all. To solve this problem we introduce the Gaussian elimination method as follows. The main idea of GEM is to perform the Gaussian elimination on the data set after extracting a new endmember. In the first step, we find the first endmember by searching the maximum position in the first band. Then we use Gaussian elimination on the hyperspectral image remove the first endmembers spectral signature from all other spectra. Next, the algorithm searches the position with maximum absolute value in the second band of the eliminated hyperspectral image, and then performs Gaussian elimination again on the eliminated hyperspectral image to remove the second endmember’s spectral values on all other spectra. During this Gaussian elimination, the first band of the eliminated hyperspectral image actually does not need to participate in the computation. And as more endmembers have been extracted, fewer bands are needed to be involved in the Gaussian elimination process, which will greatly reduce the computing time. This process iterates until all the endmembers have been found. Each step of Gaussian elimination is equivalent to conducting the linear transformation against the original hyperspectral image. The simplex will remain as a simplex after the linear transformation. Moreover, the Gaussian elimination in each iteration ensures that one endmember will be extracted only once. Therefore, all the endmembers can be obtained by GEM. Usually, the number of endmembers present in a given scene is much smaller than the number of bands. Therefore DR is necessary for many endmember extraction methods for both computational time saving and signal-to-noise ratio (SNR) improvements. Many algorithms can be used for DR, such as Principle Component Analysis (PCA), Minimum Noise Fraction (MNF), and Independent
Algorithm 1. Gaussian Elimination Method (GEM)
INPUT: p and R = [r1r2 . . . rM] (1) Z := UTR {U is obtained by PCA} (2) E := 0; {E is the p p mixing matrix}; (3) f :¼ ½1; 0; ; 0T ; {f is the projection vector} |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} p
(4) (5) (6)
G := Ipp;{G is the complete Gaussian elimination matrix} for i := 1to p do y := fTZ; k :¼ arg max jyj j; {find the projection extreme} j¼1;2;;M
(7) (8)
[E]:,i = [Z]:,k;
(9)
Let ei ¼ ½E:;i , g ¼ Gei , 2 1 0 .. .. . 6 6 0 .. . . 6. .. 6 .. 0 1 . 6 .. . .. Gi ¼ 6 g iþ1 6 .. . g . 6 i 6. .. .. 6. 4. . . 0 gp 0 0 0 g
u :¼ ½0; ; 0; 1 ; 0 ; ; 0 T 1
i iþ1 iþ2
i
p
.. . .. . .. . .. . 0
3 0 .. 7 .7 .. 7 .7 7 .. 7; fGi is the ith .7 7 7 7 05 1
elementary Gaussian elimination matrix} (10) Update G := GGi; {G is (i1)th complete Gaussian elimination matrix} (11) f := (uTG)T; endfor
213
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Compared to VCA, GEM is guaranteed to produce consistent results. The projection vector of VCA for each endmember is randomly selected in the orthogonal space spanned by endmembers already extracted, so the result of VCA is unrepeatable (Liu and Zhang, 2012). In contrast, the initial projection vector of GEM is fixed (referring to Step 3). And the rest projection vectors are equal to f=(uTG)T, where u is a fixed vector in the Algorithm 1. For a fixed data set, Gi in each iteration is fixed, so G is also fixed. Therefore, for a fixed data set, the result of GEM is constant for different runs, and even the extraction order is fixed. It should be pointed out that the endmember corresponding to one surface object type will be extracted only one time in the process of Gaussian elimination. There will be new endmember generated during each iteration. Therefore, all of the endmembers can be obtained automatically via GEM. Of course, like PPI, NFINDR, VCA and MVHT, the endmembers obtained by GEM are more mathematically than physically meaningful. Furthermore, GEM also makes the assumption that there are pure pixels for all the surface objects in the image. Since LMM is only an approximation for real imagery, the extracted endmembers are also the approximation of the real endmembers. 2.3. The time consumption of GEM Here, the computational complexity of VCA, MVHT and GEM is measured by the number of floating point operations (ignoring the lower order items). Ignoring the computation time of PCA, the time consumption of GEM during each iteration mainly includes two parts: (1) the projection of the data Z on the projection vector f (referring to the Step 5); (2) the generation of the projection vector f (referring to the Steps 10–11). For the first part, since the projection vector of GEM (noted as fGEM for clarity) has only i non-zero values in the ith iteration, the computational complexity of GEM is 2iM. However, since all of the elements of the projection vector of VCA (noted as fVCA) are non-zero, the computation of fVCA has a complexity of 2pM flops. Therefore, for all the p endmembers, the total computational complexity of the first part of GEM is p2M. And the total computational complexity of VCA is 2p2M, which is two times of that of GEM. For the second part, the calculation of fVCA needs to perform not only the matrix multiplication but also the pseudo-inversion operation, while the calculation of fGEM only needs to perform the multiplication of a lower triangular and a sparse matrix. After a little algebra, we can get the computational P 2 2 3 complexity of this part is p4 for GEM and 2 pk¼1 ðkp þ 2k p þ k Þ for VCA. Combining these two parts, we can calculate the total computational complexity tabulated in Table 1. It can be shown that computational savings using GEM are quite evident. The detailed computational analysis of MVHT can be found in Liu and Zhang (2012). 3. Experiment 3.1. Experiment with simulated data In the first experiment we compared the performance of the three methods, namely, VCA, GEM, and MVHT using simulated hyperspectral data. All of the algorithms were implemented using Matlab 6.5b; the system used was based on Intel (R) Core (TM) 2 Duo CPU T8300 @ 2.40 GHz, Hard Disk of 160 GB and Memory of
2 GB. To evaluate the performance of these methods, the metric of rmsSAE and rmsSID (Nascimento and Dias, 2005) was used to measure the similarity of the real endmembers and endmembers extracted from VCA, GEM and MVHT respectively. We used four reference spectra, Alunite, Calcite, Desert Varnish and Kaolinite from the U.S. Geological Survey (USGS) digital spectral library (Clark et al., 1993) included in ENVI (the Environment for Visualizing Images) software, as shown in Fig. 1, to generate the mixing data with 1000 pixels. The corresponding abundance fractions were generated by Dirichlet distribution (Nascimento and Dias, 2005). A simulated white Gaussian noise with SNR defined by SNR = 10 log10(E(xTx)/E(nTn)) was added to the synthetic data mentioned above, where E(.) denotes the expectation operator. Since the results of the three algorithms will all be influenced by the data set used, a total of 50 hyperspectral data sets were simulated by adding noise under different SNR levels, and the average of the 50 results was used for comparison. Before endmember extraction, the dimensionality of the simulated data was reduced to 4 by PCA. Fig. 2 shows the performance of the three methods as function of SNR. For each SNR, both rmsSAE and rmsSID were calculated from an average of 50 runs. It can be seen that GEM provides a similar result with VCA, and the performance of both GEM and VCA is much better than that of MVHT. Meanwhile, the presence of noise degrades the performance of all algorithms, whether in term of rmsSAE or rmsSID. It is noticeable that though the values of all bands of the Desert Varnish are minimum (which can be regarded as dark object), GEM can extract it under all SNR levels. The next experiment was used to evaluate run times for GEM, VCA and MVHT. It was assumed that the number of bands and endmembers were equal in this experiment. Two sets of simulated data were generated. For the first one, p vectors with p dimensions (p ranged from 3 to 52) were generated randomly as endmembers with M = 1000; and for the second one, M (M ranged from 4000 to 20,000) pixels based on the LMM were generated with p = 3. In order to make sure that every endmember exists in the data set, the top p pixels are replaced by the p endmembers. GEM, VCA and MVHT can all extract endmembers correctly from the data set. Fig. 3 shows the computing time of GEM, VCA and MVHT as function of the number of endmembers (Fig. 3a) and pixels (Fig. 3b). We can see that compared to VCA and MVHT, GEM has obvious advantages in terms of run time. This was especially true when the number of endmembers is larger. For example, when the number of endmembers was greater than 20, the run time for GEM was significantly lower than it was for VCA and MVHT. Fig. 3b presents the relationship between the computational time and the number of pixels for the three methods when the number of the endmembers was fixed at 3. It is shown that the computational time of MVHT was much more than run times for VCA and GEM, while the computational speed of GEM is slightly better than that of VCA. The last experiment was used to evaluate the effect of DR on the performance of GEM. The simulated data in the first experiment was used here. Though in Algorithm 1, DR was performed before endmember extraction, GEM can run on data without DR. Theoretically, GEM can extract 4 endmembers for any 4 linearly independent bands of the original data. In this case, GEM may produce different results when starting on different bands. Fig. 4 shows the performance of GEM running on data with DR (using the first 4 PC bands) and without DR (using the first 4 bands of original
Table 1 Computational complexity of VCA, MVHT, and GEM. Algorithm
VCA
Complexity (flops)
P 2 2 3 2p2 M þ 2 pk¼1 ðkp þ 2k p þ k Þ
MVHT P 2ð p1 4p 4k þ 3ÞM k¼2
GEM p2 M þ p4
214
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
bands. From Fig. 4 we can see that for all SNR levels, the performance of GEM running on PC bands is always better than that on original data. It should be pointed out that GEM can extract all the p endmembers from any p linearly independent bands if the LMM holds and the data are noiseless. However, real hyperspectral scenes usually cannot fully satisfy the LMM and are easily corrupted by sensor noise. Therefore, to improve endmember extraction accuracy, GEM needs to perform DR first, which can improve the data SNR and at the same time remove the redundancies between bands. 3.2. Experiment with real data
Fig. 1. The four endmembers (Alunite, Calcite, Desert Varnish and Kaolinite) from USGS library used to generate the simulated data.
data). Since the first p PC bands retain much more information than any p bands of the original data, the results of GEM running on PC bands are usually better than that of GEM running on original
From the above experiments, we found that the computational complexity of VCA is much closer to that of GEM. As a result, in this section, a more detailed comparative analysis between GEM and VCA is presented with measured hyperspectral data, collected by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS) sensor over Cuprite, Nevada. AVIRIS is a high quality low noise imaging spectrometer, which acquires data in 224 contiguous spectral bands from 370 to 2510 nm (Green et al., 1998). The selected subscene is shown in Fig. 5 with a size of 190 250 pixels. Due to water vapor absorption or low SNR, bands 1–3, 107–113, 153–
Fig. 2. The performances of GEM, VCA and MVHT as function of SNR. (a) rmsSAE and (b) rmsSID.
Fig. 3. Computing time for VCA, MVHT and GEM. (a) M = 1000 and (b) p = 3.
215
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Fig. 4. The effect of PCA DR to the performance of GEM. (a) rmsSAE and (b) rmsSID.
Table 2 Spectral angle (in degrees) between USGS reference spectra and extracted endmembers by GEM and VCA in 3 different runs (VCA1 for the 1st run, VCA2 for the 2nd run, and VCA3 for the 3rd run).
Fig. 5. The band 8 (468.71 nm) of the real hyperspectral subscene.
169, 221–224 were removed. The resampled spectra from the USGS Digital Spectral Library included in ENVI software were selected as ground references for the comparison. According to the ground truth presented in reference (Swayze et al., 1992) and endmember extraction results in references (Liu and Zhang, 2012; Nascimento and Dias, 2005, 2007), there are less than 18 materials in this subscene. Thereby, the estimated number of endmembers was set to p = 18, and the dimensionality of the image was reduced to 18 by PCA. As GEM may find more than one endmember for one mineral, they are numbered successively, such as Alunite #1 and Alunite #2. For each mineral, the library spectrum, which has the minimum spectral angle degree (SAD) with the first endmember of that mineral, was selected as the reference for all the endmembers belonging to that mineral. Taking Alunite for example, there are 6 library spectra of Alunite in USGS spectral library, the one having the minimum SAD with Alunite #1 was chosen as the reference spectrum for both Alunite #1 and 2. Since VCA may produce different results for different runs, we ran VCA randomly for 3 times, noted as VCA1-3 in Table 2. Table 2 compares the SAD between extracted endmembers and USGS spectra quantitatively. We found that only 10 endmembers could be extracted by both the GEM and the 3 VCA runs. Compared to GEM, there were three endmembers missed by VCA1-3. It should be noted that GEM also missed some of the remaining VCA endmembers. For a fair comparison, the average SAD in Table 2 was calculated for the 10 common endmembers only. Within these 10 endmembers, 2 endmembers corresponded to the same minerals, so only 8 endmembers’ spectral curves are shown in Fig. 6.
Substance
GEM
VCA1
VCA2
VCA3
Alunite #1 Alunite #2 Andradite Buddingtonite Chalcedony #1 Chalcedony #2 Dumortierite #1 Dumortierite #2 Dumortierite #3 Kaolinite #1 Kaolinite #2 Montmorillonite Muscovite #1 Muscovite #2 Nontronite #1 Nontronite #2 Nontronite #3 Pyrope Average
4.3262 16.0481 5.7873 3.7080 5.4079 5.7812 9.4846 4.5170 9.3771 4.7599 9.8903 7.6854 3.9345 7.9540 3.5958 6.2838 4.8177 3.0012 6.8701
7.3982 18.5156 – 4.4668 – 6.1915 9.4846 – 8.6622 4.6720 – 5.7613 3.9345 – 3.5958 4.5381 – – 7.2682
6.2043 18.1011 – 5.0350 – 4.8733 9.4846 – 8.6622 4.9172 11.7588 6.8411 3.9345 – 3.6867 – – 5.4229 7.1740
15.0414 7.5699 6.8183 4.4668 – 5.9371 9.4846 – 8.6622 5.6708 11.1889 6.5376 3.9345 8.4721 5.5266 – – – 7.2832
From Table 2 we can see that GEM has a better performance for 3 endmembers (Alunite #1, Buddingtonite and Nontronite #1), an equal performance for 2 endmembers (Dumortierite #1 and Muscovite #1) and a worse performance for 2 endmembers (Montmorillonite and Dumortierite #3). For the remaining 3 endmembers (Alunite #2, Chalcedony #2 and Kaolinite #1), GEM performed much better than the worst VCA. Overall GEM produced a smaller average angle than VCA. Fig. 7 shows the abundance images of the 18 endmembers by the fully constrained least squares method in (Chang and Heinz, 2000). The abundance maps for the rest of the endmembers are in good agreement with the ground truth presented in reference(Swayze, et al., 1992) and the results presented in (Liu and Zhang, 2012; Nascimento and Dias, 2005, 2007). Fig. 8 shows the computing time of GEM and VCA versus p ranging from 10 to 18 using the real hyperspectral data. Similar to Fig. 3a, the advantage of GEM in computing time is more and more obvious as p increases.
4. Discussion and conclusions Based on the pure pixel assumption, we propose an endmember extraction method using Gaussian elimination method for
216
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Fig. 6. Comparison of endmember signatures extracted by (a) GEM, (b–d) VCA1-3, and (e) the corresponding laboratory spectra.
hyperspectral imagery, named Gaussian elimination method (GEM). For simulated data with p endmembers, GEM can extract
all the p endmembers from any p linearly independent bands based on LMM. However, real hyperspectral images usually do not satisfy
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
217
Fig. 7. Abundance maps of different endmembers from GEM. (a) Muscovite #1, (b) Alunite #1, (c) Dumortierite #1, (d) Kaolinite #1, (e) Montmorillonite, (f) Kaolinite #2, (g) Buddingtonite, (h) Nontronite #1, (i) Alunite #2, (j) Andradite, (k) Pyrope, (l) Nontronite #2, (m) Chalcedony #1, (n) Dumortierite #2, (o) Nontronite #3, (p) Muscovite #2, (q) Dumortierite #3 and (r) Chalcedony #2.
LMM and are easily corrupted by noise, so we need to perform DR before endmember extraction for better results. Besides, GEM assumes the existence of pure pixels for all the endmembers in an image. However, GEM might perform slightly worse than those endmember generation methods (Chan et al., 2009; Miao and Qi, 2007) if the pure-pixel assumption is not satisfied. In addition, it has to be pointed out that the performance of GEM was evaluated
on a relative small area of hyperspectral image. Both the performance and processing speed of GEM for large images will be studied in the future. Compared to the orthogonal matrix of VCA, the introducing of the Gaussian elimination matrix, which is a lower triangular matrix, can reduce the computational complexity of GEM a lot (referring to Table 1, Figs. 3 and 8). Another difference between
218
X. Geng et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 211–218
Fig. 8. Computing time (s) for VCA and GEM versus the number of endmembers, p when using real hyperspectral image.
GEM and VCA is that GEM uses the fixed projection vector for each endmembers, i.e. GEM always searches endmembers from the first PC band to the pth PC bands. Thereby, GEM can produce repeatable results. Moreover, the performance of GEM is comparable to that of VCA (see Fig. 2 and Table 2). Acknowledgements This work was subsidized by the National High Technology Research and Development Program (863 Program) of China (Nos. 2009AA122101, 2009AA12ZX1486531). Appendix A The proof of the GEM is given as follows. Theorem. Suppose Ri = [ri1, ri2, . . ., riM] represents all the M pixels in ith (1 6 i 6 L) band of a hyperspectral image, and there exists only one maximum absolute value in Ri and denoted its index as J, that is absðriJ Þ ¼ max ðabsðrij ÞÞ, then rJ ¼ ½r 1J ; r 2J ; ; r LJ T is one of the 16j6M
endmembers of the image. Proof. If rJ is not an endmember, it must be a linear combination P of the endmembers of the image, i.e. rJ ¼ pk¼1 ck ek , where Pp k¼1 ck ¼ 1, 0 6 ck 6 1, fek ; 1 6 k 6 pg are all the endmembers of the image. P In the ith band of image, we have riJ ¼ pk¼1 ck eik , which can be Pp transformed to 1 ¼ k¼1 ck erikiJ by dividing both sides with riJ . Since riJ is the only maximum absolute value, there is at least one e ik k0 ð1 6 k0 6 pÞ which meets the condition that abs riJ0 < 1. Then Pp P P we have 1 ¼ absð k¼1 ck erikiJ Þ 6 pk¼1 absðck ÞabsðerikiJ Þ < pk¼1 ck ¼ 1 which results in the contradiction. Therefore rJ ¼ ½r1J ; r 2J ; ; rLJ T must be one of the endmembers. h References Adams, J.B., Smith, M.O., Johnson, P.E., 1986. Spectral mixture modeling: a new analysis of rock and soil types at the Viking Lander I site. Journal of Geophysical Research 91, 8098–8812. Adams, J.B., Smith, M.O., Gillespie, A.R., 1989. Simple models for complex natural surfaces: a strategy for the hyperspectral era of remote sensing. In: Paper presented at the IGARSS’89, Vancouver, BC, Canada.
Asner, G.P., Heidebrecht, K.B., 2003. Imaging spectroscopy for desertification studies: comparing AVIRIS and EO-1 Hyperion in Argentina drylands. IEEE Transactions on Geoscience and Remote Sensing 41 (6), 1283–1296. Boardman, J.W., 1993. Automated spectral unmixing of AVIRIS data using convex geometry concepts. In: Paper presented at the 4th JPL Airborne Geoscience Workshop. Chan, T.-H., Chi, C.-Y., Huang, Y.-M., Ma, W.-K., 2009. A convex analysis-based minimum-volume enclosing simplex algorithm for hyperspectral ummixing. IEEE Transaction on Signal Processing 57, 4418–4432. Chang, C.I., Heinz, D.C., 2000. Constrained subpixel target detection for remotely sensed imagery. IEEE Transactions on Geoscience and Remote Sensing 38 (3), 1144–1159. Chein, I.C., Chao-Cheng, W., Wei-min, L., Yen-Chieh, O., 2006. A new growing method for simplex-based endmember extraction algorithm. IEEE Transactions on Geoscience and Remote Sensing 44 (10), 2804–2819. Clark, R.N., Swayze, G., Gallagher, A., King, T.V., Calvin, W.M., 1993. The U.S. geological survey digital spectral library: Version 1: 0.2 to 3.0 um, U.S. Geological Survey, Washington, DC, Open File Rep. 93–592. Combe, J.P., Le Mouélic, S., Sotin, C., Gendrin, A., Mustard, J.F., Le Deit, L., et al., 2008. Analysis of OMEGA/Mars Express data hyperspectral data using a MultipleEndmember Linear Spectral Unmixing Model (MELSUM): Methodology and first results. Planetary and Space Science 56 (7), 951–975. http://dx.doi.org/10.1016/ j.pss.2007.12.007. Full, W.E., R, E., Bezdek, J.C., 1982. Fuzzy Qmodel- a new model approach for linear mixing. Mathematical Geology 14 (3), 259–270. Gill, T.K., Phinn, S.R., 2009. Improvements to ASTER-derived fractional estimates of bare ground in a Savanna Rangeland. IEEE Transactions on Geoscience and Remote Sensing 47 (2), 662–670. Gong, P., Miller, J.R., Spanner, M., 1994. Forest canopy closure from classification and spectral unmixing of scene components-multisensor evaluation of an open canopy. IEEE Transactions on Geoscience and Remote Sensing 32 (5), 1067– 1080. Green, R.O., Eastwood, M.L., Sarture, C.M., Chrien, T.G., Aronsson, M., Chippendale, B.J., et al., 1998. Imaging spectroscopy and the airborne visible/infrared imaging spectrometer (AVIRIS). Remote Sensing of Environment 65, 227–248. Liu, J., Zhang, J., 2012. A new maximum simplex volume method based on householder transformation for endmember extraction. IEEE Transactions on Geoscience and Remote Sensing 50 (1), 104–118. Martin, G., Plaza, A., 2011. Region-based spatial preprocessing for endmember extraction and spectral unmixing. IEEE Geoscience and Remote Sensing Letters 99, 745–749. Miao, L., Qi, H., 2007. Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization. IEEE Transactions on Geoscience and Remote Sensing 45 (3), 765–777. Nascimento, J.M.P., Dias, J.M.B., 2003. Vertex component analysis: a fast algorithm to extract endmembers spectra from hyperspectral data. In: Paper presented at the Pattern Recognition and Image Analysis: First Iberian Conference. Nascimento, J.M.P., Dias, J.M.B., 2005. Vertex component analysis: a fast algorithm to unmix hyperspectral data. IEEE Transactions on Geoscience and Remote Sensing 43 (4), 898–910. Nascimento, J.M.P., Dias, J.M.B., 2007. Hyperspectral Signal Subspace Estimation. In: Paper presented at the Proceedings of IGARSS. Neville, R.A., Staenz, K., Szeredi, T., Lefebvre, J., Hauff, P., 1999. Automatic endmember extraction from hyperspectral data for mineral exploration. In: Paper presented at the Proceedings of 21st Canadian Symposim Remote Sensing. Roberts, D.A., Gardmer, M., Church, R., Ustin, S., Scheer, G., Green, R.O., 1998. Mapping chaparral in the Santa Monica mountains using multiple endmember spectral mixture models. Remote Sensing of Environment 65, 267–279. Shaohui, M., Mingyi, H., Zhiyong, W., Dagan, F., 2010. Spatial purity based endmember extraction for spectral mixture analysis. IEEE Transactions on Geoscience and Remote Sensing 48 (9), 3434–3445. Smith, M.O., Ustin, S.L., Adams, J.B., Gillespie, A.R., 1990. Vegetation in deserts: I. A regional measure of abundance from multispectral images. Remote Sensing of Environment 31, 1–52. Somers, B., Asner, G.P., Tits, L., Coppin, P., 2011. Endmember variability in spectral mixture analysis: a review. Remote Sensing of Environment 115, 1603–1616. Swayze, G., Clark, R., Sutley, S., Gallagher, A., 1992. Ground-truthing AVIRIS mineral mapping at cuprite, Navada. In: Paper presented at the Proceedings of Summaries Third ANNu. JPL Airborne Geoscience Workshop. Tao, X., Wang, B., Zhang, L., 2009. Orthogonal bases approach for the decomposition of mixed pixels in hyperspectral imagery. IEEE Geoscience and Remote Sensing Letters 6 (2), 219–223. Winter, M.E., 1999. Fast autonomous spectral endmember determination in hyperspectral data. In: Paper Presented at the Proceedings of SPIE. Xia, W., Wang, B., Zhang, L., Lu, Q., 2011. Simplex volume analysis based on triangular factorization: A framework for hyperspectral for hyperspectral unmixing. In: Paper presented at the Geoscience and Remote Sensing Symposium (IGARSS), IEEE International. Zare, A., Gader, P., 2010. PCE: piecewise convex endmember detection. IEEE Transactions on Geoscience and Remote Sensing 48 (6), 2620–2632. Zhang, J.K., Rivard, B., Rogge, D.M., 2008. The successive projection algorithm (SPA), an algorithm with a spatial constraint for the automatic search of endmembers in hyperspectral data. Sensors 8, 1321–1342.