Endmember Bundle Extraction Based on Multiobjective Optimization

A number of endmember extraction methods have been developed to identify pure pixels in hyperspectral images (HSIs). The majority of them use only one spectrum to represent one kind of material, which ignores the spectral variability problem that particularly characterizes a HSI with high spatial resolution. Only a few algorithms have been developed to identify multiple endmembers representing the spectral variability within each class, called endmember bundle extraction (EBE). This article introduces multiobjective particle swarm optimization for the identification of multiple endmember spectra with variability. Unlike existing convex geometry-based EBE methods, which operate on a single geometry of the dataspace, the proposed method divides the observed data into subsets along the spectral dimension and simultaneously operates on multiple dataspaces to obtain candidate endmembers based on multiobjective particle swarm optimization. The candidate endmembers are then refined by spatial post-processing and sequential forward floating selection to produce the final result. Experiments are conducted on both synthetic and real hyperspectral data to demonstrate the effectiveness of the proposed method in comparison with several state-of-the-art methods.


Endmember Bundle Extraction Based on
Multiobjective Optimization

Rong Liu , Member, IEEE and Xiaoxiang Zhu , Senior Member, IEEE
Abstract-A number of endmember extraction methods have been developed to identify pure pixels in hyperspectral images (HSIs). The majority of them use only one spectrum to represent one kind of material, which ignores the spectral variability problem that particularly characterizes a HSI with high spatial resolution. Only a few algorithms have been developed to identify multiple endmembers representing the spectral variability within each class, called endmember bundle extraction (EBE). This article introduces multiobjective particle swarm optimization for the identification of multiple endmember spectra with variability. Unlike existing convex geometry-based EBE methods, which operate on a single geometry of the dataspace, the proposed method divides the observed data into subsets along the spectral dimension and simultaneously operates on multiple dataspaces to obtain candidate endmembers based on multiobjective particle swarm optimization. The candidate endmembers are then refined by spatial post-processing and sequential forward floating selection to produce the final result. Experiments are conducted on both synthetic and real hyperspectral data to demonstrate the effectiveness of the proposed method in comparison with several state-of-the-art methods.

I. INTRODUCTION
W ITH the ability to record abundant spectral information about materials, hyperspectral imagery has been widely used for various applications, including vegetation mapping [1], mineral exploration [2], agricultural assessment [3], and many others [4]- [6]. Hyperspectral unmixing (HU), currently a hot topic in the processing of hyperspectral images (HSIs), involves estimating quantitative abundances of pure ground components within a pixel, so as to derive quantitative information at the subpixel scale. The selection of pure ground components, referred to as endmembers, is important for the successful application of HU. Endmembers can be obtained either from observed data or from field or laboratory measurement. To reduce the amount of time and expense involved in field measurement and keep similar atmospheric effects between endmembers and the data to be unmixed, a number of methods have been proposed to automatically extract endmembers directly from the image data. Many of these approaches use a single spectrum to represent one kind of material and only extract a single endmember spectrum for each endmember class. Techniques in this category either directly extract endmembers from the image, for example, the pixel purity index (PPI) [7], N-FINDR [8], and vertex component analysis (VCA) [9], among many others [10]- [12], or they generate virtual endmembers without assuming the presence of pure signatures in the input data, such as minimum volume-based methods [13]- [15] and nonnegative matrix factorization-based methods [16]- [18]. The major drawback of these methods is that they ignore the endmember variability problem within each endmember class. However, the endmember variability problem is usually unavoidable in real HSIs [19]. For example, illumination differences in the scene can cause shape and magnitude variations within one endmember class. For a scene with large spectral variations, ignoring the endmember variability problem can lead to poor unmixing results. Some solutions developed to solve the endmember variability problem incorporate multiple endmembers within each endmember class [20], [21]. Multiple endmember spectral mixture analysis (MESMA) [22], one of the most widely used and successful methods, selects an optimal endmember combination for each pixel from a spectral library that includes spectral variability. However, there is a heavy computation burden when there is a large number of candidate endmember combinations. To alleviate this problem, MESMA selects a small number of endmember spectra to represent spectral variability within the data instead of assessing all the available endmembers; however, this solution may lead to estimation error of the abundance fractions. More efficient unmixing methods to overcome endmember variability have been presented in the recent literature [19], [23]. Yet the premise of these unmixing methods is the availability of a spectral library that contains endmembers representing the spectral This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ variability within each endmember class. To date, there are only a few methods that can extract multiple endmembers with spectral variability, also called endmember bundles, from the image, indicating the need for new efficient endmember bundle extraction (EBE) methods.
The method developed in [24] divides the global image into spatial subsets and uses a traditional endmember extraction (EE) method such as N-FINDR to extract endmembers within each spatial subset, then integrates the endmembers of all the subsets to obtain the final multiple endmembers. The automated EBE method proposed in [25] sequentially generates a number of subsets by randomly selecting pixels from the observations, extracts endmembers from each subset, then integrates all of the extracted endmembers as the final endmember bundles. The above methods can extract endmember bundles with high efficiency. However, they do not consider spatial information, and they cannot guarantee that endmembers appear in each subset, which may result in the presence of mixed pixels in the extracted endmembers. To enhance the ability of EBE, Xu et al. [26] proposed an image-based EBE method using both spatial and spectral information (SSEBE). SSEBE uses PPI to select candidate endmembers. Since it takes all pixels that have a positive PPI index as candidate endmembers, there is a high probability that the candidate set contains mixed pixels. To solve this problem, SSEBE uses a homogeneity index (HI) to retain the candidate endmembers that are spectrally similar to their spatially adjacent pixels as final endmember bundles. Another work of Xu et al. [27] also utilizes PPI to select candidate endmembers, removing the mixed and redundant endmembers by analyzing the reconstruction error between the chosen endmembers and the remaining candidate endmembers. The effectiveness of these two methods largely depends on the quality of the candidate endmember set. However, the geometric structure of a real HSI in the feature space is not a simple simplex and endmembers may locate within the boundary of the simplex; PPI may fail to extract those endmembers, resulting in an incomplete candidate endmember set.
The convex geometry-based methods proposed in [28] and [29] have the same difficulty of extracting endmembers present within the boundary of the data simplex. The spectral curve-based endmember extraction (SCEE) method [30] obtains spectral curves by processing the original observed data with wavelet transform with different scale factors, and chooses a user-defined number of pixels with maximal or minimal values in each dimension of the curves as candidate endmembers. SCEE uses connected-component labeling to remove mixed pixels from candidate endmembers: an endmember region with an area of more than eight pixels is retained and candidate endmembers that locate outside of the region are removed. SCEE can extract complete candidate endmember sets when the user-defined number of pixels is large enough, which may lead to high redundancy. Furthermore, rare pure pixels may be removed in the mixed pixel removal step.
In recent years, intelligent optimization has been successfully applied in EE. These methods consider EE to be a combinatorial optimization problem and use different strategies to optimize the designed objective functions [31]- [36]. The first work applying intelligent optimization in EE was proposed by Zhang et al. [31], who took the minimization of rootmean-square error (RMSE) between the original image and the reconstructed image as the objective function and used discrete particle swarm optimization to search the optimal endmember combination. Du et al. [37] systematically constructed a quantum behavior-driven particle swarm optimization algorithm to effectively extract endmembers from HSIs. Other strategies, such as ant colony optimization and genetic optimization, have also been employed to minimize the RMSE [33], [36]. It has been proven that intelligent optimization-based methods can obtain results with smaller RMSE than traditional EE methods. Another objective function used in the intelligent optimization-based methods is the maximization of the volume of the simplex constructed by the chosen endmembers. It has been shown that the EE results obtained by these two objective functions are different [38]. In order to get robust results for different real images, multiobjective optimization is used to simultaneously optimize the two objective functions [39], [40]. Although the existing optimization-based EE methods have achieved better results than the traditional EE methods, they do not consider the endmember variability problem, which cannot best fit the situation of real HSIs.
This article proposes a new method to enhance the performance of EBE by leveraging the outstanding optimization ability of intelligent optimization. The proposed method is called multiobjective endmember bundle extraction (MOEBE). For the real HSI, the observed data in the feature space is usually not a simple simplex. It is easy to lose some of the endmembers if the EE is operated in a single dataspace. Considering that different materials have different characteristics among different wavelength ranges, we divide the original data equally into three subsets along its spectral dimension. The feature spaces of these three subsets are different, which means that the distribution of endmembers may differ in the three dataspaces. Even so, the pixels that locate on the vertices of the data simplex for all three subsets are endmembers. The idea of the proposed work is to simultaneously operate EE on the three constructed spaces to obtain multiple endmembers with variability. The simplex volume is used to measure the positions of the endmembers, and the objective function is used to simultaneously maximize the volume of the simplex constructed by endmembers from each subset. A set of Pareto solutions will be obtained by a modified multiobjective particle swarm optimization method [41], and the Pareto solutions can be integrated to produce candidate endmembers. To remove the possible mixed pixels and redundant endmembers from the candidate set, a post-processing step inspired by the technique used in SSEBE [26] as well as the sequential forward floating selection (SFFS) method [42] are utilized to reach the final result. The main contributions of this work can be summarized as follows: 1) Unlike existing EBE methods that only operate EE in a single feature space, the proposed method jointly operates in multiple feature spaces to obtain multiple endmembers. In multiple feature spaces, more endmembers will be located in the vertices of the data simplex, which will benefit the completeness of the extracted multiple endmembers. 2) The multiobjective particle swarm optimization method is modified to better fit the problem. Specifically, the coding of the particles and searching strategy of the population are designed based on the characteristics of the HSI, which helps to find the optimal solution and accelerate the optimization process. 3) Both the endmembers located in the vertices of the data simplex and those located within the boundary of the data simplex can be extracted by taking advantage of the Pareto solutions. As long as the endmember combination is not dominated by other combinations, it is a Pareto solution and they will be taken as candidate endmembers. In fact, this kind of endmember combination can also contain endmembers located within the boundary of the simplex. Therefore, the endmembers located within the boundary of the simplex can be extracted by the proposed method. This article is structured as follows. Section II briefly introduces the linear mixture model (LMM) and multiobjective optimization. Section III gives a detailed description of the proposed method. Section IV describes the comparison experiments between MOEBE and several representative EBE algorithms, with both synthetic HSIs and real data sets. The conclusion is given in Section V.

A. Linear Mixture Model
The LMM [43] is widely used in HU. The majority of existing EE methods are based on the LMM. In the LMM, a mixed pixel is assumed to be the linear combination of all its constituent materials and the corresponding abundance coefficients. Suppose that the image contains N pixels with L spectral bands, and let y = (y 1 , y 2 , . . . , y L ) T be one of the N pixel vectors. Without considering endmember variability, the spectral signature y can be represented by the LMM as follows: where A = [a 1 , a 2 , . . . , a P ] denotes the L × P endmember matrix, with a i = (a 1 , a 2 , . . . , a L ) T being the ith endmember signature, and P is the number of endmembers. The expression s = (s 1 , s 2 , . . . , s P ) T is a P-dimensional vector associated with y, and s i denotes the abundance fraction of the ith endmember present in the pixel y. The term e represents the L × 1 additive observation noise and error vector. The LMM for all the observed pixels can be expressed by the matrix notation as where Y = [y 1 , y 2 , . . . , y N ], S = [s 1 , s 2 , . . . , s N ], and E = [e 1 , e 2 , . . . , e N ]. The abundance is subject to two constraints with the physical meaning: the abundance nonnegative constraint (ANC) and the abundance sum-to-one constraint (ASC), which can be given by s i ≥ 0, i = 1, 2, . . . , P and 1 T s = 1, respectively. If the endmember variability problem is considered, we can use endmember bundles that contain multiple endmember spectra with variability within each endmember class to substitute the corresponding endmember spectra in (1). The endmember bundle matrix is denoted by B in this article.

B. Multiobjective Optimization
Considering the maximization optimization problem, a multiobjective optimization problem can be expressed as where the decision vector z belongs to the feasible solution space and m(≥ 2) conflicting objective functions are to be maximized simultaneously. A decision vector z 2 is said to be dominated by z 1 if A vector z 1 is called Pareto optimal if it is not dominated by any other vectors. Fig. 1 shows the Pareto optimal solutions in the objective space when m = 2. It is obvious that Pareto optimal solutions are nondominated solutions. There is no single optimal solution for multiobjective optimization problems. The results of the multiobjective optimization are a set of Pareto optimal solutions. The corresponding objective vector set of all Pareto optimal solutions is called the Pareto front, and the task of multiobjective optimization is to achieve the Pareto optimal solutions.

C. Evolutionary Algorithm
In the past decades, evolutionary algorithms have attracted increasing interest for the solution of multiobjective optimization problems and a large number of multiobjective evolutionary algorithms (MOEAs) have been developed, including genetic algorithm [44]- [46], differential evolution [47], [48], particle swarm optimization [41], [49], [50], and memetic algorithm [51]- [53]. The goal of MOEA is to reach a good distribution of Pareto solutions with good convergence and diversity. An MOEA usually maintains a population consisting of a set of individuals, where an individual represents a solution to the problem. The individuals are generated by operators and the population is updated in each generation of the evolution to approach the optimal result. Most of the MOEAs mainly focus on three categories. One category contains the decomposition-based algorithms such as those proposed in [54]- [56]. The second category is the indicator-based approach such as [57]- [59]. The Pareto domination approaches [60], [61] belong to the third category.
The competition mechanism-based multiobjective particle swarm optimization (CMOPSO) algorithm [41], which belongs to the Pareto domination approach, is adopted as the basic optimization model for the multiobjective EBE problem due to the high convergence speed and simple implementation. In particle swarm optimization, a particle searches in the feasible solution space by moving along a trajectory depicted by its position and velocity until reaching the optimal. In CMOPSO, a competition mechanism-based learning strategy is used for the updating of particles. In this strategy, particles are pairwise randomly selected from the current swarm for competition. The loser in the competition is updated by learning from the winner, whereas the winner is directly passed to the swarm of next generation. The elite particles, which are used to provide candidate particles to be used in the pairwise competitions to guide the search of the swarm, are selected by the nondominated sorting and crowding distance-based ranking as adopted in nondominated sorting genetic algorithm (NSGA)-II [62].

III. EBE BASED ON MULTIOBJECTIVE PARTICLE SWARM OPTIMIZER
The proposed MOEBE method implements the task of EBE through a modified multiobjective particle swarm optimizer and seeks to find the Pareto optimal solutions of multiple objective functions in order to obtain candidate endmembers with variability. The final results are obtained after removing the mixed pixels and redundant endmembers. In the following, we will introduce the method in detail.

A. Objective Functions
In the convex geometry-based EE methods, volume maximization is quite often utilized to extract endmembers. It assumes that the vertices of the simplex with the largest volume are endmembers. In MOEBE, the original image cube is divided into three subsets along the spectral dimension, as shown in Fig. 2. The volume maximization objective function is applied for each subset as follows: where A i (i = 1, 2, 3) are endmember matrices of each subset. In order to calculate the volume, the dimensionality of endmember matrices A i (i = 1, 2, 3) is reduced to P − 1 by minimum noise fraction (MNF) [63]; A i (i = 1, 2, 3) are the dimensionality reduced matrices. The motivation of dividing the image cube into subsets is that the observed data in the feature space is no longer a simple simplex if there exists endmember variability. In this situation, endmembers may locate within the boundary of the data simplex, which makes the extraction of these endmembers difficult. The division into subsets enables us to operate in multiple spaces, making it possible to solve the variability problem. A set of synthetic pixels is generated to show the mechanism of MOEBE. The endmembers are chosen from the United States Geological Survey (USGS) spectral library, which are montmorillonite, calcite, and topaz, as shown in Fig. 3. Abundance fractions are generated according to the Dirichlet distribution. To display the pixels in feature spaces, dimensionality reduction is implemented on all the pixels by MNF. The structures of both the original data and the three subsets are shown in Fig. 4. It is clear that both the structure of the data and the locations of endmembers are different, which indicates that operating EE in multiple spaces would be useful.

B. MOEBE Procedure
A modified multiobjective particle swarm optimization algorithm based on CMOPSO is employed to optimize the designed objective functions. For the EE problem, the feasible solution space is discrete, accordingly the particle's position must be discrete. Spatial locations of endmembers are used to code the position of particles. The position of a particle can be written as X = (X r1 , X r2 , . . . , X r P , X c1 , X c2 , . . . , X cP ) where X ri (i = 1, 2, . . . , P) is the row number of the ith endmember, and X ci (i = 1, 2, . . . , P) is the column number of the ith endmember. With the row and column numbers, we can accordingly build a combination of endmembers, which is a feasible solution to the EE problem. The velocity of a particle can be written as The competition mechanism used in CMOPSO is not adopted by the proposed method. In MOEBE, we have designed two kinds of particles, one that performs experienced searching and another that performs local searching. The number of particles for the two kinds is the same, and the total number of particles is denoted by M. Assuming the velocity and position of a particle at the current time are V old and X old , and the velocity and position of a particle at the next time are V new and X new . The update of velocity based on the experienced searching is where r 1 and r 2 are random numbers in the interval (0, 1).
The term X b is the historical best solution of the particle. The round() operation rounds each element to the nearest integer. The update of velocity based on the local searching is where the operation of R is to randomly choose values from the integers in the interval [−5, 5]. This operation makes the particle search in a local window with the size 5 × 5. The position of the particle is updated by According to the definition of the particle's position, there should be a lower bound and an upper bound for each element of the position. If the spatial size of the image is r × c, then the lower band is 1 and the upper bound is r for the first P elements, and the lower band is 1 and the upper bound is c for the last P elements. To accelerate the optimization process, the endmembers extracted by VCA are used to initialize a part of the endmembers as [37] did. In order to increase the diversity of the population, the polynomial mutation operator in [64] is used after the update of particles. Nondominated sorting [60] is utilized to choose the optimal solutions. The nondominated solutions are saved in the archive gbest. When the number of solutions in the archive is greater than M, we calculate the product of three objective function values for each solution and only keep the first M solutions with larger values to avoid explosion. Multiobjective particle swarm optimization is an iterative process. The maximum number of iterations max_iter is used as the stopping condition for the optimization process. All the solutions in the final archive constitute the candidate endmember set. The optimization results of independent runs of the multiobjective particle swarm optimization method are different due to the Initialize the best position of each particle X b . 3: For each particle with experienced searching, update the velocity by (8); for each particle with local searching, update the velocity by (9). 4: Update the position of each particle by (10), and set the elements to valid range by X new = mi n(max (X new , X l ), X u ). 5: Perform the polynomial mutation. 6: Calculate the three objective function values for each particle. 7: Do nondominated sorting on the objective function values.
If X new is not dominated by X old , then gbest = gbest ∪ X new . Update X b . 8: If the size of gbest is larger than M, then only keep the first M solutions with larger prod( f i (gbest), i = 1, 2, 3). 9: t = t + 1. 10: If t ≤ max_iter, go to step 3, otherwise output gbest. 11: Find the corresponding endmembers from Y according to gbest, and obtain the endmember candidates B.
intrinsic mechanism of the method. To ensure the stability of the results, ten independent runs are implemented and the candidate endmember set comes from the integration of the results. The procedure of the multiobjective particle swarm optimization method is shown in Algorithm 1.
To show that nondominated solutions can help to extract endmembers with variability, we calculated the objective function values of all the endmember combinations among different endmember classes from Fig. 3; the results are displayed in Table I. The endmembers of montmorillonite, calcite, and topaz are denoted by m, c, and t, respectively, and index values 1, 2, 3 are used to distinguish different endmembers within an endmember class. All the nondominated combinations are displayed in bold. It can be seen that the union set of all nondominated solutions contains all of the nine endmembers, which means that we can successfully achieve the endmember candidate set containing all the endmembers if the nondominated solutions are obtained by multiobjective optimization.

C. Removal of Mixed Pixels and Redundant Endmembers
Spatial post-processing is utilized to remove mixed pixels from the candidate endmembers. MOEBE adopts the postprocessing method used in [26] but uses a slightly different approach. Its steps are: 1) calculate the spectral angle dis- tance (SAD) between each candidate endmember and its spatially neighboring pixels within a 5 × 5 window and keep the maximum SAD for each candidate endmember and 2) show the histogram of the SAD for all candidate endmembers and determine a threshold value τ according to the histogram. The candidate endmembers with maximum SAD larger than the threshold are then removed. After the spatial post-processing, the candidate endmember set is updated to B 1 . The post-processing is a rough screening and can only remove some mixed pixels. In addition to the spatial postprocessing, the SFFS method is utilized to further remove mixed pixels and redundant endmembers. Spatial postprocessing is applied before SFFS because SFFS is time consuming: it saves time if the size of candidate endmember set is reduced in advance.
The aim of the SFFS-based method is to select a subset from the endmember candidate set. Two criterion functions are required to search the optimal subset. One examines the significance of a new endmember, while the other tests if the newly selected endmember can replace some of existing endmembers and form a downsized subset with a higher performance score than that for the same size subset. The criterion function for identifying a new endmember is expressed as where B k is the set of k selected endmembers, B k is the remaining set of endmembers after removing k endmembers Algorithm 2 SFFS-Based Endmember Selection 1: Identify two initial endmembers. The first endmember b 1 is selected as the one in B 1 ∈ R L×K 1 that yields the maximum norm: b 1 = arg max b 2 2 . The second endmember b 2 is the one in B k that has the biggest distance with b 1 : b 2 = arg max b − b 1 2 2 . Update B k and B k . 2: Select the endmember corresponding to the largest J 1 (B k ) value as the candidate of (k+1)th endmember b new and set , then set k = k + 1 and return to step 2, else exclude b j , 1 < j < k from B k+1 to form a new endmember set B k = B k+1 − b j . If k = 2, then set B k = B k and return to step 2, else go to step 4.
k . Return to step 2. 5: The termination condition for step 2 ∼ 4 is: k ≥ 3 and abs(( )) < ω. 6: Output the selected endmembers B 2 that contains in the set B k .
from B 1 . The term K 1 is the number of candidate endmembers in B 1 , b i is the ith endmember in B k , andb i is the reconstructed spectrum of b i by using the endmembers in B k and the nonnegative constrained least square (NCLS) method [65]. The reconstruction error for the ith endmembers in B k is denoted by e i . The criterion function for testing the newly selected endmembers is expressed as whereb i is the reconstructed spectrum of b i by using the endmembers in B k and the NCLS method. The term b i is the ith endmember in B 1 . The procedure of the SFFS-based endmember selection method is described in Algorithm 2.

D. Overall Workflow of MOEBE
The overall workflow of the proposed MOEBE method is shown in Algorithm 3.

IV. EXPERIMENTS AND ANALYSIS A. Data for Experiments
Four data sets were used to validate the proposed method. One of them was simulated by library endmembers and synthetic abundance fractions. The other three were real images: the Samson data set, the Jasper Ridge data set, and the Urban data set.
A synthetic data set was used for experiments because it enabled methods to be precisely validated using known endmembers and abundances. The endmembers used to generate the synthetic data were chosen from the USGS spectral library, Step 1: Division of the input hyperspectral imagery along the spectral dimension The observed data Y ∈ R r×c×L is divided into three subsets: where the number of bands for each subset is the same if L is divisible by 3, otherwise the remainder are added to the third subset.
Step 2: Dimensionality reduction of the subsets Apply the MNF transformation to the input hyperspectral data Y to obtain the transformation matrix. Use the corresponding parts of the transformation matrix to perform dimensionality reduction for Y i , i = 1, 2, 3, generating the data Y i ∈ R r×c×P−1 , i = 1, 2, 3. A row with all elements equal to 1 is added to Y i , i = 1, 2, 3, constructing data Y i ∈ R r×c×P , i = 1, 2, 3 that is used to calculate the volume.
Step 3: Selection of endmembers by MOEBE Perform MOEBE in Algorithm 1 to get the candidate endmembers B ∈ R L×K . Step 4: Removal of mixed pixels and redundant spectra from candidate endmembers Perform spatial post-processing to get a downsized endmember candidate set B 1 ∈ R L×K 1 . Perform SFFS in Algorithm 2 to get the final result B 2 ∈ R L×K 2 .
including 12 endmember classes as shown in Fig. 5. The abundances with size 200 × 200 were generated by the "synthesis tools" package, which is a MATLAB toolbox available online. 1 The abundance maps for the 12 endmember classes are displayed in Fig. 6. To generate synthetic pixels, a single spectrum was randomly selected from multiple endmember spectra for each endmember class, and these selected endmembers were linear combined weighted by the corresponding abundance fractions to produce the synthetic spectral mixture. White Gaussian noise with SNR of 50 dB was added to the synthetic data. Samson is a simple data set. A region of 95 × 95 pixels was used for this experiment, as shown in Fig. 7(a). Each pixel is recorded with 156 spectral bands covering the wavelengths from 401 to 889 nm. There are three main endmember classes in this image, i.e., soil, tree, and water.
The Jasper Ridge data set used for this experiment contains 100 × 100 pixels, as shown in Fig. 7(b). The wavelength ranges from 380 to 2500 nm. After removing the water absorption and noisy bands (1-3, 108-112, 154-166, and 220-224)  from the original 224 bands, the remaining 198 bands were used for the experiment. There are four main endmember classes in this data: road, soil, water, and tree.
Urban is one of the most widely used hyperspectral data used in the HU study and is shown in Fig. 7(c). There are 307×307 pixels. The wavelength ranges from 400 to 2500 nm. After removing the water absorption and noisy bands (1-4, 76, 87, 101-111, 136-153, and 198-210) from the original 210 bands, the remaining 162 bands were used for the experiment. The number of endmember classes was set to 6 for this experiment: asphalt road, grass, tree, roof#1, roof#2, and soil.

B. Performance Metrics
Four metrics, SAD, RMSE, average deviation of the mean (ADM), and average deviation of the standard deviations (ADS), were used for quantitative validation [30]. SAD measures the SAD between each extracted endmember and the reference endmember. RMSE is the root-mean-square error between estimated abundances and true abundances. ADM is the average deviation between the mean of the reference endmembers for each endmember class and the mean of the extracted endmembers for each endmember class. ADS is the average deviation between the standard deviations of the reference endmembers for each endmember class and the standard deviations of the extracted endmembers for each endmember class. Smaller values of the four metrics indicate better results. The true endmembers and abundances of the synthetic data were known. The reference endmembers and abundances of the three real images were provided by Zhu et al. [66]- [68], and are available online. 2 Abundances were estimated by sparse unmixing via the variable splitting and augmented Lagrangian (SUnSAL) method [69] for all the experiments. The extracted multiple endmembers were not clustered into bundles by any clustering algorithm. Instead, they were compared with the reference endmembers by calculating the SAD between the extracted endmember and each reference endmember; the extracted endmember was then assigned to the endmember class with the smallest SAD. The abundance map for one endmember class was the sum of the abundances of all the endmembers within the endmember class.
All four of the metrics were used to evaluate performance for the synthetic experiments, whereas only SAD and RMSE were used for the real image experiments. This was because references of the synthetic data are endmember classes with variability, whereas there is only one reference spectrum for one class for the real data sets.

C. Methods Used for Comparison With MOEBE and Parameter Settings
The EBE [25], SSEBE [26], and SCEE [30] methods were used for comparison with MOEBE. All four methods required parameters. For EBE, the number of subsets and the percentage of the pixel number of subsets against the pixel number of the original image are required parameters. For SSEBE, the block size and the percentage of endmembers in each block are required parameters. For SCEE, the user-defined number of endmember candidates is the required parameter. For MOEBE, the number of particles, the maximum number of iterations, and the threshold for the SFFS method are required parameters. The number of particles and the maximum number of iterations were empirically set to 40 and 2000, respectively, for all the experiments, and only the threshold for the SFFS method was tested and specially determined for each experimental data set. All parameters were determined for each data set by choosing parameters that generated the best result in terms of SAD and ADS for the synthetic data set, and in terms of SAD and RMSE for the three real data sets.

1) Results of Experiment 1 (Using the Synthetic Data Set):
To validate the effectiveness of the velocity and position updating strategies used in MOEBE, the velocity and position updating strategies of the proposed method are replaced with the original ones in CMOPSO (named as MOEBE_CM), and the method has been tested on the synthetic image. The comparison results of SAD, RMSE, ADS, and ADM from MOEBE and MOEBE_CM are shown in Fig. 8. The proposed velocity and position updating strategies produced better results than that of CMOPSO, which demonstrated the effectiveness of the proposed method.
In experiment 1, the endmembers extracted by EBE, SSEBE, SCEE, and MOEBE were evaluated quantitatively with SAD, RMSE, ADS, and ADM under the condition that the true endmembers with variability and abundances were known. Ten independent runs were implemented for each method in order to test the stability of each method. The mean and standard deviations of SAD, RMSE, ADS, and ADM of the ten runs are recorded in Table II, and are also shown in the form of a bar graph in Fig. 9. Positive feedback for the stability of each method was received from the small standard Comparison results of the two velocity and position updating strategies. deviations of the ten runs. In particular, the standard deviations for SCEE were zero, which indicated that the endmembers extracted by SCEE were the same for each run with the same parameters. The average SAD of all the endmember classes for MOEBE was slightly smaller than that for SSEBE, as well as smaller than those for EBE and SCEE. Results from MOEBE had much smaller RMSE, ADS, and ADM than those from EBE, SSEBE, and SCEE. The SAD, RMSE, ADM, and ADS were calculated for each endmember class (see Figs. 10-13, respectively). The SAD, ADM, and ADS of the sixth endmember class for SCEE were empty because SCEE failed to extract any endmembers in this class, which resulted in a large RMSE for this endmember class. Mixed pixels may have been extracted by EBE, resulting in much larger SAD of the ninth endmember class than those of other methods. The ADM and ADS for EBE, SSEBE, and SCEE had big differences among different endmember classes, whereas more balanced ADM and ADS were obtained by MOEBE. MOEBE outperformed other methods in terms of RMSE for each endmember class.   Finally, the mean and standard deviations of computational times for ten runs are shown in Table III. The 64-b version of MATLAB was implemented on a 3.4-GHz Intel Core i7. The results indicate that EBE was much cheaper computationally than other methods.
2) Results of Experiment 2 (Using the Samson Data Set): There was only a single reference endmember spectrum to be compared with the multiple endmembers within an     . 14).
Although MOEBE had the smallest mean SAD among the four methods (see Table IV), SCEE estimated abundances more accurately than the other methods (see Table V). The overall performances of EBE, SSEBE, SCEE, and MOEBE were similar. From the results of computational times (see Table VI),   TABLE VI   COMPUTATIONAL TIMES REQUIRED FOR THE SAMSON DATA SET   TABLE VII   SAD BETWEEN EXTRACTED ENDMEMBERS AND REFERENCES  FOR THE JASPER RIDGE DATA SET   TABLE VIII   RMSE BETWEEN ESTIMATED ABUNDANCES AND REFERENCES  FOR THE JASPER RIDGE DATA SET   TABLE IX  COMPUTATIONAL TIMES REQUIRED FOR THE JASPER RIDGE DATA SET EBE was the most computationally efficient method, and SCEE spent more time than other methods.
3) Results of Experiment 3 (Using the Jasper Ridge Data Set): The estimated abundance maps created by using multiple endmembers extracted by EBE, SSEBE, SCEE, and MOEBE are shown in Fig. 15. The abundance map of tree produced by EBE and soil produced by SSEBE were obviously darker than the reference maps. The abundance maps of water produced by all methods lost some detail, compared with the reference map. From the results of SAD in Table VII and RMSE in Table VIII, MOEBE produced smaller SAD and RMSE than did EBE, SSEBE, and MOEBE. From the results of computational times (see Table IX), EBE was the most computationally efficient method, and MOEBE spent the most time.

4) Results of Experiment 4 (Using the Urban Data Set):
The scene in the Urban data is more complex than that of Samson and Jasper Ridge. In addition to natural ground objects, man-made objects also appear in the Urban data. The estimated abundance maps of EBE, SSEBE, SCEE, and MOEBE are shown in Fig. 16. For EBE, there were obvious differences between the produced and reference abundance maps of asphalt road and soil, and the abundance fractions of roof#2 were overestimated. For SSEBE, abundance maps  of asphalt road, roof#2, and soil were obviously different from those of the reference maps, where the abundance fractions of roof#2 were underestimated for some pixels, and the abundance fractions of asphalt road and soil were overestimated. For SCEE, the abundance map of soil was partly different from that of the reference abundance map, and the map of grass was darker than that of the reference map. For MOEBE, the abundance fractions of tree were overestimated for some pixels. The overall performance of MOEBE was the best visually. SSEBE had the smallest SAD and MOEBE had the smallest RMSE among all the methods (see Tables X and XI). The SAD and RMSE for SSEBE and MOEBE were much smaller than those of EBE and SCEE. From the results of computational times shown in Table XII, EBE was the most computationally efficient method.

E. Discussion
The results of the four experiments indicate that EBE, SSEBE, SCEE, and MOEBE performed well for the two relatively simple images in Samson and Jasper Ridge, in which only endmembers of natural materials were to be extracted and the number of endmember classes was small. When these four methods were applied to extract multiple endmembers for the synthetic and Urban images, some of the results were not satisfactory. Although the synthetic image has only a simple type of noise (no noise or error from the imaging process  and no anomaly) compared with the real images, it contains 12 classes of endmembers with 11 kinds of minerals and the vegetation class, which increases the complexity of the image. In this experiment, some extracted endmembers from EBE had a large SAD from the reference endmembers. This showed that mixed pixels may have been extracted as endmembers by EBE. EBE extracts endmembers from randomly selected subsets, so the mixed pixels could be mistakenly extracted as endmembers once the endmember of one class is absent in the subset. The ADM and ADS of some endmember classes were large for SSEBE, which indicated that some endmembers within the endmember class may have been lost by SSEBE. SCEE failed to extract one class of endmembers. Because SCEE utilizes the connected region to remove mixed pixels from candidate endmembers, true endmembers may also be removed if they are not located in the connected region. This is not good for the identification of rare endmembers. MOEBE performed well in extracting multiple endmembers for each of the endmember classes, and the variability of spectra within each class was close to the true variability. This is an indicator of MOEBE's potential. The Urban image has only six endmember classes, but it is affected by various noises or interferences and contains both natural and man-made materials, which result in a complex scene. EBE and SSEBE performed poorly in this experiment, which implies that the multiple spectra extracted by these two methods were not able to represent true endmember variability well for this image. EBE may have introduced mixed pixels in the endmember set. SSEBE failed to extract the endmembers located within the boundary of data simplex, leading to inaccurate variability within the endmember class. SCEE and MOEBE were able to provide multiple endmembers describing the variability within the endmember class more accurately. EBE had very high efficiency in computational time for the four experiments. For SSEBE and SCEE, the computational time increased as the image size increased. The PPI method used in SSEBE conducts a large number of random projections for the pixels, so it will take more time to do the projections if the number of pixels increases. Since one step of SCEE applies the wavelet transform to the input image with several different scale factors, an image of a larger size will cost more time. In MOEBE, it only takes tens of seconds to get the candidate endmembers and the image size has little effect on the time cost. The most time-consuming part for MOEBE is the mixed pixel removal by SFFS. It will spend more time to get the result if the number of extracted endmembers increases.

V. CONCLUSION
A novel method, MOEBE, which extracts multiple endmembers with variability within endmember classes has been demonstrated. MOEBE constructs subspaces of the original image, and simultaneously searches endmembers in multiple spaces by using multiobjective particle swarm optimization. The comparison of MOEBE with EBE, SSEBE, and SCEE using synthetic data showed that MOEBE obtained more accurate results than did the other methods. Multiple endmember spectra obtained by MOEBE represented the true variability well for each endmember class. The abundance maps generated with multiple endmembers extracted by MOEBE were reliable under the complex scene of a real image.