Support Vector Machine-Based Endmember ... - Semantic Scholar

2 downloads 570 Views 2MB Size Report
Sep 22, 2009 - remote sensing, support vector machines (SVMs). I. INTRODUCTION. HYPERSPECTRAL remote sensors acquire radiant flux data in typically ...
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

771

Support Vector Machine-Based Endmember Extraction Anthony M. Filippi and Rick Archibald

Abstract—Introduced in this paper is the utilization of support vector machines (SVMs) to semiautomatically perform endmember extraction from hyperspectral data. The strengths of SVM are exploited to provide a fast and accurate calculated representation of high-dimensional data sets that may consist of multiple distributions. Once this representation is computed, the number of distributions can be determined without prior knowledge. For each distribution, an optimal transform can be determined that preserves informational content while reducing the data dimensionality and, hence, the computational cost. Finally, endmember extraction for the whole data set is accomplished. Results indicate that this SVM-based endmember extraction algorithm has the capability of semiautonomously determining endmembers from multiple clusters with computational speed and accuracy while maintaining a robust tolerance to noise. Index Terms—Endmember extraction, hyperspectral imaging, remote sensing, support vector machines (SVMs).

I. I NTRODUCTION

H

YPERSPECTRAL remote sensors acquire radiant flux data in typically hundreds of narrow contiguous spectral bands, yielding a continuous signature for each pixel in the image [1]. Such data can be exploited to detect materials that cannot be regularly distinguished via multispectral cameras [2]. Endmember extraction algorithms can be applied to multispectral and hyperspectral remote sensor images to facilitate the estimation of the subpixel fractional abundance of a given material [3]. An endmember is an idealized pure signature for a given class. In general, remote-sensing image pixels tend to be mixed rather than pure [4]; sensor noise and withinclass reflectance variance translates into endmembers being a conceptual convenience in real imagery [2], [5]. An endmember may thus represent only one material (i.e., a pure endmember), or it could contain a large proportion of a material, combined with a distinctive mixture of other materials. Also, note that there are usually many more endmembers than materials within Manuscript received February 15, 2008; revised June 30, 2008. First published December 9, 2008; current version published February 19, 2009. This work was supported in part by an appointment to the U.S. Department of Energy (DOE) Higher Education Research Experiences (HERE) for Faculty at the Oak Ridge National Laboratory (ORNL) administered by the Oak Ridge Institute for Science and Education. The work of R. Archibald was supported by the Householder Fellowship that is supported under the Mathematical, Information, and Computational Sciences Division, Office of Advanced Scientific Computing Research, U.S. Department of Energy under Grant DE-AC05-00OR22725. A. M. Filippi is with the Department of Geography, Texas A&M University, College Station, TX 77843-3147 USA (e-mail: [email protected]). R. Archibald is with the Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 USA (e-mail: archibaldrk@ ornl.gov). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2008.2004708

a given scene. That is, for any given material type, exemplars may exist in an image corresponding to extremes of the condition of that material (e.g., with respect to shadowing or solar illumination, weathering, pigmentation, etc.). The net effect is a potentially large number of endmembers for a single material type [6]. This was addressed in [7], where endmember variability was integrated into spectral mixture analysis via the representation of each endmember by an endmember bundle, or set of spectra. A simulated annealing algorithm was used to derive endmembers as vertices of a simplex. Endmembers can be extracted via manually driven methods or semiautomated/automated algorithms. Manual/interactive endmember identification methods (e.g., the pixel purity index (PPI) algorithm [8], [9] and the manual endmember selection method [10]) are often slow if significant human intervention in the endmember selection process is required. In addition, the final set of identified endmembers may be relatively human analyst dependent. Therefore, a limited number of semiautonomous/autonomous endmember extraction algorithms have been proposed to date. For instance, Chang and Plaza [11] noted a variety of shortcomings of PPI (e.g., no guarantee of true endmembers, no guidance for retention of endmembers, noise sensitivity, no criteria for how to select certain parameters, requirement of human intervention, computational complexity, etc.) and thus developed an unsupervised fast iterative PPI algorithm to address such issues. Other such automated or semiautomated methods include the optical real-time adaptive spectral identification system (ORASIS) [12], [13], convex cone analysis (CCA) [14], iterative error analysis (IEA) [15], N-FINDR [16], [17], the simulated annealing algorithm of Bateson et al. [7] (noted previously), automated morphological endmember extraction (AMEE) [18], iterated constrained endmembers (ICE) [19], the sequential maximum angle convex cone (SMACC) endmember model [6], vertex component analysis (VCA) [20], and the single individual evolutionary (SIE) strategy proposed in [21]. For example, ORASIS utilizes a modified Gram–Schmidt code to factor the data matrix, and a shrink-wrapping algorithm finds an outer simplex [13]. CCA is rooted in the notion that radiance is a nonnegative quantity; vectors composed of such spectra are linear combinations of nonnegative elements that lie inside a nonnegative convex region. CCA seeks to locate the region boundary points, which can serve as endmembers [14]. IEA locates endmembers via iterative unmixing. Ultimately, IEA will iterate until a predetermined quantity of endmembers is located [15], [22]. AMEE exploits the spatial correlation between pixels and jointly considers spatial and spectral information; it applies mathematical morphology [23] to the spectral domain while maintaining the

0196-2892/$25.00 © 2008 IEEE

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

772

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

spatial characteristics. AMEE is based on the following two parameters: a minimum and a maximum spatial kernel size [18]. ICE is a statistical method of endmember determination that employs regularization, and it is an iterative solution to a nonlinear continuous parameter optimization problem. ICE will terminate after a finite, although unknown, number of iterations [19]. SMACC uses a convex cone model to represent vector data; this convex factorization method locates extreme vectors and employs them as endmembers. SMACC begins with a single endmember and sequentially/incrementally increases in dimension. Each new endmember is determined by the angle made with the existing cone [6]. VCA is an unsupervised endmember extraction method that is based upon the notion that endmembers are simplex vertices and that the affine transformation of a simplex is a simplex as well [20]. Reference [21] used an SIE strategy for endmember search, where an endmember set is considered to be an evolving population. N-FINDR was employed in this research as a baseline algorithm, as it is commonly in other endmember studies, and will be briefly described in the following discussion. Although various semiautonomous/autonomous endmember extraction methods currently exist, the variety is still quite limited [18], and new approaches with unique properties and the potential for improvement should be considered. Support vector machines (SVMs), introduced and developed by Boser et al. [24] and Vapnik [25], have proven to be a powerful method in solving supervised classification and regression problems. The simple, yet elegant, idea behind SVM is to find the optimal hyperplane that maximizes the margin between the different classes. Kernels are used as computationally efficient nonlinear transforms within SVM, thereby allowing complex decision surfaces to be generated. Thus, SVM is considered to be a kernel-based classifier method. The hyperspectral community has embraced the SVM paradigm to solve challenging problems [26]. Reference [27] has utilized SVM sensitivity analysis for data-driven knowledge discovery in order to provide a ranking of input band significance. Reference [28] developed composite kernels to integrate spectral knowledge into SVM classification. A complete theoretical and experimental comparison of SVM and other kernelbased approaches, specifically geared toward hyperspectral data classification, was performed by [29]. Reference [30] utilized the structure of SVM to derive optimal training sets. This brief summary only represents a sampling of the most recent research on this topic. Although SVMs have gained much popularity for classification problems, they have only been applied to the problem of endmember extraction in an extremely limited sense. The constrained least squares linear spectral mixture model (LSMM) has been found to be related to the linear SVM (LSVM) and to be equivalent, given certain circumstances (i.e., when the same information is employed in their designs, and where the same thresholding and normalization process are used). This property has led to the proposed use of linear SVM for “pure pixel” selection in [31] and [32] (reference [33] offered a further proof for the equivalence of LSVM and LSMM and extended the work of Brown et al. [31] and [32] by combining LSVM and LSMM to produce a double-unmixing algorithm). It was

also noted that SVMs can handle instances where pure-pixel sets for various classes overlap, as well as where nonlinear mixing regions exist [32]. The other support vector approach [34] utilized in this domain is support vector data description (SVDD), which is related to SVM and was developed for oneclass classification [35], [36]. SVDD enables nonparametric modeling of the distribution of spectra and the determination of convex boundaries in high-dimensional space. The SVDD endmember approach in [34] treats the spectra that lie on the convex boundaries, the support vectors, as endmembers. Support vector clustering (SVC) [37] is a clustering method that can be seen as an extension of SVDD, with the added capability, similar to the method proposed in this paper, of detecting multiple clearly separated distributions. To our knowledge, however, SVC has not been applied to the problem of endmember extraction. Common among these methods based on or related to SVM is their requirement to use the entire data set in a single support vector minimization problem. The computational cost of the support vector minimization problem grows as the square of the number of samples. For large data sets, applying a single SVM to the entire data set would be computationally prohibitive. In contrast, the SVM-based method presented in this paper is an iterative procedure that is based on small subdata sets, and can be used on large data sets; training multiple SVMs is among the key properties that clearly distinguish our proposed method from prior support vector approaches. Furthermore, the proposed algorithm has an architecture that is amenable to the parallel environment of high-performance computing (HPC). This paper introduces a support vector machine-based endmember extraction (SVM-BEE) method. This method uses the rich source of information contained within SVM in a variety of different modes. First, a decision surface that approximates the convex hull of each separated distribution is iteratively generated through multiple applications of SVM on small subsets of data. The support vectors that determine the decision surface are, by design, close to the convex hull and provide the necessary information to extract all endmembers. Second, the number of separate distributions or clusters is found directly from testing the connectivity of all the support vectors with the decision function. Third, a dimension-reducing transform is employed for each cluster by finding the lowest dimension hyperplane that minimizes the total distance to the corresponding support vectors. Finally, the endmembers are determined by finding the simplex that minimizes the total distance to the transformed support vectors. SVM-BEE has many advantages that are inherited from SVM. The computational cost of the method is minimal since all optimization is done on a small sample size. High accuracy is achieved because the data used in calculation, or support vectors, provide an ideal representation of the convex hull. SVM-BEE can determine distributions of various shapes by exploiting support vector information, for multiple clusters as necessary, from its decision surface algorithm. Using these interesting points, we focus upon extracting simplex topology in this paper, but SVM-BEE actually has no limitation with respect to the shape of the distribution to which it fits. This is a very advantageous property that lends the method

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

flexibility and limits the potential for neglecting to extract a given endmember. Finally, the method demonstrates consistency and robust noise tolerance. The outline of this paper is as follows. Section II provides a brief introduction to SVM in order to define the notation necessary to explain SVM-BEE. Section III is an explanation of the SVM-BEE method with specific details of each stage. Section IV presents numerical results for SVM-BEE on both synthetic and real data and compares this method to the wellknown N-FINDR algorithm [16], [17], as well as the SVDD method [34]. We end with a discussion and conclusion in Section V. Here, we define some notation that will be consistent for the remainder of this paper. The notation x is the standard distance norm for a vector. The notation pij is the ith band of the jth pixel, whereas pj is the vector of all bands for the jth pixel.

We present a brief introduction to SVMs with the intention to define the notation necessary to explain SVM-BEE. For a more detailed explanation of SVMs, we refer the interested reader to [25]. Within the context of hyperspectral imaging, suppose that NS image pixels, each containing NB spectral bands, can be separated in two classes, notated as lj ∈ {−1, 1},

j = 1, . . . , NS .

(1)

Therefore, each of the sampled pixels belongs either to what we will call in this paper the outer class Sout = {pj : lj = −1, j = 1, . . . , NS }

(2)

or the inner class Sin = {pj : lj = 1, j = 1, . . . , NS }.

(3)

Given this setup, the solution to the SVM minimization problem that maximizes the margin between the two transformed classes is found by solving ⎡ ⎤ NS NS   1 min ⎣− αj + lj lk αj αk K(pj , pk )⎦ (4) α 2 j jk

with NS 

lj αj = 0,

Once the solutions of (4) and (5) are obtained, any pixel p ∈ RNb can be classified by the sign of the decision function ⎞ ⎛ NS  lj αj K(pj , p) + b⎠ (7) f (p) = sign ⎝ j=1

where b=

NS NS 1  1  lj − lk αj K(pj , pk ). NS j=1 NS

(8)

j,k=1

Here, the support vectors are defined by {pj : αj = 0, j = 1, . . . , NS }. Additionally, the decision function defines a complex decision surface as Ω = {p : f (p) = 0}

(9)

which is the hyperplane in the range space of the kernel that separates the transformed training data set into two classes.

II. SVM S

pj ∈ RNb

773

0 ≤ αj ≤ C

(5)

j

for some kernel function K and penalty parameter C. Since the following algorithms are designed to completely separate the two classes, the penalty parameter can be considered to be infinite. It is noted that, for this paper, only the radial basis function kernel is used and is defined as 2

Kσ (x, y) = e−x−y/2σ .

(6)

III. SVM-BEE A LGORITHM The SVM-BEE algorithm consists of three parts. First, a decision surface is iteratively generated that approximates the convex hull for every cluster within the entire range of sampled data. The support vectors that define the final decision surface contain a rich source of information by the fact that they are a small subset of the data that collectively approximate each convex hull of the whole data set and are individually located near a particular cluster distribution edge. These support vectors can be used in the second part of the SVM-BEE algorithm to determine the number of clusters within the data set and the dimension-reducing transform that describes each of the clusters. Finally, endmembers can be extracted for each cluster by determining the simplex that minimizes the distance to the corresponding edge pixels (which are support vectors identified to be on the distribution edge). A. Decision Surface-Generation Algorithm Suppose that we are given a hyperspectral data set, denoted as S of NS pixels, and NB bands. A complex decision surface can be generated using SVM as follows: 1) Randomly pick a small subset Sin ⊂ S of N  NS pixels of the hyperspectral data set. 2) Set the initial decision function f (p) = r − p − c with c = pj∈Sin pj /N and r = 2 maxj∈Sin pj − c. Therefore, the initial decision surface Ω0 is the surface of the NB dimensional ball with center c and radius r. 3) Sout is a set of 2NB N pseudopoints generated from Sin and designed to always be outside the convex hull of each cluster of the hyperspectral data set. Specifically Sout = {pj + μi Df (pj )} for pj ∈ Sin , i = 1, . . . , 2NB , such that μi is the solution of the equation μi = min{μ|f (pj + μDf (pj )) = ξ}

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

774

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

with unit directional derivative Df (pj ) that is rotated by 90◦ in 2NB mutually orthogonal directions and limited in the sense μi , μi ≥ ξ/2 μi = ξ, μi < ξ/2.

4)

5)

6)

7) 8)

Note that it is always possible to algebraically determine the directional derivative from the decision function generated by SVM. Implicitly determine the next decision surface Ω as described in (9) by solving the SVM minimization problem using Sin , Sout , and Kσ as defined in (6). For every pixel in the data set, evaluate the decision function (7). If there exist pixels that have a decision function value less than zero, then redefine Sin with up to N random points with a decision function value less than zero and the support vectors identified in the previous step, and repeat step 4. Note that picking random points with function values less than zero provides a scattered distribution of points, whereas picking points with the lowest decision function value has a tendency to condense Sin and hinder convergence. Add N random points to Sin from S such that their decision function value is in the first quartile (i.e., the current decision function is used to pick points that are close to the estimated distribution boundary). 0 = Sout and update Sout as outlined in step 3) Set Sout with the decision function found in step 5). If the decision surface movement is below a certain threshold, measured by maxx∈Ω miny∈Ω0 x − y < ξ, then stop. Otherwise, go to step 4). Note that, in practice, the change in the decision surface is approximated by the change in the updated set of pseudopoints, maxx∈Ω miny∈Ω0 x − y ∼ 0 x − y. maxx∈Sout miny∈Sout

The only parameters necessary other than the lower limit of the subsample size are σ, which controls the amount of deformation in the decision surface, and ξ, which controls how tightly the decision surface will conform to the underlying sample distribution. Both of these parameters are linked and can be estimated by the average distribution density as 2σ = ξ =

  1 pj − x N NS

(10)

pj ∈Sin x∈S

for σ given in (6). The default for the lower limit of the subsample size is N = 4NB . It is important for computational efficiency to keep this parameter low, and it is further noted that this parameter does not effect the final estimation of the distribution edge since the subsample size will naturally increase with the number of support vectors to approximate complicated shapes. If the sample is a fusion of multiple nonoverlapping distributions, the decision surface is designed to isolate the convex hull of each distribution. Thus, the number of separate distributions is determined directly from the data and need not be known

a priori. However, if the distributions are overlapping, then SVM-BEE will only isolate a set of pixels that are on the edges of the combined distributions. In the future, we plan to research algorithms that can detect both the type and number of distributions from the information extracted by SVM-BEE about the distribution boundary. Currently, the method is only designed to detect distributions based on linear mixing of endmembers. Although noise may be present, at this stage, it is not necessary to mitigate its effect. The goal of this process is to find the minimal number of points that is needed to accurately describe the convex hull of each possible distribution contained in the sample. B. Cluster Approximation and Dimension-Reducing Transform The decision function (7) and associated support vectors that are a product of the generation of the decision surface are the necessary input needed to determine both the number of clusters in the hyperspectral data set and the corresponding dimension-reducing transform. The number of clusters is determined by the connectivity of the support vectors. We say that two support vectors are connected if there exists a continuous path between these two points such that the decision function is never negative. The numerical testing of connectivity is simplified because the support vectors approximate the convex hull. Therefore, it is sufficient to test only the linear path between any two support vectors and use the transitive property to completely identify all connectivity. Although the decision function is known analytically, testing connectivity requires that the decision function be evaluated along the linear path between any two support vectors. However, since the parameter ξ controls the precision at which the convex hull is estimated, it is only required that the linear path be evaluated at uniform points that are a distance of ξ/2 apart. A cluster is defined as a set of support vectors that are connected to each other. The entire data set can be clustered based upon the clustering of the support vectors. If there is only one cluster of support vectors, then it is not necessary to cluster the entire data set. Otherwise, for each cluster, split the support vectors into Sin and Sout according to membership of the particular cluster tested. Apply the resultant decision function to the data set, and identify all points with a positive decision function value to be a member of the particular cluster tested. The dimension-reducing transform is performed on each identified cluster of support vectors. There are many different types of dimension-reducing transforms that could be used, but to keep the continuity of the utilization of SVM, least squares fitting of hyperplanes is utilized. Specifically, cluster points are projected onto the lowest order dimensional hyperplane in which the projected and actual points have a least squares error less than the parameter σ (an estimate of the noise in the signal). This hyperplane is found using PCA decomposition. Suppose that the eigenvalues of the PCA decomposition for one particular cluster, sorted in descending order, are notated

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

as {λ1 , λ2 , . . . , λNB }. The total variance described by the mth eigenvector is defined by the following function: Λ(m) = λm

N B 

|λi |.

(11)

i=1

The appropriate dimension-reducing transform can be determined by Λ using the scree test, first proposed by Cattell [38]. Suppose the function Λ begins to level off at the point (mI , Λ(mI )), then the dimension-reducing transform will be the projection onto the first mI eigenvectors of the PCA decomposition. We recognize that more sophisticated dimension-reducing transforms exist (e.g., [2]), and developing a more robust dimension-reducing method for SVM-BEE will be a topic of future investigation. It is demonstrated for these experiments, however, that the SVM-BEE method performs adequately with the PCA transformation and is stable when scaling the dimension (i.e., endmember estimation is not volatile to reasonable changes in the dimension size). The use of PCA here represents a simple baseline case for the dimension-reducing transform step. Note that after dimension reduction is performed, SVM-BEE is designed to search for the number of endmembers that fits the linear mixing model for this reduced dimensional space. Although dimensionality reduction is not specifically needed for SVM-BEE, it has the advantage of reducing the computational effort, and it also fixes the number of endmembers that will be found. Part of our future work will involve adapting SVM-BEE to search for other mixing models, and we will use the existing information determined by the decision surface algorithm to provide flexibility in the number of endmembers that is sought after. C. Simplex Approximation SVM-BEE is capable of determining distributions of various shapes. Here, we describe how the information contained within the decision surface algorithm can be utilized to extract simplex topology from hyperspectral data. Assuming the well-studied linear mixture model of hyperspectral data, any given pixel within the data set is a linear combination of endmember spectra pij =



eik ckj + ε

(12)

k

where pij is the ith band of the jth pixel, eik is the ith band of the kth endmember, ckj is the mixing proportion for the jth pixel from the kth endmember, and ε is the Gaussian random error that is assumed to be relatively small. Since ckj is the mixing proportion for the jth pixel, we must have ckj ≥ 0



ckj = 1

(13)

k

for k = 1, . . . , Ne and j = 1, . . . , NS . Here, Ne is the number of endmembers and is determined by the dimension (11) of the reducing transform as Ne = mI .

775

The decision surface algorithm is a data-driven iterative procedure that can be used to find an excellent subset of the hyperspectral data for approximating endmembers. In fact, if pure endmembers are contained in the hyperspectral data set, then, for sensitive ξ (control parameter for convex hull estimation), these pure endmembers will also be contained within the support vector subset. Even if pure endmembers are not included, it is possible to estimate them from dimensionally reduced support vectors. Within the SVM-BEE algorithm, the decision surface-generation occurs multiple times. The first application of the decision surface algorithm is to determine the number of different clusters so that the appropriate dimensionreducing transform can be done within each cluster. The second application of the decision surface algorithm is done on each dimension-reduced cluster with a slight modification to determine the pixels on the edges of the cluster. As the decision surface algorithm is executed, any pixel in step 3) that satisfies μ < ξ/2, for any of the tested directions, is classified as an edge pixel. The set of all edge pixels identified during the execution of the decision surface algorithm is what is used to approximate the simplex geometry and, hence, the endmembers as explained in the following steps. 1) Determine the initial set of Ne endmembers. a) e1 is the edge pixel that is farthest from the center of mass of all the edge pixels. b) Shift the origin by subtracting the first endmember e1 from the data set. c) For j = 2 : Ne − 1 i) ej is the edge pixel that is farthest from the origin. ii) Normalize the data set by |ej |. iii) Rotate the data set with matrix multiplication by e −1 ), where columns are mutuR = (e1j e2j , . . . , eN j ally orthogonal with the first column set as e1j = ej . iv) Remove the first dimension from all the edge pixels, and set j = 2 : Ne − 1. d) eNe is the edge pixel that is farthest from the origin. It is noted that if the edge pixels follow a linear mixing model of (12) with ε = 0, then this step will exactly identify the true endmembers (assuming that the true endmembers are contained within the set of edge pixels). 2) Group all of the edge pixels according to the closest hyperplane. 3) Update each hyperplane by finding the least squares hyperplane solution to each of the associated groups of edge pixels. Note that it is necessary that there be at least Ne edge pixels in each group. This requirement can be satisfied in all cases by including all but the jth endmember to the groups of associated edge pixels. 4) Update each endmember by adding all vectors that project that endmember to all associated hyperplanes. 5) Repeat steps 2)–4) until endmember update is negligible. The simplex of the resulting endmembers of this algorithm is that which minimizes the total distance to each edge pixel. The vertices of this simplex are the estimate of the endmembers.

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

776

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

D. Order of Computational Cost The SVM-BEE algorithm almost exclusively uses the support vectors in the estimation of endmembers. The only instance where the SVM-BEE algorithm uses the entire data set is in step 5) of the decision surface-generation algorithm. Since the number of pixels in a hyperspectral image NS is significantly more than the number of bands NB or the number of support vectors NSup , step 5) of the decision surface-generation algorithm dominates the computational order of the SVM-BEE algorithm. Therefore, a reasonable estimate for the computational order of SVM-BEE is O(Nit NSup NB NS ), where Nit is the relatively small number of iterations of step 5) during the execution of the decision surface-generation algorithm. This can be compared to the computationally dominant step of N-FINDR [39], the repeated calculation of different simplex volumes, which has a computational order of O(N it NB (NB + 1)3 NS ), where N it is the number of iterations of the N-FINDR algorithm. It can be seen that both of these methods will benefit by reducing NB and NS through dimension-reducing transforms and effective sampling of the hyperspectral image. A note regarding an SVDD computational/memory problem will be provided operationally in Section IV-D and, in more detail, Section IV-E. IV. E XPERIMENTAL E VALUATION OF THE SVM-BEE A LGORITHM The SVM-BEE algorithm is first presented on a simple 2-D problem so that the details of the procedure may be better understood. The algorithm is then tested against the N-FINDR [16], [17] and SVDD [34] algorithms on synthetic data. SVMBEE is also tested against N-FINDR using a real data set. Throughout this paper, the parameters used in the SVM-BEE algorithm are defined by the data according to (10) and (11). Additionally, all data sets are normalized so that the average Euclidean norm of the whole set is unity. This preprocessing step does not change the characteristics of the data and is beneficial to the stability of SVM and, hence, SVM-BEE. A. Example: Multiple Simplexes Suppose that a data set is generated from two separated sets of endmembers given as e1 = (1, 1)

e2 = (1.5, 2)

e3 = (2.5, 1.25)

e˜1 = (3, 3)

e˜2 = (4.75, 4.25)

e˜3 = (4.75, 3.5).

(14)

Fig. 1 shows specific steps of the decision surface-generation algorithm for a data set of size NS = 3000, based on the linear combination endmember spectra model of (12) and (13) for the two sets of endmembers given in (14). The initial conditions of the decision surface-generation algorithm shown in Fig. 1(a) display the small set of points used to generate a decision surface that encloses the whole data set. As the algorithm progresses, the decision surface eventually morphs to encircle each separate distribution of endmembers. It can be seen in Fig. 1(c) that the final decision surface provides

Fig. 1. (a) Initial conditions of the decision surface-generation algorithm. Ten points are picked from the total sampled points. The initial decision surface Ω0 is the circle, and (green) each sampled point generates (yellow) four pseudopoints, two in either direction of the radial projection of the sampled point out to the decision surface and two in either direction perpendicular to the radial projection. It is noted that extreme points or points farthest from the center of mass of sample are purposefully chosen as initial conditions. (b) Decision surface and resulting pseudopoints after the first iteration of the algorithm. (c) Final iteration of the algorithm.

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

777

Fig. 2. (a) SVM decision surface with the support vectors for the first cluster defined to be Sin , with the remaining support vectors defined as Sout . Evaluation of the decision function for each data point will return a positive value for each member of the first simplex. (b) Similarly, when the procedure is repeated for the second cluster of support vectors, the membership of the second simplex can be identified.

Fig. 3. (a) Edge pixels identified during the execution of the decision surface algorithm for each cluster. (b) For each cluster, the simplex that minimizes the distance to the edge pixels identifies the endmembers.

sufficient information to determine the connectivity of the final set of support vectors. In this case, the support vectors have two clusters, and this reduced set has the necessary information to cluster the entire data set, as shown in Fig. 2. Once the data set is clustered, the hyperplane that minimizes the total distance of each cluster is determined. It is found that the dimension of each hyperplane is the same as the original data, and therefore, no dimension-reducing transform is necessary. Next, the decision surface algorithm is applied separately to each cluster in order to determine the pixels on the edges of each distribution. By construction, the support vectors of the decision surface algorithm are often at the edges of the distribution. During the execution of the algorithm, edge pixels are identified when a point has solution μ < ξ/2 for at least one search direction in step 3). The set of all identified edge pixels for each cluster is shown in Fig. 3(a), where it can be seen that these pixels approximate well the convex hull of each cluster. Finally, as shown in Fig. 3(b), finding the simplex for each cluster that minimizes the total distance to the edge pixels identifies the endmembers.

In this example, 2(NB + 1) endmembers were properly identified. For other endmember extraction algorithms, such as N-FINDR, the upper limit on the number of endmembers that can be identified is NB + 1. The SVM-BEE algorithm, as demonstrated in this example, does not share the same limitation. The upper limit for SVM-BEE is NC (NB +1), where NC is the number of separate distributions or clusters. B. Example: Noise Tolerance The robust nature of SVM-BEE is demonstrated by the problem shown in Fig. 4. Here, the first set of endmembers in (14) is used to generate a linear combination of NS = 3000 points, according to the model given (10) and (11), with a signal-to-noise ratio (SNR) of 20 dB. Fig. 4(b) shows the edge pixels identified by SVM-BEE. The outer simplex in Fig. 4(c) minimizes the total distance to the edge pixels. However, as expected from the noise, this outer simplex overestimates the true endmembers. For this case, the simplex that minimizes the total distance to the edge pixels is a biased estimator of the endmembers. The correction of this bias is determined by

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

778

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

moving the simplex inward according to the SNR. The support vectors are extreme points of the distribution, and therefore, it can be assumed that these points will lie an average distance of δ=

RS SNR/20 10

(15)

away from the true simplex. Here, the root mean square of the signal is notated as RS , and the SNR is measured in decibels (dB). The inner simplex in Fig. 4(a) depicts the bias correction and is a much better approximation to the true endmembers. C. N-FINDR Overview In this paper, the performance of SVM-BEE was compared to that of the N-FINDR algorithm (version 3.1) [16], [17] for a baseline comparison. N-FINDR “inflates” a simplex inside the data set until the volume of the simplex is maximized. N-FINDR assumes that the true endmembers are contained within the data set, and therefore, the simplex “inflates” by choosing vertices from the data set. Note that N-FINDR and other endmember algorithms (e.g., PPI) assume that, for each endmember, there exists at least one pure pixel within the image, which is an assumption that may not be valid for some scenes. Such algorithms thus determine the purest set of pixels in the image data [20]. Reference [19] found that N-FINDR and PPI endmember extraction performances were particularly problematic when pure endmembers were absent in an image. N-FINDR also essentially assumes that the number of endmembers present in an image is known, which is rarely the case [19]. Neither of these requirements applies to SVM-BEE. Two user-defined parameters that N-FINDR requires are as follows: 1) the maximum number of endmembers that N-FINDR can identify and 2) an estimate of the sensor SNR. It has been suggested that the maximum number of endmembers parameter (maxEM) be set higher than the maximum probable number of endmembers in the image. Setting this parameter too high though can trigger an N-FINDR failure or a splitting of genuine endmembers. It has also been recommended that average SNR be set high, at least at first. The number of endmembers that N-FINDR will extract increases as the SNR increases [40]. Note that since its original publication, N-FINDR has been undergoing some proprietary modifications pertaining to optimizations of the code for speed and accuracy [39], [M. Winter, University of Hawai’i at Manoa, Honolulu, HI, private communication, May 2006]. SVM-BEE, N-FINDR, and SVDD experiments were conducted using synthetic hyperspectral image data, followed by a real Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) hyperspectral image experiment, discussed in the following sections. D. Example: Synthetic Hyperspectral Image Test #1 Fig. 4. (a) Initial conditions of the decision surface-generation algorithm for a set of 3000 points based on the linear model given in (10) and (11), with endmembers given in the first set of (14), where the SNR is 20 dB. Here, the green dots are sampled points determined to be support vectors, and the yellow dots are pseudopoints. (b) Final edge pixels determined by SVM-BEE. (c) Outer simplex minimizes the total distance to the edge pixels, and the inner simplex is the bias correction. It can be seen that the inner simplex is a much better approximation to the true endmembers, which are the red dots.

Laboratory-measured spectral library data were used to generate a synthetic hyperspectral image. We refer to this first effort as synthetic hyperspectral image test #1. Sample reflectance spectra for human-made materials and soils were taken from the standard Johns Hopkins University (JHU) Spectral Library [Jet Propulsion Laboratory Spectral Library. URL:

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

Fig. 5.

779

Endmember spectra used for pseudoimage data set generation. TABLE I RELATIVE ERROR OF THE N-FINDR- AND SVDD-EXTRACTED AND SVM-BEE-ESTIMATED ENDMEMBERS TO TRUE ENDMEMBERS FROM SYNTHETIC DATA AT FOUR DIFFERENT NOISE (SNR) LEVELS. FOR EACH ENTRY, THE FIRST RELATIVE ERROR VALUE WAS COMPUTED WITH (16), WHEREAS THE SECOND VALUE [AFTER THE SLASH (/)] WAS COMPUTED WITH (17). FOR THE N-FINDR TRIAL RESULTS GIVEN HERE, ACROSS ALL SNR LEVELS, THE MAXIMUM NUMBER OF ENDMEMBERS WAS SET TO TEN. DASHES ( ) INDICATE THAT A GIVEN ENDMEMBER WAS NOT EXTRACTED

http://speclib.jpl.nasa.gov/ (URL last updated September 24, 2002; URL last accessed July 30, 2008)], and sample vegetation spectra were obtained from the United States Geological Survey (USGS-Reston) Spectral Library (e.g., [41]). All library reflectance spectra were spectrally subset and convolved, using a Gaussian filter function, to 276 bands in the 420–2500-nm wavelength range. The laboratory-based spectra/endmembers that were used were aspen leaf, black tar roofing paper, brown loamy fine sand, brown sandy loam, brown silty loam, construction concrete, construction tar, fir tree, and maple leaf. Fig. 5 shows all the endmember spectra used for generating the pseudoimage data set. The simulated data set consisted of 256 × 256 pixels, where each pixel was a random linear mixture for two clusters of endmembers. The first third of the data set was generated from the random linear mixing of the three man-made-material endmembers, and the remainder of the data set was generated from the linear mixing of the six natural-material endmembers using the model given in (12) and (13), where, again, the mixing proportions were random. Thus, there were no pure endmember

pixels in the simulated data set. Gaussian noise was imposed, resulting in synthetic images with four different SNR levels, namely, 5, 10, 20, and 30 dB. Both SVM-BEE and N-FINDR were applied to these synthetic-image data for endmember extraction. Table I compares the relative error of N-FINDR- and SVDDextracted and SVM-BEE-estimated endmembers to the true endmembers, or eest − etrue /etrue .

(16)

Note that since the spectra used in this synthetic-data experiment are in reflectance units, they are actually implicitly normalized; we therefore computed the following additional error measure: eest − etrue 

(17)

which may also be of utility in the quantitative evaluation of the extracted endmembers. Thus, relative error values computed using both (16) and (17) are summarized in Table I. It can

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

780

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

be seen from Table I that SVM-BEE accurately estimates all the endmembers even with substantial noise. Fig. 6 qualifies this statement by displaying the SVM-BEE estimation for the most extreme case, where the SNR is 5 dB. Although the final estimated endmembers suffer from the added noise, they capture the signal structure of each true endmember. SVDD produces results similar to those of SVM-BEE, with SVDD performing at a slightly lesser level, particularly in the most extreme cases of SNR. However, it is noted here and explained in greater detail at the end of the next section that SVDD could not be applied to the entire scene due to memory requirements. Thus, to address this problem and limit the amount of data with which SVDD had to contend, for this synthetic data set, as well as the synthetic data set used in the following section, we applied the SVDD method to subsets of the synthetic image that contained the true endmembers and maintained a uniform sampling of the range of fractional abundances. The problem with the memory requirements of SVDD highlights one of its major differences with SVM-BEE. Because SVM-BEE trains SVMs only on small subsets, the problems associated with training an SVM on the entire data set are avoided. It is noted that some endmember estimations have relative error that is significantly greater than others. This is a direct result of how the noise was added to the system. Random noise was added after the linear mixture of the signals at a constant power relative to the average power of all the true endmembers. Thus, endmembers such as black tar roofing paper and construction tar that have low power actually have a much lower SNR. Using the same JHU/USGS spectral library-based syntheticimage data, at the same varying noise levels, the N-FINDR algorithm was run to extract endmembers. The N-FINDR synthetic-image experiments proceeded as follows: For each SNR level [i.e., 5 dB (3.16 : 1), 10 dB (10 : 1), 20 dB (100 : 1), and 30 dB (1000 : 1)], N-FINDR was run three times such that the user-specified maximum number of endmembers varied from 10 to 12 to 15, resulting in a total of 12 N-FINDR synthetic-image trials. Note that these upper limit values are greater than the actual number of endmembers in the synthetic image (which was nine), so the number of endmembers that was returned was N-FINDR determined rather than user specified. These values were set as per the N-FINDR recommendations, as noted in Section IV-C [40]. In addition, for all synthetic N-FINDR trials, a bad-pixel detection and removal algorithm was enabled. N-FINDR compares endmember and average image spectra; bad pixels are considered to have single-band spikes in their spectra. If bad pixels are not removed, N-FINDR will identify such pixels as endmembers [40] (note that badpixel removal is not a preprocessing step of SVM-BEE, as discussed in the following discussion). Also, in accordance with [11], for example, endmember averaging, where groups of pixels with similar spectral features are averaged, which is sometimes performed to reduce noise, was not performed here with any of the endmember extraction algorithms in order to facilitate a direct comparison of algorithm performances. This also applies to all other experiments conducted in this paper.

Fig. 6. SVM-BEE estimation of endmembers based on the synthetic data set for (a) fir tree and (b) black tar roofing paper for the case when the SNR is 5 dB, as well as at other noise levels. (c) Results of the scree test for each cluster. It is clear that the number of endmembers will be six and three for the natural- and man-made-material clusters, respectively.

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

781

Fig. 7. Number of endmembers extracted by N-FINDR by synthetic-data trial. For each N-FINDR synthetic data trial, the SNR, maxEM, and NEE values are as follows: (a) Trial 1: SNR = 5 dB, maxEM = 12, and NEE = 2; (b) Trial 2: SNR = 5 dB, maxEM = 10, and NEE = 2; (c) Trial 3: SNR = 5 dB, maxEM = 15, and NEE = 2; (d) Trial 4: SNR = 10 dB, maxEM = 12, and NEE = 2; (e) Trial 5: SNR = 10 dB, maxEM = 10, and NEE = 2; (f) Trial 6: SNR = 10 dB, maxEM = 15, and NEE = 3; (g) Trial 7: SNR = 20 dB, maxEM = 12, and NEE = 4; (h) Trial 8: SNR = 20 dB, maxEM = 10, and NEE = 5; (i) Trial 9: SNR = 20 dB, maxEM = 15, and NEE = 5; (j) Trial 10: SNR = 30 dB, maxEM = 12, and NEE = 7; (k) Trial 11: SNR = 30 dB, maxEM = 10, and NEE = 8; (l) Trial 12: SNR = 30 dB, maxEM = 15, and NEE = 5.

Note that for synthetic N-FINDR Trials 1–3, SNR = 5 dB; for Trials 4–6, SNR = 10 dB; for Trials 7–9, SNR = 20 dB; and for Trials 10–12, SNR = 30 dB. For each set of three trials at each SNR level, the maximum number of endmembers was first set to 12, followed by 10, and then 15 endmembers. In terms of the number of endmembers actually extracted via N-FINDR, the number varied between two and eight endmembers across the 12 trials, and in Trials 1–5, only two endmembers were found. The largest number of N-FINDR-extracted endmembers (eight) was in Trial 11, where SNR = 30 dB, and

the maximum number of endmembers was set to ten. Thus, the best N-FINDR performance in terms of the number of endmembers extracted was realized at the lowest noise level tested, and when the maximum number of endmembers parameter was set just above the true number of endmembers, which agrees with the assessment of Berman et al. [19] regarding N-FINDR. A summary of these results by trial is shown in Fig. 7, which demonstrates that there is wide variability in the N-FINDRextracted endmembers as a function of the maximum number of endmembers and the SNR parameters. As indicated earlier, in

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

782

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

TABLE II FRACTIONAL ABUNDANCE, IN PERCENTAGE, OF THE N-FINDR- AND SVDD-EXTRACTED AND SVM-BEE-ESTIMATED ENDMEMBERS FROM THE CS1 MAN-MADE-/NATURAL-MATERIAL SYNTHETIC IMAGE AT SNR = 30 : 1 AND 1000 : 1. FOR EACH PIXEL, THE FRACTIONAL ABUNDANCE OF ONLY THE MOST PREDOMINANT SPECTRUM IS GIVEN AS AN INDICATOR OF SIGNATURE PURITY. MULTIPLE N-FINDR TRIALS ARE DESIGNATED (T1, T2, ETC.). DASHES ( ) INDICATE THAT A GIVEN ENDMEMBER WAS NOT EXTRACTED

no case did N-FINDR find all nine endmembers in the synthetic data set. Regarding N-FINDR-extracted endmember accuracy, since the result of synthetic N-FINDR Trial 11 most closely approached the true number of endmembers, these endmembers were evaluated for error. Table I thus compares the relative error of these N-FINDR-extracted endmembers to the true endmembers, in addition to the SVM-BEE and SVDD results. We also consider in Table I all other synthetic N-FINDR trial results, across all SNR levels, where the maximum number of endmembers was set to ten (as was the case in Trial 11), which, again, is just above the actual number of true endmembers. Therefore, N-FINDR-extracted endmembers from Trials 2, 5, and 8 are evaluated in Table I, in addition to those from Trial 11 (see Fig. 7 caption for details). These results demonstrate that N-FINDR had some difficulty in extracting endmembers in the absence of pure endmember pixels and that the endmember set was too small in the presence of noise; these N-FINDR findings are in agreement with those in [19], for example, noted previously. Further synthetic-image experimentation was conducted, as described in the following paragraph. E. Example: Synthetic Hyperspectral Image Test #2 To better approximate a real hyperspectral image, our second test of the algorithms was patterned after the synthetic image generation logic established in [42]. Synthetic hyperspectral image test #2 specifically involved generating an image essentially in accordance with the first computer-simulated image in [42] (i.e., what is termed image “CS1”), except that the man-made- and natural-material endmembers from synthetic hyperspectral image test #1 were used (Fig. 5) rather than the mineral spectra set utilized in [42]. Note that a “CS2”-type image [42] could not be generated using this set of man-made and natural spectra since not enough other exemplar spectra for these materials were available in these databases to provide the necessary material variability, as was done in [42] using mineralogical spectra. However, the use of this set of manmade- and natural-material spectra maintains consistency with synthetic hyperspectral image test #1, and it presents a case

where such materials are mixed and where overlapping clusters exist. We generated a square 100 × 100-pixel image, with the four corner pixels being comprised of pure spectral signatures for aspen leaf, black tar roofing paper, brown loamy fine sand, and brown sandy loam, and the center pixel was generated with a pure spectrum of maple leaf (i.e., five pure pixels in the image). The signature fractional abundance of each pure/endmember spectrum was made to linearly decrease from the corner and center 100%-abundance pure pixels. Abundance fractions progressively and regularly decreased in a radial pattern emanating from these pure pixels until an abundance value of 0% was reached outside of an imaginary circular perimeter. The imaginary circles centered on the corner pixels had radii of 60 pixels, whereas that centered on the center image pixel had a radius of 25 pixels. Four spectral mixtures were generated and inserted midway between pairs of corner pixels [42]. These composite spectra were as follows: 71.48% brown silty loam/28.52% maple leaf; 71.48% construction concrete/28.52% maple leaf; 72.53% construction tar/27.47% maple leaf; and 72.53% fir tree/27.47% maple leaf. The radius of each imaginary circle centered on these impure spectra was 49 pixels, and the signal purity decreased linearly from slightly greater than 70% for these spectra. All fractional abundances were fully constrained [42]. As in [42], Gaussian noise with an SNR of 30 : 1 was added; we also conducted additional trials where the same CS1 image was used but with an SNR of 1000 : 1. See [42] for more details on the process used for generating these synthetic image data. SVM-BEE was then used to estimate the endmembers in the CS1 man-made-/natural-material synthetic image at SNR values of both 30 : 1 and 1000 : 1. SVM-BEE signal purity results are given in Table II, where the fractional abundance of just the dominant signature in each pixel of interest is given. Overall, SVM-BEE yielded excellent fractional abundance estimates. N-FINDR was also applied to the man-made-/naturalmaterial CS1 synthetic image. Four N-FINDR trials were run on this CS1 image for each of the two SNR levels, where we varied the N-FINDR maxEM parameter since the proper value

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

was unknown a priori. Specifically, maxEM = 9, 12, 15, and 20 in the respective trials; the N-FINDR SNR parameter was set at 30 : 1 and 1000 : 1 in separate sets of trials, which matches the actual SNR of the respective images. N-FINDR results for CS1 trials are given in Table II. For all trials at SNR = 30 : 1, N-FINDR consistently extracted only three endmembers. However, when SNR = 1000 : 1, N-FINDR extracted nine or more endmembers across all maxEM steps. This demonstrates N-FINDR’s sensitivity to noise, which was shown to be greater than that of SVM-BEE, as well as SVDD. As for the results of applying SVDD to the CS1 man-made-/natural-material synthetic image (Table II), signal purity/fractional abundance values were similar to those of SVM-BEE, although, where the results varied, SVM-BEE usually performed somewhat better than SVDD. Both SVMBEE and SVDD almost always posted results superior to those of N-FINDR, where the computed fractional abundances of the dominant signatures were typically closer to the respective actual maximum partial fractions. For certain endmember materials, there was a marked difference in support-vectorversus N-FINDR-computed signal purity values. Finally, it is observed that SVDD had problems in identifying the composite spectra of construction tar/maple leaf and fir tree/maple leaf as compared to SVM-BEE. For these composite spectra, the endmembers were identified as support vectors with SVM-BEE, but not for SVDD. Both SVM-BEE and SVDD identify a set of support vectors from which the endmembers are extracted. However, for these composite spectra—even with varying the radial basis parameter and penalty parameter over reasonable ranges—SVDD did not identify these as endmembers. In addition to these synthetic-data experiments, the SVMBEE and N-FINDR endmember extraction algorithms were also tested with real airborne hyperspectral imagery, as discussed in the next section. It should be noted that we did not apply the SVDD method to this real airborne hyperspectral image due to the following limitations. Given that the number of pixels to be analyzed is NS , the Lagrangian minimization problem in SVDD requires a double-precision array of the size Ns × Ns to be used in computing the solution. For an image that is 100 × 100, the memory requirement for the Lagrangian minimization problem is on the order   O 8Ns2 /10243 ≈ 0.75 GB. With 4 GB of RAM, the Matlab programming environment, which we exclusively utilized for all calculations with the exception of the commercially available N-FINDR, has a maximum variable size of 0.677 GB. Although increased performance and memory can be obtained by computing outside Matlab, the problem of using the entire data set in computation becomes an issue very fast. The SVDD Lagrangian minimization problem for a 165 × 165-pixel image contains an array on the order of 5.5 GB. Computing with an array of this size with 4 GB of RAM requires computationally expensive input–output (IO) operations. In this paper, we applied SVDD on both synthetic data sets using only Ns = 5000 because it is possible to generate subsets

783

Fig. 8. Band 29 (645.96-nm) subscene extracted from the AVIRIS image acquired on June 19, 1997, over the Cuprite mining district, NV, USA.

that contain the true endmembers and maintain a uniform sampling of the range of fractional abundances. In the case of the real airborne hyperspectral image, which is discussed in the next section, we cannot generate a subsample that we can guarantee as an accurate representation of the entire scene. For this reason, no comparison involving SVDD was performed on the real hyperspectral image data discussed next. F. Example: AVIRIS Hyperspectral Imagery 1) Study Area and Data Description: The Cuprite mining district, NV, USA, served as the study area and was selected due to its common use as a remote-sensing validation site (e.g., [43]). The geologic and geomorphologic characteristics of the Cuprite site have been discussed in detail in [44]–[49]. The Cuprite district is located in the southwestern Great Basin, approximately 15 km south of Goldfield, NV, in Esmeralda and Nye Counties, and it covers an area of ∼12 km2 . Only very sparse vegetation exists at the site, and three degrees of hydrothermal alteration exist in the area, with alteration occurring in two elongate north–south areas on either side of US-95. The most altered rocks to the least altered are referred to as silicified, opalized, and argillized, respectively [50]. The hyperspectral image data used in the experiments were obtained from AVIRIS [51]. AVIRIS is a hyperspectral whiskbroom scanner that entails 224 contiguous spectral channels in the 370–2500-nm nominal wavelength region, measuring upwelling radiance at ∼10-nm intervals. When flown in high-altitude mode at 20 000 m above ground level, given the 1.0-mrad instantaneous field-of-view (IFOV), the nominal ground-projected IFOV is 20 m [51]. The same AVIRIS scene employed in [11] and a variety of other endmember studies, which was acquired over the Cuprite mining district on June 19, 1997, was also analyzed in this paper (Fig. 8). The Cuprite AVIRIS data are available in reflectance units from the AVIRIS Web site [Jet Propulsion Laboratory.

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

784

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

Fig. 9. SVM-BEE-estimated endmember set extracted from the 1997 Cuprite AVIRIS image for (a) the five-endmember and (b) 22-endmember cases, respectively.

“Ordering free AVIRIS standard data products,” http:// aviris.jpl.nasa.gov/html/aviris.freedata.html (URL last updated October 30, 2007; URL last accessed July 30, 2008)]. According to Rowan et al. [52], this calibrated 1997 AVIRIS scene was atmospherically corrected using the MODTRAN radiative transfer model [53] at the Jet Propulsion Laboratory, facilitating the comparison of image-extracted endmembers with a laboratory-based spectral library. Similar to that in [11], bands 1–3, 105–115, and 150–170 were removed from this analysis due to atmospheric water vapor absorption and low SNR. In addition, in this paper, bands 221–224 were discarded due to noise and lack of spatial coherency, and the remaining AVIRIS spectral overlap bands (bands 33–34 and 97–98) were also removed prior to the analysis, resulting in a 181-channel data set. Approximately the same spatial subset used in [11] was used here (355 samples × 356 lines) (Fig. 8). No attempt was made to geometrically rectify the image, consistent with other studies (e.g., [50] and [54]). These AVIRIS reflectance data were used as input to the SVM-BEE and N-FINDR algorithms. All AVIRIS-image-based SVM-BEE-estimated and N-FINDR-extracted endmembers were compared against the laboratory-based USGS mineral spectral library described in [41], and subsequent unmixing results were compared against reference mineral maps [47], [55] and AVIRIS-derived USGS Tetracorder results [56], where the latter is also available at http://speclab.cr.usgs.gov/. 2) AVIRIS Endmember Experiments: The SVM-BEE method for the AVIRIS experiment was conducted precisely as outlined in the algorithm and previous examples. Parameters were estimated using (10), and the scree point of Λ was determined by visual inspection to be mI = 5. That is, most of the variance in the image can be attributed to five endmembers; it is possible for other endmembers explaining less variance in the scene to be included in the endmember set, if desired. Note that a problem with a PCA-based method of determining the number of endmembers is that PCA deals with the statistical significance of the spectra rather than their uniqueness. Consequently, small objects or materials that are limited in areal extent, which may be potentially important for a particular application, may be relegated to the last principal components [57]. Thus, other methods for

determining the number of endmembers in a given scene, such as those based on estimating virtual dimensionality (VD) [2], [58], should be investigated in the future. Kernel PCA plus a linear discriminant, or a kernel Fisher discriminant, would also likely yield an improvement over PCA, as would an independent-component-analysis (ICA)-based approach, perhaps combined with VD estimation [59], [60]. Spatial ICA has been used in conjunction with spectral Bayesian positive source separation (BPSS), where spatial ICA generates a relevant pixel subset, and BPSS estimates source spectra based on this reduced pixel set [61]. The minimum noise fraction (MNF) transform [62] could also be employed as an alternative to PCA. Thus, estimating the number of endmembers/signal sources in a given image can very likely be better accomplished by using methods other than linear PCA, but this issue is not the central focus of what we are proposing here with SVM-BEE. In the current investigation, convergence of the decision surface occurred in ∼40 iterations with only one cluster identified and therefore minimizing a single simplex to fit the edge points, identifying five endmembers. The SVM-BEE method was run three different times on the AVIRIS data with consistent endmembers estimated across trials within similar computational time. Shown in Fig. 9 are the estimated endmembers from one of these five-endmember SVM-BEE trials. In addition, based on the findings of Chang and Plaza [11], where VD of an AVIRIS spatial subset nearly identical to the one analyzed in this paper was estimated using the Neyman–Pearson-detection-theory-based Harsanyi–Farrand–Chang eigenthresholding method [2], [63], an SVM-BEE trial was also run where the total number of estimated endmembers was set to 22. These 22 extracted endmembers are also shown in Fig. 9. Finally, note that unlike N-FINDR [40], it was not necessary to remove “bad pixels” for SVM-BEE to correctly identify endmembers. This advantage can be attributed to SVM-BEE’s ability to tolerate noise. Regarding N-FINDR, several trials were performed, where bad pixels were detected and removed. As with one of the SVM-BEE trials, for N-FINDR AVIRIS Trial 1, the number of endmembers for N-FINDR to extract was explicitly set based on the findings of Chang and Plaza [11] (i.e., 22 endmembers). The reasonability of utilizing 22 endmembers was also

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

785

Fig. 10. Example SVM-BEE- and N-FINDR-extracted AVIRIS endmembers from the five-endmember trials and the corresponding USGS library reference spectra for (a) alunite and (b) hematite. The SAM angle (in radians) between each extracted endmember and the respective USGS reference spectrum is noted. SVM-BEE and N-FINDR endmembers have been scaled.

independently assessed via computation of an MNF transform [62] of the Cuprite AVIRIS image, with a subsequent joint inspection of the MNF scree plot and the spatial coherency of the eigenimages. N-FINDR AVIRIS Trial 2 represented the least restrictive case. In Trial 2, the approximate SNR was set to 1280 : 1, which was the mean across all bands of the maximum SNR value for each image band, calculated directly from the image. The SNR was computed as the mean/standard deviation in a 3 × 3 moving window [64]. Given the recommendation for the N-FINDR SNR parameter [40] (discussed previously), setting the SNR to this value for this trial is less restrictive, relative to Trial 3 (see in the following discussion), in the sense that N-FINDR has the ability to identify more endmembers as the SNR increases. In Trial 2, triplicate runs were performed, and the number of extracted endmembers (NEEs) in the three runs varied (13, 15, and 17 endmembers, respectively), even though all input parameters were constant across the runs. For both N-FINDR AVIRIS Trails 2 and 3, the maximum number of endmembers was set to 30. N-FINDR AVIRIS Trial 3 was a more restrictive case than Trial 2 in terms of the SNR parameter. For Trial 3, the

Fig. 11. SVM-BEE- and N-FINDR-extracted AVIRIS endmembers from the 22-endmember trials and the corresponding USGS library reference spectra for (a) alunite, (b) calcite, and (c) kaolinite. The SAM angle (in radians) between each extracted endmember and the respective USGS reference spectrum is noted. SVM-BEE and N-FINDR endmembers have been scaled.

SNR parameter value was set to 125 : 1, which was the mean across all bands of the mean SNR value for each image band; this was also computed directly from the image [64]. Given the lower SNR parameter value used, N-FINDR extracted fewer endmembers (i.e., five endmembers), which is equivalent in quantity to the number of AVIRIS endmembers

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

786

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

TABLE III SAM-BASED SPECTRAL SIMILARITY VALUES (IN RADIANS) GENERATED BY COMPARING 22 SVM-BEE-ESTIMATED ENDMEMBERS AGAINST THREE USGS LIBRARY SPECTRA IN THE 1.979–2.467-µm WAVELENGTH RANGE. USGS MINERAL NAMES ARE GIVEN IN THE COLUMN HEADERS, WHEREAS SVM-BEE ENDMEMBER NUMBERS ARE GIVEN IN THE ROWS (i.e., S1 . . . S22). SAM VALUES IN BOLD INDICATE THE BEST ENDMEMBER MATCHES FOR A PARTICULAR MINERAL TYPE

estimated by SVM-BEE when the number of endmembers was PCA determined. For this reason, the N-FINDR result from Trial 3 was evaluated here. In addition, the result from the 22-endmember N-FINDR AVIRIS trial (Trial 1) was also assessed in this paper to facilitate comparison with other studies. The spectral angle mapper (SAM) algorithm [65] was used as a spectral similarity measure to compare endmember spectra, extracted via SVM-BEE and N-FINDR, against USGS Spectral Library reference spectra. SAM determines the angle between spectral signatures as a measure of spectral similarity, and it is relatively insensitive to variation in illumination/gain factors [65], [66]. This is an attractive property in the context of this evaluation, as there are scale and illumination differences between the laboratory reference spectra and image-derived endmembers [18]. Two pairs of example endmember spectra that were common to both the SVM-BEE and N-FINDR five-endmember trials are discussed here. Alunite (alunite GDS82 Na82 [i.e., alunite #3)] and hematite (hematite GDS27 [i.e., hematite #2)] were identified via SAM values computed using SVMBEE- and N-FINDR-extracted endmembers, and both alunite and hematite have previously been detected by AVIRIS at the Cuprite site, as well as identified via X-ray diffraction [67]. For each mineral type, Fig. 10 shows the SVMBEE-estimated and N-FINDR-extracted endmembers, as well as the corresponding USGS Spectral Library reference spectrum. Following a variety of mineralogical studies at the Cuprite site (e.g., [16] and [68]), the spectra were plot-

ted in the 1.978–2.467-μm wavelength range, and the SAM values were also computed based on the 50 bands in this range. That is, 181 bands were used in the endmember extraction step, and endmembers were evaluated over this 50-band range, as this SWIR wavelength region often contains diagnostic absorption features for various minerals [16], [56], [68]. The 1.978–2.467-μm wavelength range was also utilized in the assessment of spectra in the 22-endmember set (see in the following discussion). For the SVM-BEE-generated alunite endmember (from the five-endmember set), the SAM angle was 0.034089 rad; N-FINDR posted a SAM angle of 0.041220 rad for its alunite endmember. For hematite, the SVM-BEE endmember SAM angle was 0.017448 rad, whereas it was 0.022540 rad for the corresponding N-FINDR endmember (Fig. 10). Thus, based on the five-endmember trials, in terms of the SAM spectral similarity measure, SVM-BEE yielded AVIRIS-derived endmembers with closer matches to the USGS library spectra than N-FINDR. For the 22-endmember SVM-BEE and N-FINDR trials, extracted-endmember evaluations for the following minerals were conducted: alunite, calcite, and kaolinite [18], [68], as these minerals are known to exist within the Cuprite study area (e.g., [16], [47], [54], [56], and [67]). The SVM-BEE results for the 22-endmember case are given in Table III, which consists of a matrix of SAM-based spectral similarity values computed using SVM-BEE-estimated endmembers and three USGS reference spectra for alunite, calcite, and kaolinite, respectively. The corresponding N-FINDR results are given in Table IV. Fig. 11 provides a visual comparison between AVIRIS-image-derived endmembers and corresponding USGS laboratory spectra. To minimize the mean square error with respect to the laboratory spectra, SVM-BEE- and N-FINDRextracted endmembers were scaled [20]; endmember spectral shape was not altered (Figs. 10 and 11). Joint inspection of Tables III and IV and Fig. 11 indicates that the smallest SVMBEE- and N-FINDR-based SAM angles, which indicate the closest matches to the reference spectrum, are nearly the same for both alunite and kaolinite. In terms of the SAM angles, N-FINDR performed somewhat better with calcite though. In Fig. 11, some differences in the reflectances and spectral shapes are evident, which may be due to atmospheric transmission effects [18], [68]. However, SVM-BEE was generally able to obtain the shapes of the spectral features quite well, relative to the USGS reference spectra. For instance, SVM-BEE seemed to approximate the kaolinite doublet feature around 2.2 μm quite closely, relative to N-FINDR. As a caveat, note that other minerals also exhibit absorption features near 2.2 μm that overlap the kaolinite feature [56], so caution should be exercised when interpreting the results. As mentioned earlier, N-FINDR does seem to yield a better calcite endmember, based on these respective 22-endmember sets. However, the SVM-BEE endmember labeled as calcite represents only the best SAM match from the 22-endmember set with the USGS calcite laboratory spectrum (i.e., it may not actually be a calcite endmember). It is possible that if the number of endmembers extracted was increased, then SVM-BEE may approximate a better calcite endmember. For instance, for approximately the same AVIRIS subset, the number of endmembers could

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

TABLE IV SAM-BASED SPECTRAL SIMILARITY VALUES (IN RADIANS) GENERATED BY C OMPARING 22 N-FINDR-E XTRACTED E NDMEMBERS A GAINST THREE USGS LIBRARY SPECTRA IN THE 1.979–2.467-µm WAVELENGTH RANGE. USGS MINERAL NAMES ARE GIVEN IN THE COLUMN HEADERS, WHEREAS N-FINDR ENDMEMBER NUMBERS ARE GIVEN IN THE ROWS (i.e., N1 . . . N22). SAM VALUES IN BOLD INDICATE THE BEST ENDMEMBER MATCHES FOR A PARTICULAR MINERAL TYPE

787

is evident that the SVM-BEE- and N-FINDR-based abundance distribution patterns are consistent with the reference maps, although, as expected, based on the SAM angular match with the USGS calcite spectrum and the endmember plots, the SVMBEE-derived calcite abundance map appears less than optimal. The CEM fractional abundance images (Fig. 12) support the findings given in Tables III and IV and in Fig. 11. V. D ISCUSSION AND C ONCLUSION

potentially reasonably be increased from 22, based on VD and the selection of different false-alarm probabilities PF [11], [42]. The same notion also applies to other SVM-BEE and N-FINDR endmembers. In general, as suggested by their similar SAM values (Tables III and IV) and endmember plots (Figs. 10 and 11), SVM-BEE and N-FINDR both extract reasonable AVIRIS endmembers. As an additional means of endmember evaluation, we used constrained energy minimization (CEM) [50], [69] to map the endmember fractional abundances. Unlike some other methods, CEM requires only the spectrum of the target material of interest as input (i.e., not all image endmembers need be known). This advantageous property enables only the several endmembers for the materials of interest to be evaluated individually. CEM maximizes the response of the target spectrum while suppressing the response of the unknown background signatures. The CEM operator has a unity constraint, although derived abundance values can go negative or superunitary. Negative values can occur due to existence of independent identically distributed noise with zero mean in the image data. If the target spectrum is not completely pure or of the same composition as the image pixels, then the computed abundances can be greater than one [50], [69], [70]. The AVIRIS-image-derived mineral abundance maps in Fig. 12 give the approximate subpixel CEM fractional abundances of target materials as grayscale images. Note that no abundances are superunitary, although some are negative. From a joint visual assessment of the CEM-based abundances and existing mineral map sources (e.g., [47], [54], [55], and [56]), it

In general, SVM-BEE performed much better than N-FINDR on synthetic data, in terms of both the number of endmembers extracted and the error associated with those endmembers, although there were several instances where the relative endmember extraction errors for SVM-BEE and N-FINDR were the same or very nearly the same. For real-image AVIRIS data, SVM-BEE and N-FINDR performed comparably on the mineralogical scene analyzed here, although evaluation is more difficult in the case of real data [56]. As noted previously, SVM-BEE was able to estimate some endmembers more accurately than N-FINDR (e.g., kaolinite and its doublet feature near 2.2 μm) and vice versa. Since the details of the current commercial N-FINDR algorithm are proprietary and thus unknown to the user [39], it is difficult to surmise all reasons underlying its performance in this paper. However, N-FINDR has previously been extensively tested on Cuprite AVIRIS image data (e.g., [16] and [17]), and based upon its known properties, it would be expected for N-FINDR to perform well in identifying endmembers such as those investigated here since relatively pure pixels may exist within the Cuprite scene for these minerals [42], [54]. If the synthetic-data experiment in this paper, as well as in other studies (e.g., [19]), is any indication, however, N-FINDR results may not be as favorable when pure endmember pixels for given materials do not exist within a given scene. Further comparative validation efforts for such cases using real remote-sensor images are needed, as is further testing of SVM-BEE compared against various other endmember extraction algorithms, including SVDD, in varied environments. Also, with respect to the real-image experiment, note that SAM may not be an optimal spectral similarity measure [71]; for instance, Chang et al. [72] found that, for real hyperspectral images, correlation matched-filter-based distance performed significantly better than SAM in discrimination and identification of subpixel targets and mixed pixels, and this method was also employed in [54]. A variety of other spectral similarity measures or shape-matching algorithms have also been used in the context of endmember evaluation (e.g., [2] and [56]). However, various studies (e.g., [42]) have previously utilized SAM as a spectral similarity measure in this context, and its use here represents a basic common baseline measure. Other metrics can be deployed in further studies. Regarding fractional abundance mapping, more accurate results—for both SVM-BEE and N-FINDR—may be accrued via other unmixing algorithms, such as mixture-tuned matched filtering (MTMF) [9], [73] which additionally provides a means of minimizing false positives, the evaluation of which is an important dimension of verification in imaging spectroscopy, although it is often difficult to

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

788

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

Fig. 12. CEM-generated fractional abundance maps for alunite GDS82 Na82 (a–d), calcite CO2004 (e–h), and kaolinite KGa-2 (pxyl) (i–l), based on SVM-BEEand N-FINDR-extracted endmembers, as noted in the column headers. The images in the two left columns have no stretch applied, whereas the images in the two right columns have been stretched over specific ranges. SVM-BEE- and N-FINDR-based alunite CEM fractional abundance images [in (c) and (d), respectively] are each stretched between 0 and 0.43, calcite abundance images in (g) and (h) are each stretched between 0 and 0.5, and kaolinite abundance images in (k) and (l) are each stretched between 0 and 0.39.

practically assess [56]. MTMF analysis yields an “infeasibility” image that indicates which pixels are likely to be matchedfilter false positives [9]. Numerous other abundance estimation methods could also be used [2]. We have identified a variety of unique features and/or advantages of SVM-BEE. As mentioned in Section III-A, SVMBEE has the ability to identify multiple separated distributions. It is emphasized that the number of separate distributions does not need to be known a priori but rather can be identified by the final decision function. In the setting of endmember extraction, the identification of separated distributions provides additional speed and accuracy. We have observed the computational benefits of SVM-BEE compared with N-FINDR and SVDD. Specifically regarding SVDD, in the experimentation performed, it was necessary to reduce the data volume/size of the synthetic images so that SVDD could process them without expensive IO operations. As for SVC, which has not been used for endmember extraction to date, SVC performs some postprocessing to SVDD to detect multiple separated

distributions; as such, it suffers from the same computational problems as SVDD, as discussed previously. Thus, even if SVC were to be implemented as an endmember extraction method, SVM-BEE would still hold a comparative advantage over an SVC approach. Although the intent of SVM-BEE is to extract endmembers, it is recognized that determining the number of distributions and the corresponding extreme points can provide valuable information for a broad range of analyses—a topic for future investigation. SVM-BEE is semiautonomous in the sense that the algorithm does not self-select all parameters that will necessarily work in absolutely every circumstance with which it is presented; however, we can say that for the cases investigated in this paper, the correct parameters were determined by SVMBEE, and in our experiments, the decision surface algorithm of SVM-BEE has operated in an effective essentially autonomous manner. We have also empirically found SVM-BEE to be consistent in its solutions; for instance, the synthetic experiments demonstrated that SVM-BEE’s NEEs were constant as a function of SNR, and the estimated endmembers were tolerant

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

to noise. In addition, when triplicate runs were performed with the number of endmembers set to 5 and 22, respectively (data not shown), very nearly the same set of SVM-BEE-estimated endmembers was extracted from the Cuprite AVIRIS image. Furthermore, SVM-BEE-estimated endmembers from the fiveendmember trial were reproduced in the 22-endmember trials (i.e., they closely matched a subset of the endmembers in the 22-endmember trial). Regarding other future work, since SVM-BEE is not constrained to the NB + 1 upper limit on endmembers, our algorithm holds much potential in not only hyperspectral endmember estimation but also with multispectral applications, where the number of bands and, thus, the number of endmembers that can be extracted via conventional means are limited. Future work will thus be conducted to assess SVM-BEE’s utility with low-dimensional data. The previously described properties of SVM-BEE suggest a distinct comparative advantage for SVM-BEE relative to SVDD when operating upon not only high-dimensional/hyperspectral data but also lowdimensional/multispectral data sets, and such a specific multispectral comparison will be further tested in future research. The application of SVDD to hyperspectral imaging, as given in [34], does not search for split distributions or very irregularly shaped distributions, a potential that is contained within the original design of SVDD, as indicated by the development of SVC. A comparative advantage of SVM-BEE relative to the SVDD endmember method of Banerjee et al. [34] follows from this characteristic. A major strength of the SVM-BEE algorithm is its ability to determine sample spectra that lie on the edges of multi/hyperdimensional distributions. This paper used the underlying assumption that all distributions came from linear mixing of endmembers. To the extent that intimate nonlinear spectral mixing exists in the Cuprite scene, this could have adversely affected the SVM-BEE results. By extending the SVM-BEE algorithm to operate on a broad range of possible mixing models, the accuracy of the endmember extraction of real data can be improved. SVM-BEE does an excellent job in approximating points that are on the edge of the distribution. If an adequate sampling of edge points exists, SVM-BEE could employ some model to fit the distribution—a model that is appropriate for a nonlinear mixing scenario. It should be noted, however, that the effects of nonlinear mixing are complex and case specific, whereas widely used linear mixture models are simpler and more generalized [74]; thus, most endmember extraction algorithms employ a linear mixing assumption as a standard for general-purpose uses [75]. Also, as noted in [32], nonlinear models may be used if convincing evidence exists for a given margin shape. Additionally, we plan to explore the use of these spectra generated by SVM-BEE as input for previously established endmember extraction methods (e.g., N-FINDR) to determine what computational savings and accuracy are gained by using SVM-BEE as a preprocessing step. Finally, the explosive amount of geospatial data being generated has created a strong need for the application of HPC algorithms that can accurately and semiautonomously or autonomously analyze and classify this continuous stream of data in a timely fashion. The SVM-BEE algorithm can be adapted

789

to HPC and has the potential to cut through massive volumes of data and extract meaningful samples for further analysis. SVMBEE is thus very open to parallelization, which is in contrast to many of the other endmember extraction algorithms, including the prior support-vector-based approaches. This is because the SVM-BEE workload can be divided upon multiple clusters with limited communication requirements, offering a dramatic potential advantage in terms of speed. Such a property may be particularly advantageous in time-sensitive applications, such as disaster response, among many others. ACKNOWLEDGMENT The authors would like to thank the anonymous reviewers for the useful comments that helped improve the quality of this paper. A. M. Filippi would also like to thank B. L. Bhaduri and E. A. Bright, Computational Sciences and Engineering Division, ORNL, for their support of this paper. M. Winter graciously provided an evaluation version of the N-FINDR algorithm. ASTER Spectral Library data were made available through the courtesy of the Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA. Copyright 1999, California Institute of Technology. All rights reserved. The AVIRIS data were also provided by the Jet Propulsion Laboratory. R EFERENCES [1] J. R. Jensen and A. M. Filippi, “Thematic information extraction: Hyperspectral image analysis,” in Introductory Digital Image Processing: A Remote Sensing Perspective, 3rd ed. J. R. Jensen, Ed. Upper Saddle River, NJ: Prentice-Hall, 2005, pp. 431–465. [2] C.-I. Chang, Hyperspectral Imaging: Techniques for Spectral Detection and Classification. New York: Kluwer, 2003. [3] D. M. Rogge, B. Rivard, J. Zhang, and J. Feng, “Iterative spectral unmixing for optimizing per-pixel endmember sets,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 12, pp. 3725–3736, Dec. 2006. [4] J. J. Settle and N. A. Drake, “Linear mixing and the estimation of ground cover proportions,” Int. J. Remote Sens., vol. 14, no. 6, pp. 1159–1177, Apr. 1993. [5] R. A. Schowengerdt, Remote Sensing: Models and Methods for Image Processing, 2nd ed. San Diego, CA: Academic, 1997. [6] J. Gruninger, A. J. Ratkowski, and M. L. Hoke, “The sequential maximum angle convex cone (SMACC) endmember model,” in Proc. SPIE—Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery X, S. S. Shen and P. E. Lewis, Eds. Bellingham, WA: SPIE, 2004, vol. 5425, pp. 1–14. [7] C. A. Bateson, G. P. Asner, and C. A. Wessman, “Endmember bundles: A new approach to incorporating endmember variability into spectral mixture analysis,” IEEE Trans. Geosci. Remote Sens., vol. 38, no. 2, pp. 1083–1094, Mar. 2000. [8] J. W. Boardman, F. A. Kruse, and R. O. Green, “Mapping target signatures via partial unmixing of AVIRIS data,” in Proc. Summaries JPL Airborne Earth Sci. Workshop, Pasadena, CA, 1995, pp. 23–26. [9] ENVI User’s Guide, RSI (Research Systems, Inc.), Boulder, CO, 2004. ENVI Vers. 4.1. [10] A. Bateson and B. Curtiss, “A method for manual endmember selection and spectral unmixing,” Remote Sens. Environ., vol. 55, no. 3, pp. 229– 243, Mar. 1996. [11] C.-I. Chang and A. Plaza, “A fast iterative algorithm for implementation of pixel purity index,” IEEE Geosci. Remote Sens. Lett., vol. 3, no. 1, pp. 63–67, Jan. 2006. [12] J. H. Bowles, P. J. Palmadesso, J. A. Antoniades, M. M. Baumback, and L. J. Rickard, “Use of filter vectors in hyperspectral data analysis,” in Proc. SPIE—Infrared Spaceborne Remote Sensing III, 1995, vol. 2553, pp. 148–157. [13] J. H. Bowles, M. Daniel, J. M. Grossman, J. A. Antoniades, M. M. Baumback, and P. J. Palmadesso, “Comparison of output from

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

790

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 3, MARCH 2009

[14] [15]

[16]

[17] [18] [19]

[20] [21] [22] [23] [24] [25] [26]

[27] [28] [29] [30]

[31] [32] [33]

[34] [35] [36]

ORASIS and pixel purity calculations,” in Proc. SPIE—Imaging Spectroscopy IV, 1998, vol. 3438, pp. 148–156. A. Ifarraguerri and C.-I. Chang, “Multispectral and hyperspectral image analysis with convex cones,” IEEE Trans. Geosci. Remote Sens., vol. 37, no. 2, pp. 756–770, Mar. 1999. R. A. Neville, K. Staenz, T. Szeredi, J. Lefebvre, and P. Hauff, “Automatic endmember extraction from hyperspectral data for mineral exploration,” in Proc. 4th Int. Airborne Remote Sens. Conf. Exhib./21st Can. Symp. Remote Sens., Ottawa, ON, Canada, Jun. 21–24, 1999, vol. II, pp. 891–896. M. E. Winter, “N-FINDR: An algorithm for fast autonomous spectral endmember determination in hyperspectral data,” in Proc. SPIE—Imaging Spectrometry V, M. R. Descour and S. S. Shen, Eds, 1999, vol. 3753, pp. 266–275. M. E. Winter, “Fast autonomous spectral end-member determination in hyperspectral data,” in Proc. 13th Int. Conf. Appl. Geologic Remote Sens., Vancouver, BC, Canada, 1999, vol. II, pp. 337–344. A. Plaza, P. Martínez, R. Pérez, and J. Plaza, “Spatial/spectral endmember extraction by multidimensional morphological operations,” IEEE Trans. Geosci. Remote Sens., vol. 40, no. 9, pp. 2025–2041, Sep. 2002. M. Berman, H. Kiiveri, R. Lagerstrom, A. Ernst, R. Dunne, and J. F. Huntington, “ICE: A statistical approach to identifying endmembers in hyperspectral images,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 10, pp. 2085–2095, Oct. 2004. J. M. P. Nascimento and J. M. B. Dias, “Vertex component analysis: A fast algorithm to unmix hyperspectral data,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 4, pp. 898–910, Apr. 2005. M. Graña, C. Hernandez, and J. Gallego, “A single individual evolutionary strategy for endmember search in hyperspectral images,” Inf. Sci., vol. 161, no. 3/4, pp. 181–197, Apr. 2004. K. Staenz, T. Szeredi, and J. Schwarz, “ISDAS—A system for processing/ analyzing hyperspectral data,” Can. J. Remote Sens., vol. 24, no. 2, pp. 99–113, 1998. J. Chanussot, J. A. Benediktsson, and M. Fauvel, “Classification of remote sensing images from urban areas using a fuzzy possibilistic model,” IEEE Geosci. Remote Sens. Lett., vol. 3, no. 1, pp. 40–44, Jan. 2006. B. E. Boser, I. M. Guyon, and V. N. Vapnik, “A training algorithm for optimal margin classifiers,” in Proc. 5th Annu. Workshop Comput. Learning Theory, 1992, pp. 144–152. V. N. Vapnik, The Nature of Statistical Learning Theory. New York: Springer-Verlag, 1995. M. Fauvel, J. Chanussot, and J. A. Benediktsson, “Evaluation of kernels for multiclass classification of hyperspectral remote sensing data,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., Toulouse, France, May 14–19, 2006, pp. II-813–II-816. G. Camps-Valls et al., “Robust support vector method for hyperspectral data classification and knowledge discovery,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 7, pp. 1530–1542, Jul. 2004. F. Melgani and L. Bruzzone, “Classification of hyperspectral remote sensing images with support vector machines,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 8, pp. 1778–1790, Aug. 2004. G. Camps-Valls and L. Bruzzone, “Kernel-based methods for hyperspectral image classification,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 6, pp. 1351–1362, Jun. 2005. G. M. Foody and A. Mathur, “The use of small training sets containing mixed pixels for accurate hard image classification: Training on mixed spectral responses for classification by a SVM,” Remote Sens. Environ., vol. 103, no. 2, pp. 179–189, Jul. 2006. M. Brown, S. R. Gunn, and H. G. Lewis, “Support vector machines for optimal classification and spectral unmixing,” Ecol. Model., vol. 120, no. 2/3, pp. 167–179, Aug. 1999. M. Brown, H. G. Lewis, and S. R. Gunn, “Linear mixture models and support vector machines for remote sensing,” IEEE Trans. Geosci. Remote Sens., vol. 38, no. 5, pp. 2346–2360, Sep. 2000. L. Wang, Y. Zhang, and C. Zhao, “Combination of linear support vector machines and linear spectral mixed model for spectral unmixing,” in Intelligent Computing in Signal Processing and Pattern Recognition, vol. 345, D.-S. Huang, K. Li, and G. W. Irwin, Eds. Berlin, Germany: SpringerVerlag, Sep. 2006, pp. 767–772. A. Banerjee, P. Burlina, and J. Broadwater, “A machine learning approach for finding hyperspectral endmembers,” in Proc. IEEE IGARSS, Jul. 23–27, 2007, pp. 3817–3820. D. M. J. Tax and R. P. W. Duin, “Data domain description using support vectors,” in Proc. 7th ESANN, Bruges, Belgium, Apr. 21–23, 1999, pp. 251–256. D. M. J. Tax and R. P. W. Duin, “Support vector domain description,” Pattern Recognit. Lett., vol. 20, no. 11–13, pp. 1191–1199, Nov. 1999.

[37] A. Ben-Hur, D. Horn, H. T. Siegelmann, and V. Vapnik, “Support vector clustering,” J. Mach. Learn. Res., vol. 2, pp. 125–137, Dec. 2001. [38] R. B. Cattell, “The scree test for the number of factors,” Multivariate Behav. Res., vol. 1, pp. 245–276, 1966. [39] C.-I. Chang, C.-C. Wu, W. Liu, and Y.-C. Ouyang, “A new growing method for simplex-based endmember extraction algorithm,” IEEE Trans. Geosci. Remote Sens., vol. 4, no. 10, pp. 2804–2819, Oct. 2006. [40] N-FINDR Visualization Package Documentation, Tech. Res. Assoc., Inc., Honolulu, HI, May 31, 2006. [41] R. N. Clark, G. A. Swayze, A. J. Gallagher, T. V. V. King, and W. M. Calvin, “The U. S. Geological Survey, Digital Spectral Library: Version 1: 0.2 to 3.0 microns,” U.S. Geological Survey, Denver, CO, Open File Rep. 93-592, 1993. [42] A. Plaza and C.-I. Chang, “Impact of initialization on design of endmember extraction algorithms,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 11, pp. 3397–3407, Nov. 2006. [43] A. Goetz and V. Strivastava, “Mineralogical mapping in the Cuprite mining district,” in Proc. AIS Data Anal. Workshop, 1985, pp. 22–29. JPL Pub. 85-41. [44] J. P. Albers and J. H. Stewart, “Geology and mineral deposits of Esmeralda County, Nevada,” Nevada Bureau Mines Geology Bull., vol. 78, pp. 1–80, 1972. [45] M. J. Abrams, R. P. Ashley, L. C. Rowan, and A. F. H. Goetz, “Mapping of hydrothermal alteration in the Cuprite mining district, Nevada, using aircraft scanner images for the spectral region 0.46 to 2.36 µm,” Geology, vol. 5, no. 12, pp. 713–718, Dec. 1977. [46] M. J. Abrams, R. P. Ashley, L. C. Rowan, A. F. H. Goetz, and A. B. Kahle, “Use of imaging in the 0.46–2.36 µm spectral region for alteration mapping in the Cuprite mining district, Nevada,” U.S. Geological Survey, Denver, CO, Open-File Rep. 77-585, 1977. [47] R. P. Ashley and M. J. Abrams, “Alteration mapping using multispectral images, Cuprite mining district, Esmeralda County, Nevada,” U.S. Geological Survey, Denver, CO, Open-File Rep. 80-367, 1980. [48] H. Shipman and J. B. Adams, “Detectability of minerals on desert alluvial fans using reflectance spectra,” J. Geophys. Res., vol. 92, pp. 10 391– 10 402, Sep. 1987. [49] F. A. Kruse, K. S. Kierein-Young, and J. W. Boardman, “Mineral mapping at Cuprite, Nevada with a 63-channel imaging spectrometer,” Photogramm. Eng. Remote Sens., vol. 56, no. 1, pp. 83–92, Jan. 1990. [50] R. G. Resmini, M. E. Kappus, W. S. Aldrich, J. C. Harsanyi, and M. Anderson, “Mineral mapping with HYperspectral Digital Imagery Collection Experiment (HYDICE) sensor data at Cuprite, Nevada, U.S.A.,” Int. J. Remote Sens., vol. 18, no. 7, pp. 1553–1570, May 1997. [51] R. O. Green et al., “Imaging spectroscopy and the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS),” Remote Sens. Environ., vol. 65, no. 3, pp. 227–248, Sep. 1998. [52] L. C. Rowan, S. J. Hook, M. J. Abrams, and J. C. Mars, “Mapping hydrothermally altered rocks at Cuprite, Nevada, using the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), a new satellite-imaging system,” Econ. Geol., vol. 98, no. 5, pp. 1019–1027, Aug. 2003. [53] R. O. Green, “Retrieval of reflectance from AVIRIS using a radiative transfer code,” in Proc. 3rd AVIRIS Workshop, R. O. Green, Ed., Pasadena, CA, May 20–21, 1991, pp. 200–210. JPL Pub. 91-28. [54] C.-I. Chang and B. Ji, “Fisher’s linear spectral mixture analysis,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 8, pp. 2292–2304, Aug. 2006. [55] G. Swayze, “The hydrothermal and structural history of the Cuprite Mining District, Southwestern Nevada: An integrated geological and geophysical approach,” Ph.D. dissertation, Univ. Colorado, Boulder, CO, 1997. [Online]. Available: http://speclab.cr.usgs.gov and ftp://ftpext.cr.usgs.gov/pub/cr/co/denver/speclab/pub/cuprite/gregg. thesis.images/mineralzone.thesis.gif [56] R. N. Clark et al., “Imaging spectroscopy: Earth and planetary remote sensing with the USGS Tetracorder and expert systems,” J. Geophys. Res., vol. 108, no. E12, p. 5131, 2003. DOI:10.1029/2002JE001847. [57] A. Plaza, D. Valencia, J. Plaza, and P. Martinez, “Commodity clusterbased parallel processing of hyperspectral imagery,” J. Parallel Distrib. Comput., vol. 66, no. 3, pp. 345–358, Mar. 2006. [58] C.-I. Chang and Q. Du, “Estimation of number of spectrally distinct signal sources in hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 3, pp. 608–619, Mar. 2004. [59] J. Wang and C.-I. Chang, “Independent component analysis-based dimensionality reduction with applications in hyperspectral image analysis,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 6, pp. 1586–1600, Jun. 2006.

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.

FILIPPI AND ARCHIBALD: SUPPORT VECTOR MACHINE-BASED ENDMEMBER EXTRACTION

[60] J. Wang and C.-I. Chang, “Applications of independent component analysis in endmember extraction and abundance quantification for hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 9, pp. 2601–2616, Sep. 2006. [61] S. Moussaoui, H. Hauksdóttir, F. Schmidt, C. Jutten, J. Chanussot, D. Brie, S. Douté, and J. A. Benediktsson, “On the decomposition of Mars hyperspectral data by ICA and Bayesian positive source separation,” Neurocomputing, vol. 71, no. 10–12, pp. 2194–2208, Jun. 2008. [62] A. A. Green, M. Berman, P. Switzer, and M. D. Craig, “A transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE Trans. Geosci. Remote Sens., vol. 26, no. 1, pp. 65–74, Jan. 1988. [63] J. Harsanyi, W. Farrand, and C.-I. Chang, “Determining the number and identity of spectral endmembers: An integrated approach using Neyman–Pearson eigenthresholding and iterative constrained rms error minimization,” in Proc. 9th Thematic Conf. Geologic Remote Sens., Feb. 1993. [64] ERDAS Field Guide, Leica Geosystems GIS Mapping, LLC, Atlanta, GA, 2003. [65] F. A. Kruse et al., “The Spectral Image Processing System (SIPS)— Interactive visualization and analysis of imaging spectrometer data,” Remote Sens. Environ., vol. 44, no. 2/3, pp. 145–163, May/Jun. 1993. [66] C. Kwan et al., “A novel approach for spectral unmixing, classification, and concentration estimation of chemical and biological agents,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 2, pp. 409–419, Feb. 2006. [67] G. A. Swayze, R. N. Clark, S. Sutley, and A. Gallagher, “Ground-truthing AVIRIS mineral mapping at Cuprite, Nevada,” in Proc. Summaries 3rd Annu. JPL Airborne Geosci. Workshop, Pasadena, CA, 1992, vol. 1, pp. 47–49. JPL Pub. 92-14. [68] M. E. Winter and E. M. Winter, “Comparison of approaches for determining end-members in hyperspectral data,” in Proc. IEEE Aerosp. Conf., Big Sky, MT, Mar. 18–25, 2000, vol. 3, pp. 305–313. [69] J. C. Harsanyi, “Detection and classification of subpixel spectral signatures in hyperspectral image sequences,” Ph.D. dissertation, Univ. Maryland, Baltimore County, Baltimore, MD, 1993. [70] W. H. Farrand and J. C. Harsanyi, “Mapping the distribution of mine tailings in the Coeur d’Alene River Valley, Idaho, through the use of a constrained energy minimization technique,” Remote Sens. Environ., vol. 59, no. 1, pp. 64–76, Jan. 1997. [71] O. A. de Carvalho, Jr. and P. R. Meneses, “Spectral correlation mapper (SCM): An improvement on the spectral angle mapper (SAM),” in Proc. 9th Airborne Earth Sci. Workshop, 2000. JPL Pub. 00-18. [72] C.-I. Chang, W. Liu, and C.-C. Chang, “Discrimination and identification for subpixel targets in hyperspectral imagery,” in Proc. IEEE Int. Conf. Image Process., Singapore, Oct. 24–27, 2004, pp. 3339–3342. [73] J. W. Boardman, “Leveraging the high dimensionality of AVIRIS data for improved sub-pixel target unmixing and rejection of false positive: Mixture tuned matched filtering,” in Proc. Summaries 7th Annu. JPL Sci. Workshop, Greenbelt, MD, 1998, vol. 1.

791

[74] X. Miao et al., “Estimation of yellow starthistle abundance through CASI2 hyperspectral imagery using linear spectral mixture models,” Remote Sens. Environ., vol. 101, no. 3, pp. 329–341, Apr. 2006. [75] J. Plaza, A. Plaza, P. Martínez, and R. Pérez, “H-COMP: A tool for quantitative and comparative analysis of endmember identification algorithms,” in Proc. IGARSS, Jul. 21–25, 2003, pp. 291–293.

Anthony M. Filippi received the B.A. (summa cum laude) degree in geography from Kansas State University, Manhattan, in 1995 and the M.S. and Ph.D. degrees in geography from the University of South Carolina, Columbia, in 1998 and 2003, respectively. He is currently an Assistant Professor with the Department of Geography, Texas A&M University, College Station. He was an Office of Naval Research/NASA Summer Fellow at Friday Harbor Laboratories, University of Washington, in 1998, and he has been a Faculty Fellow at Oak Ridge National Laboratory during the summers of 2005–2008. His research interests include hyperspectral remote sensing, environmental and ocean optics, machine learning, and GIS-based modeling. Dr. Filippi is a Fellow of the American Geographical Society, and he is a member of the American Geophysical Union, American Society of Limnology and Oceanography, American Society for Photogrammetry and Remote Sensing, Association of American Geographers, and The International Society for Optical Engineers. He is also a member of Phi Beta Kappa, Phi Kappa Phi, Golden Key National Honor Society, and Gamma Theta Upsilon.

Rick Archibald received the B.Sc. degree in honors physics and the M.Sc. degree in applied mathematics from the University of Alberta, Edmonton, AB, Canada, in 1996 and 1998, respectively, and the Ph.D. degree in mathematics from Arizona State University, Tempe, in 2002. He has completed postdoctorate work at the Center for System Science and Engineering Research, Arizona State University, and the Neuroscience Department, Brown University, Providence, RI. Currently, he is the Householder Fellow at the Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN. His current research interests include image reconstruction and analysis, biological modeling, and boundary detection.

Authorized licensed use limited to: Oak Ridge National Laboratory. Downloaded on September 22, 2009 at 14:46 from IEEE Xplore. Restrictions apply.