Multipopulation Particle Swarm Optimization for Evolutionary Multitasking Sparse Unmixing

Abstract

Recently, the multiobjective evolutionary algorithms (MOEAs) have been designed to cope with the sparse unmixing problem. Due to the excellent performance of MOEAs in solving the NP hard optimization problems, they have also achieved good results for the sparse unmixing problems. However, most of these MOEA-based methods only deal with a single pixel for unmixing and are subjected to low efficiency and are time-consuming. In fact, sparse unmixing can naturally be seen as a multitasking problem when the hyperspectral imagery is clustered into several homogeneous regions, so that evolutionary multitasking can be employed to take advantage of the implicit parallelism from different regions. In this paper, a novel evolutionary multitasking multipopulation particle swarm optimization framework is proposed to solve the hyperspectral sparse unmixing problem. First, we resort to evolutionary multitasking optimization to cluster the hyperspectral image into multiple homogeneous regions, and directly process the entire spectral matrix in multiple regions to avoid dimensional disasters. In addition, we design a novel multipopulation particle swarm optimization method for major evolutionary exploration. Furthermore, an intra-task and inter-task transfer and a local exploration strategy are designed for balancing the exchange of useful information in the multitasking evolutionary process. Experimental results on two benchmark hyperspectral datasets demonstrate the effectiveness of the proposed method compared with the state-of-the-art sparse unmixing algorithms.

1. Introduction

With the progress of remote sensing technology, hyperspectral imagery, which can obtain hundreds of sequential spectrum bands, has been widely applied in both civilian and military scenarios, for example, land-cover classification [1,2,3], environmental monitoring [4,5,6] and target detection [7,8], and so forth. However, there remains the problem of mixed pixels due to the low spatial resolution of sensors and the mixture of the surface features [9,10]. Therefore, spectral unmixing aims at extracting the collection of constituent spectra (called endmembers) from the mixed pixels and calculating the fractional abundances of these endmembers [11,12]. Accordingly, different spectral unmixing methods can be divided into three categories, that is, the geometrical-based, statistical-based and sparse-regression-based approaches. Traditional geometrical-based and statistical-based methods are extensively used as they can be utilized easily and flexibly, but they also suffer from the weakness of poor performance on highly mixed scenes spectra and the limitedness of time consumption, respectively [13]. Sparse unmixing, as an emerging spectral unmixing technology in recent years, is devised to find out the optimal solution that can represent each pixel of the hyperspectral image the most from a spectral library known in advance. Among these algorithms, the sparse unmixing via variable splitting and augmented Lagrangian (SUnSAL) based on the alternating direction method of multipliers has been proposed to relax the $l_0$ norm [14]. To overcome the disadvantage of SUnSAL that only utilizes spectral information without considering the spatial-contextual information, Iordache et al. proposed the collaborative SUnSAL (CLSUnSAL) which improves the unmixing results by solving a joint sparse regression problem, where the sparsity is simultaneously imposed to all pixels in the dataset [15,16].

Mathematically, sparse unmixing is an NP-hard problem. Multiobjective evolutionary algorithms (MOEAs), which are able to optimize some contradictory objectives and acquire a set of nondominated solutions called the Pareto-optimal front, are suitable for solving the NP-hard problems and overcoming the aforementioned difficulty in sparse unmixing [17]. A multiobjective sparse unmixing (MOSU) model was first proposed by Gong et al. [18] to deal with the sparse unmixing for hyperspectral imagery. Xu et al. [19] developed a multiobjective optimization based sparse unmixing (SMoSU) to take full advantage of the spectral characteristics of hyperspectral images under the framework of the multiobjective evolutionary algorithm based on decomposition (MOEA/D). In [20], the SMoSU was further improved and a classification-based model called CM-MoSU was designed. The estimation of distribution algorithms is modified to pay more attention to the feasible space with high quality.

However, the existing sparse unmixing algorithms based on MOEAs are limited to the pixel-based unmixing, which leads to the disadvantage of the low efficiency and the lack of the spatial structure information [21]. In some recent studies [22,23,24], a hyperspectral image is clustered into multiple homogeneous regions based on the assumption that the probability of the active endmember set in the homogeneous region is likely to be the same, which not only reduces the complexity of unmixing, but also further enhances the spatial correlation of pixels in the same category. Interestingly, this coincides with the idea of evolutionary multitasking framework emerging in recent years. The evolutionary multitasking [25] aims to solve different multiobjective optimization problems simultaneously to take advantage of the implicit parallelism from different tasks. Therefore, it is promising to employ the evolutionary multitasking multiobjective framework to efficiently solve the sparse unmixing problem. Besides, the particle swarm optimization (PSO) algorithm, which simulates the regularity of bird cluster activities, has proved to be effective in solving multiobjective endmember extraction problems [26,27,28]. From this, the current multitasking paradigm can be further explored and applied to sparse unmixing problems.

In this paper, we propose a novel evolutionary multitasking multipopulation particle swarm optimization (EMMPSO) framework for sparse unmixing. In the proposed method, a hyperspectral image is clustered into multiple homogeneous regions first, then the multipopulation particle swarm optimization is employed to explore each sparsity. Finally, the multiobjective optimization is applied to each task simultaneously to obtain a compromise between the reconstruction error and the endmember sparsity. Significantly, it is different from the traditional MOEA-based algorithms that EMMPSO can process the entire matrix due to the decomposition strategy of evolutionary multitasking, aiming at pixel-based unmixing only. In addition, we design a novel intra-task and inter-task transfer strategy to overcome the impact of negative transfer in multitasking. It can not only utilize the effective information in the same task to speed up the convergence of each sub-particle swarm, but also explore the similarities between different tasks to improve the overall convergence performance. Finally, the Pareto optimal solution in each task can be obtained to reverse the final endmember abundance.

The contributions of the proposed EMMPSO algorithm are summarized as follows:

(1)

A novel evolutionary multitasking multipopulation particle swarm optimization framework is proposed to solve the sparse unmixing problem. With the decomposition of the evolutionary multitasking, multiple homogeneous regions of a hyperspectral image can be processed simultaneously, which can accelerate the convergence by exploring the relevance of all the tasks. In addition, the Pareto optimal solution between the reconstruction error and the endmember sparsity can be obtained with the multiobjective optimization.

(2)

A multipopulation particle swarm optimization is designed in the multitasking framework for the major evolution. In addition, the intra-task and inter-task transfer strategy are proposed to balance the evolutionary process of exploration and utilization. An efficient local exploration strategy with MOEA is designed to facilitate the search process to obtain the optimal points.

(3)

The superiority of EMMPSO on the convergence speed, global optimization performance and unmixing accuracy is substantiated compared with the classical mathemati cal-based and MOEA-basd sparse unmixing algorithms.

The remainder of this paper is structured as follows: Section 2 briefly reviews some related work on sparse unmixing. In Section 3, our method is introduced in detail. Section 4 gives the experimental settings and the analysis of the experimental results. Finally, the conclusions and future works are described in Section 5.

2. Related Work

Generally, the mixed pixels are usually unmixed in the linear mixing model. For a single mixed pixel $\mathit{y}\in {\mathbb{R}}^{\mathit{L}\times 1}$ with L spectral bands, which can be expressed as:

$$\begin{array}{c}\hfill \mathit{y}=\mathit{Ax}+\mathit{n},\end{array}$$

(1)

where $\mathit{A}\in {\mathbb{R}}^{\mathit{L}\times \mathit{D}}$ is the spectral library. It is worth noting that all the spectral information is known in advance in the spectral library. In addition, $\mathit{x}\in {\mathbb{R}}^{\mathit{D}\times 1}$ is the corresponding fractional abundance vector, that is, the proportion of each endmember, and $\mathit{n}\in {\mathbb{R}}^{\mathit{L}\times 1}$ represents the noise term for the mixed pixel. In normal circumstances, a hyperspectral image $\mathit{Y}\in {\mathbb{R}}^{\mathit{L}\times \mathit{n}}$ contains $\mathit{n}$ pixels, the matrix form of (1) can be formulated as:

$$\begin{array}{c}\hfill \mathit{Y}=\mathit{AX}+\mathit{N}.\end{array}$$

(2)

Therefore, the purpose of sparse unmixing is to obtain the most suitable set of endmembers for the reconstructing remote sensing image from the huge spectral library. Mathematically, this is an NP hard optimization problem, which can be expressed as:

$$\begin{array}{c}\underset{\mathit{x}}{min}{∥\mathit{x}∥}_0,\mathit{s}.\mathit{t}.{∥\mathit{y}-\mathit{Ax}∥}_2^2\le \delta .\hfill \end{array}$$

(3)

Many studies employed the relaxation methods to solve the ${\mathit{l}}_0$-norm problem. SUnSAL [14] resorted to the ${\mathit{l}}_0$-norm to match ${\mathit{l}}_0$-norm, and the mathematical optimization formula is as follows:

$$\begin{array}{c}\underset{\mathit{x}}{min}(1/2){∥\mathit{y}-\mathit{Ax}∥}_2^2+\lambda {∥\mathit{x}∥}_1+{\iota }_{\mathit{R}+}\left(\mathit{x}\right)+{\iota }_{\left\{1\right\}}\left({\mathbf{1}}^{\mathbf{T}}\mathit{x}\right),\hfill \end{array}$$

(4)

where $\lambda$ stands for a regularization parameter that controls the relative weight between the sparse term and the error term. In [15], the CLSUnSAL takes spatial information into account and directly processes the whole matrix, which is shown as follows:

$$\begin{array}{c}\underset{\mathit{X}}{min}{∥\mathit{Y}-\mathit{AX}∥}_F^2+\lambda {∥\mathit{X}∥}_{2,1}+{\iota }_{\mathit{R}+}\left(\mathit{X}\right).\hfill \end{array}$$

(5)

Considering the excellent performance of MOEAs in solving NP-hard optimization problems, many studies have turned their attention to MOEAs to solve the sparse unmixing problem in recent year. Gong [18] proposed a novel multiobjective cooperative coevolutionary algorithm to optimize the reconstruction term, the sparsity term and the total variation regularization term simultaneously, which can be expressed as:

$$\begin{array}{c}\underset{\mathit{x}}{min}({∥\mathit{y}-\mathit{Ax}∥}_2^2,{∥\mathit{x}∥}_0,\sum _{\mathit{j}\in \epsilon }{∥\mathit{x}-{\mathit{x}}_j∥}_1),\hfill \end{array}$$

(6)

where $\epsilon$ stands for the set of the horizontal and vertical neighbors in $\mathit{X}$. Jiang [29] decomposed the sparse unmixing problem into two stages and employed the MOEAs to solve them separately. In the first phase, it is mainly aimed at the endmember extraction, the optimized formula is as follows: $\underset{\mathit{M}}{min}(RSE1,SP1)$, where RSE1 is the residual of the measured hyperspectral image, SP1 represents the size of the measured estimated endmembers (M). In the second phase, the extracted abundance estimation becomes the focus, which can be expressed as: $\underset{\mathit{M}}{min}(RSE2,SP2)$, where RSE2 is the residuals of the hyperspectral unmixing, SP2 represents the favorable abundance matrix obtained by incorporating the spatial–contextual information. In addition, Jiang [30] improved the Tp-MoSU to settle the problems of the limited performance in identifying real endmembers from high-noise data in the first phase, and cannot effectively use the spatial context information in the second phase due to the similarity metric used. Besides, many sparse unmixing algorithms based on evolutionary multiobjective decomposition [19,20,31] have also been explored.

Recently, evolutionary multitasking optimization [25,32] has become a new favorite in the field of evolutionary computing. In a nutshell, evolutionary multitasking aims to deal with multiple optimization problems at the same time, and promote the optimization of each task by exploring the hidden relationship between these optimization problems. It is worth noting that many evolutionary multitasking optimization related algorithms have been explored and applied to many fields, such as feature selection [33], reinforcement learning [34] and sparse regression [22] and so forth. In sparse unmixing, a hyperspectral image can be clustered into multiple homogeneous regions according to spatial information, so this coincides with the concept of evolutionary multitasking. It is very promising to model each homogeneous region as an optimization task, though the decomposition of multiple tasks can effectively reduce the impact of dimensional disasters.

3. EMMPSO Framework

The pseudo code of EMMPSO is shown in Algorithm 1. In this section, the proposed framework is introduced in detail from initialization, multipopulation particle swarm optimization and the decision making with MOEA.

Algorithm 1 The EMMPSO Framework

1: %Initialization   2: Set t = 0, G = ⌀.   3: forj = 1 to K do   4:    for s = 1 to S do   5:      for i = 1 to N/KS do   6:         ${\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}$ = $\left\{({\mathit{x}}_1,{\mathit{x}}_2,...,{\mathit{x}}_{\mathit{n}})\right|{\mathit{x}}_{\mathit{i}}\in \{0,1\},\left|\right|{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}|{|}_0=\mathit{s}\}$.   7:      end for   8:    end for   9:    Evaluate the fitness of each particle in task ${\mathit{T}}_{\mathit{j}}$. 10:    Assign the skill factor ${\tau }_{\mathit{i}}$. 11:    Initialize the ${\left\{{\mathit{pbest}}_{\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}$ and ${\left\{{\mathit{gbest}}_{\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}$. 12:    ${\left\{{\mathit{V}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}=({\left\{{\mathit{pbest}}_{\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}-{\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}})+({\left\{{\mathit{gbest}}_{\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}-{\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}})$. 13:    ${\mathit{G}}_{\mathit{j}}^{\mathit{t}}={\sum }_{\mathit{s}=1}^{\mathit{S}}{\left\{{\mathit{gbest}}_{\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}$. 14: end for 15: %Evolution 16: whilet < Maxt do 17:    for j = 1 to K do 18:      Update the ${\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}+1}\right\}}_{\mathit{j}}$ and ${\left\{{\mathit{V}}_{\mathit{i},\mathit{s}}^{\mathit{t}+1}\right\}}_{\mathit{j}}$ based on (9). 19:    end for 20:    Update the particle according to Algorithm 2. 21:    Evaluate the fitness of each particle. 22:    Update the ${\left\{{\mathit{pbest}}_{\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}$, ${\left\{{\mathit{gbest}}_{\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}$, and ${\mathit{G}}_{\mathit{j}}^{\mathit{t}}$ with the Local Exploration Strategy. 23:    t = t + 1. 24: end while 25: %Decision Making 26: Obtain the optimal point in each task from ${\mathit{G}}_{\mathit{j}}$.

3.1. Initialization and Representation

In sparse unmixing, the spectral library known in advance and a hyperspectral image are input for processing, and the endmember set selected from the library and the corresponding abundance map are output. In the proposed EMMPSO, a hyperspectral image is first clustered into $\mathit{K}$ homogeneous regions, and each homogeneous region is processed as a task [22], which is shown in Figure 1. The spectra of the entire spectrum library are coded into each particle in order, that is, the length of particle is equal to the number of spectra. Considering that the sparsity of particles remains unchanged in the evolution for most current discrete particle swarm optimization algorithms, the population in each task is divided into multiple subpopulations according to the sparsity to ensure that there are particles to explore in each sparsity. For the $\mathit{s}$-th subpopulation in the $\mathit{j}$-th task, the position of each particle is initialized as follows:

$${\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}=\{({\mathit{x}}_1,{\mathit{x}}_2,...,{\mathit{x}}_{\mathit{n}})|{\mathit{x}}_{\mathit{i}}\in \{0,1\},\left|\right|{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}{\left|\right|}_0=\mathit{s}\},$$

(7)

where ${\mathit{x}}_{\mathit{i}}$ is composed of two elements, 0 or 1. If the ${\mathit{x}}_{\mathit{i}}$ is equal to 1, it means that the spectrum at the corresponding position in the spectral library is selected, and vice versa.

Then, each particle is evaluated with the reconstruction error ($\left|\right|\mathit{Y}-{\mathit{A}}_{\mathit{v}}{\mathit{X}}_{\mathit{v}}{\left|\right|}_{\mathit{F}}$) in the corresponding task, where the $\mathit{Y}$ is the hyperspectral image, $\mathit{v}$ represents the endmember set from the particle ${\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}$, ${\mathit{A}}_{\mathit{v}}$ and ${\mathit{X}}_{\mathit{v}}$ are the subset of spectral library $\mathit{A}$ and the abundances of endmembers, respectively. After the evaluation is completed, the skill factor ${\tau }_{\mathit{i},\mathit{s}}$, defined as the task with the best performance of the subpopulations with sparsity $\mathit{s}$ in all the tasks, is assigned to each particle. Besides, the ${\left\{{\mathit{pbest}}_{\mathit{s}}\right\}}_{\mathit{j}}$ and ${\left\{{\mathit{gbest}}_{\mathit{s}}\right\}}_{\mathit{j}}$ for the subpopulation with sparsity $\mathit{s}$ in the $\mathit{j}$-th task can be obtained. The velocity of particle is initialized as:

$$\begin{array}{c}{\left\{{\mathit{V}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}=({\left\{{\mathit{pbest}}_{\mathit{s}}\right\}}_{\mathit{j}}-{\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}})+({\left\{{\mathit{gbest}}_{\mathit{s}}\right\}}_{\mathit{j}}-{\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}).\hfill \end{array}$$

(8)

3.2. Multipopulation Particle Swarm Optimization for Knowledge Transfer

Considering the discreteness of decision variables in sparse unmixing, the population in each task is divided into multiple subpopulations according to the sparsity during initialization. In the process of particle swarm evolution, the position and velocity of the particles in the j-th task with the sparsity s are updated as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}+1}\right\}}_{\mathit{j}}={\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}+{\left\{{\mathit{V}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}},\hfill \\ & {\left\{{\mathit{V}}_{\mathit{i},\mathit{s}}^{\mathit{t}+1}\right\}}_{\mathit{j}}=\left\{\begin{array}{c}\mathit{T}\left({\left\{{\mathit{V}}_{\mathit{i},\mathit{s}}^{\mathit{t}+1}\right\}}_{\mathit{j}}\right),\mathrm{if}(\mathrm{any}\left({\left\{{\mathit{V}}_{\mathit{i},\mathit{s}}^{\mathit{t}}\right\}}_{\mathit{j}}\right)\ge 0)\hfill \\ \mathit{R}\left({\left\{{\mathit{X}}_{\mathit{i},\mathit{s}}^{\mathit{t}+1}\right\}}_{\mathit{j}}\right),\mathrm{otherwise},\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(9)

where T and R are both the selection functions [35].

After updating the positions and velocities of all the particles, we designed an efficient knowledge transfer of intra-task and inter-task to explore the useful information, which is shown in Algorithm 2. Firstly, two particles are randomly selected from the current generation of particles. In the intra-task transfer, the same positions of the particles focus on exploitation. $\cap ({\mathit{p}}_{\mathit{a}},{\mathit{p}}_{\mathit{b}})$ represents the positions where the elements in ${\mathit{p}}_{\mathit{a}}$ and ${\mathit{p}}_{\mathit{b}}$ are both 1. Then, the new particles directly inherit positions in $\cap ({\mathit{p}}_{\mathit{a}},{\mathit{p}}_{\mathit{b}})$, and the remaining randomly inherit the position on the original particle. On the contrary, the exploration of randomness focuses on the inter-task transfer. $\cup ({\mathit{p}}_{\mathit{a}},{\mathit{p}}_{\mathit{b}})$ represents the positions where the elements are equal to 1 in ${\mathit{p}}_{\mathit{a}}$ or ${\mathit{p}}_{\mathit{b}}$. For the new particles ${\mathit{p}}_{{\mathit{a}}^{{}^{\prime }}}$ and ${\mathit{p}}_{{\mathit{b}}^{{}^{\prime }}}$, $\left|\right|{\mathit{p}}_{\mathit{a}}{\left|\right|}_0$ and $\left|\right|{\mathit{p}}_{\mathit{b}}{\left|\right|}_0$ positions are directly selected from $\cup ({\mathit{p}}_{\mathit{a}},{\mathit{p}}_{\mathit{b}})$, respectively. Then the ${\mathit{p}}_{\mathit{a}}$ and ${\mathit{p}}_{\mathit{b}}$ are updated with the better fitness particles. In order to more intuitively illustrate the essence of Algorithm 2, Figure 2 shows a simple example for the genetic knowledge transfer. Two particles with sparsity of 3 and 4 are selected from the current generation first, in the intra-task transfer, the new particles are updated by inheriting all positions in their same positions which refer to the positions in ${\mathit{p}}_{\mathit{c}}$, then randomly set the rest of positions to 1 to ensure that the sparsity of the new particles is the same as the previous particles. Similar operations are also performed in the inter-task transfer, but the difference is that the new particles are updated by selecting form the positions with 1 in the previous particle, and randomly set them to 1 with the same sparsity.

Algorithm 2 Genetic Knowledge Transfer

Input:${\mathit{P}}^{\mathit{t}}$: the current generation of particles.   1: forg = 1 to N/2 do   2:    Randomly select two particles ${\mathit{p}}_{\mathit{a}}$ and ${\mathit{p}}_{\mathit{b}}$ in ${\mathit{P}}^{\mathit{t}}$.   3:    if ${\tau }_{\mathit{a}}={\tau }_{\mathit{b}}$ then   4:      %Intra-task Transfer   5:      ${\mathit{p}}_{\mathit{c}}\leftarrow \cap ({\mathit{p}}_{\mathit{a}},{\mathit{p}}_{\mathit{b}})$.   6:      ${\mathit{p}}_{{\mathit{a}}^{{}^{\prime }}}\leftarrow$ Inherit all positions of ${\mathit{p}}_{\mathit{c}}$ and randomly set ($\left|\right|{\mathit{p}}_{\mathit{a}}{\left|\right|}_0-\left|\right|{\mathit{p}}_{\mathit{c}}{\left|\right|}_0$) positions to 1.   7:      ${\mathit{p}}_{{\mathit{b}}^{{}^{\prime }}}\leftarrow$ Inherit all positions of ${\mathit{p}}_{\mathit{c}}$ and randomly set ($\left|\right|{\mathit{p}}_{\mathit{b}}{\left|\right|}_0-\left|\right|{\mathit{p}}_{\mathit{c}}{\left|\right|}_0$) positions to 1.   8:    else   9:      %Inter-task Transfer 10:      ${\mathit{p}}_{\mathit{c}}\leftarrow \cup ({\mathit{p}}_{\mathit{a}},{\mathit{p}}_{\mathit{b}})$. 11:      ${\mathit{p}}_{{\mathit{a}}^{{}^{\prime }}}\leftarrow$ Randomly select $\left|\right|{\mathit{p}}_{\mathit{a}}{\left|\right|}_0$ positions in ${\mathit{p}}_{\mathit{c}}$. 12:      ${\mathit{p}}_{{\mathit{b}}^{{}^{\prime }}}\leftarrow$ Randomly select $\left|\right|{\mathit{p}}_{\mathit{b}}{\left|\right|}_0$ positions in ${\mathit{p}}_{\mathit{c}}$. 13:    end if 14:    Evaluate the fitness of ${\mathit{p}}_{{\mathit{a}}^{{}^{\prime }}}$ and ${\mathit{p}}_{{\mathit{b}}^{{}^{\prime }}}$. 15:    Update the ${\mathit{p}}_{\mathit{a}}$ and ${\mathit{p}}_{\mathit{b}}$. 16:    g = g + 1. 17: end for

3.3. An Efficient Local Exploration Strategy with MOEA

After the optimization of multipopulation particle swarms, the set of globally optimal particles with all the sparsity levels on each task $(\mathit{G}=\{{\sum }_{\mathit{j}=1}^{\mathit{K}}\left\{{\sum }_{\mathit{s}=1}^{\mathit{S}}{\left\{{\mathit{gbest}}_{\mathit{s}}\right\}}_{\mathit{j}}\right)$ can be obtained. Two conflicting parameters are included in each particle, that is, the endmember sparsity and the reconstruction error. Therefore, we employ the multiobjective optimization algorithm to facilitate the search process to obtain the optimal points in each task. In the evolutionary multitasking multiobjective framework, the optimized function is expressed as follows:

$$\left\{\begin{array}{c}\{{\mathit{X}}_1^{*},{\mathit{X}}_2^{*},...,{\mathit{X}}_{\mathit{K}}^{*}\}=argmin\left\{\mathit{F}\left({\mathit{X}}_1\right),\mathit{F}\left({\mathit{X}}_2\right),...,\mathit{F}\left({\mathit{X}}_{\mathit{K}}\right)\}\right.,\hfill \\ \mathit{F}\left({\mathit{X}}_{\mathit{j}}\right)=\underset{{\mathit{X}}_{\mathit{j}}}{min}\left(\right||{\mathit{X}}_{\mathit{j}}{\left|\right|}_0,\left|\right|{\mathit{Y}}_{\mathit{j}}-{\mathit{AX}}_{\mathit{j}}{\left|\right|}_{\mathit{F}}),\hfill \end{array}\right.$$

(10)

where the ${\mathit{Y}}_{\mathit{j}}$ and ${\mathit{X}}_{\mathit{j}}$ represent the original image and inversion abundance in the $\mathit{j}$-th task, respectively.

The local exploration strategy processes are in Figure 3. First, the globally optimal particles are transcoded to the first generation of the evolutionary algorithm for the NSGA-II framework. The roulette selection, single-point crossover and bitwise mutation operators are employed to participate in the evolution of multiobjective optimization. Then, the generated offspring are evaluated to update the Pareto front in each task according to the nondominated sorting and crowding distance, and the nondominated solutions are transcoded back to the globally optimal particles.

With the above design, the optimal point in each task can finally be obtained by: ${\mathit{X}}_{\mathit{jv}}^{*}=argmin\left|\right|{\mathit{Y}}_{\mathit{j}}-{\mathit{A}}_{\mathit{v}}{\mathit{X}}_{\mathit{jv}}{\left|\right|}_{\mathit{F}},$ which can be solved simply with the least squares method. Finally, the optimal abundance map obtained from each task constitutes the final inverted abundance map.

4. Experimental Results

4.1. Data Sets

Data 1 provided by Iordache et al. [36] is an image which contains 100 × 100 pixels and 224 bands in each pixel, and the related abundance map of nine endmembers is shown Figure 4. It contains nine randomly selected signatures from a sublibrary of 230 spectral signals, and the fractional abundances are piecewise smooth. Data 2 provided by Tang et al. [37] is an image which contains 64 × 64 pixels and 224 bands in each pixel, the related abundance map of five endmembers is shown Figure 5. It includes five endmembers from a sublibrary of 498 spectral signals, and the fractional abundances are also homogeneous. These two benchmark datasets were tested at different levels of white noise, that is, SNR = 20, 30 and 40 dB. The number of tasks was set to three on these two datasets as recommended in [22]. In order to maintain the fairness of the experiments, all experimental results were taken from the average results of 20 experiments, which is the same as in the comparative method paper.

4.2. Performance Analysis of EMMPSO

In this section, the ablation experiments were performed to demonstrate the effectiveness of the knowledge transfer and the local exploration strategy. The hypervolume indicator was used to compare the evolution process and the convergence procedure of the EMMPSO and the EMMPSO without transfer. Hypervolume was calculated using a reference point 1% larger in every component than the corresponding nadir point [38]. As an important indicator to measure the Pareto-optimal front (PF), the larger the value of the hypervolume, or the faster the convergence speed of the hypervolume, the better the PF obtained by the algorithms. The evolution of the hypervolume indicator is shown in Figure 6. It is clear that, after a few iterations, our method can obtain the higher hypervolume values with the help of the intra-task and inter-task transfer strategy. When several related tasks are optimized simultaneously under the framework of evolutionary multitasking, the convergence rate is improved significantly.

Secondly, to test the efficiency of the local exploration, the performance of EMMPSO and EMMPSO without local exploration was compared. Usually, signal to reconstruction error (SRE) is used to measure the quality of the reconstruction of a signal.

Table 1. Comparison of EMMPSO and EMMPSO without Local Exploration on data 1 and data 2.

Data 1	$\mathbf{SRE} \left(\mathbf{dB}\right)$
	20	30	40
EMMPSO without LE	7.9435	13.3654	22.7536
EMMPSO	8.2783	15.8039	25.2174
Data 2	$\mathbf{SRE} \left(\mathbf{dB}\right)$
	20	30	40
EMMPSO without LE	10.7025	14.7089	17.0224
EMMPSO	12.3572	20.5891	25.7204

shows the SRE (dB) with different noise levels of our proposed method and the EMMPSO without the local exploration on the simulated data. It can be observed that our method can achieve values of SRE (dB) higher than the EMMPSO without local exploration on both simulated datas. It is obvious to see that the local exploration is useful for facilitating the search process to obtain the optimal points.

4.3. Comparing with State-of-Art Algorithms

In order to reflect the superiority of our proposed algorithm, EMMPSO compares with the state-of-art algorithms, including SUnSAL, CLSUnSAL, two-phase multiobjective sparse unmixing (Tp-MOSU) and evolutionary multitasking sparse reconstruction (MTSR). Among them, SUnSAL and CLSUnSAL are the traditional pixel-based and matrix-based processing algorithms. Tp-MOSU and MTSR are two algorithms based on the multiobjective optimization and multitasking optimization, respectively. In order to reflect the advantage of the proposed method, Figure 7 and Figure 8 depict the estimated abundance maps for the endmember 2, 5, 8 on data 1 and the endmember 1, 3, 5 data 2, respectively. The rightmost column represents the abundance map of the real endemembers. The closer the inverted abundance map is to the real abundance map, the better the unmixing performance of the modified algorithm is. It can be seen that the Tp-MOSU, MTSR and EMMPSO exhibit better performances than the other two methods in the similarity with the original abundance map. Although the abundance maps obtained by the Tp-MOSU, MTSR and EMMPSO are similar, the abundance map of EMMPSO has much less noise.

Table 2. Comparison of EMMPSO and other methods on data 1 and data 2.

Data 1	$\mathbf{SRE} \left(\mathbf{dB}\right)$
	20	30	40
SunSAL	4.5568	8.5833	12.9890
CLSUnSAL	5.5164	11.4842	18.7935
Tp-MOSU	8.4083	14.4070	22.5478
MTSR	7.0496	13.7802	22.7329
EMMPSO	8.5783	15.2039	24.2174
Data 2	$\mathbf{SRE} \left(\mathbf{dB}\right)$
	20	30	40
SunSAL	3.5823	8.0323	12.9896
CLSUnSAL	8.2382	13.0988	14.3502
Tp-MOSU	11.3578	15.7132	18.0457
MTSR	10.7254	14.6143	17.6775
EMMPSO	12.3572	20.5891	25.7204

shows the results of SRE (dB) obtained by the five methods on data 1 and data 2. At different levels of noise, the proposed EMMPSO can always achieve the highest values of SRE (dB) on both simulated datasets. The experimental results on two datasets have proved that our proposed EMMPSO is able to achieve a competitive performance by evolutionary multitasking and local exploration strategy.

5. Conclusions

In this paper, we propose a novel evolutionary multitasking multiobjective particle swarm optimization framework called EMMPSO to solve the sparse unmixing problem. With processing multiple homogeneous regions of a hyperspectral image simultaneously, the evolution convergence is accelerated. The local exploration strategy with MOEA is also designed to obtain the optimal solution. For the case study, the proposed EMMPSO is compared with some state-of-the-art methods on benchmark simulated datasets. The results demonstrate the superiority of the EMMPSO.

In future work, we will focus on reducing the time complexity of EMMPSO, and design an efficient multiobjective particle swarm optimization paradigm for the sparse unmixing problem.