# Hyperspectral Image Reconstruction Based on Reference Point Nondominated Sorting Genetic Algorithm

Authors: Li Wang; Wei Wang
Journal: Mobile Information Systems
 (2022)

---

## Abstract

Spatial and spectral features of hyperspectral imagery reconstruction have gained increasing attention in the latest years. Based on the study of orthogonal matching pursuit (OMP) idea, a hyperspectral image reconstruction algorithm based on reference point nondominated sorting genetic algorithm (NSGA) is proposed. Instead of directly reconstructing the entire hyperspectral data as a traditional OMP reconstruction algorithm, the proposed algorithm explores the idea of the evolution process in the reconstruction. The Gabor redundancy dictionary is established as the sparse basis of hyperspectral images, and the reconstruction model of multiobjective optimization is constructed. In the reconstruction process, the NSGA-III algorithm is used to find the optimal atoms to represent the original signal, and Hermitian inversion lemma is also used to realize the recursive update of the residuals. The initial solution generation, the definition of reference points, the association and niche-preservation operation, and the crossover and mutation operation in NSGA-III are presented in detail. Experimental results on hyperspectral data demonstrate that the proposed algorithm could maintain the reconstruction accuracy, as well as the computational efficiency, and are superior to the state-of-the-art reconstruction algorithms. The proposed algorithm could be applied in the classification and unmixing in hyperspectral images.


---

## Body

## 1. Introduction

Hyperspectral images (HSIs) contain rich spatial geometric information and spectral feature information, suitable for target detection, image classification, sparse representation, image denoising, image fusion, biomedicine, and other fields [1]. The progress of hyperspectral imaging spectrometer, the analysis and processing of hyperspectral image data, and the improvement of image interpretation have become frontier research fields that have attracted much attention. The continuous expansion of application fields requires hyperspectral images to provide complete information on the target, to improve the spatial resolution and spectral resolution of the imaging system to achieve more complete spatial and interspectral information [2]. The greater the spatial and spectral resolution, the more difficult it is for the data storage and transmission capacities of the system. With the rapid increase of people’s demand for information [3], data dimensionality reduction caught people’s attention. Furthermore, finding effective data compression and reconstruction algorithms is very important for the growth of hyperspectral remote sensing technology [4].According to the Nyquist sampling theorem, the data acquisition process is to first calculate the highest frequency of the signal, then perform uniform sampling at more than twice the Nyquist rate, and then restore the original data information [5]. The consequence of this is that a huge amount of data is collected and discarded in the compression stage. This method required more storage and transmission resources [6].Compressed sensing (CS) theory was proposed in 2006 [7]. The theory assumes that the information flow of continuous-time signals should be significantly lower than the signal’s maximum frequency and that the degree of freedom of discrete-time signals should be significantly lower than the signal’s length [8]. Therefore, the compressed sensing theory believes that the information characteristics of the signal can be directly collected, completing the reconstruction of the signal [9].According to compressed sensing theory, if the signal has sparsity, the original signal can be remodeled from the measured value [10, 11]. The signal reconstruction process is accomplished by solving convex optimization problems [12], including l1 minimization [13] or gradient projection for sparse reconstruction (GPSR) algorithm [14]. Some other greedy algorithms, such as total variation-based algorithms [15–17], have also been proposed to reconstruct the original signal. Compared to Nyquist-based sampling methods, using CS to acquire and process signals has the advantage of saving transmission and storage resources. It has been widely used in SAR [18–20], airborne [21, 22], classification [23, 24], shadow detection [25], and other fields [26].Analyzing the imaging process of hyperspectral images, we can find that the spatial information and interspectral information of hyperspectral images are redundant, indicating that the images can be compressed. Using compressed sensing theory to reconstruct, classify, and detect hyperspectral images [27–31] can greatly decrease the amount of data. The point of this article is to think about how to restructure the original image from a lesser number of measurements.When using compressed sensing theory to reconstruct hyperspectral images, spatial correlation and interspectral correlation are need to be considered. The same measurement matrix is employed to measure each band image, and after the measurement data is obtained, each band image is reconstructed using prior knowledge [32, 33], such as the dictionary sparse representation algorithm proposed by Geng [34] and the deep residual attention network algorithm proposed by Kohei [35]. Block compression sampling can reduce the storage space required by the measurement matrix at the sampling end, and better reconstructed images can be obtained by using the smoothed projected Landweber algorithm (BCS_SPL) [36, 37]. The same measurement matrix is used to measure each pixel vector [38, 39], and after the measurement data is obtained, principal component analysis [40, 41] or multihypothesis prediction [42] is used to complete the signal reconstruction. Bacca [43] proposed recursive spectral band reconstruction from single pixel hyperspectral measurements.Hyperspectral images also have spectral mixing properties [44, 45], which are described using a linear mixture model. This property suggests that hyperspectral images can be decomposed into endmember matrices and abundance matrices [46, 47], and unmixing using spectral mixing models facilitates subpixel mapping [48] and resolution enhancement [49].A compressed sensing reconstruction algorithm with reference point nondominated sorting genetics is proposed for hyperspectral images. The method samples spatially during the sampling process to acquire available information for reconstruction. In the reconstruction process, a Gabor sparse basis redundancy dictionary and a multiobjective optimization reconstruction model are constructed. The algorithm uses the reference point nondominated sorting genetic algorithm (NSGA-III) [50] to optimize the matching process and finds the optimal atoms. Furthermore, the algorithm uses Hermitian inversion lemma to complete the residual updating in a recursive way. The proposed algorithm is established on the idea of orthogonal matching pursuit (OMP) [51], using the Hermitian matrix characteristics and NSGA-III to obtain the multiobjective optimization reconstruction model; therefore, the algorithm is named HN_OMP. The reconstruction results would be further applied to image classification and interpretation.The remainder of the manuscript is organized into five sections. In Section2, the mathematical model includes the Gabor dictionary basis construction and the multiobjective reconstruction model construction. And then in Section 3, the NSGA-III method is introduced firstly, with its selection operator, crossover operator, mutation operator, and its steps are given. The Hermitian inversion process used to accelerate the residual updating is described. The reconstruction process of algorithm HN_OMP and the diagram are described. In Section 4, experimental datasets and evaluation metrics are presented. The results and comparisons of hyperspectral data are analyzed in Section 5, and Section 6 is about conclusion.


## 2. Mathematical Model

Different from the Nyquist sampling theorem, CS theory does not sample each point of the signal but uses a measurement matrixΦ, to measure sparse signals x linearly, in which M<<N, to obtain the observation value y. The process is represented as
(1)y=Φx=ΦΨs=Θs,whereΘ=ΦΨ is a M×N dimensional matrix, called the CS information operator and Ψ represents the sparse basis of the signal. The sampling rate is defined as SR=M/N. The reconstruction process is to recover the sparse signal s with the known observation value y and information operator Θ, so as to acquire the original signal x. The sampling and reconstruction procedure of compressed sensing theory is given in Figure 1.Figure 1
Sampling and reconstruction process of CS theory.According to the compressed sensing theory that if the sensing matrix fulfills the constraint isometric property, then the measurement is performed in the above manner. The measured value already carries enough information of the original signal. By solving the inverse process of the above measurement process, the original signal can be restored. Because ofM<<N, it is hard to obtain an effective solution x directly from the measurement y . When the signal is sparse, the above problem can be converted into solving the minimum l0 norm problem [52]:
(2)x^=argminx∈RNΨ−1s0,s.t.y=ΦΨs,where⋅−1 denotes the inverse of matrix.However,l0 optimization is still an NP-hard problem. Under certain conditions, the solutions of l1 norm and l0 norm are equivalent, namely:
(3)x^=argminx∈RNΨ−1s1,s.t.y=ΦΨs.After the above analysis, it is found that the key point to reconstruct the signal is how to find the optimal atom which can signify the original sparse signal from the sparse basis. This process consists of two steps: one is to construct a sparse basis, and the second is the method to find the best atom. In the following, we elaborate on the composition of sparse basis and the establishment of reconstruction model.


### 2.1. Sparse Basis Construction

The time-frequency atom analysis approach involves performing a time-frequency joint analysis on a super-complete redundant dictionary made up of time-frequency atoms. The time-frequency atoms can effectively express the localized structure of the signal, which is beneficial to extract useful information from the complex signal. The key of this method is how to select the time-frequency atom library suitable for the local structural characteristics of the signal and how to search for the most characteristic time-frequency atom from the library. By performing time translation, scale transformation, and frequency modulation on standard Gaussian function atoms, a Gabor redundant dictionary with time-frequency characteristics is obtained. Using the Gabor function to generate a redundant dictionary for sparse representation of the signal can ensure the global information of the signal and the intensity of the signal change in any local time. The generating function of Gabor atoms [53] can be expressed using Gaussian window function:
(4)gn=1swinn−uscosυn+ω,winn=e−πn2,wheren=1,2,⋯,N, N indicates the signal length, γ=s,u,υ,ω is time-frequency parameter, and Γ shows the set of parameters γ.The time-frequency parameters are discretized according to the following methods:(5)γ=aα1,α2aα1Δu,α3a−α1Δυ,α4Δω,wherea=2, Δu=1/2, Δυ=π, Δω=π/6, 0<α1≤logN2, 0≤α2≤2−α1+1N, 0≤α3<2α1+1, and0≤α4≤12. The Gabor dictionary is highly redundant, and the number of atoms is Natom=52Nlog2N+N−1.The hyperspectral image is composed of multiple images of different bands. Each band image reflects the spatial distribution information of ground objects. It is the same as the ordinary two-dimensional natural image. There is information redundancy between adjacent pixels, which is sparse and can be used for redundancy. The time-frequency atoms contained in the constructed Gabor redundancy dictionary are various transformations of Gaussian function atoms. Therefore, it could be constructed as the sparse basis to perform sparse approximation for hyperspectral image reconstruction.


### 2.2. Multiobjective Reconstruction Model

Using the basic idea of the orthogonal matching pursuit algorithm, through continuous iteration, we find some optimal atoms that can characterize the sparse signals from the information operator Θ. The residual between the signal and the measured value is minimized, thereby completing the reconstruction of the signal.The redundant dictionary isΨ=gγγ∈Γ, and gγ is the atom defined by the parameter group γ. The process of using OMP to achieve reconstruction is described in Algorithm 1.Algorithm 1: OMP.
Inputs: the measurement value y, the measurement matrix Φ, the sparse basis Ψ, the atom number K.Output: the reconstructed signal x^.1: Initialization parameters, residuals R0=y, optimal atom set Ψ0=, number of iterations k=1.2: Whilek≤Kdo3: Find the atom gk that best matches the residuals from the redundant dictionary, which satisfies gk=maxγ∈ΓRk−1,Φgγ.4: Update the optimal atoms set with gk to obtain Ψk=Ψk−1,gk.5: Update the residual with Rk=y−ΘkΘkTΘk−1ΘkTy, in which Θk=ΦΨk.6: Compute k=k+1.7: end8: Compute the reconstructed signal according x^=ΨKΘKTΘK−1ΘKTyAnalyzing the reconstruction process of OMP, we find that the selection of the optimal atom only considers the correlation between the atoms in the redundant dictionary and the residuals whereas ignores the correlation between the atoms in the redundant dictionary and the selected optimal atom. From the perspective of signal analysis, only when the selected atoms are uncorrelated, the smallest set of uncorrelated atoms representing the signal can be found. Therefore, we choose to build a multiobjective reconstruction model to enhance the accuracy of the algorithm.Afterk−1 iteration, the residual is Rk−1, and the optimal atom set is Θk−1, and then in the kth reconstruction process, the first objective function that needs to be optimized is the maximum inner product of matching atoms with residuals. In other words, the first objective function is computed as
(6)obj1g=minγ∈Γ1Rk−1,Φgγ.At the same time, we require that the currently selected atom is not related to the selected atom. In other words, the inner product between the currently selected atom and the selected atoms should be minimized. Therefore, the second objective function is computed as(7)obj2g=minγ∈ΓΨk−1,gγ.In the multiobjective reconstruction model, we choose the maximum inner product of the matching atom with the image residual and the minimum correlation between the atoms, as the two optimization goal. The multiobjective reconstruction model is expressed as(8)obj=minγ∈Γ1Rk−1,Φgγ,minγ∈ΓΨk−1,gγ.Finding the optimal atom is actually a search in a set of optimal solutions; we can try to simulate this process with a genetic algorithm based on reference point nondominated sorting to explore the feasibility of improving computational efficiency. The detailed steps are described in Section3.In the residual update process of OMP, as the number of atoms increases, the inversion operation of a large matrix appears, and the computational complexity will increase. Through analysis, we found that the matrixΘkTΘk is a positive definite Hermitian matrix. Thence, in Section 3, we will introduce the Hermitian inversion lemma to accelerate its inversion process, thereby improving the calculation efficiency of the residual update process.


## 2.1. Sparse Basis Construction

The time-frequency atom analysis approach involves performing a time-frequency joint analysis on a super-complete redundant dictionary made up of time-frequency atoms. The time-frequency atoms can effectively express the localized structure of the signal, which is beneficial to extract useful information from the complex signal. The key of this method is how to select the time-frequency atom library suitable for the local structural characteristics of the signal and how to search for the most characteristic time-frequency atom from the library. By performing time translation, scale transformation, and frequency modulation on standard Gaussian function atoms, a Gabor redundant dictionary with time-frequency characteristics is obtained. Using the Gabor function to generate a redundant dictionary for sparse representation of the signal can ensure the global information of the signal and the intensity of the signal change in any local time. The generating function of Gabor atoms [53] can be expressed using Gaussian window function:
(4)gn=1swinn−uscosυn+ω,winn=e−πn2,wheren=1,2,⋯,N, N indicates the signal length, γ=s,u,υ,ω is time-frequency parameter, and Γ shows the set of parameters γ.The time-frequency parameters are discretized according to the following methods:(5)γ=aα1,α2aα1Δu,α3a−α1Δυ,α4Δω,wherea=2, Δu=1/2, Δυ=π, Δω=π/6, 0<α1≤logN2, 0≤α2≤2−α1+1N, 0≤α3<2α1+1, and0≤α4≤12. The Gabor dictionary is highly redundant, and the number of atoms is Natom=52Nlog2N+N−1.The hyperspectral image is composed of multiple images of different bands. Each band image reflects the spatial distribution information of ground objects. It is the same as the ordinary two-dimensional natural image. There is information redundancy between adjacent pixels, which is sparse and can be used for redundancy. The time-frequency atoms contained in the constructed Gabor redundancy dictionary are various transformations of Gaussian function atoms. Therefore, it could be constructed as the sparse basis to perform sparse approximation for hyperspectral image reconstruction.


## 2.2. Multiobjective Reconstruction Model

Using the basic idea of the orthogonal matching pursuit algorithm, through continuous iteration, we find some optimal atoms that can characterize the sparse signals from the information operator Θ. The residual between the signal and the measured value is minimized, thereby completing the reconstruction of the signal.The redundant dictionary isΨ=gγγ∈Γ, and gγ is the atom defined by the parameter group γ. The process of using OMP to achieve reconstruction is described in Algorithm 1.Algorithm 1: OMP.
Inputs: the measurement value y, the measurement matrix Φ, the sparse basis Ψ, the atom number K.Output: the reconstructed signal x^.1: Initialization parameters, residuals R0=y, optimal atom set Ψ0=, number of iterations k=1.2: Whilek≤Kdo3: Find the atom gk that best matches the residuals from the redundant dictionary, which satisfies gk=maxγ∈ΓRk−1,Φgγ.4: Update the optimal atoms set with gk to obtain Ψk=Ψk−1,gk.5: Update the residual with Rk=y−ΘkΘkTΘk−1ΘkTy, in which Θk=ΦΨk.6: Compute k=k+1.7: end8: Compute the reconstructed signal according x^=ΨKΘKTΘK−1ΘKTyAnalyzing the reconstruction process of OMP, we find that the selection of the optimal atom only considers the correlation between the atoms in the redundant dictionary and the residuals whereas ignores the correlation between the atoms in the redundant dictionary and the selected optimal atom. From the perspective of signal analysis, only when the selected atoms are uncorrelated, the smallest set of uncorrelated atoms representing the signal can be found. Therefore, we choose to build a multiobjective reconstruction model to enhance the accuracy of the algorithm.Afterk−1 iteration, the residual is Rk−1, and the optimal atom set is Θk−1, and then in the kth reconstruction process, the first objective function that needs to be optimized is the maximum inner product of matching atoms with residuals. In other words, the first objective function is computed as
(6)obj1g=minγ∈Γ1Rk−1,Φgγ.At the same time, we require that the currently selected atom is not related to the selected atom. In other words, the inner product between the currently selected atom and the selected atoms should be minimized. Therefore, the second objective function is computed as(7)obj2g=minγ∈ΓΨk−1,gγ.In the multiobjective reconstruction model, we choose the maximum inner product of the matching atom with the image residual and the minimum correlation between the atoms, as the two optimization goal. The multiobjective reconstruction model is expressed as(8)obj=minγ∈Γ1Rk−1,Φgγ,minγ∈ΓΨk−1,gγ.Finding the optimal atom is actually a search in a set of optimal solutions; we can try to simulate this process with a genetic algorithm based on reference point nondominated sorting to explore the feasibility of improving computational efficiency. The detailed steps are described in Section3.In the residual update process of OMP, as the number of atoms increases, the inversion operation of a large matrix appears, and the computational complexity will increase. Through analysis, we found that the matrixΘkTΘk is a positive definite Hermitian matrix. Thence, in Section 3, we will introduce the Hermitian inversion lemma to accelerate its inversion process, thereby improving the calculation efficiency of the residual update process.


## 3. Proposed Reconstruction Algorithm

Nondominated sorting genetic algorithm is a multiobjective optimization method based on genetic algorithm, which optimizes multiobjectives based on the Pareto frontier. Considering the accuracy and efficiency, NSGA-II [54] using the elite strategy and nondominated quick sort is a good choice. To enhance the performance of the algorithm, the literature [50] uses the target spatial reference point to replace the crowded distance in NSGA-II and proposes NSGA-III. NSGA-III, unlike NSGA-II, helps maintain population diversity by supplying and adaptively updating a number of widely distributed reference points.The main advantage of the algorithm NSGA-III is that, by defining a hyperplane and multiple reference points [55, 56], it selects the individuals which can best adapt the environment, in order to maintain the diversity of the population and ensure that decision-makers can find the optimal solution.Firstly, we describe the NSGA-III algorithm in detail, including the initial solution generation, the definition of reference points, the association and niche-preservation operation, and the crossover and mutation operation. Following that, the Hermitian inversion process is presented, and the implementation of HN_OMP is at the last.


### 3.1. NSGA-III Method

The constructed optimization model is a multiobjective system with many constraints. The traditional NSGA-III algorithm puts forward the concept of constraint control, and formulate that the definite feasible solution must be better than the infeasible solution. However, its completely random generation rules will lead to cause the algorithm to converge slowly. To improve the convergence speed, the method of generating initial solutions is improved in generating initial solutions.


#### 3.1.1. Initial Solution Generation

In light of the evolutionary algorithm’s properties, to solve the multiobjective reconstruction model, there is no need to generate a Gabor dictionary in advance, and a real-valued encoding method is used to construct chromosomes, and each chromosome represents an atom. According to the generating formula of Gabor atom, the length of chromosome isD=4, and its value is the vector represented by s,u,υ,ω. Suppose the population size is P, the maximum evolution generation is G, and the number of objective functions is T.According to the discretization method of thes,u,υ,ω, all atoms are generated to form the initial feasible solution set Zint, the size of Zint is N×Natom, and the element in Zint is expressed in formula (9), where s,u,υ,ω is the natomth time-frequency parameter in set Γ.
(9)Zintn,natom=1se−πn−us2cosυn+ω.We select all alternative atoms as the initial feasible solutionZint and randomly select P solutions as the initial solution in our algorithm, which can meet the constraints and increase the diversity of the solutions. The specific method to obtain the initial solution Z0Par is described in Algorithm 2.Algorithm 2: Initialization.
Inputs: the initial feasible solution Zint, the population size P.Output: the initial solution Z0Par.1: Set the current individual p=1.2: Whilep≤Pdo3: Randomly generate one number Nrand between 1,Natom, set the Nrand th atom in Zint as the p-th initial solution in Z0Par, Z0Par:,p=Zint:,Nrand.4: Compute p=p+1.5: end


#### 3.1.2. Classification of Population into Nondominated Levels

Assume the parent population at this generation isZgPar with size P, and ZgSpr is the offspring population with size P. The first phase is to take the top P members from the combined population ZgCom=ZgPar∪ZgSpr, thus allowing preserve best members of the parent population. To complete this, the combined population ZgCom is sorted based on different nondomination levels using the usual domination principle.The individuals of the population are divided intoL layers according to the nondominated rules [50]. The total population members from nondominated front level F1 to level Fl are denoted as Sg, and the number is Sg. If Sg=P, no further operations are needed, and the next generation is started with Zg+1Par=Sg, or else you need to select PL=P−Sg population members from the Lth level FL.


#### 3.1.3. Identification of Reference Points on a Hyperplane

The user can either predefine the reference points in a structured way or supply them preferentially. In this study, we employ an efficient technique that spaces points on a normalization hyperplane with an intercept of one on each axis and is equally inclined to all objective axes. The number of objective functions isT, the length of chromosome is D, then the number of reference points is J=CT+D−1D, j=1,2,⋯,J, and the reference points can be set on the normalized hyperplane of dimension T−1 [50].We highlight population members who are in some way related with each of these reference sites, in addition to nondominated solutions. The derived solutions are expected to be dispersed on or near the Pareto-optimal front, as the reference points are widely scattered across the complete normalized hyperplane.


#### 3.1.4. Normalized Objective Function

Initially, the ideal point of the populationSg is identified by recognizing the minimum value btmin, for each objective function t=1,2,⋯,T and by building the ideal point b1min,b2min,⋯,bTmin(10)btmin=minobjtz,z∈Sg.Each objective value ofSg is then interpreted by deducting objective objt by ztmin, so that the ideal point of decoded Sg turns out to be a zero vector. This translated objective is computed as
(11)objt′z=objtz−btmin.Thereafter, the extreme point in each objective axis is identified by finding the solutionz∈Sg that makes the following achievement scalarzing function minimum with weight vector ω being the axis direction:
(12)ASFz,ω=maxt=1Tobjt′zωt,forz∈Sg.The interceptat of the tth objective axis and the linear hyperplane can then be computed, and the objective functions can be normalized as follows.
(13)objtnz=objt′zat−btmin=objtz−btminat−btmin,fort=1,2,⋯,T.The structured reference points are already located on the normalized hyperplane. The NSGA-III process maintains a variety in the space covered bySg, because the normalizing procedure and the building of the hyperplane are done at each generation using terminal points. The procedure is described in Algorithm 3.Algorithm 3: Normalize.
Inputs: population Sg, objective function objt,t=1,2,⋯,T, structured reference points Hs.Output: normalized objective function objtn,t=1,2,⋯,T, reference points on normalized hyper-plane Hn.1: fort=1toTdo2: Compute ideal point using formula (10).3: Translate objectives using formula (11).4: Compute extreme point using formula (12).5: endfor6: Compute intercept at for t=1,2,⋯,T.7: Normalize objective functions using formula (13).8: Obtain Hn=Hs and save the points in Hn.


#### 3.1.5. Association Operation

After adaptively normalizing each objective function, we must correlate each population member ofSg with a reference point. By linking each reference point with the origin of the hyperplane, we may create a reference line that relates to each reference point on the hyperplane. Next, the perpendicular distance of each population member of Sg from each of the reference lines is computed. In the normalized objective space, the reference point whose reference line is closest to a population member is considered related with the population member. The procedure is presented in Algorithm 4.


#### 3.1.6. Niche-Preservation Operation

Through association operation, we associate individuals with reference points. It may appear that a reference point is associated with one or more population individuals, or some reference points are not associated with any population individual. At this point, we need to perform niche-preservation operation.Determine the numberρj of first l−1 layer individuals associated with each reference point j, and let the reference points with a smaller number of associations form a set Jmin. Randomly select a reference point from j1 among them, if ρj1=0 and j1 is associated with multiple individuals in lth layer, and then select the individual closest to the reference line j1O to enter the next generation.Ifρj1=0 and j1 is not associated with multiple individuals in lth layer, then randomly select an individual to enter the next generation. Repeatedly select reference points from Jmin, until the individuals are supplemented and the parent population is completed. The procedure is presented in Algorithm 5.Algorithm 4: Associate.
Inputs: the populationSg, the reference points HnOutput: the closest reference pointπz, the distance dz between πz and z1: for each reference pointh∈Hn do2: Compute reference line:k=h3: endfor4: for each population memberz∈Sg do5: for each reference linek∈Hn do6:  Computed⊥z,k=z−kTz/k7: endfor8: Assignπz=k:argmind⊥z,k,k∈Hn9: Assigndz=d⊥z,πz10: endforAlgorithm 5: Niching.
Inputs: the number needed to selectPL, the number associated with the reference point ρj, the closest reference point πz, the distance dz, the reference points Hn, the L-th level FL.Output: the next generation populationZg+1Par1:m=12: Whilem≤PL do3:Jmin=j:argminρj,j∈Hn4:j¯=randomJmin5:Ij¯=z:πz=j¯,z∈FL6: ifIj¯≠φ then7:  ifρj¯=0 then8:Zg+1Par=Zg+1Par∪z:argmindz,z∈Ij¯9:  else10:Zg+1Par=Zg+1Par∪randomIj¯11:  endif12:ρj¯=ρj¯+1,FL=FL\z13:m=m+114: endif15: endwhile


#### 3.1.7. Crossover and Mutation Operation

AfterZg+1Par is formed, crossover and mutation operation are used to create the offspring population Zg+1Spr. According to the crossover probability, it is judged whether the individual has crossover and the gene position where the crossover occurs. The analog binary crossover operator is used to generate the progeny population. According to the mutation probability, determine whether the individual has mutation and determine the gene location where the mutation occurred and adopt the basic mutation method to determine the individual value after mutation.


#### 3.1.8. Steps of NSGA-III

Combined with the description of the above initial solution generation and genetic operator, the implementation steps of NSGA-III is summarized in Algorithm6.After obtaining the optimal solution setZg+1Spr, according to the contribution degree of the individual objective function value in the total objective function value, the individual Zoo=1,2,⋯,P with the smallest value of objo is selected as the optimal atom found in this evolution.
(14)objo=∑t=1TobjtZo−objtZpminobjtZpmin,whereZpp=1,2,⋯,P∈Zg+1Spr.Algorithm 6: NSGA-III.
Inputs: the evolution generation G, the initial parent population Z0Par, the objective function objt,t=1,2,⋯,T, the structured reference points Hs.Output: the optimal solution set Zg+1Spr.1: g=02: Perform crossover and mutation operations to generate the initial offspring population Z0Spr.3: Whileg≤Gdo4: Sg=φ, l=15: ZgCom=ZgPar∪ZgSpr6: F1,F2,⋯,FL=NondominatedsortZgCom7: repeat8: Sg=Sg∪Flandl=l+19: untill=L−110: ifSg=Pthen11: Zg+1Par=Sg,break12: else13:  Zg+1Par=Sg14:  Points to be choosen from FL: PL=P−Sg15:  Normalize the objective functions and create the reference points using Algorithm 3.16:  Associate each member of Sg with a reference point using Algorithm 4.17:  Compute the niche count of the reference point using Algorithm 5 to obtain Zg+1Par.18: endif19: Perform crossover and mutation operations to generate the offspring population Zg+1Spr.20: g=g+121: endWhile


### 3.2. Hermitian Inversion Process

The Hermitian inversion lemma mainly uses the characteristics of positive definite Hermitian matrices to improve the speed of large-scale matrix inversion operations. The block form of the Hermitian matrix is(15)Hk=Hk‐1hkhkTρk,whereHk−1 and Hk‐1−1 exists, if let bk=−Hk‐1−1hk, then
(16)Hk−1=Hk‐1−10k0kT0+1ρk+hkTbkbkbkTbkbkT1.During the iteration process, the update of the optimal atoms set is expressed asΨk=Ψk−1,gk, where Ψk−1 is the set of atoms formed after the first k‐1 iteration and gk is the optimal atom found in the kth iteration, and then Θk=Θk‐1,dk, where Θk‐1=ΦΨk−1 and dk=Φgk.LetHk−1=Θk−1TΘk−1, then
(17)Hk=ΘkTΘk=Θk−1TdkTΘk−1dk=Θk−1TΘk−1Θk−1TdkdkTΘk−1dkTdk.Compared with formula (15)
(18)hk=Θk−1Tdk,ρk=dkTdk.SubstitutingHk−1 and Hk into equation (16), the recursive expression of matrix inversion in the residual update process is obtained.
(19)Hk−1=Θk−1TΘk−1‐10k0kT0+1dkTdk−dkTΘk−1Θk−1TΘk−1‐1Θk−1Tdk⋅Θk−1TΘk−1‐1Θk−1TdkdkTΘk−1Θk−1TΘk−1‐T−Θk−1TΘk−1‐1Θk−1Tdk−dkTΘk−1Θk−1TΘk−1‐T1.


### 3.3. Implementation of HN_OMP

The implementation HN_OMP algorithm is summarized in Algorithm7.Algorithm 7: HN_OMP.
Inputs: the measurement valuey, the measurement matrix Φ, the atom number K, the parameters of NSGA (including the population size P, the maximum evolution generation G).Output: the reconstructed signalx^.1: Initialize residualR0=y, the optimal atoms set Ψ0=, and the current atom number k=1.2: Whilek≤K do3: According to Algorithm6: NSGA-III, with the parameters P and G, the optimal individual Zo is obtained as the optimal atom gk.3: Usegk to update the optimal atoms set Ψk=Ψk−1,gk.4: Update residual with Hermitian inversion, using formula (19) to obtainRk=y−ΘkΘkTΘk−1ΘkTy.5: Computek=k+1.6: end7: Compute reconstructed signalx^ according to x^=ΨKΘKTΘK−1ΘKTy


## 3.1. NSGA-III Method

The constructed optimization model is a multiobjective system with many constraints. The traditional NSGA-III algorithm puts forward the concept of constraint control, and formulate that the definite feasible solution must be better than the infeasible solution. However, its completely random generation rules will lead to cause the algorithm to converge slowly. To improve the convergence speed, the method of generating initial solutions is improved in generating initial solutions.


### 3.1.1. Initial Solution Generation

In light of the evolutionary algorithm’s properties, to solve the multiobjective reconstruction model, there is no need to generate a Gabor dictionary in advance, and a real-valued encoding method is used to construct chromosomes, and each chromosome represents an atom. According to the generating formula of Gabor atom, the length of chromosome isD=4, and its value is the vector represented by s,u,υ,ω. Suppose the population size is P, the maximum evolution generation is G, and the number of objective functions is T.According to the discretization method of thes,u,υ,ω, all atoms are generated to form the initial feasible solution set Zint, the size of Zint is N×Natom, and the element in Zint is expressed in formula (9), where s,u,υ,ω is the natomth time-frequency parameter in set Γ.
(9)Zintn,natom=1se−πn−us2cosυn+ω.We select all alternative atoms as the initial feasible solutionZint and randomly select P solutions as the initial solution in our algorithm, which can meet the constraints and increase the diversity of the solutions. The specific method to obtain the initial solution Z0Par is described in Algorithm 2.Algorithm 2: Initialization.
Inputs: the initial feasible solution Zint, the population size P.Output: the initial solution Z0Par.1: Set the current individual p=1.2: Whilep≤Pdo3: Randomly generate one number Nrand between 1,Natom, set the Nrand th atom in Zint as the p-th initial solution in Z0Par, Z0Par:,p=Zint:,Nrand.4: Compute p=p+1.5: end


### 3.1.2. Classification of Population into Nondominated Levels

Assume the parent population at this generation isZgPar with size P, and ZgSpr is the offspring population with size P. The first phase is to take the top P members from the combined population ZgCom=ZgPar∪ZgSpr, thus allowing preserve best members of the parent population. To complete this, the combined population ZgCom is sorted based on different nondomination levels using the usual domination principle.The individuals of the population are divided intoL layers according to the nondominated rules [50]. The total population members from nondominated front level F1 to level Fl are denoted as Sg, and the number is Sg. If Sg=P, no further operations are needed, and the next generation is started with Zg+1Par=Sg, or else you need to select PL=P−Sg population members from the Lth level FL.


### 3.1.3. Identification of Reference Points on a Hyperplane

The user can either predefine the reference points in a structured way or supply them preferentially. In this study, we employ an efficient technique that spaces points on a normalization hyperplane with an intercept of one on each axis and is equally inclined to all objective axes. The number of objective functions isT, the length of chromosome is D, then the number of reference points is J=CT+D−1D, j=1,2,⋯,J, and the reference points can be set on the normalized hyperplane of dimension T−1 [50].We highlight population members who are in some way related with each of these reference sites, in addition to nondominated solutions. The derived solutions are expected to be dispersed on or near the Pareto-optimal front, as the reference points are widely scattered across the complete normalized hyperplane.


### 3.1.4. Normalized Objective Function

Initially, the ideal point of the populationSg is identified by recognizing the minimum value btmin, for each objective function t=1,2,⋯,T and by building the ideal point b1min,b2min,⋯,bTmin(10)btmin=minobjtz,z∈Sg.Each objective value ofSg is then interpreted by deducting objective objt by ztmin, so that the ideal point of decoded Sg turns out to be a zero vector. This translated objective is computed as
(11)objt′z=objtz−btmin.Thereafter, the extreme point in each objective axis is identified by finding the solutionz∈Sg that makes the following achievement scalarzing function minimum with weight vector ω being the axis direction:
(12)ASFz,ω=maxt=1Tobjt′zωt,forz∈Sg.The interceptat of the tth objective axis and the linear hyperplane can then be computed, and the objective functions can be normalized as follows.
(13)objtnz=objt′zat−btmin=objtz−btminat−btmin,fort=1,2,⋯,T.The structured reference points are already located on the normalized hyperplane. The NSGA-III process maintains a variety in the space covered bySg, because the normalizing procedure and the building of the hyperplane are done at each generation using terminal points. The procedure is described in Algorithm 3.Algorithm 3: Normalize.
Inputs: population Sg, objective function objt,t=1,2,⋯,T, structured reference points Hs.Output: normalized objective function objtn,t=1,2,⋯,T, reference points on normalized hyper-plane Hn.1: fort=1toTdo2: Compute ideal point using formula (10).3: Translate objectives using formula (11).4: Compute extreme point using formula (12).5: endfor6: Compute intercept at for t=1,2,⋯,T.7: Normalize objective functions using formula (13).8: Obtain Hn=Hs and save the points in Hn.


### 3.1.5. Association Operation

After adaptively normalizing each objective function, we must correlate each population member ofSg with a reference point. By linking each reference point with the origin of the hyperplane, we may create a reference line that relates to each reference point on the hyperplane. Next, the perpendicular distance of each population member of Sg from each of the reference lines is computed. In the normalized objective space, the reference point whose reference line is closest to a population member is considered related with the population member. The procedure is presented in Algorithm 4.


### 3.1.6. Niche-Preservation Operation

Through association operation, we associate individuals with reference points. It may appear that a reference point is associated with one or more population individuals, or some reference points are not associated with any population individual. At this point, we need to perform niche-preservation operation.Determine the numberρj of first l−1 layer individuals associated with each reference point j, and let the reference points with a smaller number of associations form a set Jmin. Randomly select a reference point from j1 among them, if ρj1=0 and j1 is associated with multiple individuals in lth layer, and then select the individual closest to the reference line j1O to enter the next generation.Ifρj1=0 and j1 is not associated with multiple individuals in lth layer, then randomly select an individual to enter the next generation. Repeatedly select reference points from Jmin, until the individuals are supplemented and the parent population is completed. The procedure is presented in Algorithm 5.Algorithm 4: Associate.
Inputs: the populationSg, the reference points HnOutput: the closest reference pointπz, the distance dz between πz and z1: for each reference pointh∈Hn do2: Compute reference line:k=h3: endfor4: for each population memberz∈Sg do5: for each reference linek∈Hn do6:  Computed⊥z,k=z−kTz/k7: endfor8: Assignπz=k:argmind⊥z,k,k∈Hn9: Assigndz=d⊥z,πz10: endforAlgorithm 5: Niching.
Inputs: the number needed to selectPL, the number associated with the reference point ρj, the closest reference point πz, the distance dz, the reference points Hn, the L-th level FL.Output: the next generation populationZg+1Par1:m=12: Whilem≤PL do3:Jmin=j:argminρj,j∈Hn4:j¯=randomJmin5:Ij¯=z:πz=j¯,z∈FL6: ifIj¯≠φ then7:  ifρj¯=0 then8:Zg+1Par=Zg+1Par∪z:argmindz,z∈Ij¯9:  else10:Zg+1Par=Zg+1Par∪randomIj¯11:  endif12:ρj¯=ρj¯+1,FL=FL\z13:m=m+114: endif15: endwhile


### 3.1.7. Crossover and Mutation Operation

AfterZg+1Par is formed, crossover and mutation operation are used to create the offspring population Zg+1Spr. According to the crossover probability, it is judged whether the individual has crossover and the gene position where the crossover occurs. The analog binary crossover operator is used to generate the progeny population. According to the mutation probability, determine whether the individual has mutation and determine the gene location where the mutation occurred and adopt the basic mutation method to determine the individual value after mutation.


### 3.1.8. Steps of NSGA-III

Combined with the description of the above initial solution generation and genetic operator, the implementation steps of NSGA-III is summarized in Algorithm6.After obtaining the optimal solution setZg+1Spr, according to the contribution degree of the individual objective function value in the total objective function value, the individual Zoo=1,2,⋯,P with the smallest value of objo is selected as the optimal atom found in this evolution.
(14)objo=∑t=1TobjtZo−objtZpminobjtZpmin,whereZpp=1,2,⋯,P∈Zg+1Spr.Algorithm 6: NSGA-III.
Inputs: the evolution generation G, the initial parent population Z0Par, the objective function objt,t=1,2,⋯,T, the structured reference points Hs.Output: the optimal solution set Zg+1Spr.1: g=02: Perform crossover and mutation operations to generate the initial offspring population Z0Spr.3: Whileg≤Gdo4: Sg=φ, l=15: ZgCom=ZgPar∪ZgSpr6: F1,F2,⋯,FL=NondominatedsortZgCom7: repeat8: Sg=Sg∪Flandl=l+19: untill=L−110: ifSg=Pthen11: Zg+1Par=Sg,break12: else13:  Zg+1Par=Sg14:  Points to be choosen from FL: PL=P−Sg15:  Normalize the objective functions and create the reference points using Algorithm 3.16:  Associate each member of Sg with a reference point using Algorithm 4.17:  Compute the niche count of the reference point using Algorithm 5 to obtain Zg+1Par.18: endif19: Perform crossover and mutation operations to generate the offspring population Zg+1Spr.20: g=g+121: endWhile


## 3.1.1. Initial Solution Generation

In light of the evolutionary algorithm’s properties, to solve the multiobjective reconstruction model, there is no need to generate a Gabor dictionary in advance, and a real-valued encoding method is used to construct chromosomes, and each chromosome represents an atom. According to the generating formula of Gabor atom, the length of chromosome isD=4, and its value is the vector represented by s,u,υ,ω. Suppose the population size is P, the maximum evolution generation is G, and the number of objective functions is T.According to the discretization method of thes,u,υ,ω, all atoms are generated to form the initial feasible solution set Zint, the size of Zint is N×Natom, and the element in Zint is expressed in formula (9), where s,u,υ,ω is the natomth time-frequency parameter in set Γ.
(9)Zintn,natom=1se−πn−us2cosυn+ω.We select all alternative atoms as the initial feasible solutionZint and randomly select P solutions as the initial solution in our algorithm, which can meet the constraints and increase the diversity of the solutions. The specific method to obtain the initial solution Z0Par is described in Algorithm 2.Algorithm 2: Initialization.
Inputs: the initial feasible solution Zint, the population size P.Output: the initial solution Z0Par.1: Set the current individual p=1.2: Whilep≤Pdo3: Randomly generate one number Nrand between 1,Natom, set the Nrand th atom in Zint as the p-th initial solution in Z0Par, Z0Par:,p=Zint:,Nrand.4: Compute p=p+1.5: end


## 3.1.2. Classification of Population into Nondominated Levels

Assume the parent population at this generation isZgPar with size P, and ZgSpr is the offspring population with size P. The first phase is to take the top P members from the combined population ZgCom=ZgPar∪ZgSpr, thus allowing preserve best members of the parent population. To complete this, the combined population ZgCom is sorted based on different nondomination levels using the usual domination principle.The individuals of the population are divided intoL layers according to the nondominated rules [50]. The total population members from nondominated front level F1 to level Fl are denoted as Sg, and the number is Sg. If Sg=P, no further operations are needed, and the next generation is started with Zg+1Par=Sg, or else you need to select PL=P−Sg population members from the Lth level FL.


## 3.1.3. Identification of Reference Points on a Hyperplane

The user can either predefine the reference points in a structured way or supply them preferentially. In this study, we employ an efficient technique that spaces points on a normalization hyperplane with an intercept of one on each axis and is equally inclined to all objective axes. The number of objective functions isT, the length of chromosome is D, then the number of reference points is J=CT+D−1D, j=1,2,⋯,J, and the reference points can be set on the normalized hyperplane of dimension T−1 [50].We highlight population members who are in some way related with each of these reference sites, in addition to nondominated solutions. The derived solutions are expected to be dispersed on or near the Pareto-optimal front, as the reference points are widely scattered across the complete normalized hyperplane.


## 3.1.4. Normalized Objective Function

Initially, the ideal point of the populationSg is identified by recognizing the minimum value btmin, for each objective function t=1,2,⋯,T and by building the ideal point b1min,b2min,⋯,bTmin(10)btmin=minobjtz,z∈Sg.Each objective value ofSg is then interpreted by deducting objective objt by ztmin, so that the ideal point of decoded Sg turns out to be a zero vector. This translated objective is computed as
(11)objt′z=objtz−btmin.Thereafter, the extreme point in each objective axis is identified by finding the solutionz∈Sg that makes the following achievement scalarzing function minimum with weight vector ω being the axis direction:
(12)ASFz,ω=maxt=1Tobjt′zωt,forz∈Sg.The interceptat of the tth objective axis and the linear hyperplane can then be computed, and the objective functions can be normalized as follows.
(13)objtnz=objt′zat−btmin=objtz−btminat−btmin,fort=1,2,⋯,T.The structured reference points are already located on the normalized hyperplane. The NSGA-III process maintains a variety in the space covered bySg, because the normalizing procedure and the building of the hyperplane are done at each generation using terminal points. The procedure is described in Algorithm 3.Algorithm 3: Normalize.
Inputs: population Sg, objective function objt,t=1,2,⋯,T, structured reference points Hs.Output: normalized objective function objtn,t=1,2,⋯,T, reference points on normalized hyper-plane Hn.1: fort=1toTdo2: Compute ideal point using formula (10).3: Translate objectives using formula (11).4: Compute extreme point using formula (12).5: endfor6: Compute intercept at for t=1,2,⋯,T.7: Normalize objective functions using formula (13).8: Obtain Hn=Hs and save the points in Hn.


## 3.1.5. Association Operation

After adaptively normalizing each objective function, we must correlate each population member ofSg with a reference point. By linking each reference point with the origin of the hyperplane, we may create a reference line that relates to each reference point on the hyperplane. Next, the perpendicular distance of each population member of Sg from each of the reference lines is computed. In the normalized objective space, the reference point whose reference line is closest to a population member is considered related with the population member. The procedure is presented in Algorithm 4.


## 3.1.6. Niche-Preservation Operation

Through association operation, we associate individuals with reference points. It may appear that a reference point is associated with one or more population individuals, or some reference points are not associated with any population individual. At this point, we need to perform niche-preservation operation.Determine the numberρj of first l−1 layer individuals associated with each reference point j, and let the reference points with a smaller number of associations form a set Jmin. Randomly select a reference point from j1 among them, if ρj1=0 and j1 is associated with multiple individuals in lth layer, and then select the individual closest to the reference line j1O to enter the next generation.Ifρj1=0 and j1 is not associated with multiple individuals in lth layer, then randomly select an individual to enter the next generation. Repeatedly select reference points from Jmin, until the individuals are supplemented and the parent population is completed. The procedure is presented in Algorithm 5.Algorithm 4: Associate.
Inputs: the populationSg, the reference points HnOutput: the closest reference pointπz, the distance dz between πz and z1: for each reference pointh∈Hn do2: Compute reference line:k=h3: endfor4: for each population memberz∈Sg do5: for each reference linek∈Hn do6:  Computed⊥z,k=z−kTz/k7: endfor8: Assignπz=k:argmind⊥z,k,k∈Hn9: Assigndz=d⊥z,πz10: endforAlgorithm 5: Niching.
Inputs: the number needed to selectPL, the number associated with the reference point ρj, the closest reference point πz, the distance dz, the reference points Hn, the L-th level FL.Output: the next generation populationZg+1Par1:m=12: Whilem≤PL do3:Jmin=j:argminρj,j∈Hn4:j¯=randomJmin5:Ij¯=z:πz=j¯,z∈FL6: ifIj¯≠φ then7:  ifρj¯=0 then8:Zg+1Par=Zg+1Par∪z:argmindz,z∈Ij¯9:  else10:Zg+1Par=Zg+1Par∪randomIj¯11:  endif12:ρj¯=ρj¯+1,FL=FL\z13:m=m+114: endif15: endwhile


## 3.1.7. Crossover and Mutation Operation

AfterZg+1Par is formed, crossover and mutation operation are used to create the offspring population Zg+1Spr. According to the crossover probability, it is judged whether the individual has crossover and the gene position where the crossover occurs. The analog binary crossover operator is used to generate the progeny population. According to the mutation probability, determine whether the individual has mutation and determine the gene location where the mutation occurred and adopt the basic mutation method to determine the individual value after mutation.


## 3.1.8. Steps of NSGA-III

Combined with the description of the above initial solution generation and genetic operator, the implementation steps of NSGA-III is summarized in Algorithm6.After obtaining the optimal solution setZg+1Spr, according to the contribution degree of the individual objective function value in the total objective function value, the individual Zoo=1,2,⋯,P with the smallest value of objo is selected as the optimal atom found in this evolution.
(14)objo=∑t=1TobjtZo−objtZpminobjtZpmin,whereZpp=1,2,⋯,P∈Zg+1Spr.Algorithm 6: NSGA-III.
Inputs: the evolution generation G, the initial parent population Z0Par, the objective function objt,t=1,2,⋯,T, the structured reference points Hs.Output: the optimal solution set Zg+1Spr.1: g=02: Perform crossover and mutation operations to generate the initial offspring population Z0Spr.3: Whileg≤Gdo4: Sg=φ, l=15: ZgCom=ZgPar∪ZgSpr6: F1,F2,⋯,FL=NondominatedsortZgCom7: repeat8: Sg=Sg∪Flandl=l+19: untill=L−110: ifSg=Pthen11: Zg+1Par=Sg,break12: else13:  Zg+1Par=Sg14:  Points to be choosen from FL: PL=P−Sg15:  Normalize the objective functions and create the reference points using Algorithm 3.16:  Associate each member of Sg with a reference point using Algorithm 4.17:  Compute the niche count of the reference point using Algorithm 5 to obtain Zg+1Par.18: endif19: Perform crossover and mutation operations to generate the offspring population Zg+1Spr.20: g=g+121: endWhile


## 3.2. Hermitian Inversion Process

The Hermitian inversion lemma mainly uses the characteristics of positive definite Hermitian matrices to improve the speed of large-scale matrix inversion operations. The block form of the Hermitian matrix is(15)Hk=Hk‐1hkhkTρk,whereHk−1 and Hk‐1−1 exists, if let bk=−Hk‐1−1hk, then
(16)Hk−1=Hk‐1−10k0kT0+1ρk+hkTbkbkbkTbkbkT1.During the iteration process, the update of the optimal atoms set is expressed asΨk=Ψk−1,gk, where Ψk−1 is the set of atoms formed after the first k‐1 iteration and gk is the optimal atom found in the kth iteration, and then Θk=Θk‐1,dk, where Θk‐1=ΦΨk−1 and dk=Φgk.LetHk−1=Θk−1TΘk−1, then
(17)Hk=ΘkTΘk=Θk−1TdkTΘk−1dk=Θk−1TΘk−1Θk−1TdkdkTΘk−1dkTdk.Compared with formula (15)
(18)hk=Θk−1Tdk,ρk=dkTdk.SubstitutingHk−1 and Hk into equation (16), the recursive expression of matrix inversion in the residual update process is obtained.
(19)Hk−1=Θk−1TΘk−1‐10k0kT0+1dkTdk−dkTΘk−1Θk−1TΘk−1‐1Θk−1Tdk⋅Θk−1TΘk−1‐1Θk−1TdkdkTΘk−1Θk−1TΘk−1‐T−Θk−1TΘk−1‐1Θk−1Tdk−dkTΘk−1Θk−1TΘk−1‐T1.


## 3.3. Implementation of HN_OMP

The implementation HN_OMP algorithm is summarized in Algorithm7.Algorithm 7: HN_OMP.
Inputs: the measurement valuey, the measurement matrix Φ, the atom number K, the parameters of NSGA (including the population size P, the maximum evolution generation G).Output: the reconstructed signalx^.1: Initialize residualR0=y, the optimal atoms set Ψ0=, and the current atom number k=1.2: Whilek≤K do3: According to Algorithm6: NSGA-III, with the parameters P and G, the optimal individual Zo is obtained as the optimal atom gk.3: Usegk to update the optimal atoms set Ψk=Ψk−1,gk.4: Update residual with Hermitian inversion, using formula (19) to obtainRk=y−ΘkΘkTΘk−1ΘkTy.5: Computek=k+1.6: end7: Compute reconstructed signalx^ according to x^=ΨKΘKTΘK−1ΘKTy


## 4. Experimental Datasets and Evaluation Metrices

This section describes the datasets used in the tests, and the metrics are used to compare the proposed algorithm’s performance.


### 4.1. Hyperspectral Datasets

Four hyperspectral datasets [57, 58] with diverse signatures are included in the investigations, allowing for a complete quantitative and quantitative comparative evaluation of the proposed scheme. The Airborne Visible Infrared Imaging Spectrometer collected the Cuprite1 and Cuprite2 images, which correspond to the first two scenes in the online data (AVIRIS). The sensor collects 224 spectral bands between 0.4 and 2.5 m, with a half-maximum of 10 nm and spatial resolution of 20 m per pixel. Due to water absorption or noisy, several bands were removed from the study, leaving a total of 188 bands for the experiments [46].The Pines photos are divided into 128 by 128 pixels to avoid impacting the reconstruction algorithm’s performance, whereas the other datasets are represented in 256 by 256 pixels. The basic situation of the four datasets is shown in Table1. Extract three band images; the pseudo-color images of four datasets are shown in Figure 2. In the next discussions, we use X=X1,X2,⋯Xl,⋯,XL∈RN×L to represent the hyperspectral images, N is the pixels of each band image, and L is the number of bands, Xl is the vector representation of one band image, and Xn is the vector representation of one pixel.Table 1
The basic situation of four datasets.
DatasetOriginal bandsOriginal image sizeAvailable bands (remove the water absorption and noisy bands)Cropped image size (spatially cropping for computation convenience)Cuprite1224614∗512188256∗256Cuprite2224614∗512188256∗256Indian Pines220145∗145200128∗128Pavia University115610∗340103256∗256Figure 2
Pseudo-color images of four datasets: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)AVIRIS in northwest Indiana acquired the Indian Pines photos, which have a spatial resolution of 145 by 145 pixels and 220 bands with a wavelength range of 0.4 to 2.5 m. Similarly, the bands including the water absorption area have been removed, leaving 200 bands open. The images obtained by the Reflective Optics System Imaging Spectrometer (ROSIS) at Pavia University comprise 115 bands and a spatial resolution of 1.3 meter per pixel. The 12 noisiest bands in this dataset should be eliminated before analysis, leaving 103 bands in the trials [46].


### 4.2. Evaluation Metrics

Different metrics are adopted to evaluate the performance of the proposed algorithm, the average peak signal-to-noise ratio (PSNR) [46], structural similarity (SSIM) [59], vector-based SNR [42], and spectral-angle distortion (SAD) [42].The band-based PSNR measured in dB is defined as(20)PSNRXl,X^l=20log10maxXlMSEXl,X^l,whereXl and X^l are the original and reconstructed band images, maxXl is the peak value of Xl, and MSEXl,X^l is the mean squared error [46] and is computed as
(21)MSEXl,X^l=1NXl−X^l22.The band-based SSIM betweenXl and X^l is defined as [46]
(22)SSIMXl,X^l=2μ1μ2+C12σ12+C2μ12+μ22+C1σ12+σ22+C2,whereμ1 and μ2 are the mean values of Xl and X^l, σ1 and σ2 are the standard deviation values of Xl and X^l, σ12 represents the correlation coefficient between Xl and X^l, and C1 and C2 are constants related to the dynamic range of the pixel values. The details for these parameters can refer to [59].The average PSNR and average SSIM are calculated by averaging the band-based PSNR and band-based SSIM over all bands of the hyperspectral dataset.Vector-based SNR measured in dB is defined as(23)SNRXn,X^n=10log10varXnMSEXn,X^n,whereXn and X^n shows the original and rebuilt spectral vectors, varXn shows the variance of Xn, and the mean squared error is MSEXn,X^n.The vector-based SAD measured in degrees is presented as(24)ξn=∠Xn−X^n.The average SNR and average SAD are calculated by averaging the vector-based SNR and vector-based SAD over all spectral vectors in the hyperspectral dataset.


## 4.1. Hyperspectral Datasets

Four hyperspectral datasets [57, 58] with diverse signatures are included in the investigations, allowing for a complete quantitative and quantitative comparative evaluation of the proposed scheme. The Airborne Visible Infrared Imaging Spectrometer collected the Cuprite1 and Cuprite2 images, which correspond to the first two scenes in the online data (AVIRIS). The sensor collects 224 spectral bands between 0.4 and 2.5 m, with a half-maximum of 10 nm and spatial resolution of 20 m per pixel. Due to water absorption or noisy, several bands were removed from the study, leaving a total of 188 bands for the experiments [46].The Pines photos are divided into 128 by 128 pixels to avoid impacting the reconstruction algorithm’s performance, whereas the other datasets are represented in 256 by 256 pixels. The basic situation of the four datasets is shown in Table1. Extract three band images; the pseudo-color images of four datasets are shown in Figure 2. In the next discussions, we use X=X1,X2,⋯Xl,⋯,XL∈RN×L to represent the hyperspectral images, N is the pixels of each band image, and L is the number of bands, Xl is the vector representation of one band image, and Xn is the vector representation of one pixel.Table 1
The basic situation of four datasets.
DatasetOriginal bandsOriginal image sizeAvailable bands (remove the water absorption and noisy bands)Cropped image size (spatially cropping for computation convenience)Cuprite1224614∗512188256∗256Cuprite2224614∗512188256∗256Indian Pines220145∗145200128∗128Pavia University115610∗340103256∗256Figure 2
Pseudo-color images of four datasets: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)AVIRIS in northwest Indiana acquired the Indian Pines photos, which have a spatial resolution of 145 by 145 pixels and 220 bands with a wavelength range of 0.4 to 2.5 m. Similarly, the bands including the water absorption area have been removed, leaving 200 bands open. The images obtained by the Reflective Optics System Imaging Spectrometer (ROSIS) at Pavia University comprise 115 bands and a spatial resolution of 1.3 meter per pixel. The 12 noisiest bands in this dataset should be eliminated before analysis, leaving 103 bands in the trials [46].


## 4.2. Evaluation Metrics

Different metrics are adopted to evaluate the performance of the proposed algorithm, the average peak signal-to-noise ratio (PSNR) [46], structural similarity (SSIM) [59], vector-based SNR [42], and spectral-angle distortion (SAD) [42].The band-based PSNR measured in dB is defined as(20)PSNRXl,X^l=20log10maxXlMSEXl,X^l,whereXl and X^l are the original and reconstructed band images, maxXl is the peak value of Xl, and MSEXl,X^l is the mean squared error [46] and is computed as
(21)MSEXl,X^l=1NXl−X^l22.The band-based SSIM betweenXl and X^l is defined as [46]
(22)SSIMXl,X^l=2μ1μ2+C12σ12+C2μ12+μ22+C1σ12+σ22+C2,whereμ1 and μ2 are the mean values of Xl and X^l, σ1 and σ2 are the standard deviation values of Xl and X^l, σ12 represents the correlation coefficient between Xl and X^l, and C1 and C2 are constants related to the dynamic range of the pixel values. The details for these parameters can refer to [59].The average PSNR and average SSIM are calculated by averaging the band-based PSNR and band-based SSIM over all bands of the hyperspectral dataset.Vector-based SNR measured in dB is defined as(23)SNRXn,X^n=10log10varXnMSEXn,X^n,whereXn and X^n shows the original and rebuilt spectral vectors, varXn shows the variance of Xn, and the mean squared error is MSEXn,X^n.The vector-based SAD measured in degrees is presented as(24)ξn=∠Xn−X^n.The average SNR and average SAD are calculated by averaging the vector-based SNR and vector-based SAD over all spectral vectors in the hyperspectral dataset.


## 5. Experimental Results and Comparisons

This section describes the experimental results and findings used to assess the proposed algorithm’s performance. The parameters are firstly determined, including population size and evolution generation. Following that, a comparison is given between the proposed algorithm HN_OMP and existing reconstruction algorithms. AMD quad-core CPU (3.80GHz), 16G memory, and Matlab2012b are the hardware and software environments employed for these studies.


### 5.1. Parameter Selection for Reconstruction Algorithm

The parameter selection would undoubtedly have a major impact on the recon-construction performance. The results of the proposed algorithm HN_OMP on four hyperspectral pictures with varying evolution parameters are reported in this section.


#### 5.1.1. Parameters Selection of Population Size and Evolution Generation

Firstly, the algorithm HN_OMP is used to compress reconstruct the 40th band image of the 4 groups of hyperspectral data, and the influence of the maximum evolution generation, population size, and atom number on the performance of the algorithm is analyzed. The variation range of population sizeP is 10~50 at the interval of 5, the variation range of the maximum evolution generation G is 10~50 at the interval of 5, and the variation range of the atom number K is 10~100 at the interval of 10. The sampling rates are set as 0.1 to 0.5 at the interval of 0.1.Figure3 shows the variation of the average PSNR of the reconstructed image of Cuprite1 with the population size and evolution generation, when the atom number is 40. The sampling rate in Figure 3(a) is 0.3, while the sampling rate in Figure 3(b) is 0.4. From the two subfigures, we can find that under the same evolution generation, after the population size increases to 20, the PSNR tends to be stable no matter what the sampling rate is. Under the same population size, the PSNR has a small oscillation with the increase of evolution generation. It demonstrates that the reconstructed PSNR is more sensitive to population size than that to evolution generation. Without loss of generality, we set the population size as 20 to discuss the influence of other parameters on the reconstruction results.Figure 3
Influence of evolution generation and population size on reconstructed PSNR, when the atom number is 40: (a) sampling rate is 0.3; (b) sampling rate is 0.4.
(a)(b)When the population size is 20, the influence of the evolution generation and atom number on the reconstructed PSNR of dataset Cuprite1 is shown in Figure4. The sampling rate in Figure 4(a) is 0.2, while the sampling rate in Figure 4(b) is 0.3. From the two subfigures, we can clearly see that the influence of the atom number on the reconstruction result is far greater than the influence of the evolution generation on the reconstruction result. As the number of atoms increases, the reconstructed PSNR steadily increases.Figure 4
Influence of evolution generation and atom number on reconstructed PSNR, when the population size is 20: (a) sampling rate is 0.2; (b) sampling rate is 0.3.
(a)(b)However, when the number of atoms further increases to more than 50, the reconstruction performance no longer increases significantly, but gradually stabilizes. The reason for this situation may be that particularly large number of atoms will affect the selecting of elite atoms. Therefore, the algorithm performance cannot be further improved. Comparing Figure3(a) and Figure 4(a), we find that the atom number and the population size have a relatively large impact on the results, and the evolution generation has a small impact. Without loss of generality, we set the evolution generation to 10 to compare the effects of the atom number and the population size on the results.When the maximum evolution generation is 10, the influence of population size and atom number on reconstruction performance is shown in Figure5. The sampling rate in Figure 5(a) is 0.3, while the sampling rate in Figure 5(b) is 0.4. Figure 5(a) shows that the atom number has a greater influence on the reconstruction result than the population size. As the number of atoms increases, PSNR can grow rapidly and stabilize after the number of atoms reaches 50. This conclusion is consistent with that of Figure 4(b).Figure 5
Influence of atom number and population size on reconstructed PSNR, when the evolution generation is 10: (a) sampling rate is 0.3; (b) sampling rate is 0.4.
(a)(b)As shown in Figure5(b), we can see that when the population increases to 20, the reconstructed PSNR reaches the maximum value and is very stable, which is consistent with the conclusion of Figure 3. This also fully shows that it is reasonable for us to set the number of populations to 20.The experimental results of the other three groups of hyperspectral images are similar to Cuprite1; considering the reconstruction accuracy, the selected population size in HN_OMP isP=20, and the maximum evolution generation is G=10.


#### 5.1.2. Parameters Selection of Atom Number

From the above experimental results, we have found that after the number of atoms reaches 50, the algorithm performance may have stabilized. Now we employ OMP algorithm and HN_OMP algorithm to compress and reconstruct the hyperspectral image to accurately determine the number of atoms. The HN_OMP algorithm and the OMP algorithm are used to reconstruct the 40th band image of the hyperspectral datasets, and the algorithm is set to terminate when the atom number reaches 50. In the experiments, the variation range of the atom number is 1-50 with the interval of 1, the sampling rate range is 0.1 to 0.5, and the interval is 0.1.The changes of the reconstructed PSNR of Cuprite2 obtained by the two algorithms are shown in Figure6. Whether it is the OMP algorithm or the HN_OMP algorithm, except for the sampling rate of 0.1, the PSNR will gradually increase as the atom number increases and finally towards stability. The reason is, the higher the sampling rate, the more information is available in the reconstruction process and the reconstruction accuracy is higher.Figure 6
Reconstructed PSNR of Cuprite2 vs atom number at different sampling rates: (a) OMP; (b) HN_OMP.
(a)(b)When the sampling rate is 0.1, after the atom number increases to 27, the performance of the two algorithm drops sharply, and accurate reconstruction is no longer possible. This is because the accurate reconstruction using compressed sensing theory needs to meet the condition that the number of measurements is greater than or equal to the sparsity of the signal. In the two reconstruction algorithms, the signal is sparsely represented by several atoms in the redundant dictionary, so the sparsity is equal to the number of atoms. Therefore, for accurate reconstruction, the number of measurements should be greater than or equal to the number of atoms. When the signal length is 256 and the sampling rate is 0.1, the number of measurements is 26. Therefore, when the number of atoms increases to 27, the conditions for accurate reconstruction are no longer met, and the algorithm performance will drop sharply.The PSNR comparison with two algorithms is shown in Figure7; the sampling rates are 0.2 in (a), 0.3 in (b), 0.4 in (c), and 0.5 in (d), respectively. We found that the reconstruction PSNR of HN_OMP algorithm rises less slowly than the OMP algorithm. This is because the HN_OMP algorithm has a certain degree of randomness when looking for the optimal atom, and it cannot find the optimal atoms step by step like the OMP algorithm. Fortunately, as the number of atoms increases, the HN_OMP algorithm adopts an elite strategy to ensure that the newly added optimal atom can restore the original signal to the greatest extent, thereby ultimately ensuring the accuracy of reconstruction.Figure 7
Reconstructed PSNR of Cuprite2 vs atom number at the same sampling rate: (a) sampling rate is 0.2; (b) sampling rate is 0.3; (c) sampling rate is 0.4; (d) sampling rate is 0.5.
(a)(b)(c)(d)Figure8 shows the changes of the reconstructed SSIM of Cuprite2 obtained by the two algorithms. Whether it is the OMP algorithm or the HN_OMP algorithm, except for the sampling rate of 0.1, the SSIM will gradually increase as the atom number increases and finally towards stability. Specially, the SSIM of the HN_OMP algorithm continues to rise and tends to stabilize, while the OMP algorithm has a slight downward trend when the number of atoms is relatively large. Similar to the result of PSNR, when the sampling rate is 0.1, the result is not good. Under other sampling rate conditions, the result of HN_OMP algorithm is better than OMP.Figure 8
Reconstructed SSIM of Cuprite2 vs atom number at different sampling rates: (a) OMP; (b) HN_OMP.
(a)(b)Not surprisingly, we also noticed that the conclusions obtained here are also consistent with section5.1.1. When the number of atoms reaches 50, the performance of the algorithm will stabilize. Comprehensive analysis of the above PSNR and SSIM results, in the following experiments, the optimal atom number of the two algorithms is set to K=25 when the sampling rate is 0.1, and other sampling rates are set to K=50 to compare reconstruction performance.


### 5.2. Reconstruction Performance Comparison

Four exemplary reconstruction techniques are compared in order to better understand the possibilities of the proposed HN_OMP. The first is the minimization total variation (min-TV) algorithm with quadratic constraints [60]. This approach uses a random Gaussian matrix as the spatial measurement matrix. It performs the reconstruction by maximizing the total variation minimization in the reconstruction process, which makes use of the image’s gradient sparseness. The second is BCS_SPL, the block compressed sensing measures the image in blocks at the sampling end and uses SPL to efficiently and accurately restore the original image [36]. The third one is GPSR [14]; this algorithm uses random Gaussian matrix to sample each band image and uses gradient projection to solve the convex reconstruct problem. The fourth one is OMP [51]; the hyperspectral images are sampled using random Gaussian matrix to sample each band image block. In the reconstruction process, each image block is reconstructed by looking for optimal atoms to linearly represent the original signal, and the specific steps are described in Section 2.2.Similar with the fourth algorithm, the proposed algorithm HN_OMP also uses random Gaussian matrix to sample each band image block. But its reconstruction process is different from all the above algorithms. In the reconstruction process, utilizing the concept of evolutionary algorithm, algorithm HN_OMP transformed the reconstruction process into the process of population evolution. Through analyzing the process of OMP algorithm, a multiobjective reconstruction model is established, and the optimal atoms are found through population evolution, thereby completing the reconstruction of the original signal. Furthermore, the NSGA-III algorithm is used to complete the matching process of the optimal atom. For the residual update process, the Hermitian inversion lemma is introduced to speed up the execution of the algorithm.The sample rates in the simulations range from 0.1 to 0.5 with a 0.1 interval. With reference, the parameters for algorithm min-TV and BCS SPL are set to default. The parameters of OMP is set asK=50 at sampling rate 0.2~0.5, and K=25 at sampling rate 0.1. The parameters of HN_OMP algorithm are set as P=20, G=10, and the atom number is set as K=50 at sampling rate 0.2~0.5, and K=25 at sampling rate 0.1.The band-based PSNR of four datasets is shown in Figure9. It can be seen from the figure that as the sampling rate increases, the amount of measurement data becomes larger and larger, and the more information available for reconstruction, the higher the reconstruction accuracy.Figure 9
Band-based PSNR for hyperspectral datasets at different sampling rates: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)Figures10(a) and 10(b) show the average reconstructed PSNR of Cuprite1 and Cuprite2, respectively. Because the two photos were captured by the same sensor and their imaging objects were the Cuprite region, the simulation results are equivalent. When the sampling rate is greater than 0.2 for the Cuprite1 image, the proposed HN_OMP has the highest PSNR, followed by GPSR. The algorithm BCS_SPL which explores smoothed projected Landweber in the reconstruction could gain almost the same performance as GPSR. The OMP algorithm using optimal atoms to recover the signal has better PSNR than that of min-TV algorithm. Because the reconstruction procedure solely uses the spatial correlation of hyperspectral pictures, the PSNR of min-TV is the lowest. The results of Cuprite2 could get the same conclusion with the results of Cuprite1 images.Figure 10
Average reconstructed PSNR for hyperspectral datasets at different sampling rates: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)Figures10(c) and 10(d) illustrate the results of Indian Pines and Pavia University, respectively. The performance of min-TV algorithm, OMP algorithm, and GPSR algorithm are almost the same. The proposed HN_OMP algorithm has the same performance as the BCS_SPL algorithm. The difference is that the GPSR algorithm cannot get superior results for these two datasets. The reason might lie in that, compared with the other two datasets, the spatial distribution of these two datasets is quite different, and the image gradient can be used to better represent the image characteristics, to obtain better reconstruction results.When the sample rate is reduced to 0.1, the proposed HN_OMP method can successfully rebuild Indian Pines and Pavia University images, whereas it cannot get good results for Cuprite images. Overall, compared with min-TV algorithm, OMP algorithm or GPSR algorithm, the proposed algorithm HN_OMP has a certain degree of improvement in terms of PSNR, which proves the efficacy of the proposed method.Table2 shows the SSIM comparison results for four hyperspectral datasets. For Cuprite1 and Cuprite2 images, the SSIM of OMP algorithm is the lowest, which is not different from the PSNR comparison in Figures 10(a) and 10(b). The SSIM of the proposed HN_OMP algorithm has a weak advantage over that of min-TV and GPSR. The SSIM of the algorithm HN_OMP is lower than BCS_SPL at lower sampling rates; fortunately, the SSIM of HN_OMP finally approaches BCS_SPL. For Indian Pines and Pavia University, the situation will have a little change. At this moment, the SSIM of min-TV would be the lowest, while the SSIM of OMP and GPSR are very close to each other. Unfortunately, at this time, the proposed HN_OMP does not have an advantage. In general, the HN_OMP algorithm can restore the structure of the original image to a certain extent, which is a compelling proof of the proposed algorithm’s usefulness.Table 2
Average SSIM for hyperspectral datasets under different sampling rates.
AlgorithmSR=0.1SR=0.2SR=0.3SR=0.4SR=0.5Cuprite1min-TV0.79270.84270.87050.89810.9223OMP0.37280.54170.69500.79340.8535GPSR0.62580.82840.86710.90020.9267BCS_SPL0.88650.90970.92420.93400.9424HN_OMP0.62870.81450.88380.92700.9508Cuprite2min-TV0.63020.70210.75400.80270.8482OMP0.31580.47680.57580.68270.7528GPSR0.43480.67140.74230.79580.8488BCS_SPL0.81890.86180.88890.90920.9246HN_OMP0.61850.79200.82690.87760.9278Indian Pinesmin-TV0.13850.22790.33060.42860.5271OMP0.14990.25380.34210.44780.5421GPSR0.11750.30040.39670.50380.5981BCS_SPL0.57550.67480.74560.79100.8293HN_OMP0.37700.48480.61170.65450.7136Pavia Universitymin-TV0.11990.21230.30720.40790.5157OMP0.25830.37750.47380.57260.6671GPSR0.18460.39690.50530.61570.7166BCS_SPL0.67530.76470.81960.85640.8854HN_OMP0.44850.58560.75330.81300.8863The following comparison with several other state-of-the-art algorithms is presented to further showcase the algorithm’s performance. Table3 lists Pavia University’s results in comparison to the algorithms CPPCA-KMEANS [41] and CPPCA-MH [42] employing vector-based SNR and SAD. To achieve the reconstruction, the two methods, both based on CPPCA, combined the advantages of classification and multihypothesis prediction. Seen from the table, the SNR of HN_OMP are better than that of the other two algorithms, when the sampling rate is 0.4 and 0.5. The SAD result is not good as the SNR; the reason might be in the reconstruction process; we did not consider the correlation between spectral vectors.Table 3
The average SNR (dB) and SAD (°) comparisons of Pavia University.
AlgorithmSR=0.2SR=0.3SR=0.4SR=0.5Average vector-based SNRCPPCA-KMEANS [41]17.5320.2822.1923.86CPPCA_MH [42]19.1623.1325.4427.95HN_OMP15.6720.7325.8728.03Average vector-based SADCPPCA-KMEANS [41]2.822.151.541.27CPPCA_MH [42]2.171.391.070.81HN_OMP5.293.161.741.36Secondly, the results of Indian Pines, comparing with the algorithms CPPCA-MH [42] and HNP [12] using vector-based and band-based SNR, are listed in Table 4. The idea of HNP is to combine the advantages of a fast l0 norm and a reliable l1 norm. Unfortunately, the performance of HN_OMP is worse than that of CPPCA-MH or HNP. The reason may lie in two aspects. One is that the reconstruction parameters are not ideal, resulting in unsatisfactory reconstruction outcomes. The other is that the optimization methods of the three algorithms differ, making it inappropriate to compare them.Table 4
The average SNR (dB) comparisons of Indian Pines.
AlgorithmSR=0.2SR=0.3SR=0.4SR=0.5Average vector-based SNRCPPCA_MH [42]36.2138.7841.0542.93HN_OMP11.1513.7215.2818.03Average band-based SNRHNP [12]35.940.042.844.9HN_OMP13.314.716.318.0Thirdly, the results of Pavia University and Indian Pines, comparing with SSAHTV [15], SSTV [16], and HSSTV [17] using PSNR and SSIM, are shown in Table 5. The PSNR and SSIM values of these algorithms are extracted from the original literature. Concerning Indian Pines, the proposed HN_OMP has less advantage; it could only be compared to SSTV. Regarding Pavia University, the PSNR of the proposed HN_OMP is higher than that of the other three algorithms in most instances. Meanwhile, the SSIM of the proposed HN_OMP can keep the difference with other algorithms within 0.1.Table 5
The average PSNR (dB) and SSIM comparisons of Indian Pines and Pavia University.
AlgorithmPSNRSSIMSR=0.2SR=0.4SR=0.2SR=0.4Indian PinesSSAHTV [15]29.9930.550.76580.7777SSTV [16]25.1127.650.34880.5066HSSTV [17]30.7131.360.75890.7806HN_OMP24.3629.020.48480.6545Pavia UniversitySSAHTV [15]24.9726.530.62420.6956SSTV [16]27.3530.020.66230.7830HSSTV [17]31.1332.880.85570.8950HN_OMP28.7434.130.58560.8130In conclusion, numerous comparisons between the proposed algorithm HN_OMP and some state-of-the-art algorithms have been demonstrated in order to assess reconstruction performance. Especially in comparison with the OMP algorithm, the results and analysis fully indicate that it is feasible to use the evolutionary algorithm to simulate the optimal atom search process. The optimal atom searched by NSGA-III could approximate the original signal, which can ensure the accuracy of image reconstruction and fully demonstrate the reliability of the algorithm. Although the results of some images are slightly worse than those of other algorithms, it can still show that the use of evolutionary algorithms can be used as a way of image reconstruction.


### 5.3. Runtime Comparison

The runtime of five reconstruction algorithms is shown in Figure11. It can be found from the histogram that the algorithm OMP needs to traverse all the atoms in the dictionary, thus leading the lowest running speed. The computational complexity of min-TV and BCS_SPL are equivalent, while the time consuming by GPSR could reduce to a half. Fortunately, our proposed algorithm HN_OMP, introducing the evolution concept, could greatly decrease computational complexity of reconstruction process. Compared with OMP, the proposed HN_OMP can increase the calculation speed by an order of magnitude. Combing with previous performance comparison, although the PSNR or SSIM of the proposed HN_OMP is in a bad situation, the computing efficiency will become the advantage of the algorithm.Figure 11
Runtime for hyperspectral datasets using different reconstruction algorithms: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)


### 5.4. Reconstruction Image Comparison

As compared to PSNR, a viewer’s preference for pleasant visual quality is more relevant. Figures12 and 13 show the reconstructed pictures of Cuprite2 and Pavia University with a sampling rate of 0.5. The images are the 50th band of the dataset, and the corresponding reconstructed PSNR of each image is also given. Seen from these reconstruction images, these five reconstruction algorithms can restore image features more accurately. The reconstructed image can well describe the detailed features of the original image, which fully demonstrates that the HN_OMP algorithm can find the optimal atoms by using the evolution process of the chromosome, could get visual effect with a cleaner image, and fully demonstrate the effectiveness of the algorithm.Figure 12
Comparison between five reconstructed Cuprite2 images and original image: (a) original image; (b) min-TV,PSNR=30.25 dB; (c) GPSR, PSNR=34.29 dB; (d) BCS_SPL, PSNR=35.84 dB; (e) OMP, PSNR=31.75 dB; (f) HN_OMP, PSNR=36.27dB.
(a)(b)(c)(d)(e)(f)Figure 13
Comparison between five reconstructed Pavia University images and original image: (a) original image; (b) min-TV,PSNR=27.31 dB; (c) GPSR,PSNR=28.26 dB; (d) BCS_SPL, PSNR=35.00dB; (e) OMP, PSNR=26.58dB; (f) HN_OMP, PSNR=36.49dB.
(a)(b)(c)(d)(e)(f)Although the algorithm has a certain degree of blocking effect, the visual effect is not as clear as the BCS_SPL algorithm, but the time-consuming of proposed HN_OMP is the lowest among these algorithms. Overall, we can conclude that the algorithm HN_OMP can improve the reconstruction efficiency through evolutionary algorithm NSGA-III and Hermitian inversion, while maintaining the reconstruction accuracy.


## 5.1. Parameter Selection for Reconstruction Algorithm

The parameter selection would undoubtedly have a major impact on the recon-construction performance. The results of the proposed algorithm HN_OMP on four hyperspectral pictures with varying evolution parameters are reported in this section.


### 5.1.1. Parameters Selection of Population Size and Evolution Generation

Firstly, the algorithm HN_OMP is used to compress reconstruct the 40th band image of the 4 groups of hyperspectral data, and the influence of the maximum evolution generation, population size, and atom number on the performance of the algorithm is analyzed. The variation range of population sizeP is 10~50 at the interval of 5, the variation range of the maximum evolution generation G is 10~50 at the interval of 5, and the variation range of the atom number K is 10~100 at the interval of 10. The sampling rates are set as 0.1 to 0.5 at the interval of 0.1.Figure3 shows the variation of the average PSNR of the reconstructed image of Cuprite1 with the population size and evolution generation, when the atom number is 40. The sampling rate in Figure 3(a) is 0.3, while the sampling rate in Figure 3(b) is 0.4. From the two subfigures, we can find that under the same evolution generation, after the population size increases to 20, the PSNR tends to be stable no matter what the sampling rate is. Under the same population size, the PSNR has a small oscillation with the increase of evolution generation. It demonstrates that the reconstructed PSNR is more sensitive to population size than that to evolution generation. Without loss of generality, we set the population size as 20 to discuss the influence of other parameters on the reconstruction results.Figure 3
Influence of evolution generation and population size on reconstructed PSNR, when the atom number is 40: (a) sampling rate is 0.3; (b) sampling rate is 0.4.
(a)(b)When the population size is 20, the influence of the evolution generation and atom number on the reconstructed PSNR of dataset Cuprite1 is shown in Figure4. The sampling rate in Figure 4(a) is 0.2, while the sampling rate in Figure 4(b) is 0.3. From the two subfigures, we can clearly see that the influence of the atom number on the reconstruction result is far greater than the influence of the evolution generation on the reconstruction result. As the number of atoms increases, the reconstructed PSNR steadily increases.Figure 4
Influence of evolution generation and atom number on reconstructed PSNR, when the population size is 20: (a) sampling rate is 0.2; (b) sampling rate is 0.3.
(a)(b)However, when the number of atoms further increases to more than 50, the reconstruction performance no longer increases significantly, but gradually stabilizes. The reason for this situation may be that particularly large number of atoms will affect the selecting of elite atoms. Therefore, the algorithm performance cannot be further improved. Comparing Figure3(a) and Figure 4(a), we find that the atom number and the population size have a relatively large impact on the results, and the evolution generation has a small impact. Without loss of generality, we set the evolution generation to 10 to compare the effects of the atom number and the population size on the results.When the maximum evolution generation is 10, the influence of population size and atom number on reconstruction performance is shown in Figure5. The sampling rate in Figure 5(a) is 0.3, while the sampling rate in Figure 5(b) is 0.4. Figure 5(a) shows that the atom number has a greater influence on the reconstruction result than the population size. As the number of atoms increases, PSNR can grow rapidly and stabilize after the number of atoms reaches 50. This conclusion is consistent with that of Figure 4(b).Figure 5
Influence of atom number and population size on reconstructed PSNR, when the evolution generation is 10: (a) sampling rate is 0.3; (b) sampling rate is 0.4.
(a)(b)As shown in Figure5(b), we can see that when the population increases to 20, the reconstructed PSNR reaches the maximum value and is very stable, which is consistent with the conclusion of Figure 3. This also fully shows that it is reasonable for us to set the number of populations to 20.The experimental results of the other three groups of hyperspectral images are similar to Cuprite1; considering the reconstruction accuracy, the selected population size in HN_OMP isP=20, and the maximum evolution generation is G=10.


### 5.1.2. Parameters Selection of Atom Number

From the above experimental results, we have found that after the number of atoms reaches 50, the algorithm performance may have stabilized. Now we employ OMP algorithm and HN_OMP algorithm to compress and reconstruct the hyperspectral image to accurately determine the number of atoms. The HN_OMP algorithm and the OMP algorithm are used to reconstruct the 40th band image of the hyperspectral datasets, and the algorithm is set to terminate when the atom number reaches 50. In the experiments, the variation range of the atom number is 1-50 with the interval of 1, the sampling rate range is 0.1 to 0.5, and the interval is 0.1.The changes of the reconstructed PSNR of Cuprite2 obtained by the two algorithms are shown in Figure6. Whether it is the OMP algorithm or the HN_OMP algorithm, except for the sampling rate of 0.1, the PSNR will gradually increase as the atom number increases and finally towards stability. The reason is, the higher the sampling rate, the more information is available in the reconstruction process and the reconstruction accuracy is higher.Figure 6
Reconstructed PSNR of Cuprite2 vs atom number at different sampling rates: (a) OMP; (b) HN_OMP.
(a)(b)When the sampling rate is 0.1, after the atom number increases to 27, the performance of the two algorithm drops sharply, and accurate reconstruction is no longer possible. This is because the accurate reconstruction using compressed sensing theory needs to meet the condition that the number of measurements is greater than or equal to the sparsity of the signal. In the two reconstruction algorithms, the signal is sparsely represented by several atoms in the redundant dictionary, so the sparsity is equal to the number of atoms. Therefore, for accurate reconstruction, the number of measurements should be greater than or equal to the number of atoms. When the signal length is 256 and the sampling rate is 0.1, the number of measurements is 26. Therefore, when the number of atoms increases to 27, the conditions for accurate reconstruction are no longer met, and the algorithm performance will drop sharply.The PSNR comparison with two algorithms is shown in Figure7; the sampling rates are 0.2 in (a), 0.3 in (b), 0.4 in (c), and 0.5 in (d), respectively. We found that the reconstruction PSNR of HN_OMP algorithm rises less slowly than the OMP algorithm. This is because the HN_OMP algorithm has a certain degree of randomness when looking for the optimal atom, and it cannot find the optimal atoms step by step like the OMP algorithm. Fortunately, as the number of atoms increases, the HN_OMP algorithm adopts an elite strategy to ensure that the newly added optimal atom can restore the original signal to the greatest extent, thereby ultimately ensuring the accuracy of reconstruction.Figure 7
Reconstructed PSNR of Cuprite2 vs atom number at the same sampling rate: (a) sampling rate is 0.2; (b) sampling rate is 0.3; (c) sampling rate is 0.4; (d) sampling rate is 0.5.
(a)(b)(c)(d)Figure8 shows the changes of the reconstructed SSIM of Cuprite2 obtained by the two algorithms. Whether it is the OMP algorithm or the HN_OMP algorithm, except for the sampling rate of 0.1, the SSIM will gradually increase as the atom number increases and finally towards stability. Specially, the SSIM of the HN_OMP algorithm continues to rise and tends to stabilize, while the OMP algorithm has a slight downward trend when the number of atoms is relatively large. Similar to the result of PSNR, when the sampling rate is 0.1, the result is not good. Under other sampling rate conditions, the result of HN_OMP algorithm is better than OMP.Figure 8
Reconstructed SSIM of Cuprite2 vs atom number at different sampling rates: (a) OMP; (b) HN_OMP.
(a)(b)Not surprisingly, we also noticed that the conclusions obtained here are also consistent with section5.1.1. When the number of atoms reaches 50, the performance of the algorithm will stabilize. Comprehensive analysis of the above PSNR and SSIM results, in the following experiments, the optimal atom number of the two algorithms is set to K=25 when the sampling rate is 0.1, and other sampling rates are set to K=50 to compare reconstruction performance.


## 5.1.1. Parameters Selection of Population Size and Evolution Generation

Firstly, the algorithm HN_OMP is used to compress reconstruct the 40th band image of the 4 groups of hyperspectral data, and the influence of the maximum evolution generation, population size, and atom number on the performance of the algorithm is analyzed. The variation range of population sizeP is 10~50 at the interval of 5, the variation range of the maximum evolution generation G is 10~50 at the interval of 5, and the variation range of the atom number K is 10~100 at the interval of 10. The sampling rates are set as 0.1 to 0.5 at the interval of 0.1.Figure3 shows the variation of the average PSNR of the reconstructed image of Cuprite1 with the population size and evolution generation, when the atom number is 40. The sampling rate in Figure 3(a) is 0.3, while the sampling rate in Figure 3(b) is 0.4. From the two subfigures, we can find that under the same evolution generation, after the population size increases to 20, the PSNR tends to be stable no matter what the sampling rate is. Under the same population size, the PSNR has a small oscillation with the increase of evolution generation. It demonstrates that the reconstructed PSNR is more sensitive to population size than that to evolution generation. Without loss of generality, we set the population size as 20 to discuss the influence of other parameters on the reconstruction results.Figure 3
Influence of evolution generation and population size on reconstructed PSNR, when the atom number is 40: (a) sampling rate is 0.3; (b) sampling rate is 0.4.
(a)(b)When the population size is 20, the influence of the evolution generation and atom number on the reconstructed PSNR of dataset Cuprite1 is shown in Figure4. The sampling rate in Figure 4(a) is 0.2, while the sampling rate in Figure 4(b) is 0.3. From the two subfigures, we can clearly see that the influence of the atom number on the reconstruction result is far greater than the influence of the evolution generation on the reconstruction result. As the number of atoms increases, the reconstructed PSNR steadily increases.Figure 4
Influence of evolution generation and atom number on reconstructed PSNR, when the population size is 20: (a) sampling rate is 0.2; (b) sampling rate is 0.3.
(a)(b)However, when the number of atoms further increases to more than 50, the reconstruction performance no longer increases significantly, but gradually stabilizes. The reason for this situation may be that particularly large number of atoms will affect the selecting of elite atoms. Therefore, the algorithm performance cannot be further improved. Comparing Figure3(a) and Figure 4(a), we find that the atom number and the population size have a relatively large impact on the results, and the evolution generation has a small impact. Without loss of generality, we set the evolution generation to 10 to compare the effects of the atom number and the population size on the results.When the maximum evolution generation is 10, the influence of population size and atom number on reconstruction performance is shown in Figure5. The sampling rate in Figure 5(a) is 0.3, while the sampling rate in Figure 5(b) is 0.4. Figure 5(a) shows that the atom number has a greater influence on the reconstruction result than the population size. As the number of atoms increases, PSNR can grow rapidly and stabilize after the number of atoms reaches 50. This conclusion is consistent with that of Figure 4(b).Figure 5
Influence of atom number and population size on reconstructed PSNR, when the evolution generation is 10: (a) sampling rate is 0.3; (b) sampling rate is 0.4.
(a)(b)As shown in Figure5(b), we can see that when the population increases to 20, the reconstructed PSNR reaches the maximum value and is very stable, which is consistent with the conclusion of Figure 3. This also fully shows that it is reasonable for us to set the number of populations to 20.The experimental results of the other three groups of hyperspectral images are similar to Cuprite1; considering the reconstruction accuracy, the selected population size in HN_OMP isP=20, and the maximum evolution generation is G=10.


## 5.1.2. Parameters Selection of Atom Number

From the above experimental results, we have found that after the number of atoms reaches 50, the algorithm performance may have stabilized. Now we employ OMP algorithm and HN_OMP algorithm to compress and reconstruct the hyperspectral image to accurately determine the number of atoms. The HN_OMP algorithm and the OMP algorithm are used to reconstruct the 40th band image of the hyperspectral datasets, and the algorithm is set to terminate when the atom number reaches 50. In the experiments, the variation range of the atom number is 1-50 with the interval of 1, the sampling rate range is 0.1 to 0.5, and the interval is 0.1.The changes of the reconstructed PSNR of Cuprite2 obtained by the two algorithms are shown in Figure6. Whether it is the OMP algorithm or the HN_OMP algorithm, except for the sampling rate of 0.1, the PSNR will gradually increase as the atom number increases and finally towards stability. The reason is, the higher the sampling rate, the more information is available in the reconstruction process and the reconstruction accuracy is higher.Figure 6
Reconstructed PSNR of Cuprite2 vs atom number at different sampling rates: (a) OMP; (b) HN_OMP.
(a)(b)When the sampling rate is 0.1, after the atom number increases to 27, the performance of the two algorithm drops sharply, and accurate reconstruction is no longer possible. This is because the accurate reconstruction using compressed sensing theory needs to meet the condition that the number of measurements is greater than or equal to the sparsity of the signal. In the two reconstruction algorithms, the signal is sparsely represented by several atoms in the redundant dictionary, so the sparsity is equal to the number of atoms. Therefore, for accurate reconstruction, the number of measurements should be greater than or equal to the number of atoms. When the signal length is 256 and the sampling rate is 0.1, the number of measurements is 26. Therefore, when the number of atoms increases to 27, the conditions for accurate reconstruction are no longer met, and the algorithm performance will drop sharply.The PSNR comparison with two algorithms is shown in Figure7; the sampling rates are 0.2 in (a), 0.3 in (b), 0.4 in (c), and 0.5 in (d), respectively. We found that the reconstruction PSNR of HN_OMP algorithm rises less slowly than the OMP algorithm. This is because the HN_OMP algorithm has a certain degree of randomness when looking for the optimal atom, and it cannot find the optimal atoms step by step like the OMP algorithm. Fortunately, as the number of atoms increases, the HN_OMP algorithm adopts an elite strategy to ensure that the newly added optimal atom can restore the original signal to the greatest extent, thereby ultimately ensuring the accuracy of reconstruction.Figure 7
Reconstructed PSNR of Cuprite2 vs atom number at the same sampling rate: (a) sampling rate is 0.2; (b) sampling rate is 0.3; (c) sampling rate is 0.4; (d) sampling rate is 0.5.
(a)(b)(c)(d)Figure8 shows the changes of the reconstructed SSIM of Cuprite2 obtained by the two algorithms. Whether it is the OMP algorithm or the HN_OMP algorithm, except for the sampling rate of 0.1, the SSIM will gradually increase as the atom number increases and finally towards stability. Specially, the SSIM of the HN_OMP algorithm continues to rise and tends to stabilize, while the OMP algorithm has a slight downward trend when the number of atoms is relatively large. Similar to the result of PSNR, when the sampling rate is 0.1, the result is not good. Under other sampling rate conditions, the result of HN_OMP algorithm is better than OMP.Figure 8
Reconstructed SSIM of Cuprite2 vs atom number at different sampling rates: (a) OMP; (b) HN_OMP.
(a)(b)Not surprisingly, we also noticed that the conclusions obtained here are also consistent with section5.1.1. When the number of atoms reaches 50, the performance of the algorithm will stabilize. Comprehensive analysis of the above PSNR and SSIM results, in the following experiments, the optimal atom number of the two algorithms is set to K=25 when the sampling rate is 0.1, and other sampling rates are set to K=50 to compare reconstruction performance.


## 5.2. Reconstruction Performance Comparison

Four exemplary reconstruction techniques are compared in order to better understand the possibilities of the proposed HN_OMP. The first is the minimization total variation (min-TV) algorithm with quadratic constraints [60]. This approach uses a random Gaussian matrix as the spatial measurement matrix. It performs the reconstruction by maximizing the total variation minimization in the reconstruction process, which makes use of the image’s gradient sparseness. The second is BCS_SPL, the block compressed sensing measures the image in blocks at the sampling end and uses SPL to efficiently and accurately restore the original image [36]. The third one is GPSR [14]; this algorithm uses random Gaussian matrix to sample each band image and uses gradient projection to solve the convex reconstruct problem. The fourth one is OMP [51]; the hyperspectral images are sampled using random Gaussian matrix to sample each band image block. In the reconstruction process, each image block is reconstructed by looking for optimal atoms to linearly represent the original signal, and the specific steps are described in Section 2.2.Similar with the fourth algorithm, the proposed algorithm HN_OMP also uses random Gaussian matrix to sample each band image block. But its reconstruction process is different from all the above algorithms. In the reconstruction process, utilizing the concept of evolutionary algorithm, algorithm HN_OMP transformed the reconstruction process into the process of population evolution. Through analyzing the process of OMP algorithm, a multiobjective reconstruction model is established, and the optimal atoms are found through population evolution, thereby completing the reconstruction of the original signal. Furthermore, the NSGA-III algorithm is used to complete the matching process of the optimal atom. For the residual update process, the Hermitian inversion lemma is introduced to speed up the execution of the algorithm.The sample rates in the simulations range from 0.1 to 0.5 with a 0.1 interval. With reference, the parameters for algorithm min-TV and BCS SPL are set to default. The parameters of OMP is set asK=50 at sampling rate 0.2~0.5, and K=25 at sampling rate 0.1. The parameters of HN_OMP algorithm are set as P=20, G=10, and the atom number is set as K=50 at sampling rate 0.2~0.5, and K=25 at sampling rate 0.1.The band-based PSNR of four datasets is shown in Figure9. It can be seen from the figure that as the sampling rate increases, the amount of measurement data becomes larger and larger, and the more information available for reconstruction, the higher the reconstruction accuracy.Figure 9
Band-based PSNR for hyperspectral datasets at different sampling rates: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)Figures10(a) and 10(b) show the average reconstructed PSNR of Cuprite1 and Cuprite2, respectively. Because the two photos were captured by the same sensor and their imaging objects were the Cuprite region, the simulation results are equivalent. When the sampling rate is greater than 0.2 for the Cuprite1 image, the proposed HN_OMP has the highest PSNR, followed by GPSR. The algorithm BCS_SPL which explores smoothed projected Landweber in the reconstruction could gain almost the same performance as GPSR. The OMP algorithm using optimal atoms to recover the signal has better PSNR than that of min-TV algorithm. Because the reconstruction procedure solely uses the spatial correlation of hyperspectral pictures, the PSNR of min-TV is the lowest. The results of Cuprite2 could get the same conclusion with the results of Cuprite1 images.Figure 10
Average reconstructed PSNR for hyperspectral datasets at different sampling rates: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)Figures10(c) and 10(d) illustrate the results of Indian Pines and Pavia University, respectively. The performance of min-TV algorithm, OMP algorithm, and GPSR algorithm are almost the same. The proposed HN_OMP algorithm has the same performance as the BCS_SPL algorithm. The difference is that the GPSR algorithm cannot get superior results for these two datasets. The reason might lie in that, compared with the other two datasets, the spatial distribution of these two datasets is quite different, and the image gradient can be used to better represent the image characteristics, to obtain better reconstruction results.When the sample rate is reduced to 0.1, the proposed HN_OMP method can successfully rebuild Indian Pines and Pavia University images, whereas it cannot get good results for Cuprite images. Overall, compared with min-TV algorithm, OMP algorithm or GPSR algorithm, the proposed algorithm HN_OMP has a certain degree of improvement in terms of PSNR, which proves the efficacy of the proposed method.Table2 shows the SSIM comparison results for four hyperspectral datasets. For Cuprite1 and Cuprite2 images, the SSIM of OMP algorithm is the lowest, which is not different from the PSNR comparison in Figures 10(a) and 10(b). The SSIM of the proposed HN_OMP algorithm has a weak advantage over that of min-TV and GPSR. The SSIM of the algorithm HN_OMP is lower than BCS_SPL at lower sampling rates; fortunately, the SSIM of HN_OMP finally approaches BCS_SPL. For Indian Pines and Pavia University, the situation will have a little change. At this moment, the SSIM of min-TV would be the lowest, while the SSIM of OMP and GPSR are very close to each other. Unfortunately, at this time, the proposed HN_OMP does not have an advantage. In general, the HN_OMP algorithm can restore the structure of the original image to a certain extent, which is a compelling proof of the proposed algorithm’s usefulness.Table 2
Average SSIM for hyperspectral datasets under different sampling rates.
AlgorithmSR=0.1SR=0.2SR=0.3SR=0.4SR=0.5Cuprite1min-TV0.79270.84270.87050.89810.9223OMP0.37280.54170.69500.79340.8535GPSR0.62580.82840.86710.90020.9267BCS_SPL0.88650.90970.92420.93400.9424HN_OMP0.62870.81450.88380.92700.9508Cuprite2min-TV0.63020.70210.75400.80270.8482OMP0.31580.47680.57580.68270.7528GPSR0.43480.67140.74230.79580.8488BCS_SPL0.81890.86180.88890.90920.9246HN_OMP0.61850.79200.82690.87760.9278Indian Pinesmin-TV0.13850.22790.33060.42860.5271OMP0.14990.25380.34210.44780.5421GPSR0.11750.30040.39670.50380.5981BCS_SPL0.57550.67480.74560.79100.8293HN_OMP0.37700.48480.61170.65450.7136Pavia Universitymin-TV0.11990.21230.30720.40790.5157OMP0.25830.37750.47380.57260.6671GPSR0.18460.39690.50530.61570.7166BCS_SPL0.67530.76470.81960.85640.8854HN_OMP0.44850.58560.75330.81300.8863The following comparison with several other state-of-the-art algorithms is presented to further showcase the algorithm’s performance. Table3 lists Pavia University’s results in comparison to the algorithms CPPCA-KMEANS [41] and CPPCA-MH [42] employing vector-based SNR and SAD. To achieve the reconstruction, the two methods, both based on CPPCA, combined the advantages of classification and multihypothesis prediction. Seen from the table, the SNR of HN_OMP are better than that of the other two algorithms, when the sampling rate is 0.4 and 0.5. The SAD result is not good as the SNR; the reason might be in the reconstruction process; we did not consider the correlation between spectral vectors.Table 3
The average SNR (dB) and SAD (°) comparisons of Pavia University.
AlgorithmSR=0.2SR=0.3SR=0.4SR=0.5Average vector-based SNRCPPCA-KMEANS [41]17.5320.2822.1923.86CPPCA_MH [42]19.1623.1325.4427.95HN_OMP15.6720.7325.8728.03Average vector-based SADCPPCA-KMEANS [41]2.822.151.541.27CPPCA_MH [42]2.171.391.070.81HN_OMP5.293.161.741.36Secondly, the results of Indian Pines, comparing with the algorithms CPPCA-MH [42] and HNP [12] using vector-based and band-based SNR, are listed in Table 4. The idea of HNP is to combine the advantages of a fast l0 norm and a reliable l1 norm. Unfortunately, the performance of HN_OMP is worse than that of CPPCA-MH or HNP. The reason may lie in two aspects. One is that the reconstruction parameters are not ideal, resulting in unsatisfactory reconstruction outcomes. The other is that the optimization methods of the three algorithms differ, making it inappropriate to compare them.Table 4
The average SNR (dB) comparisons of Indian Pines.
AlgorithmSR=0.2SR=0.3SR=0.4SR=0.5Average vector-based SNRCPPCA_MH [42]36.2138.7841.0542.93HN_OMP11.1513.7215.2818.03Average band-based SNRHNP [12]35.940.042.844.9HN_OMP13.314.716.318.0Thirdly, the results of Pavia University and Indian Pines, comparing with SSAHTV [15], SSTV [16], and HSSTV [17] using PSNR and SSIM, are shown in Table 5. The PSNR and SSIM values of these algorithms are extracted from the original literature. Concerning Indian Pines, the proposed HN_OMP has less advantage; it could only be compared to SSTV. Regarding Pavia University, the PSNR of the proposed HN_OMP is higher than that of the other three algorithms in most instances. Meanwhile, the SSIM of the proposed HN_OMP can keep the difference with other algorithms within 0.1.Table 5
The average PSNR (dB) and SSIM comparisons of Indian Pines and Pavia University.
AlgorithmPSNRSSIMSR=0.2SR=0.4SR=0.2SR=0.4Indian PinesSSAHTV [15]29.9930.550.76580.7777SSTV [16]25.1127.650.34880.5066HSSTV [17]30.7131.360.75890.7806HN_OMP24.3629.020.48480.6545Pavia UniversitySSAHTV [15]24.9726.530.62420.6956SSTV [16]27.3530.020.66230.7830HSSTV [17]31.1332.880.85570.8950HN_OMP28.7434.130.58560.8130In conclusion, numerous comparisons between the proposed algorithm HN_OMP and some state-of-the-art algorithms have been demonstrated in order to assess reconstruction performance. Especially in comparison with the OMP algorithm, the results and analysis fully indicate that it is feasible to use the evolutionary algorithm to simulate the optimal atom search process. The optimal atom searched by NSGA-III could approximate the original signal, which can ensure the accuracy of image reconstruction and fully demonstrate the reliability of the algorithm. Although the results of some images are slightly worse than those of other algorithms, it can still show that the use of evolutionary algorithms can be used as a way of image reconstruction.


## 5.3. Runtime Comparison

The runtime of five reconstruction algorithms is shown in Figure11. It can be found from the histogram that the algorithm OMP needs to traverse all the atoms in the dictionary, thus leading the lowest running speed. The computational complexity of min-TV and BCS_SPL are equivalent, while the time consuming by GPSR could reduce to a half. Fortunately, our proposed algorithm HN_OMP, introducing the evolution concept, could greatly decrease computational complexity of reconstruction process. Compared with OMP, the proposed HN_OMP can increase the calculation speed by an order of magnitude. Combing with previous performance comparison, although the PSNR or SSIM of the proposed HN_OMP is in a bad situation, the computing efficiency will become the advantage of the algorithm.Figure 11
Runtime for hyperspectral datasets using different reconstruction algorithms: (a) Cuprite1; (b) Cuprite2; (c) Indian Pines; (d) Pavia University.
(a)(b)(c)(d)


## 5.4. Reconstruction Image Comparison

As compared to PSNR, a viewer’s preference for pleasant visual quality is more relevant. Figures12 and 13 show the reconstructed pictures of Cuprite2 and Pavia University with a sampling rate of 0.5. The images are the 50th band of the dataset, and the corresponding reconstructed PSNR of each image is also given. Seen from these reconstruction images, these five reconstruction algorithms can restore image features more accurately. The reconstructed image can well describe the detailed features of the original image, which fully demonstrates that the HN_OMP algorithm can find the optimal atoms by using the evolution process of the chromosome, could get visual effect with a cleaner image, and fully demonstrate the effectiveness of the algorithm.Figure 12
Comparison between five reconstructed Cuprite2 images and original image: (a) original image; (b) min-TV,PSNR=30.25 dB; (c) GPSR, PSNR=34.29 dB; (d) BCS_SPL, PSNR=35.84 dB; (e) OMP, PSNR=31.75 dB; (f) HN_OMP, PSNR=36.27dB.
(a)(b)(c)(d)(e)(f)Figure 13
Comparison between five reconstructed Pavia University images and original image: (a) original image; (b) min-TV,PSNR=27.31 dB; (c) GPSR,PSNR=28.26 dB; (d) BCS_SPL, PSNR=35.00dB; (e) OMP, PSNR=26.58dB; (f) HN_OMP, PSNR=36.49dB.
(a)(b)(c)(d)(e)(f)Although the algorithm has a certain degree of blocking effect, the visual effect is not as clear as the BCS_SPL algorithm, but the time-consuming of proposed HN_OMP is the lowest among these algorithms. Overall, we can conclude that the algorithm HN_OMP can improve the reconstruction efficiency through evolutionary algorithm NSGA-III and Hermitian inversion, while maintaining the reconstruction accuracy.


## 6. Conclusion

Analyzing the characteristics of hyperspectral images, a hyperspectral image reconstruction algorithm based on HN_OMP is proposed in this study. The algorithm utilizes reference point nondominated sorting genetic method and Hermitian inversion to solve the multiobjective reconstruction model. The algorithm uses the chromosomes of the genetic algorithm to simulate the atoms of the redundant dictionary and simulates the crossover, selection, mutation, and other operators as the optimal atom search process. Furthermore, the algorithm uses the Hermitian inversion lemma to complete the update operation of the residual matrix in a recursive way. The proposed HN_OMP algorithm is used to reconstruct four groups of hyperspectral images, including the influence of the population size, evolution generation, and atom number on the reconstruction performance is discussed, and reasonable parameters are set to ensure the reliability of the algorithm. The proposed algorithm HN_OMP provides benefits over advanced reconstruction algorithms in terms of high reconstruction accuracy and low computing complexity, as demonstrated by results on real hyperspectral data. Since the proposed algorithm’s reconstruction accuracy is lower than that of some classic reconstruction algorithms, future research will focus on the possibility of local convergence and randomness concerns in evolutionary processes to conduct high-accuracy reconstruction.



---

*Source: 8455150-2022-04-05.xml*