Computer Science
Computer Science-2021.gz:143279
Computer Science Engineering
9122be0dc716fe1c17800a67275caa8b690e1103
Multiobjective Bilevel Evolutionary Approach for Off-Grid Direction-of-Arrival Estimation

The source number identification is an essential step in direction-of-arrival (DOA) estimation. Existing methods may provide a wrong source number due to inferior statistical properties (in low SNR or limited snapshots) or modeling errors (caused by relaxing sparse penalties), especially in impulsive noise. To address this issue, we propose a novel idea of simultaneous source number identification and DOA estimation. We formulate a multiobjective off-grid DOA estimation model to realize this idea, by which the source number can be automatically identified together with DOA estimation. In particular, the source number is properly exploited by the $l_0$ norm of impinging signals without relaxations, guaranteeing accuracy. Furthermore, we design a multiobjective bilevel evolutionary algorithm to solve the proposed model. The source number identification and sparse recovery are simultaneously optimized at the on-grid (lower) level. A forward search strategy is developed to further refine the grid at the off-grid (upper) level. This strategy does not need linear approximations and can eliminate the off-grid gap with low computational complexity. Simulation results demonstrate the outperformance of our method in terms of source number and root mean square error.


Introduction
Direction-of-arrival (DOA) estimation refers to finding the direction information of electromagnetic waves/sources according to the outputs of receiving antennas/sensors that form an array, as shown in Fig. 1. It is a crucial subject in array signal processing and is ubiquitous in radar, sonar, and wireless communications [1].
Since the number of sources is normally unknown a priori, methods for DOA estimation need to identify the source number as well as estimate the direction of each source. Existing methods can be divided into two categories. The first category is subspace-based methods [2], which exploits the eigenvalues decomposition of the sample covariance matrix of the array outputs. They use intuitive thresholding or information criterion to remove spurious sources from the obtained ones [3][4] [5], such that the source number is identified a posteriori.
The second category is sparse methods [6][7][8] [9], which can be further divided into on-grid, gridless, and off-grid types. The on-grid type separates the continuous angular range into finite grid points and assumes sources fall exactly on the grid [10], but this assumption is not always the case in practice. The gridless type directly operates in the continuous angular domain. This type needs to solve a time-consuming semidefinite programming problem and only works well in cases with well-separated sources [11] [12]. In contrast, the off-grid type adjusts the initial grid points to approach the sources [8] [13] [14], which alleviates the gap between grid points and sources.
One major challenge in the above methods is the identification of source number. For subspace-based methods, the a posteriori processing tends to provide a wrong source number in limited snapshots or low signal-to-noise ratio (SNR) scenarios because of inferior statistical properties. For sparse methods, apart from postprocessing, l p -norm (0 < p 1) sparse penalties, Laplace prior [8], and Gaussian scale mixtures [15] have also been employed to capture the sparsity of impinging signals. These penalties/priors are relaxations of the l 0 norm of impinging signal matrix, which bring modeling error and bias the sparse distribution, subsequently overestimating the source number [16]. To address the challenge in source number identification, in this paper, we propose a novel idea of simultaneous source number identification and DOA estimation. We design a multiobjective bilevel evolutionary approach (MoBEA) with two significant innovations to realize this idea. Firstly, we bulid a multiobjective DOA estimation model, in which the source number and a robust fitting error are taken as two conflicting objectives. Unlike existing DOA estimation models, our model can automatically identify the source number together with DOA estimation. Moreover, our model perfectly exploits the source number via the original l 0 -norm of impinging signals without relaxation, providing a more precise source number. Besides, different from existing models explicitly or implicitly working with the assumption of Gaussian noise, we incorporate the robustness metric "correntropy" [17] into the model, greatly improving the performance in impulsive noise.
Secondly, we design a multiobjective bilevel evolutionary algorithm to solve the proposed model. Our algorithm includes on-grid (lower) level and off-grid (upper) level optimization. At the on-grid level, we develop a population-based genetic algorithm to simultaneously identify the source number and recover the impinging signal matrix on a coarse grid. The diversified solutions in the population evolve various source numbers and communicate with each other, providing multiple pathways to the optima. The knee solution is finally extracted from the obtained solutions, which represents the identified source number and active grid points. At the off-grid level, we introduce a forward search strategy to refine the active grid points. This strategy does not need linear approximations and can efficiently eliminate the off-grid gap with low computational complexity. Overall, MoBEA's main contributions are: • Multiobjective DOA estimation model. The proposed model can automatically identify the source number together with DOA estimation. Besides, it perfectly exploits the source number via the l 0 -norm of impinging signals without relaxation, guaranteeing its accuracy.
• Multiobjective bilevel evolutionary DOA estimation algorithm. It solves the proposed model via alternate execution of two levels' optimization. A forward search strategy is introduced in the off-grid level, which efficiently eliminates the off-grid gap with low computational complexity.
• Empirical validation of MoBEA's performance. Simulation results confirm the superiority of MoBEA in both source number identification and DOA estimation over state-of-the-art methods.
In this paper, bold-face letters represent vectors and matrices. R and C are the real and complex domain. E denotes the mathematical expectation. T , * , and H denote transpose, conjugate, and conjugate transpose of a vector or matrix, respectively. A i,: and A i,j denote the i-th row and (i, j)-th element of matrix A, respectively. diag(w) is a diagonal matrix with w as the diagonal elements. In particular, s| e stands for the sub-vector of s with entries indexed by the set I = {i|e i = 1}. Similarly, S| e denotes the sub-matrix of S with rows indexed by the set I = {i|e i = 1}.
The rest of the paper includes related work in Section II, the sampling signal model in Section III, the proposed DOA estimation model and algorithm in Sections IV and V, respectively, simulation results in Section VI, and conclusions in Section VII.

Related Work
DOA estimation aims at finding the direction information of electromagnetic sources according to the outputs of receiving antennas/sensors that form an array. Methods for this purpose are mainly divided into subspace-based ones and sparse ones. Subspace-based methods exploit the eigenvalues decomposition of the sample covariance matrix of the array output [2]. They employ intuitive thresholding or information criterion as postprocessing to determine the source number [3][4] [5]. In scenarios with limited snapshots or low SNR, this postprocessing step may provide a wrong source number because of inferior statistical properties.
Sparse methods are inspired by the compressed sensing theory. Representatives include greedy algorithms [6], convex optimization approaches [7], and sparse Bayesian learning (SBL)-based methods [8] [9]. These methods show many prominent properties over subspace-based methods, e.g., enhanced robustness to lower SNR, fewer snapshots, and highly correlated sources [18] [19]. Sparse methods can be divided into on-grid, gridless, and off-grid types. The on-grid type separates the continuous angular range into finite grid points and assumes all sources fall exactly on the grid [10]. However, some sources may deviate from the initial grid in practice (i.e., the off-grid gap), degrading on-grid methods' performance.
Gridless methods are based on the atomic norm theory [11] [12] or covariance fitting criterion [20]. They operate in the continuous angular domain, thereby eliminating the off-grid gap. However, they need to solve a time-consuming semidefinite programming problem and only work well in cases with sufficiently separated sources.
Off-grid methods are another type for alleviating the off-grid gap [13][8] [14]. They parameterize the gap into the models, which modifies the initial grid to approach sources. Several strategies have been proposed for off-grid gap estimation, such as linear approximations [13] [8], orthogonality-based methods [14], and root-finding strategies [15] [21]. Their performance depends on the trade-off between estimation accuracy and computational complexity.
Apart from the aforementioned intuitive thresholding or information criterionbased methods [3][4] [5], l p -norm (0 < p 1) sparse penalties, Laplace prior [8], and Gaussian scale mixtures [15] have also been used in sparse methods to identify the source number. These penalties/priors are relaxations of the l 0 norm of impinging signal matrix, which bring modeling error and may not perfectly represent the sparse distributions, leading to overestimation of the source number [16]. Identifying the source number is still an open challenge in DOA estimation.
In addition, DOA estimation methods are normally developed for Gaussian noise. However, the noise may exhibit non-Gaussian characteristics, e.g., impul-sive noise (also known as burst noise) [22], in practice. To deal with impulsive noise, maximum likelihood subspace methods [23] [24] have been developed for Gaussian mixture model-based noise environments. Other robust subspacebased methods include fractional lower-order moment-based methods [25] [26], and the zero-memory nonlinear functions based methods [27]. Different from them, the l p -MUSIC considers the sample covariance matrix with a robust metric, i.e., minimizing the l p -norm (p ∈ [1, 2)) of the residual fitting error matrix to handle outliers [28]. In [29] [30], an effective correntropy induced estimator is integrated into MUSIC to reduce the influence of outliers. Several robust sparse methods have also been proposed. In [31], an additional l 2 -norm constraint on the parameter of interest is introduced to improve the estimation robustness. In [32] and [33], robust relaxation methods by imposing a l p fidelity (with p = 1 and 1.2 respectively) for denoising are presented. Two robust SBL-based methods [15] model the additive noise as the mix of dense Gaussian noise and sparse noise outliers to get rid of the adverse effect of noise outliers. However, the assumption that the impulsive noise has a sparse structure is not reasonable, leading to a performance loss [34].

Signal Model
Consider K narrowband far-field sourcess 1 (t),s 2 (t), ...,s K (t) impinging on a uniform linear array of M omnidirectional sensors from directions ofθ = [θ 1 ,θ 2 , ...,θ K ] T . The array output at time t is modeled as where y(t) = [y 1 (t), y 2 (t), ..., y M (t)] T is the array output,s(t) = [s 1 (t),s 2 (t), ..., s K (t)] T are the signal waveforms, A(θ) = [a(θ 1 ), a(θ 2 ), ..., a(θ K )] is the array manifold matrix, a(θ k ) = [1, e −j2πd/λ sin(θ k ) , ..., e −j2π(M−1)d/λ sin(θ k ) ] T contains the time delay of the k-th signal received at each sensor, λ is the wavelength of sources, d is the spacing distance between adjacent sensors, n(t) = [n 1 (t), n 2 (t), ..., n M (t)] T is an unknown noise vector. When T snapshots are collected, the array output is To accomplish this goal, existing grid-based sparse methods cast DOA estimation as a sparse recovery problem [8]. The angular range [−π/2, π/2] is divided into an equi-spaced grid θ 0 = [θ 1 , ..., θ N ] T with the assumption that N ≫ K, where N denotes the total number of grid points, and r = θ i+1 − θ i is the grid interval. If all the sources fall exactly on this grid, it would be an on-grid case. The array output Y for the on-grid model is where the impinging signal matrix S is the extension ofS fromθ to θ 0 , with S being row-sparse, i.e., all columns of S are sparse and share the same support. However, some sources may deviate from the predefined grid (known as the offgrid case), as shown in Fig. 2. In this case, the sparse coefficient S is not sparse any more. Thus, the off-grid gap may lead to estimation degradation.
To alleviate the off-grid gap, the gap between true DOAs and the initial grid (denoted as grid mismatch ζ) is usually parameterized into the on-grid model (3). As a result, the array output for the off-grid model can be formulated as The grid mismatch and impinging signal matrix are estimated by solving the following regularization optimization problem: where ρ is a balancing parameter. Y − A(θ 0 + ζ)S 2 F is a fitting error term, denoted as "measurement error"; S is a row-sparse matrix whose i-th row corresponds to a possible source impinging on the array at (θ 0 + ζ) i ; S 2,0 stands for the source number, which equals to the number of nonzero rows of S: S 2,0 is a l 0 -norm penalty. It is the proper and exact sparsity-enforcing penalty, which plays a prominent role in identifying the source number.

Proposed Multiobjective DOA Estimation Model
Problem (P1) is NP-hard because of the l 0 norm [35]. Existing off-grid methods [6][8] [9] relaxed the l 0 norm to other suboptimal penalties. Such relaxation brings modeling error, which may provide some spurious DOAs and degrade the localization performance [16]. Also, the performance is susceptible to the balancing parameter (i.e., ρ in (P1)), which is hard to tune.
To address the above issues, we naturally model the regularization optimization problem (P1) as a multiobjective optimization problem (MOP). Our model holds two significant advantages: 1) it allows automatically identifying the source number together with DOA estimation, and 2) it perfectly exploits the source number via the l 0 -norm of impinging signals without relaxation, guaranteeing accuracy. In detail, our model simultaneously minimizes the two conflicting objectives: where S 2,0 and Y − A(θ 0 + ζ)S 2 F are the source number and measurement error, respectively.
The MOP (P2) inherently involves two problems. 0ne is sparse recovery of the impinging signal matrix S with given grid, where S 2,0 denotes the source number. The other is grid refinement with given S. It should be noted that, compared to source number identification by sparse recovery, perceiving true DOAs is more critical. Therefore, the decision maker of the grid refinement problem should have the complete knowledge about sparse recovery, while the sparse recovery problem only observes the decisions of grid refinement. Therefore, we cast the MOP (P2) as a bilevel MOP: where the first and the third rows are the off-grid (upper) level and the on-grid (lower) level problems, respectively. The on-grid level simultaneously identifies the source number and recovers S with given grid. Based on the decisions of the on-grid level, the off-grid level minimizes the measurement error to refine the grid and update S. The bilevel MOP (P3) explicitly works in Gaussian noise. In practice, the impulsive noise may exist, bringing great challenges for accurate estimation. To cope with this, we incorporate a robust metric "correntropy" [17] to reduce the detrimental effect of the impulsive noise. Correntropy is a local and nonlinear similarity measure in a feature space [17]. Given two arbitrary vectors x = [x 1 , ..., x M ] T and z = [z 1 , ..., z M ] T , the correntropy is approximated as [17] V

Variable Description
G the outer generation number lt the inner generation number at the on-grid level ut the inner generation number in forward search θ 0 the initial grid ζ the grid mismatch the active set of the knee solution at G-th generation S G the impinging signal matrix of the knee solution at G-th generation where κ σ (γ) is the Gaussian kernel function with a kernel size σ. Based on correntropy, correntropy-based loss function (CLF) [36] is defined as This function is related to the Welsch's cost function [17]. Compared to the l 2 -norm error term of the bilevel MOP (P3), V CLF increases much slower and is bounded, thus large noise outliers have a limited effect on V CLF . Hence, it has better robustness to impulsive noise. With incorporating V CLF , we give the final bilevel multiobjective optimization model where the first and the third rows are the off-grid level and on-grid level problems, respectively. The second row is the grid mismatch-related constraint. Compared to the regularization problem (P1), no balancing parameter is needed. The sparsity of impinging signals is perfectly captured by the l 0 norm term without evoking relaxations. Hence the source number can be accurately identified.

Proposed Multiobjective BiLevel Evolutionary DOA Estimation Algorithm
We design a multiobjective bilevel evolutionary algorithm to solve the proposed bilevel MOP (P4). The designed algorithm has three key features: 1) it simultaneously identifies the source number and recovers the impinging signal matrix on a coarse grid at the on-grid level; 2) A forward search strategy is developed at the off-grid level, which eliminates the off-grid gap with lower computational complexity than other grid refinement strategies [13][14] [21]; and 3)  the population-based evolutionary algorithm evolves diversified search pathways to the optima, promoting the effectiveness and efficiency of the algorithm.
The workflow of the proposed algorithm is shown in Fig. 3, with the pseudocode and variable definitions exhibited in Algorithm 1 and Table 1, respectively. It starts with initialization, following by a bilevel optimization. At the on-grid level, an improved framework based on NSGA-II [37] is developed to simultaneously identify the source number and recover the impinging signal matrix S. Instead of searching the complex-valued S directly, we encode S into a binary vector for simpleness. After evolutionary search and selection, the knee solution is extracted to represent the active grid points, where "active" means a source located at that grid point. After that, at the off-grid level, a forward search strategy is developed and executed for refining the grid represented by the knee solution. The population is then updated. When the stopping criterion is satisfied, the active grid points are the estimated DOAs. Key components of the proposed algorithm are detailed below.

Chromosome Encoding and Decoding
(S, ζ) represents the solution of the proposed model (P4), but it is non-trivial to execute genetic operators on the complex-valued impinging signal matrix S. Therefore, we search the locations of the nonzero rows of S instead and then recover these nonzero values. In detail, we encode (S, ζ) as (e, ζ), where e is Algorithm 1 Pseudocode of the proposed algorithm Input: A, Y, θ 0 , ζ 0 ← 0 N ×1 Output:θ 1: G = 0; 2: (P G , S G ,ë G ,S G ) =Initialization(A, Y, ζ 0 ); 3: while "the algorithm does not teminate" do 4: /*Population updates*/ 7: a binary vector with "1" and "0" denoting the corresponding row of S being nonzero and zero entries, respectively. For unity, we denote e as "active set".
To recover S from e, we develop a decoding method improved from the correntropy matching pursuit (CMP) algorithm [36]. For clarity, we briefly introduce CMP before describing the decoding process.

Correntropy matching pursuit
CMP is a matching pursuit method for sparse recovery in impulsive noise [36]. In compressive sensing, let the measurement being y = As + n with the assumption that the dictionary A and measurements y are both known, CMP solves a correntropy loss minimization problem to recover s s = arg min where z = y − As is the residual. With the half-quadratic theory [38], problem (7) can be reformulated as a weighted least square problem and solved by [36]: where s l+1 | e is the sub-vector of s l+1 with entries indexed by I = {i|e i = 1}, g σ (p) = exp(− pp * 2σ 2 ), and supp(s) denotes the nonzero entries of s. w is a weight vector. It indicates the importance of measurements, i.e., the small coefficients in w suppress the severely contaminated measurements, and the

Algorithm 2 Decoding
Input: A, Y, E,S, ζ, σ, S = ∅ Output: put S into S; 8: end for large coefficients enable to preserve the clean measurements. w can improve CMP to identify the true atoms, thus boosting the recovery performance [36].

Decoding
We extend the CMP to the case of multiple measurement vectors, as the decoding process of the proposed algorithm (see Algorithm 2). In Algorithm 2, W is a weight matrix with each column being the weight vector at each snapshot. To decode the active set, we first compute the weight matrix based on the impinging signal matrix of the knee solutionS, where the knee solution will be detailed in Section 5.2. Then, the impinging signal matrix is recovered by solving a series of reweighted least square problems.
In the decoding process, the choice of kernel size σ is important. To balance the recovery accuracy and convergence speed, we employ the kernel annealing method [17] instead of (8) due to its greater performance enhancement found in [39]. Specifically, the kernel size is defined by where G is the generation counter of the outer cycle, σ min =0.03, σ max is calculated by 0.5(abs(Y) 0.875 − abs(Y) 0.125 ) − σ min with abs(Y) l being the l-th quantile of {|Y i,j |} M×T , and ν = 2 × 10 −4 is the decay rate.

Initialization
We leverage the statistical correlation knowledge to provide a preliminary choice for the active set to enhance the search efficiency. Since a uniform linear array can identify up to M − 1 sources, the condition K ≤ M − 1 can be considered as a priori [40]. In noise-free cases, the sources can be roughly located by selecting the K atoms that are most correlated with the measurements. Motivated by this and consider the atoms' rank distortion caused by impulsive noise, we give priority to the atoms that have larger correlations, e.g., those who rank the top 2M . In this way, some proper atoms are involved with great probability, thus improving the search efficiency. The process is shown in Line Algorithm 3 Initialization Input: A, Y, ζ 0 Output: P 0 , S 0 ,ë 0 ,S 0 1: τ = T t=1 | < y(t), A :,j > |, j = 1, ..., N ; 2: Λ ← save the indices of 2M highest correlation in τ ; 3: E 0 ← randomly choose no more than M − 1 indices from Λ to constructN active sets; 4: P 0 = (E 0 , ζ 0 ); 5: S 0 =Decoding(E 0 , 0 N ×T , ζ 0 , σ); // Algorithm 2 6:ë 0 ,S 0 ←Knee Identification(P 0 , S 0 ).

1-5 of Algorithm 3. For the missing atoms, they can be further searched in the mutation step of the on-grid level.
After that, the knee solution is selected from the current population for preparing the next decoding. The reason for the above operation is that in each decoding process, a solution is required to estimate the weight matrix. Here we choose the knee solution based on the following considerations. The knee solution is commonly an interesting point along the Pareto front (PF), and it provides a promising trade-off between the two objectives [41]. Therefore, it is expected to use the knee solution to obtain reasonable weights and bring performance enhancement. Even when the knee solution fails to be the most satisfactory, it is still a Pareto nondominate solution from the perspective of multiobjective optimization. Here the kink method [42] is employed to identify the knee solution, where the solution with the largest slope variance over the PF is the knee solution.

On-grid Level Optimization
At the on-grid level, we develop a framework based on NSGA-II [37] to simultaneously identify the source number and recover the signals, where the grid mismatch acts as parameters. The procedure is exhibited in Algorithm 4. The active set E lt firstly undergoes the one-point crossover, and the bitwise mutation operations [43] to create offspring solution setÊ lt .Ê lt is then decoded according to Algorithm 2 to recover the corresponding signals S lt . After that, the environmental selection operator of NSGA-II [37] is executed to selectN elite solutions for the next generation. Lastly, the knee solution is identified. This solution will be used to estimate the weight matrix in the next decoding and to identify the active grid points for off-grid level refinement.

Forward Search
We propose a straightforward forward search strategy at the off-grid level for solving the grid mismatch. Unlike existing grid refinement strategies, this strategy does not require modeling approximation, hence achieving higher DOA estimation quality. To save computational cost, the grid points corresponding to the knee solution's nonzero entries are viewed as "active" and allowed to be P lt+1 , S lt+1 ←Environmental Selection(P lt ∪ Q lt , S lt ∪Ŝ lt ); 6:ë lt+1 ,S lt+1 ←Knee Identification(P lt+1 , S lt+1 ); 7: end for refined. The perturbation direction of each active grid point where to decrease the correntropy-based loss function is obtained separately. Then all the active grid points are perturbed together along with the obtained directions.
The procedure of the forward search is shown in Algorithm 5. The indices of active grid points are saved into I according to the active set of the knee solutionë. Subsequently, the perturbation directions of the active grid points are detected and saved into β (Lines 3 to 10 of Algorithm 5): for each active grid point, a random perturbation direction is given, and it is accepted only if the corresponding correntropy-based loss function decreases; otherwise, the perturbation direction is switched to the opposite direction or zero. After that, the active grid points are perturbed together along with β with a stepsize µ (Line 13), where µ is set to r/100 to provide a high resolution. The perturbation is repeated several times until ut reaches its maximum value ut max or the correntropy-based loss function no longer decreases.

Remark 1.
For any active grid point, it is finally accepted if it is within the set of −r/2 < ζ i r/2. Otherwise, it is rejected, and the corresponding grid point remains unchanged.

Computational Complexity
The main computational complexity of the proposed algorithm lies in the decoding process, i.e., Line 5 of algorithm 3, Line 3 of Algorithm 4 and Line 22 of Algorithm 5. For each solution, the complexity of decoding is T × O(N 3 ) in the worst case, where T and N are the number of snapshots and the number of grid points, respectively. Thus, the algorithm's complexity is (lt where lt is the number of inner iterations of the on-grid level,N is the population size. It can be seen that the number of grid points has a significant impact on the computational complexity. To reduce the execution time, the decoding is suggested to be implemented in parallel in practical applications.
• l p -MUSIC: A representative robust subspace-based algorithm that adopts the l p -norm of the residual fitting error matrix for subspace decomposition.
• MCC-MUSIC: A state-of-the-art robust subspace-based algorithm that estimates the signal subspace by solving an optimization problem under the maximum correntropy criterion.
• Bayes-optimal: A state-of-the-art off-grid SBL-based algorithm that models the measurement noise as the mix of Gaussian noise and outliers.
• Fast-alternating: A state-of-the-art off-grid SBL-based algorithm which is a fast execution version of the Bayes-optimal algorithm.

Simulation Settings
The comparison algorithms need the source number K as the input. Here we set K = M − 1 1 . For l p -MUSIC, we set p = 1.1. The grid interval is set to 0.1 • for l p -MUSIC and MCC-MUSIC, and 2 • for the remaining methods if not stated. For MoBEA, the crossover probability and the mutation probability are set to 0.9 and 1/N , respectively, following the practice of [37]. The population size N is 50. The on-grid level optimization terminates when its inner generation reaches 50 or the knee solution remains unchanged in five consecutive generations.
In all simulations, we consider a small number of uncorrelated sources impinging on a uniform linear array of sensors with an inter-sensor spacing of d = λ/2. The number of receiving sensors is M = 8. For a fair comparison, all methods stop running when the variance change of the impinging signal S G is less than 10 −6 in five consecutive generations, or the outer-cycle generation reaches 200; and the total number of Monte Carlo trials is 100 for all algorithms. All the experiments are implemented in MATLAB R2018b on a laptop with Intel i5-8265U CPU and 8GB RAM.
Two widely-used PDF models, i.e., Gaussian mixture model (GMM) [23], and symmetric α-stable distribution (SαS) [44], are considered to model the impulsive noise. All sources are assumed to possess the equal power η 2 s .

GMM
The PDF of a two-term Gaussian mixture noise n(t) can be described as where 0 ≤ c i ≤ 1 and η 2 i are the probability and variance of the i-th term with c 1 + c 2 = 1, respectively. Assume η 2 2 = 100η 2 1 and c1 > c2 > 0, the SNR noise model can be viewed as the large noise outliers of variance η 2 2 with a smaller probability c 2 embedded into Gaussian noise of variance η 2 1 with a larger probability c 1 . According to [15], the SNR is simplified as SNR = η 2 s /η 2 1 .

α-Stable distribution
The symmetric SαS distribution with zero-location is used to model the noise. Its characteristic function is defined as where 0 < α ≤ 2 denotes the characteristic exponent that indicates the tail of the distribution, and γ is the scale. When α = 2 and α = 1, the SαS reduces to the Gaussian distribution and the Cauchy distribution, respectively. The smaller value of α, the more impulsive noise is. Since the closed-form PDF of SαS does not exist when α = 2 and α = 1 [44], the SNR becomes meaningless. Instead, the generalized SNR (GSNR) [28] is reformulated as GSNR = η 2 s /γ α .

Performance metrics
Since the Bayesian-based algorithms and the proposed MoBEA do not output the spatial spectrum, the spatial spectrum will not be used for comparison in this paper. We employ two statistical measures, i.e., the root mean square error (RMSE), and the average estimated source number.
Let (θ k ) i stands for the estimate ofθ k at the i-th Monte Carlo trial, the RMSE can be formulated as where K is the true source number, and υ = 100 is the total number of Monte Carlo trials. For a given simulation point, the RMSE is obtained by averaging only the trials in which the estimated source number is greater than or equals to K. The assignment of estimated DOAs to the true one is executed based on the Hungarian algorithm [45]. The average estimated source number is obtained by averaging the empirical source numbers of all Monte Carlo runs. These two metrics well reveal the identification capacity of the DOAs' positions and number.

Detailed Analysis of MoBEA
In this subsection, the efficacy of the knee solution and the forward search scheme are investigated to demonstrate the superior performance of MoBEA.

Efficacy of the knee solution
At the on-grid level, the knee solution extracted from the PF is employed to represent the active grid points for grid refinement. It provides a promising trade-off between the source number and the estimation accuracy.
We first investigate the relationship between the knee solution and the source number. Fig. 4 shows the typical final PF obtained by the on-grid level optimization of 100 trails under different GSNRs in SαS noise with α = 1.4. The corresponding settings are as follows: three uncorrelated sources are from −9.7 • , 6.8 • , and 12.7 • . The number of snapshots is 20. Since a uniform linear array can identify up to M − 1 sources [40], it is expected to identify the knee solution from the non-dominated solutions whose source number is less than M (where M = 8). It can be observed that the knee solution in the red circle corresponds to the true source number K = 3.
To further exploit the impact of the knee solution on the accuracy of DOA estimate, the performance of MoBEA with the knee solution and a randomly selected nondominated solution entering the off-grid level optimization are compared in SαS noise environment. The experimental parameters are the same as those in Fig. 4 except that the GSNR varies from −10 ∼ 15dB. Fig. 5 exhibits the RMSE and average estimated source number versus GSNR. It can be seen that the version with the knee solution achieves a smaller DOAs error compared to the random nondominated solution case, thanks to its more precise estimate of the source number. This result indicates that the knee solution is more valuable than other nondominated solutions in terms of the source number and DOAs.

Efficacy of forward search
To validate the efficacy of forward search, the performance of MoBEA are compared to two versions, i.e., without the off-grid level optimization (denoted as the "on-grid" version) and with Taylor expansion. In the last version, the manifold matrix is approximated using the first order Taylor approximation [13]: where θ 0 and ζ are the initial grid and grid mismatch, respectively, A ′ (θ 0 ) is the first derivative of A(θ 0 ) with respective to θ 0 ; similar to [13], ζ can be obtained by letting the derivative of V CLF (Y, (A(θ 0 ) + A ′ (θ 0 )diag(ζ))S) with respect to ζ be zero.
In this simulation, the DOAs of two signals are from 1.6 • and 13.2 • , respectively. The number of snapshots is T = 20. Fig. 6 shows the results of DOA estimate versus grid interval in GMM noise with c 2 = 0.1 and SNR = 10dB. It is observed that the three strategies have similar results to the source number. The localization accuracy of the "on-grid" version and "Taylor expansion" version deteriorates with the coarser grid due to the ignorance of off-grid mismatch and the large modeling error caused by Taylor expansion, respectively. As expected, the MoBEA with forward search achieves the best performance with different grid intervals, which validates the forward search's superiority in evaluating the grid mismatch.

Performance under various SNRs
This subsection investigates algorithms' performance under different SNRs in GMM noise. Three uncorrelated sources from −2.7 • , 5.8 • , and 20.2 • are considered. The number of snapshots is 20. Fig. 7 shows the RMSE and estimated source number in GMM noise versus SNR with c 2 = 0.1. It can be seen that the performance of all algorithms degrades obviously. l p -MUSIC and  MCC-MUSIC perform the worst under lower SNRs, which may be due to their high sensitivity to lower SNR and imprecise information of source number. The two Bayesian methods behave better in most cases, but they are inferior to MoBEA in localization accuracy. The reason is MoBEA estimating the source number more precisely, which avoids energy leakage on spurious DOAs. Fig. 8 reports the performance in SαS noise versus GSNR with α = 1.4. The two MUSIC-based methods still behave badly in impulsive noise. Compared to the SBL-based methods, MoBEA achieves higher or comparable DOA estimation accuracy but always predicts the source number more accurately.

Performance under various grid intervals
In this simulation, the performance of MoBEA under various grid intervals is investigated. Three uncorrelated sources from −19.7 • , 6.8 • , and 32.7 • are considered. The number of snapshots and SNR/GSNR are 20 and 10dB, respectively. The grid interval varies from 2 • to 10 • . Fig. 9 exhibits the RMSE and average estimated source number versus grid interval in SαS noise (α = 1.4). MoBEA achieves comparable or smaller RMSE results than other algorithms under lower SNRs and performs much better with higher SNRs. The enhancement of MoBEA under lower SNRs is non-significant, probably because the localization accuracy is severely affected by the ignorance of the noise variance in the proposed model. But it is worth noting that, MoBEA shows an absolute predominance in predicting the source number for all grid interval cases.

Performance under various angular separations
This simulation examines algorithms' ability to identify two closely-located sources. We consider two sources. The first is from −10.8 • and the second varies from −8.8 • to −0.8 • . The number of snapshots and SNR are set to be T = 30 and 10dB, respectively. Fig. 10 displays the performance versus angular separations increasing from 2 • to 10 • in GMM noise. The RMSE values of MoBEA retain the lowest in most cases. As the two sources separate from each other, MoBEA predicts the source number more precisely, which performs remarkably better than other algorithms. These results indicate that MoBEA has higher localization resolution when the sources are closely located.

Performance under various numbers of snapshots
In this simulation, the off-grid DOA estimation performance in SαS noise with the number of snapshots at α = 1.4 and GSNR = 10dB is studied. Three sources are randomly chosen from −2.7 • , 5.7 • and 20.2 • impinging on the uniform linear array. Fig. 11 shows the RMSE and estimated source number in SαS noise versus the number of snapshots. It can be observed that the l p -MUSIC and MCC-MUSIC algorithms gain the worst performance under fewer data samples. Fastalternating and Bayes-optimal methods obtain great enhancement in RMSE with fewer snapshots. As expected, MoBEA achieves the highest localization accuracy due to the accurate identification of the source number.   6.4.5. Performance on running time Section 5.5 demonstrates that the number of grid points N has a remarkable impact on the computational complexity of MoBEA. To further quantify the impact, this simulation provides the average running time of MoBEA versus the grid interval compared to other algorithms, as shown in Fig. 12. The experimental settings are the same as those in Fig. 9. It can be observed that the computational time of all algorithms shows a decreasing trend with increasing grid intervals (i.e., fewer grid points). The l p -MUSIC and MCC-MUSIC run much faster than other algorithms. MoBEA and Bayes-based algorithms are slower because they need to execute the time-consuming inversions. If there is a high demand for the faster running speed of MoBEA in practical applications, the decoding is suggested to be executed in parallel. Despite this, MoBEA achieves great improvements in RMSE and estimated source number over the comparing algorithms in most scenarios.

Conclusion
This paper has proposed the MoBEA for DOA estimation. MoBEA involves two innovations. The first is the multiobjective DOA estimation model, in which the source number and a robust correntropy-based fitting error function are taken as two objectives. Unlike existing DOA estimation models, the proposed model automatically identifies the source number together with DOA estimation. Besides, the original l 0 -norm penalty is used to capture the sparsity of signal-of-interest, which avoids relaxing sparse-inducing penalties and enables source number identification more accurately. The second innovation is the multiobjective bilevel evolutionary DOA estimation algorithm for solving the proposed model. The on-grid level is for joint source number identification and sparse recovery, while the off-grid level works on grid refinement via the proposed forward search strategy. This strategy avoids linear approximation and enhances the localization accuracy. Thanks to the population-based search of the proposed algorithm, the solutions own different source numbers communicate to each other during the evolutionary search, which provides diverse search pathways to the optima. Experimental results have demonstrated the effectiveness and efficiency of the proposed MoBEA in identifying the source number and DOAs in impulsive noise.
In the future, parallel processing will be considered to accelerate our algorithm. A future study on reducing the computational complexity of MoBEA will be also explored.