Multiobjective Optimized Endmember Extraction for Hyperspectral Image

Endmember extraction (EE) is one of the most important issues in hyperspectral mixture analysis. It is also a challenging task due to the intrinsic complexity of remote sensing images and the lack of priori knowledge. In recent years, a number of EE methods have been developed, where several different optimization objectives have been proposed from different perspectives. In all of these methods, only one objective function has to be optimized, which represents a specific characteristic of endmembers. However, one single-objective function may not be able to express all the characteristics of endmembers from various aspects, which would not be powerful enough to provide satisfactory unmixing results because of the complexity of remote sensing images. In this paper, a multiobjective discrete particle swarm optimization algorithm (MODPSO) is utilized to tackle the problem of EE, where two objective functions, namely, volume maximization (VM) and root-mean-square error (RMSE) minimization are simultaneously optimized. Experimental results on two real hyperspectral images show the superiority of the proposed MODPSO with respect to the single objective D-PSO method, and MODPSO still needs further improvement on the optimization of the VM with respect to other approaches.


Introduction
Each pixel of hyperspectral image (HSI) has tens or hundreds of values corresponding to its spectral bands, which can effectively represent the unique ground objects [1,2].Hyperspectral images have been successfully applied to a wide range of fields [3].However, mixed pixels, constituting more than one distinct material, may widely exist in the HSI due to the limited spatial resolution, which makes one single pixel not pure and brings troubles to accurate precision analysis of HSIs [4][5][6].Spectral unmixing (SU) is an effective technique to resolve the mixed pixels problem, which decomposes the mixed pixels into a collection of pure materials, named endmembers, as well as the corresponding abundances [7].SU has two tasks: EE and abundance estimation.It is usually assumed that there are some pixels that contain only one kind of ground object in the image, and EE is to find out such pure pixel for basic ground objects [8].Abundance estimation is the process to estimate different proportion of each endmember in a mixed pixel.This paper mainly focuses on the task of EE.
The studies of mixed pixels are mostly based on the linear mixture model (LMM) in which each observed pixel in the image can be represented as the linear combination of a set of spectrally pure constituent endmembers, weighted by the corresponding abundance coefficients that establish the proportion of each endmember in the pixel [9].Under the LMM, assuming that the image scene is dominated by P kinds of distinct materials with L bands, mathematically, a pixel vector y = [y 1 , y 2 , • • • , y L ] T can be written as: where A = [a 1 , a 2 , • • • , a P ] is an L × P endmember matrix, with each column being an endmember signature vector.The number of endmember P is a pre-defined parameter, which can be estimated by existing methods, and the commonly used ones are the virtual dimensionality (VD) estimation method [10] and the hyperspectral subspace identification (HySime) method [11].s = [s 1 , s 2 , • • • , s P ] T is a P-dimensional column vector composed of abundance coefficients of the corresponding endmembers for the pixel, and e represents the L × 1 additive observation noise and error vector.Generally, there are various kinds of noise in HSIs, and this work assumed that the error is represented by the additive white Gaussian noise [12].The LMM for all N observed pixels can be expressed by the matrix notation: where , and E = [e 1 , e 2 , • • • , e N ].Due to physical constraints, the abundance vector is subject to the nonnegative constraint (ANC, s i ≥ 0, i = 1, 2, • • • , P) and the abundance sum-to-one constraint (ASC, 1 T s = 1).With the LMM, the geometrical interpretation of the HSI is that if e = 0 and there are pure pixels of all kinds of materials in the image (pure pixel assumption), all the pixels are contained in a simplex whose vertices are corresponding to the endmembers [13].Based on the convex geometry theory, the EE problem can be converted into finding the simplex vertices.Typical methods include the pixel purity index (PPI) [14], N-FINDR [15], the simplex growing algorithm (SGA) [16], vertex component analysis (VCA) [17], as well as some new algorithms proposed in recent years [18][19][20][21][22].However, the classic algorithms such as N-FINDR and VCA have been shown easily affected by noise and outliers [23].One progress in recent years lies on the intelligent optimization methods to enhance the EE results in real HSIs [8,[23][24][25][26][27].Most of these algorithms [8,[23][24][25] consider the EE problem as a combination optimization problem, and seek the optimal endmember combination that minimizes the root-mean-square error (RMSE) between the original image and its remixed image.It is showed that intelligent optimization methods such as D-PSO can get a smaller RMSE compared to N-FINDR and VCA [24].Different from the above methods, the MOAQPSO method in [26] takes the VM as the objective function, and the experimental results showed the conflicts between the RMSE minimization and the VM objective functions.Specifically, the two objective functions did not achieve their best values for the same endmember combination.If one method got the optimal endmember combination in terms of the volume value, then there would be another method superior to it in terms of the RMSE value.From the previous studies [23,26,28], although the RMSE minimization-based methods can get superior results than the VM-based methods in terms of the RMSE value (or the VM-based methods can get superior results than the RMSE minimization-based methods in terms of the volume value), neither of them can prove completely superior to the other when comparing each one of those endmember spectra with the reference.The VM-based methods have an obvious advantage over the RMSE minimization-based methods when extracting rare endmembers, while in VM-based methods the noises and outliers located within the bounds of the data simplex may be identified incorrectly as endmembers; the RMSE minimization-based methods are more robust to noises and outliers, while they usually ignore the rare endmembers.Effective EE results can be achieved if there is a good match between the characteristic expressed by the objective function and the characteristic of the real image.However, no a priori knowledge is provided in practice, and different complex hyperspectral remote sensing images usually have different characteristics.It is concluded that the generalization ability of one single objective function is poor, and it may not be enough to provide satisfactory EE results for various images.Hence, it is natural to simultaneously optimize several objective functions so as to capture the different data characteristics.In this article, the two widely used objective functions, the VM [13] and the RMSE minimization [24] objective functions, are integrated to be simultaneously optimized.In this way, the problem of endmember extraction is transformed into multiobjective optimization (MOO) problem.Some MOO methods have been suggested to solve various multiobjective problems [29][30][31][32], in which the PSO-based MOO methods have attracted a lot of attention, and this kind of method is chosen as the optimization method of this paper due to the simplicity and good search ability of PSO.Although the existing MOO methods have provided us some ideas on how to solve MOO problems, to our knowledge, no previous works are reported for EE problem, and the difficulty lies on that the distribution characteristics of search spaces and solution spaces of different problems are usually different, so existing methods that are effective for other problems may not work while solving the EE problem.In this paper, a multiobjective discrete particle swarm optimization algorithm (MODPSO) is proposed to perform the task of EE for hyperspectral images.The work includes three aspects: (1) The update strategy of particles' velocity and position in D-PSO method [24] is selected as the basic searching strategy for the proposed MOO method.Since the search space and solution space of the EE problem are both discrete, the particle's position and the velocity must also be discrete to ensure the validity of the solution, so the update strategy of the particles should be modified to make it suitable to the EE problem.(2) For the proposed MOO method, the two objective functions often conflict with each other during the process of optimization, which means that finding a solution that optimizes both objective functions at the same time is almost impossible during the process of optimization [33].This brings a trouble for the acquisition of the particle's personal best position (pbest) and the population's global best position (gbest) in the multiobjective searching space.The nondominated sorting algorithm [34] is used to determine pbest and gbest according to the multiobjective function values.
(3) Different from the single objective optimization, there is more than one gbest for the population in the MOO, and all of the non-dominated solutions are gbests.This brings the problem to determine which gbest should be chosen when updating the velocity of each particle.To solve the problem, the Sigma method is utilized to find best local guides for each particle of the population [35].With all of the above works, the multiobjective discrete particle swarm optimization algorithm (MODPSO) is finally formed to perform the task of endmember extraction for hyperspectral images.Like common EE methods, MODPSO is based on the pure pixel assumption, and needs the number of endmembers as a priori parameter.
As far as we know, this is the first attempt to use MOO for the purpose of EE.The remainder of this paper is organized as follows.Section 2 gives a detailed description of the proposed MODPSO algorithm for EE.Section 3 reports the experimental results of the MODPSO method and several representative single objective optimization EE algorithms.Conclusions are drawn in Section 4.

MODPSO
The proposed MODPSO method implements the task of EE through a MOO technique, and it aims at finding the Pareto-optimal solutions for simultaneously optimizing multiple objective functions.Hence, the establishment of the objective functions and the optimization strategy for the multiple objective functions are two key elements of the MODPSO method.In the following, we will introduce them in detail.

Objective Functions for MODPSO
Two kinds of objective functions are elaborately chosen for the proposed algorithm.One is the maximum volume objective function, which is based on the convex geometry theory, and the other is to minimize the RMSE obtained after reconstructing the hyperspectral scene by only assuming the presence of the additive white Gaussian noise [12], like almost all the unmixing methods do [36].
We have transformed the VM into minimizing the volume inverse, so the two objective functions are both minimization problems.The two objectives are listed below: where L is the spectral dimensionality of the HSI, N is the total number of pixels, P is the number of endmembers, A is the endmember matrix, and Y and Ŷ are the original image and the remixed image, respectively.The abundances used to calculate Ŷ are estimated by Equation ( 5) rather than the fully constrained least squares method (FCLS) for the sake of efficiency: In most cases, these two objective functions will not obtain their optimal solution for the same combination of endmembers, for considering that there usually exists noise or outlier in real hyperspectral images.The strategy for optimizing the multiple objective functions in MODPSO will help to find a number of endmember combinations, and none of the obtained solutions can be further improved on the objective value without degrading another.

The Updating Strategy of the Particle's Velocity and Position in MODPSO
MODPSO use particles to search in the feasible solution space.Each particle has two properties: the position and the velocity.A particle moves along a trajectory depicted by its position and velocity in the search space, to find an optimal solution.For the EE problem, the search space is discrete, the particle's position and the velocity must also be discrete to ensure the validity of the solution.The Binary coding method used in the D-PSO method is employed here to make particles be able to search in the discrete feasible solution space.The position of the ith particle at iteration time t can be written as: where x j = 1 if y i ∈ A and x j = 0 if not.Explicitly, for the position of the ith particle X t i , all the elements of it are composed of 0 and 1, and each element x j (j = 1, . . ., N) in it represents the attribute of the corresponding pixel y j (j = 1, . . ., N), if the value of x j is 1, the pixel y j is selected as an endmember; otherwise, the pixel y j is not selected as an endmember.Hence, P elements in each particle's position are 1, and the remaining elements are 0.
V t i is used to specify the ith particle's velocity at time t.pbest t i and gbest t are used to specify the ith particle's personal best position and all population's global best position in history before time t.The updating functions of position and velocity are: where T and R are both random selection functions.The velocity obtained by T is based on self-experience and social experience, while R generates velocity without considering past experiences.Both T(X) and R(X) are vectors with the same dimension of X, and the calculation for them can be divided into three steps: (1) Predefine a random selection probability p, and randomly generate a number between 0 and 1. (2) If the generated number is greater than or equal to p, select T to obtain the velocity.First, we find the positive elements and the negative elements of X, respectively.Then, randomly select one element from all the positive elements of X, and set the element of T(X) with the same position of this randomly select one to 1. Next, randomly select another element from all the negative elements of X, and set the element of T(X) with the same position of this randomly select one to −1.The final velocity is obtained by setting the rest of the elements of T(X) to 0. (3) If the generated number is less than p, select R to obtain the velocity.First, we find the zero elements and the positive elements of X respectively.Then, randomly select one element from all the zero elements of X, and set the element of R(X) with the same position of this randomly select one to 1. Next, randomly select another element from all the positive elements of X, and set the element of R(X) with the same position of this randomly select one to -1.The final velocity is obtained by setting the rest of the elements of R(X) to 0. The acquisition of pbest t i and gbest t will be introduced in the following part.

Strategy for Updating pbest and gbest for Optimizing the Multiple Objective Functions
Considering the minimization optimization problem, a MOO problem is of the form: where the decision vectors z belong to the feasible space formed by some constraint functions.m(≥ 2) conflicting objective functions are to be minimized simultaneously.A decision vector z 1 is said to dominate A vector z 1 is called Pareto-optimal if another z 2 that dominates it does not exist.Figure 1 shows the Pareto-optimal solutions when m = 2.There is no single optimal solution in MOO, but a set of optimal solutions.The set containing all the optimal solutions is known as the Pareto front, and the task of MOO is to achieve the Pareto front.It is obvious that the solutions in the Pareto front are non-dominated solutions.
Remote Sens. 2017, 9, 558 5 of 17 of X, respectively.Then, randomly select one element from all the positive elements of X, and set the element of ( ) with the same position of this randomly select one to 1. Next, randomly select another element from all the negative elements of X, and set the element of ( ) with the same position of this randomly select one to −1.The final velocity is obtained by setting the rest of the elements of ( ) to 0. (3) If the generated number is less than p, select to obtain the velocity.First, we find the zero elements and the positive elements of X respectively.Then, randomly select one element from all the zero elements of X, and set the element of ( ) with the same position of this randomly select one to 1. Next, randomly select another element from all the positive elements of X, and set the element of ( ) with the same position of this randomly select one to -1.The final velocity is obtained by setting the rest of the elements of ( ) to 0. The acquisition of and will be introduced in the following part.

Strategy for Updating and for Optimizing the Multiple Objective Functions
Considering the minimization optimization problem, a MOO problem is of the form: where the decision vectors belong to the feasible space formed by some constraint functions.
(≥ 2) conflicting objective functions are to be minimized simultaneously.A decision vector is said to dominate if: 1, 2, , , A vector is called Pareto-optimal if another that dominates it does not exist.Figure 1 shows the Pareto-optimal solutions when = 2.There is no single optimal solution in MOO, but a set of optimal solutions.The set containing all the optimal solutions is known as the Pareto front, and the task of MOO is to achieve the Pareto front.It is obvious that the solutions in the Pareto front are non-dominated solutions.One main step in MODPSO is to determine the personal and global best positions.They are easy to be determined in single objective optimization by selecting the position best fits the objective function.However, in MOO problems, it is hard to determine which position is better if the solutions represented by two positions are not dominated by each other.To handle this problem, the nondominated sorting algorithm [21] is used to update and .Among the population, different particles are compared by the concept of Pareto domination.If the solution of one particle is not dominated by that of all the other particles, then it is a Pareto-optimal solution.
It should be noted that there is not only one in MOO, all non-dominated solutions in the optimization process are taken as .For the update of , all the pairwise comparisons of the solutions are conducted by Pareto domination after each iteration, and all the non-dominated One main step in MODPSO is to determine the personal and global best positions.They are easy to be determined in single objective optimization by selecting the position best fits the objective function.However, in MOO problems, it is hard to determine which position is better if the solutions represented by two positions are not dominated by each other.To handle this problem, the nondominated sorting algorithm [21] is used to update pbest and gbest.Among the population, different particles are compared by the concept of Pareto domination.If the solution of one particle is not dominated by that of all the other particles, then it is a Pareto-optimal solution.
It should be noted that there is not only one gbest in MOO, all non-dominated solutions in the optimization process are taken as gbest.For the update of gbest, all the pairwise comparisons of the solutions are conducted by Pareto domination after each iteration, and all the non-dominated solutions are kept as gbest.A set named global best archive (GBA) is used to store all these non-dominated solutions (gbest).
For the update of pbest, the newly generated particle's position , as shown in Figure 2. If the particle appears in the area where the purple point located, then the pbest of the particle should remain unchanged; if the particle appears in the area where the red point located, then the pbest of the particle will be updated by the red point; and if the particle appears in the area where the cyan points located, then randomly select one point as the pbest.

Choose the Best Local Guide for Each Particle
In single objective optimization, there is only one for the population, so all of the particles use the same to generate the new velocity.However, we have stated that there is not only one in MOO, all the solutions in GBA are taken as .This brings an additional problem of which solution in GBA should be used to generate the velocity for each particle.To solve this problem, The Sigma method [35] is utilized to select one best local guide from GBA for the th particle.In the Sigma method, a value σ is assigned to each point , , , , and the σ value is defined as: According to Equation (10), all the points on the line = have the same σ values.By considering the objective space, finding the best local guide among GBA for the particle i at iteration time t is as follows: in the first step, the values of each position in GBA is assigned.In the second step, for particle i is calculated.Then the distances between and all the values of GBA are calculated.Finally, the kth position in GBA which has the minimum value distance with particle i is selected as the best local guide .Figure 3 shows this method for a two objective example.If the particle appears in the area where the purple point located, then the pbest of the particle should remain unchanged; if the particle appears in the area where the red point located, then the pbest of the particle will be updated by the red point; and if the particle appears in the area where the cyan points located, then randomly select one point as the pbest.

Choose the Best Local Guide for Each Particle
In single objective optimization, there is only one gbest for the population, so all of the particles use the same gbest to generate the new velocity.However, we have stated that there is not only one gbest in MOO, all the solutions in GBA are taken as gbest.This brings an additional problem of which gbest solution in GBA should be used to generate the velocity for each particle.To solve this problem, The Sigma method [35] is utilized to select one best local guide gbest t i from GBA for the ith particle.In the Sigma method, a value σ i is assigned to each point ( f 1,i , f 2,i ), and the σ value is defined as: According to Equation ( 10), all the points on the line f 2 = a f 1 have the same σ values.By considering the objective space, finding the best local guide gbest t i among GBA for the particle i at iteration time t is as follows: in the first step, the σ values of each position in GBA is assigned.In the second step, σ i for particle i is calculated.Then the distances between σ i and all the σ values of GBA are calculated.Finally, the kth position in GBA which has the minimum σ value distance with particle i is selected as the best local guide gbest t i .Figure 3 shows this method for a two objective example.The sigma values of all the GBA members and particles are calculated and compared.For one particle, the GBA member that has the closest sigma value with it is chosen as the best local guide for it.

The Framework of MODPSO for EE
The overall process of the proposed MODPSO for EE is shown in Figure 4.

Experiments
Two real HSIs are used to test the performance of the proposed method.N-FINDR [15], VCA [17] and D-PSO [24] are comparison algorithms.There are two reasons for selecting these three algorithms as comparing algorithms.One reason is that N-FINDR and VCA are two of the most The squares stand for the GBA members, and circles stand for all the particles.The sigma values of all the GBA members and particles are calculated and compared.For one particle, the GBA member that has the closest sigma value with it is chosen as the best local guide for it.

The Framework of MODPSO for EE
The overall process of the proposed MODPSO for EE is shown in Figure 4.
Figure 3. Selection of the best local guide among the global best archive (GBA) for each particle.The squares stand for the GBA members, and circles stand for all the particles.The sigma values of all the GBA members and particles are calculated and compared.For one particle, the GBA member that has the closest sigma value with it is chosen as the best local guide for it.

The Framework of MODPSO for EE
The overall process of the proposed MODPSO for EE is shown in Figure 4.

Experiments
Two real HSIs are used to test the performance of the proposed method.N-FINDR [15], VCA [17] and D-PSO [24] are comparison algorithms.There are two reasons for selecting these three algorithms as comparing algorithms.One reason is that N-FINDR and VCA are two of the most

Experiments
Two real HSIs are used to test the performance of the proposed method.N-FINDR [15], VCA [17] and D-PSO [24] are comparison algorithms.There are two reasons for selecting these three algorithms as comparing algorithms.One reason is that N-FINDR and VCA are two of the most popular EE methods, and D-PSO is a representative method of the intelligent optimization based methods.Furthermore, the objective functions used in MODPSO have been used in these three methods, so the validity of the proposed method can be checked by comparing the objective values of these methods.For both D-PSO and MODPSO, the maximum iteration number was set to 300, the number of particles was set to 20, the random selection probability was 0.2, and the particles were randomly initialized.

HYDICE Washington DC Dataset
The first real image dataset collected by the Digital Imagery Collection Experiment (HYDICE) sensor over Washington DC, and a subset of 150 × 150 was extracted from the original image for this experiment.In the Washington DC dataset, there are 210 bands, which cover the range of 0.4-2.5 um.Low-SNR and water-vapor absorption bands were removed before unmixing, leaving 187 bands for the experiment.Figure 5 shows the false-color image composed of R-band 64, G-band 52, and B-band 36.There are six distinct materials in the image [37], so the endmember number is set to six.
popular EE methods, and D-PSO is a representative method of the intelligent optimization based methods.Furthermore, the objective functions used in MODPSO have been used in these three methods, so the validity of the proposed method can be checked by comparing the objective values of these methods.For both D-PSO and MODPSO, the maximum iteration number was set to 300, the number of particles was set to 20, the random selection probability was 0.2, and the particles were randomly initialized.

HYDICE Washington DC Dataset
The first real image dataset was collected by the Hyperspectral Digital Imagery Collection Experiment (HYDICE) sensor over Washington DC, and a subset of 150 × 150 was extracted from the original image for this experiment.In the Washington DC dataset, there are 210 bands, which cover the range of 0.4-2.5 um.Low-SNR and water-vapor absorption bands were removed before unmixing, leaving 187 bands for the experiment.Figure 5 shows the false-color image composed of R-band 64, G-band 52, and B-band 36.There are six distinct materials in the image [37], so the endmember number is set to six. Figure 6 shows the objective function value as a function of the number of iteration times of MODPSO for the Washington DC dataset.It can be seen that the proposed method can converge to a stationary point when reached the maximum iteration.Figure 6 shows the objective function value as a function of the number of iteration times of MODPSO for the Washington DC dataset.It can be seen that the proposed method can converge to a stationary point when reached the maximum iteration.Figure 7a shows the obtained GBA by MODPSO, ten non-dominated solutions are finally remained.The number of solutions in GBA is less than the number of particles, which indicates that some different particles converge to the same solution.We can see in these results that no solution has the minimum and values simultaneously.The non-dominated solution with the minimum has the largest and vice versa.The ten results are uneven distributed in the objective function value space, they can be easily divided into three parts: three solutions have relatively bigger and smaller , four solutions have relatively bigger and smaller , the remaining three solutions have the best tradeoff between the two objective functions.We have also calculated and values of the other three methods according to their extracted endmembers.The results are put together with that of MODPSO in Figure 7b, and the numerical results are shown in Table 1 as well as the computation time of them.It can be seen that solutions of MODPSO dominate the result of D-PSO, so the search ability of MODPSO is better than that of D-PSO.VCA and N-FIDNR achieved smaller and larger than MODPSO, which indicates that the results with bigger volume are obtained by VCA and N-FINDR, while the RMSE generated by them is larger.Since the results by VCA and N-FINDR are non-dominated solutions compared with the results by MODPSO, it tells us that the Pareto front found by MODPSO is not completed.In terms of the computation time, the two conventional EE methods are more efficient than the two intelligent optimization based methods, especially for the VCA method.
(a) Figure 7a shows the obtained GBA by MODPSO, ten non-dominated solutions are finally remained.The number of solutions in GBA is less than the number of particles, which indicates that some different particles converge to the same solution.We can see in these results that no solution has the minimum f 1 and f 2 values simultaneously.The non-dominated solution with the minimum f 1 has the largest f 2 and vice versa.The ten results are uneven distributed in the objective function value space, they can be easily divided into three parts: three solutions have relatively bigger f 2 and smaller f 1 , four solutions have relatively bigger f 1 and smaller f 2 , the remaining three solutions have the best tradeoff between the two objective functions.We have also calculated f 1 and f 2 values of the other three methods according to their extracted endmembers.The results are put together with that of MODPSO in Figure 7b, and the numerical results are shown in Table 1 as well as the computation time of them.It can be seen that solutions of MODPSO dominate the result of D-PSO, so the search ability of MODPSO is better than that of D-PSO.VCA and N-FIDNR achieved smaller f 1 and larger f 2 than MODPSO, which indicates that the results with bigger volume are obtained by VCA and N-FINDR, while the RMSE generated by them is larger.Since the results by VCA and N-FINDR are non-dominated solutions compared with the results by MODPSO, it tells us that the Pareto front found by MODPSO is not completed.In terms of the computation time, the two conventional EE methods are more efficient than the two intelligent optimization based methods, especially for the VCA method.Figure 7a shows the obtained GBA by MODPSO, ten non-dominated solutions are finally remained.The number of solutions in GBA is less than the number of particles, which indicates that some different particles converge to the same solution.We can see in these results that no solution has the minimum and values simultaneously.The non-dominated solution with the minimum has the largest and vice versa.The ten results are uneven distributed in the objective function value space, they can be easily divided into three parts: three solutions have relatively bigger and smaller , four solutions have relatively bigger and smaller , the remaining three solutions have the best tradeoff between the two objective functions.We have also calculated and values of the other three methods according to their extracted endmembers.The results are put together with that of MODPSO in Figure 7b, and the numerical results are shown in Table 1 as well as the computation time of them.It can be seen that solutions of MODPSO dominate the result of D-PSO, so the search ability of MODPSO is better than that of D-PSO.VCA and N-FIDNR achieved smaller and larger than MODPSO, which indicates that the results with bigger volume are obtained by VCA and N-FINDR, while the RMSE generated by them is larger.Since the results by VCA and N-FINDR are non-dominated solutions compared with the results by MODPSO, it tells us that the Pareto front found by MODPSO is not completed.In terms of the computation time, the two conventional EE methods are more efficient than the two intelligent optimization based methods, especially for the VCA method.
(a)  For the HYDICE Washington DC dataset, the ground features are easy to distinguish by visual interpretation; we manually select the endmembers of the six kinds of materials from the image by referring to [37].These spectra are taken as a rough reference to be shown together with the extracted spectra.The extracted endmembers by the four algorithms and manually selected reference spectra are shown in Figure 8, where the shown endmembers by MODPSO are the union set of the endmembers in GBA.Among the ten sets of results in GBA, there are twelve different endmember spectra.We can see from Figure 8 that N-FINDR and VCA missed the street's spectra and extracted two different paths' spectra.The spectral shapes of the extracted endmembers by the four methods are similar to the shapes of the reference spectra, while there are some differences in the scale, the endmember spectra by D-PSO and MODPSO are more close to the manually selected reference spectra than that of N-FINDR and VCA.For the HYDICE Washington DC dataset, the ground features are easy to distinguish by visual interpretation; we manually select the endmembers of the six kinds of materials from the image by referring to [37].These spectra are taken as a rough reference to be shown together with the extracted spectra.The extracted endmembers by the four algorithms and manually selected reference spectra are shown in Figure 8, where the shown endmembers by MODPSO are the union set of the endmembers in GBA.Among the ten sets of results in GBA, there are twelve different endmember spectra.We can see from Figure 8 that N-FINDR and VCA missed the street's spectra and extracted two different paths' spectra.The spectral shapes of the extracted endmembers by the four methods are similar to the shapes of the reference spectra, while there are some differences in the scale, the endmember spectra by D-PSO and MODPSO are more close to the manually selected reference spectra than that of N-FINDR and VCA.

HYDICE Urban Dataset
The second real dataset was the Urban HYDICE HSI, as shown in Figure 9

HYDICE Urban Dataset
The second real dataset was the Urban HYDICE HSI, as shown in Figure 9 by R-band 64, G-band 52, and B-band 36.This image is of size 307 × 307 and has 210 spectral bands in the range of 0.4-2.5 um.A total of 162 bands remained after removing bands 1-4, band 76, band 87, band 111, bands 101-111, bands 136-153 and bands 198-210.The number of endmembers is set to six [38].The Pareto front by MODPSO is displayed in Figure 11a.Eleven non-dominated solutions finally remained, which indicates that some different particles converge to the same solution.We can see that a more uniform distribution of the non-dominated solutions is obtained by the Urban image than the Washington DC image.The comparison results of four methods in Figure 11b and Table 2 are similar to that of the Washington DC image: most of solutions of MODPSO dominate the result of D-PSO; MODPSO and D-PSO have results with smaller RMSE than VCA and N-FINDR; and VCA and N-FIDNR obtained bigger volume than MODPSO and D-PSO, which demonstrate the validity of the MODPSO method.Meanwhile, the MOO result can be further improved.Seen from the computation time, VCA is the most efficient method, while MODPSO and D-PSO are both time consuming.The Pareto front by MODPSO is displayed in Figure 11a.Eleven non-dominated solutions finally remained, which indicates that some different particles converge to the same solution.We can see that a more uniform distribution of the non-dominated solutions is obtained by the Urban image than the Washington DC image.The comparison results of four methods in Figure 11b and Table 2 are similar to that of the Washington DC image: most of solutions of MODPSO dominate the result of D-PSO; MODPSO and D-PSO have results with smaller RMSE than VCA and N-FINDR; and VCA and N-FIDNR obtained bigger volume than MODPSO and D-PSO, which demonstrate the validity of the MODPSO method.Meanwhile, the MOO result can be further improved.Seen from the computation time, VCA is the most efficient method, while MODPSO and D-PSO are both time consuming.The Pareto front by MODPSO is displayed in Figure 11a.Eleven non-dominated solutions finally remained, which indicates that some different particles converge to the same solution.We can see that a more uniform distribution of the non-dominated solutions is obtained by the Urban image than the Washington DC image.The comparison results of four methods in Figure 11b and Table 2 are similar to that of the Washington DC image: most of solutions of MODPSO dominate the result of D-PSO; MODPSO and D-PSO have results with smaller RMSE than VCA and N-FINDR; and VCA and N-FIDNR obtained bigger volume than MODPSO and D-PSO, which demonstrate the validity of the MODPSO method.Meanwhile, the MOO result can be further improved.Seen from the computation time, VCA is the most efficient method, while MODPSO and D-PSO are both time consuming.
The extracted endmembers and manually selected reference spectra of the Urban image are shown in Figure 12.Half of the endmembers extracted by N-FINDR and VCA and one endmember extracted by MODPSO cannot be matched with the manually selected reference spectra.Only the N-FINDR algorithm extracted the sixth endmember, and the spectrum of the endmember is not so close to that of the reference endmember.In general, the endmembers extracted by D-PSO and MODPSO are better matched than that of N-FINDR and VCA.

Review of Experimental Results
Experimental results of the Washington dataset showed that N-FINDR and VCA failed to extract the fourth endmember, which resulted in a larger RMSE than the other two methods, while the volumes obtained by N-FINDR and VCA were much larger than the other two methods.Considering the fact that one of the objective functions of MODPSO is VM, we can infer that the searching ability of MODPSO needs further improvement.In term of RMSE, MODPSO was superior to the other methods.Experimental results of the Urban dataset showed that N-FINDR and VCA extracted several outliers, which indicated that they were easy to be affected by outliers.MODPSO and D-PSO were more robust to these interferences.We can infer that the RMSE objective function can play a key role when there are interferences in the image.In both experiments, MODPSO can find better solution than D-PSO.This may because the MOO mechanism of MODPSO (several gbest in MODPSO compared to one in D-PSO) increased the diversity of particles and alleviated the premature convergence problem of D-PSO, thus leading to a better optimization result.Time costs of the methods showed that the two intelligent optimization based methods were time consuming, which mainly resulted from the calculation of the RMSE objective function.

Generalization of MODPSO
In this work, MODPSO assumed that the error is represented by the additive white Gaussian noise.In fact, there may have mixed noise in the HSI such as impulse noise, multiplicative noise or vertical line strips [36,39].It should be noted that the MODPSO method can also be applied when considering other types of noise, as long as an objective function is built according to a certain type of noise or mixed noise, the RMSE function can be replaced by the newly built one.

Conclusions and Future Work
This paper proposed a multiobjective optimization method MODPSO for endmember extraction.In MODPSO, the volume maximization and RMSE minimization objective functions are simultaneously optimized, and the multiobjective optimization framework is especially designed to solve the multiobjective endmember extraction problem.Instead of obtaining one unique solution for one implementation like other endmember extraction methods, the result by MODPSO is a set of non-dominated solutions, and they can be regarded as solutions with different tradeoffs between two objective functions.The experimental results show that the search ability of MODPSO method is superior to that of the D-PSO method, and it can obtain result with smaller RMSE than N-FINDR and VCA.However, the results of N-FINDR and VCA are not dominated by that of MODPSO for the reason that the volume obtained by them is bigger than that of MODPSO, which indicates that the Pareto front obtained by MODPSO is not complete, a part of non-dominated solutions are not founded by them, which revealed the limitation of the MODPSO's search ability.
Considering the future work, in our opinion, two contents are worthy of study.One is that there exist other characteristics of the hyperspectral image that are not considered in this work, so the objective functions can be replaced by others to study the effect of different combinations of objective functions on the endmember extraction result.The other is that the multiobjective optimization method with better search ability can be studied to achieve a Pareto front with higher quality.

Figure 1 .
Figure 1.Feasible solutions for minimization optimization.Blue points stand for common solutions, and red points stand for Pareto-optimal solutions.

Figure 1 .
Figure 1.Feasible solutions for minimization optimization.Blue points stand for common solutions, and red points stand for Pareto-optimal solutions.

Figure 2 .
Figure 2. The update of the particle's personal best position.The blue point stands for the current pbest of one particle, and other points are possible locations of the particle at the next time.The plane can be divided into four parts centered on the pbest.If the particle appears in the area where the purple point located, then the pbest of the particle should remain unchanged; if the particle appears in the area where the red point located, then the pbest of the particle will be updated by the red point; and if the particle appears in the area where the cyan points located, then randomly select one point as the pbest.

Figure 2 .
Figure 2. The update of the particle's personal best position.The blue point stands for the current pbest of one particle, and other points are possible locations of the particle at the next time.The plane can be divided into four parts centered on the pbest.If the particle appears in the area where the purple point located, then the pbest of the particle should remain unchanged; if the particle appears in the area where the red point located, then the pbest of the particle will be updated by the red point; and if the particle appears in the area where the cyan points located, then randomly select one point as the pbest.

Figure 3 .
Figure3.Selection of the best local guide among the global best archive (GBA) for each particle.The squares stand for the GBA members, and circles stand for all the particles.The sigma values of all the GBA members and particles are calculated and compared.For one particle, the GBA member that has the closest sigma value with it is chosen as the best local guide for it.

Figure 3 .
Figure 3. Selection of the best local guide among the global best archive (GBA) for each particle.The squares stand for the GBA members, and circles stand for all the particles.The sigma values of all the GBA members and particles are calculated and compared.For one particle, the GBA member that has the closest sigma value with it is chosen as the best local guide for it.

Figure 5 .
Figure 5. Sub-scene extracted from the Washington DC dataset.

Figure 5 .
Figure 5. Sub-scene extracted from the Washington DC dataset.

Figure 6 .
Figure 6.The objective function value as a function of the number of iteration times for the Washington DC dataset: (a) the volume inverse; and (b) root-mean-square error (RMSE).

Figure 6 .
Figure 6.The objective function value as a function of the number of iteration times for the Washington DC dataset: (a) the volume inverse; and (b) root-mean-square error (RMSE).

Figure 6 .
Figure 6.The objective function value as a function of the number of iteration times for the Washington DC dataset: (a) the volume inverse; and (b) root-mean-square error (RMSE).

Figure 7 .
Figure 7.The results of the Washington DC image: (a) the Pareto front obtained by MODPSO; and (b) comparison of the results by four methods.

Figure 8 .
Figure 8. Endmember spectra manually selected from the image and automatically extracted by the four methods for the Washington DC dataset: (a) grass; (b) path; (c) roof; (d) street; (e) tree; and (f) water.

Figure 8 .
Figure 8. Endmember spectra manually selected from the image and automatically extracted by the four methods for the Washington DC dataset: (a) grass; (b) path; (c) roof; (d) street; (e) tree; and (f) water.

Figure 10
Figure10shows the objective function value as a function of the number of iteration times of MODPSO for the Urban dataset.It can be seen that the proposed method can converge to a stationary point when it reached the maximum iteration, and the volume value reached the stationary point earlier than the RMSE value.

Figure 10 .
Figure 10.The objective function value as a function of the number of iteration times for the Urban dataset: (a) the volume inverse; and (b) RMSE.

Figure 10
Figure10shows the objective function value as a function of the number of iteration times of MODPSO for the Urban dataset.It can be seen that the proposed method can converge to a stationary point when it reached the maximum iteration, and the volume value reached the stationary point earlier than the RMSE value.

Figure 10 Figure 10 .
Figure10shows the objective function value as a function of the number of iteration times of MODPSO for the Urban dataset.It can be seen that the proposed method can converge to a stationary point when it reached the maximum iteration, and the volume value reached the stationary point earlier than the RMSE value.

Figure 10 .
Figure 10.The objective function value as a function of the number of iteration times for the Urban dataset: (a) the volume inverse; and (b) RMSE.

Figure 11 .
Figure 11.The results of the Urban image: (a) the Pareto front obtained by MODPSO; and (b) comparison of the results by four methods.

Figure 11 .
Figure 11.The results of the Urban image: (a) the Pareto front obtained by MODPSO; and (b) comparison of the results by four methods.

Table 1 .
The results of objective function values and computation time for the Washington DC Dataset.

Table 1 .
The results of objective function values and computation time for the Washington DC Dataset.

Table 2 .
The results of objective function values and computation time for the Urban Dataset.

Table 2 .
The results of objective function values and computation time for the Urban Dataset.