1. Introduction
Dozens or even hundreds of narrow, adjacent spectral bands can effectively represent the unique materials in hyperspectral remote sensing images. However, due to the spatial resolution limitation of imaging sensors, the observation of one pixel may contain several different substances; this is called a “mixed pixel” [
1,
2]. Spectral unmixing is a technique that decomposes the mixed-pixel spectra into a collection of pure spectra, named endmembers, and the corresponding abundance fractions of each endmember [
3,
4,
5,
6,
7,
8]. Endmembers play an important role in quantitatively exploring the fine spectral information of hyperspectral images, and the appropriate extraction of endmembers from imagery is of great significance in hyperspectral image processing [
9,
10,
11,
12,
13,
14,
15].
Simplex theory is one of the most important theories on which many endmember extraction algorithms are based. Popular pure pixel based algorithms include the Pixel Purity Index (PPI) [
16], N-FINDR [
17], the Simplex Growing Algorithm (SGA) [
18], and Vertex Component Analysis (VCA) [
19], as well as some new swarm intelligence algorithms proposed in recent years, such as Discrete Particle Swarm Optimization (D-PSO), which is a relatively new technique that has been empirically shown to perform well in endmember extraction optimization problems as it searches for a global optimal combination of the endmember set in the feasible space [
20].
In essence, endmember extraction involves selecting a subset with size
k from a total set of
n variables to optimize a particular criterion under the hypothesis that pure pixels exist. Over the past few decades, population-based random optimization techniques, such as evolutionary algorithms and swarm intelligence optimization, have been widely employed to solve optimization problems [
21,
22]. The Particle Swarm Optimization (PSO) method, which was originally proposed by Kennedy [
23], is a widely-used swarm intelligence method for solving global optimization problems, while Zhang et al. have adapted discrete PSO (D-PSO) for endmember extraction. However, each particle in PSO moves in the search space with a certain velocity, which is dynamically adjusted according to its own optimum and the global optimum of all the particles. Therefore, PSO may easily get trapped in a local optimum. Although many attempts have been made to improve the performances of PSO [
24,
25], it does have a fatal flaw: There is no guarantee that it can achieve globally optimal results, which was demonstrated by Van den Bergh [
26]. In addition, PSO contains at least three parameters, which have a great influence on the algorithm efficiency.
This situation remained unchanged until Sun et al. [
27] introduced quantum theory into PSO and proposed a Quantum-Behaved PSO (QPSO) algorithm, which can be theoretically guaranteed to find the globally optimal solution in the search space. In this paper, an enhanced version of PSO, the Mutation Operator Accelerated Quantum-Behaved Particle Swarm Optimization (MOAQPSO) algorithm, is proposed, which is expected to guarantee convergence on globally optimal results, as well as having a fast convergence speed. Among all the quantum physics based theories, the QPSO algorithm has proved to be a promising intelligent optimization algorithm [
27]. All the particles in QPSO have a quantum behavior instead of the classical Newton’s laws of motion that are followed in classical PSO. This quantum motion ensures that the particles can appear in any position with a certain probability. Therefore, QPSO can theoretically guarantee the efficiency and global convergence of the algorithm. Another attractive feature of the QPSO algorithm is that only a few parameters need to be adjusted. Above all, QPSO is easy to implement, has very few parameters, and is more robust than PSO. In this paper, the MOAQPSO algorithm is proposed for hyperspectral endmember extraction. Furthermore, a high-dimensional particle definition is proposed, which is very different from the existing methods of particle definition.
In this paper we introduced QPSO to endmember extraction, which to the best of our knowledge has not previously been used for this purpose. In order to ensure the effectiveness of QPSO for endmember extraction, some strategies were used:
(1) The mutation operation. Inspired by the crossover operator and mutation operator in genetic algorithm, we used mutation operation in QPSO to accelerate convergence. The maintenance of population diversity is undoubtedly a key issue in the performance of an evolutionary algorithm, by which the proposed algorithm can enhance its ability to jump out of local optima. The mutation operation, which gives up the current position at a certain probability, is used to effectively maintain the population diversity and accelerate convergence.
(2) The novel particle definition for endmember extraction, which differs from the existing methods of particle definition. Particle definition is another important problem for QPSO application. In this paper, a high-dimensional particle definition is proposed, which ensures that each particle has a physical significance for the endmember extraction.
The remainder of the paper is organized as follows.
Section 2 describes some basic details of PSO and QPSO. The proposed MOAQPSO algorithm is presented in
Section 3.
Section 4 presents an experimental comparison of the proposed method with both simulated and real data sets.
Section 5 concludes the paper and introduces the future research.
3. The MOAQPSO Algorithm for Hyperspectral Endmember Extraction
3.1. Endmember Extraction Problem
The linear unmixing of hyperspectral images is a popular approach used to determine and quantify materials in remotely sensed images. The Linear Mixing Model (LMM) assumes that the spectrum of each pixel within the scene is a linear combination of the material constituents within that pixel [
1]:
where
represents one pixel of a hyperspectral image (
L is the number of spectral bands);
is the endmember set (
D is the number of endmembers, which is known already);
is the abundance of the
jth endmember in the pixel
; and
is a noise vector. Equation (21), the Abundance Sum-to-One Constraint (ASC), and the Abundance Nonnegative Constraint (ANC) represent the physical meaning of the abundance.
The Root-Mean-Square Error (RMSE) between the original image and its remixed image is defined as
where
and
mean the estimated value by unmixing methods.
If the number of endmembers is known, the smaller the RMSE, the better the result of the endmember extraction.
According to the LMM, all the pixels in a hyperspectral image form a simplex, and the vertices of the convex simplex will be the endmembers of the image. The volume of the simplex is defined as
If the number of endmembers is known, a greater simplex volume suggests a better result of the endmember extraction [
17,
18,
19].
Therefore, the problem of endmember extraction can be regarded as a combinatorial optimization problem by minimizing the RMSE or maximizing the volume. In order to calculate Equation (25), unconstrained least squares or constrained least squares are used according to whether the constraint conditions are considered or not. The iterative technique used for solving the nonnegative constraint problem needs more time, and the unconstrained approach or the equality constraint cannot achieve comparative results. Therefore, the volume of the simplex is chosen as the fitness function in our research [
17,
18,
19,
32]. The fitness function can be replaced by any appropriate objective function. In this paper, we discuss the endmember extraction method only, under the hypothesis that pure pixels exist. The optimization model for a maximum-volume method is
3.2. The MOAQPSO Algorithm for Endmember Extraction
MOAQPSO also uses particles to search the feasible solution space of the optimization problem. We assume that the whole image pixel set is
, where
,
n is the number of pixels, and
L is the number of bands. The feasible solution space of the endmember extraction at time
t can be defined as
where
,
represents the
jth endmember vector in the
ith particle,
M is the number of particles, and
D is the number of endmembers.
Although there are no velocity vectors in MOAQPSO, is the best position of particle i, and is the best position of the entire swarm. The individual best position of each particle and the globally best position are updated in the same way as in PSO.
The position of the particle updates according to the following equation:
where:
Some research has proved that a linearly decreasing parameter
outperforms a constant, both in global search ability and local search ability [
30]. Therefore, parameter
is linearly decreasing in this paper.
According to Definition (26), . In order to meet the above conditions, the position of the particle should be adjusted after updating according to Equation (27). One strategy is finding the most similar pixel in set R with , calculated by Equation (27). The Euclidean distance and Spectral Angle Distance (SAD) are the commonly used criteria for similarity measure. As the Euclidean distance calculation has a lower computation complexity than that of SAD does, the Euclidean distance is used in this paper. It should be pointed out that any other strategy with the same function could replace the Euclidean distance.
Finally, the position of the particle updates according to the following equations as follows:
where Equation (31) represents the distance between the
jth endmember in the
ith particle and the
kth pixel.
With the increase in the number of iterations, the diversity of the population is reduced, and the algorithm can easily fall into a local optimum. Therefore, the algorithm would achieve a poor solution at a certain probability (pp), and this probability (pp) is called the mutation rate.
3.3. The Flow Chart of the MOAQPSO Algorithm
(1) Input Image
Due to the feasible solution space being defined as shown in Definition (26), the number of pixels in the input image has little impact on the search speed, which is different from the D-PSO. Therefore, a preprocessing for candidate endmember extraction is unnecessary.
(2) Initialization of the Population
At first, each individual in the population randomly selects D pixels from the image, as defined in Equation (25). Each vector in each individual represents one pixel spectrum, which has a high dimension. In general, one particle is one possible endmember set.
(3) Calculating the Individual Best Position of Each Particle and the Global Optimum Position
Each particle, which is an endmember set, is applied with the simplex volume calculation using Equation (24).
The individual best position of each particle is updated as follows:
The globally best position is defined as
(4) Particle Evolution
Not all the particles evolve following the evolution Equation (27). Each particle is mutated with the mutation rate pp to obtain its associated mutated position by randomly selecting D pixels from the image. Otherwise, the position of the particle updates according to Equation (27), and the updated position should be mapped to the image according to Equation (30).
(5) Stopping Condition
If the result converges to a fixed value, or the iteration number reaches the maximum, the best endmember extraction solution is the output.
The flowchart of MOAQPSO is shown in
Figure 1.
4. Experiments and Analysis
In this section, the extensive experiments conducted on one simulated and two real hyperspectral data are described. The algorithms used for comparison include N-FINDR, VCA, D-PSO and MOAQPSO. These typical algorithms are all based on the LMM and adopt the pure pixel assumption.
The SAD metric, which is widely used in the study of hyperspectral unmixing, was used to evaluate the endmember extraction results. The well-known SAD metric between pixel
and
is given by
where
and
are the norms of
and
, respectively.
The RMSE was the other metric used to evaluate the results of the experiments. In the experiments on the synthetic data set, the abundance maps were available, therefore errors in the abundance maps could be calculated, and the RMSEs between the reference abundances and the calculated abundances were determined. In the experiments on the real data set, the RMSE between the original image and the reconstruction image was used because the abundance maps were not available.
N-FINDR, VCA, and D-PSO were selected for comparative purposes. In D-PSO, the random selection probability
pp was set as 0.2, the number of particles was set to 10, and the VCA algorithm was used to obtain the 50 candidate endmembers, which could constitute a feasible solution space. In MOAQPSO, parameter
was linearly decreased from 1 to 0.5, according to previous studies [
33]. The mutation rate was set to
pp = 0.4 in the experiments, and the number of particles was set to 20 (synthetic data) and 30 (real data). The sensitivity of the parameter mutation rate
pp was analyzed using the synthetic data.
4.1. Synthetic Hyperspectral Image
The simulated data cube was generated by six randomly selected spectral signatures (as shown in
Figure 2). The simulated data contained 256 × 256 pixels covering 216 bands. The data was generated by the LMM. The fractional abundances satisfied the nonnegative constraint and the sum-to-one constraint. At last, the fractional abundances were piecewise smooth (i.e., they were smooth with sharp transitions), as shown in
Figure 3. The simulated data were generated using the hyperspectral imagery synthesis toolbox for Matrix Laboratory (MATLAB) [
34]. After the simulated data was generated, the scene was contaminated with Gaussian noise (SNR = 40 dB).
The quantitative evaluation metric for the endmember extraction results was the SAD value between each extracted endmember and the set of ground-truth spectral signatures. The lowest spectral angle represents the best performance.
Table 1 shows the SAD results obtained with N-FINDR, VCA, D-PSO, and MOAQPSO, respectively. The numbers in bold correspond to the best performances. From
Table 1, we can see that the proposed MOAQPSO achieves a better result (a lower mean SAD) than VCA, N-FINDR, and D-PSO. D-PSO does not obtain a good result. The reason for this may be that the synthetic data used in the experiments only contain a small number of pure pixels, which have affected the results of the D-PSO algorithm.
Table 2 shows the volume of the endmember set obtained with N-FINDR, VCA, D-PSO, and MOAQPSO, respectively. The numbers in bold correspond to the best performances. From
Table 2, we can see that the proposed MOAQPSO achieved greater simplex volume than any other method. This is because the volume was the objective function of the proposed method.
Table 3 shows a comparison of the mean RMSE between the reference and the abundance maps calculated by the endmember sets obtained by N-FINDR, VCA, D-PSO, and MOAQPSO, with the best results marked in bold. Although the mean RMSE of MOAQPSO is slightly higher than that of N-FINDR, it is still lower than that of VCA and D-PSO. This is because the objective function of MOAQPSO in the experiments was maximum-volume, and it achieved the greatest volume as shown in
Table 2.
In this part, the sensitivity of the parameter mutation rate
pp is analyzed. The population sizes were set to 20 and 30 respectively. The values of
pp were set to 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, and 0.7, respectively. The experimental results are shown in
Figure 4, where the horizontal axis represents the number of iterations, the vertical axis represents the volume, and the colors represent the
pp value. From
Figure 4, we can see that significant changes in the fitness function values occur over the first 150 iterations. After 200 iterations, the fitness function may increase slowly. A large population size reaches the optimal value with fewer iterations. This is because with a greater population size, the search scope increases, and a better solution is more likely to be found with fewer iterations. From
Figure 4, we can also see that the effect of the mutation operation is obvious: MOAQPSO can obtain much better results with the mutation operation than without the mutation operation. However, the regularity is not so obvious: (1) all the mutation operations from 0.1 to 0.7 are clearly effective; and (2) a low mutation rate results in a fast convergence speed, and a high mutation rate results in a high precision.
4.2. AVIRIS Real Data Set
In this section, we perform MOAQPSO for endmember extraction from a real hyperspectral image. This real data set covering the Cuprite scene was captured by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS) sensor in 1997, and a sub-scene is used in the experiments, which is shown in
Figure 5. The size of the sub-scene is 350 × 400, covering the wavelength range of 2.0–2.48 μm. The data set has become a popular benchmark data set for algorithm evaluation, due to the extensive ground-truth data available for the scene from the United States Geological Survey (USGS). Swayze and Clark [
35] have also produced a report about the ground truth of the area (shown in
Figure 6).
It should be pointed out that: (1) The number of endmembers in this image is difficult to determine and the number of endmembers to be extracted was set to 10 for all the algorithms, according to
Figure 6, and the reference spectra chosen from the USGS library are shown in
Figure 7; (2) Due to the large number of pixels, which influences the speed of the mapping into image step in the MOAQPSO algorithm, another operation for acceleration is the selective mapping. In other words, mapping into image step does not work in each iteration. For example, if the length of the interval is set to 10, mapping into the image step will work once in every 10 iterations. In the real data experiments, the length of the interval was set to 10, and the number of iterations was set to 2000.
According to
Figure 5 and
Figure 6, the major minerals in the experiment region include Alunite, Buddingtonite, Calcite, Chlorite, Jarosite, Halloysite, Kaolinite, Montmorillonite, Muscovite, and Nontronite. The SADs were calculated between the reference spectra from the USGS library and the endmembers were extracted by those algorithms. The reference spectra chosen from the USGS library are shown in
Figure 7. The performances of the algorithms were compared in terms of: (1) the number, type and average spectral angle of each mineral; (2) the volumes and RMSE of the endmember set. The results are shown in
Table 4,
Table 5 and
Table 6.
Table 4 shows the SAD results obtained by N-FINDR, VCA, D-PSO, and MOAQPSO, respectively, with the best performance for each endmember marked in bold. From
Table 4 we can observe that the endmembers extracted by the various algorithms were not the same. The proposed MOAQPSO and N-FINDR extracted seven different minerals, but failed to correctly extract Buddingtonite. VCA extracted six different minerals, but failed to correctly extract Halloysite and Montmorillonite. D-PSO extracted four different minerals, but failed to correctly extract Buddingtonite, Kaolinite and Montmorillonite. Although D-PSO achieved the minimum mean SAD, more classes’ spectra were ignored or repeated. Above all, the proposed MOAQPSO extracts a comprehensive endmember set with a higher precision than that of N-FINDR and VCA.
Table 5 shows the volume of the endmember set obtained by N-FINDR, VCA, D-PSO, and MOAQPSO, respectively. The numbers in bold correspond to the best performances. Because the volume is computed by the objective function of the proposed method, the MOAQPSO algorithm achieves the maximum simplex volume, followed by N-FINDR, VCA, and D-PSO.
Table 6 shows a comparison of the mean RMSE between the original image and the reconstruction image calculated by the endmember set obtained by N-FINDR, VCA, D-PSO, and MOAQPSO. The numbers in bold correspond to the best performances. Although the mean RMSE of MOAQPSO is slightly higher than that of VCA, it is lower than that of N-FINDR and D-PSO. This is because the objective function of MOAQPSO in the experiments is maximum-volume. If the objective function of MOAQPSO is the RMSE (named as MOAQPSO-2), the mean RMSE of MOAQPSO-2 is 2.1396. Compared with
Table 6, we can see that the mean RMSE of MOAQPSO-2 is the lowest compared with that of N-FINDR, VCA, D-PSO, and MOAQPSO.
4.3. HYDICE Data Set
In this section, we perform MOAQPSO for endmember extraction from another real hyperspectral image. The hyperspectral data set, covering an area in a Washington DC Mall, was collected by the Hyperspectral Digital Imagery Collection Experiment (HYDICE) sensor. This image had 210 spectral channels ranging between 0.4 and 2.4 μm. Bands where the atmosphere was opaque or with a high noise level were omitted from the data set, leaving 180 bands. In this experiment, a subset of 80 × 130 was extracted from the original image.
Figure 8 shows the false-color image composed of bands 60, 27, and 17 for the red, green, and blue channels, respectively.
The sub-scene includes seven types of ground objects: water, tree, grass, roof 1, roof 2, street, and path. Seven kinds of pure pixels are selected manually, corresponding to the seven ground objects. In
Figure 8, the areas in the yellow circles are the selected pure pixels. The averages of the selected pixels in each class were computed as the reference endmember signatures [
36]. The SADs were calculated between the reference spectra manually selected from the image and the endmembers were extracted by the algorithms. The performances of the algorithms were compared in terms of: (1) the average spectral angle of each mineral; (2) the volumes and RMSE of the endmember set. The results are shown in
Table 7,
Table 8 and
Table 9.
As the number of pixels will influence the speed of mapping into image step in the MOAQPSO algorithm, another operation for acceleration is selectively mapping. The length of the interval was also set to 10 in this real data experiment as the first real image.
Table 7 shows the SAD results obtained with N-FINDR, VCA, D-PSO, and MOAQPSO, respectively, with the best performance marked in bold. From
Table 7, we can see that all the methods extracted seven different materials and the proposed MOAQPSO achieves the minimum mean SAD value.
Table 8 shows the volume of the endmember sets obtained with N-FINDR, VCA, D-PSO, and MOAQPSO, respectively. The numbers in bold correspond to the best performances. From
Table 8, we can see that the MOAQPSO algorithm did not achieve a greater simplex volume than other methods, although the volume was computed as the objective function. The reason for this may be that there are many pure pixels in this data set, and it is very hard to get the optimal solution and the number of iterations is not enough in this case.
Table 9 shows a comparison of the mean RMSE between the original image and the reconstructed image calculated by the endmember sets obtained with N-FINDR, VCA, D-PSO, and MOAQPSO. The numbers in bold correspond to the best performances. Although the mean RMSE of MOAQPSO is slightly higher than that of VCA and D-PSO, it is lower than that of N-FINDR. If the objective function of MOAQPSO is the RMSE (named as MOAQPSO-2), the mean RMSE of MOAQPSO-2 is 31.392. Compared with
Table 6, we can see that the mean RMSE of MOAQPSO-2 is the lowest compared with that of N-FINDR, VCA, D-PSO, and MOAQPSO.
5. Discussion
In this paper, a novel Mutation Operation Accelerated Quantum-Behaved Particle Swarm Optimization (MOAQPSO) algorithm is proposed for hyperspectral endmember extraction. In order to ensure the effectiveness of QPSO for endmember extraction, some strategies are used, including: (1) The mutation operation. The mutation operation, which gives up the current position at a certain probability, is used to effectively maintain the population diversity and accelerate convergence. The effect of the mutation operation is obvious from the experiments results; (2) In this paper, a high-dimensional particle definition is proposed, and it ensures that each particle contains a physical significance for the endmember extraction.
The endmember spectra extracted by the proposed MOAQPSO method fit well with the references. From above experiments, we can see that: (1) the proposed MOAQPSO method can extract a more complete endmembers set with a lower mean SAD than that of N-FINDR, VCA, and D-PSO; (2) the MOAQPSO algorithm with volume as the objective function can achieve the maximum simplex volume, followed by N-FINDR, VCA, and D-PSO, and the MOAQPSO-2 algorithm with RMSE as the objective function achieves the minimum mean RMSE. However, there are still some unresolved issues and research on the use of QPSO for endmember extraction is still in its infancy. First, the proposed approach assumes that the image contains pure pixels. It cannot obtain the desired extraction results for highly mixed data and data with considerable noise. Second, the volume of the simplex is chosen as the fitness function in this paper. The fitness function can be replaced by any appropriate objective function. In our future work, we will try different kinds of fitness functions for endmember searching and endmember generation.
6. Conclusions
In this paper, a novel Mutation Operator Accelerated Quantum-Behaved Particle Swarm Optimization (namely MOAQPSO) approach has been proposed. All the particles in MOAQPSO have a quantum behavior instead of the classical Newton’s laws of motion that are followed in classical PSO. This quantum motion ensures that the particles can appear in any positions at a certain probability. The proposed algorithm introduces QPSO to endmember extraction, and a mutation operation is used to increase the population diversity, which has a significant effect, as proved by the experimental results. What is more, only two parameters need to be adjusted: and pp. The default value is linearly decreased from 1 to 0.5 and 0.4, respectively, and the recommended range for these parameters is [0.5–1.5] and [0.3–0.7], respectively. The proposed MOAQPSO algorithm was evaluated by both synthetic and real hyperspectral data sets. The experimental results indicated that the proposed method obtained better results than the algorithms of VCA, N-FINDR, and D-PSO, with a lower mean SAD accuracy.