Covariance Matrix Reconstruction for Direction Finding with Nested Arrays Using Iterative Reweighted Nuclear Norm Minimization

In this paper, we address the direction finding problem in the background of unknown nonuniform noise with nested array. A novel gridless direction finding method is proposed via the low-rank covariance matrix approximation, which is based on a reweighted nuclear norm optimization. In the proposed method, we first eliminate the noise variance variable by linear transform and utilize the covariance fitting criteria to determine the regularization parameter for insuring robustness. And then we reconstruct the low-rank covariance matrix by iteratively reweighted nuclear norm optimization that imposes the nonconvex penalty. Finally, we exploit the search-free DoA estimation method to perform the parameter estimation. Numerical simulations are carried out to verify the effectiveness of the proposed method. Moreover, results indicate that the proposed method has more accurate DoA estimation in the nonuniform noise and off-grid cases compared with the state-of-the-art DoA estimation algorithm.


Introduction
Source localization is always a significant research direction in the past decades and today which is widely used in various fields including radar, sonar, wireless communication, and acoustics [1][2][3][4].However, achieving superresolution is a great challenge for more sources with the least elements.Compared with the uniform linear arrays (ULAs), the sparse linear arrays (SLAs) can provide much more degree of freedoms (DOFs) and resolve more sources with the same number of physical elements.Meanwhile, the SLAs can efficiently reduce the cost and power consumption.Therefore, the study of SLAs has attracted more and more researchers' great attention [5,6].
So far, there are some works about SLA configurations.The minimum redundancy arrays (MRAs) are class optimum lag SLAs [5], which provide a complete set of spatial lags between pairs of elements with minimum redundancies.However, it is not easy to attain the optimum design of MRAs, especially for a large array; it is time-consuming to search the optimum structure.In the past decades, two potential array configurations, i.e., coprime arrays (CPA) and nested arrays (NA), have drawn the researchers' attentions since they have exactly closed-form expressions for sensor locations and it is easy to predict the attainable DOFs [7,8].Coprime arrays can resolve up to o MN sources with only M + N − 1 elements; however, the CPAs do not generate a filled coarray since CPA has holes in the difference coarray [7].So it cannot directly apply the augmentation techniques.By contrast, nested arrays can generate a filled difference coarray except the mentioned advantages [8].Although the NAs and CPAs are not optimum lag arrays like MRAs, they are the most attractive array configuration because they are easy to build in the last decade.
The processing of SLA mainly has array interpolation (AI) [9][10][11] and the Toeplitz completion method [12][13][14].The classical AI technique [10] can effectively obtain the element data of virtual ULAs, which imposes a linear interpolation process on the element data of a real SLA and selects interpolator coefficients by minimizing the interpolation error for a source coming from a certain angular sector.The drawback is that this technique needs to know the angular sector.In AI techniques, Wiener array interpolation (WAI) [15] is a practically attractive method, which exploits a maximum likelihood method to estimate the power of signal and noise and use the calibration angles to recover the array steering matrix; hence, it can approximate the MSE optimum solution.However, WAI requires the initial DoA estimation.In the classical Toeplitz completion method, the direct augmentation approach (DAA) is a widely used method for improving the covariance matrix of the SLA [12].DAA constructs a Toeplitz matrix using the sample covariance matrix; since there is one-to-one correspondence between covariance lags and spatial lags, the diagonal elements could be obtained by redundant averaging.Unfortunately, Toeplitz completion method does not guarantee positive semidefinite augmented covariance matrix.In order to construct a positive definite augmented matrix, the authors in [13,14] proposed an iterative DAA algorithm; however, the complicated iteration procedure cannot guarantee the global convergence.Besides, the coarray MUSIC algorithm is used to estimated angle parameters [16] which utilizes the Toeplitz properties of matrices.
Sparse signal representation (SSR) framework has attracted a great interest in direction finding [17][18][19][20][21][22][23], which exploits the spatial sparsity of the signal arriving angles.Despite these aforementioned SSR-based methods have some attractive features, there are still two problems need to be considered.The first problem is nonuniform noise [23,24]; many of SSR-based methods assume that the element noises are spatially uniform, and the reason for this is the uniform noise assumption benefits to choose an appropriate regularization parameter in a sparse representation model.However, in some practical applications, due to the nonidealistic of the practical arrays, such as the nonideality of the receiving channel, the nonuniformity of the element response, and the mutual coupling between elements, the uniform white noise assumption among all elements may not be satisfied.Therefore, the diagonal elements in noise covariance matrix should be considered as arbitrary values which represent the noise levels.Once the assumption misfits the true noise levels, the performance of conventional SSR-based DoA estimation approaches may thereby degrade severely.The second problem is the off-grid problem (also called basis mismatch) [17,25]; in fact, a great deal of SSR-based DoA estimation methods assume that the true target directions exactly lie on the prespecified angular grid.In practical application, the angle parameters are continuous variables, the on-grid assumption rarely holds, and the basis mismatch always exists, which leads to the accuracy degradation of DoA estimation.
In this paper, we propose a gridless direction finding method under the unknown nonuniform noise over with nested array, which is based on low-rank covariance matrix approximation by an iteratively reweighted nuclear norm minimization.The advantages of the proposed method are in three aspects: Firstly, we mitigate the noise variance variable by linear transform to reduce the effect of unknown nonuniform noise.Secondly, we utilize the covariance fitting criteria to determine the regularization parameter for the accuracy error fitting.Last but not least, we develop a novel reweighted nonconvex penalty objective function for exactly approximating the low-rank covariance matrix.After attaining the covariance matrix, we exploit the search-free DoA method to estimate parameters.Numerical simulation is carried out to verify the effectiveness of the proposed method.
The rest of this paper is organized as follows: In Section 2, we present the signal model with nested array.In Section 3, the DoA estimation with iteratively covariance matrix reconstruction (ICMR) is introduced to the nested array, which is based on an iteratively reweighted nuclear norm minimization algorithm.Section 4 validates the proposed method by numerical simulations and comparing the algorithm with state-of-the-art DoA estimation methods based on nested array.Section 5 concludes this paper.
Notation.We define a vector x with a boldface lowercase and a matrix X with a boldface uppercase.The symbol • * denotes the complex conjugate of vector x or matrix X, and the transpose and the Hermitian transpose of vector or matirx are expressed by • T and • H , respectively.E • denote the statistical expectation, vec • denotes stacking element-wise in a column, rank • and • * denote the rank and nuclear norm of the matrix, respectively, ⊗ denotes the Kronecker product, ⊙ denotes the element-wise multiplication between two matrices with the same dimensional, diag • denotes the diagonal operator that takes a vector to a matrix with vector on the diagonal, sign is signum function, σ X denote the eigenvalue of matrix X, σ i X denote the i-th largest eigenvalue value, I k denotes a K × K dimensional identity matrix, • 2 denotes the l 2 norm, and Pr • denotes probability density function.

Array Model for Nested Array
We consider a two-level nested array consisting of two ULA subarrays, where the first ULA subarray has N 1 elements with the interelement spacing d 1 = d, while the second ULA subarray has N 2 elements with the interelement spacing N 1 + 1 d 1 ; the total number of element is M = N 1 + N 2 , as shown in Figure 1.
Let the first element of the array as the phase reference point, the position set of element in ascending is denoted as Assume that K uncorrelated narrow-band far-field sources from θ = θ 1 , ⋯, θ K T impinge on the nested array.The array output at time t can be represented as where a θ k = e jπp 1 sin θ k , e jπp 2 /sin θ k , ⋯, e jπp M /sin θ k is the array steering vector and the element position is defined in (1).We note that each position p i is measured with half wavelength of indent signal source.Let the M denote the maximum interelement distance (i.e., the array aperture).
Then M is an integral multiple of the half signal wavelength.Nested array possesses the same set of interelement distances as the ULA which composed of M elements separated by half signal wavelength.
From the above discussion, we can derive that nested array can be regarded as spatial sparse sampling on the ULA, i.e., we can take the position set P of nested array from the position subset P = 0, 1,⋯,M − 1 of the elements of the ULA, as shown in Figure 1.Thus, nested array can be represented by its element index set P ⊂ P .For the k-th signal source, the steering vector a P θ k of the ULA can be given as a P θ k = e j2πP 1 d/λ sin θ k , ⋯, e j2πP M d/λ sin θ k .According to the index set P of the NA, we design a selection matrix Γ ∈ 0, 1 M×M such that the i-th row element contains 0s except a single 1 at the P i -th location.From the spatial sparse sampling viewpoint, we have a θ k = Γa P θ k .Then, the nested array model can be expressed as Under ideal condition, the array covariance matrix of the output x t of nested array can be calculated as where N is the number of collected snapshot, R s = diag σ 2 s,1 , ⋯, σ 2 s,K , and Toep P u = A P R s A H P for some u P ∈ ℂ M×1 .Furthermore, the noise-free covariance matrix Toep P u is a Hermitian Toeplitz matrix, which is determined by M complex numbers.
It is evident that rank Toep P u = K ≤ M − 1, and Toep P u is positive definite, i.e., Toep P u ⪰ 0. For nested array, Toep u holds that For the NA, D is a redundancy array which can be regarded as a coarray defines on the NA.The matrix Toep u contains all elements u P of explicitly by (6).It implies that the relation A coarray of nested array with six elements is shown in Figure 2. When the positions of nested array are P = 0, 1, 2, 3, 6, 10 , the coarray of nested array composes of M + 1 = 11 3 International Journal of Antennas and Propagation elements with the position P = 0, 1,⋯,10 .From Figure 2, we can observed that the coarray is a redundancy array, which guaranteed the reconstruction of the covariance matrix [17].
One of the aims of this paper is to recover the noise-free covariance matrix Toep P u of the nested array output under nonuniform noise environment.In the next section, we introduce the reweighed nuclear norm minimization to approximate the low-rank noise-free covariance matrix Toep P u .

Iterative Covariance Matrix
Reconstruction for Direction Finding 3.1.Nuclear Norm Minimization Optimization.As mentioned above, the key problem is recovery of the noise-free covariance matrix Toep P u .On one hand, we note that the covariance matrix Toep P u has three key features: Toeplitz-structured, low rank, and positive definite (PSD); on the other hand, in application, the sample covariance matrix of NA is calculated by R = ∑ N n=1 x t n x H t n /N.Due to the finite samples, there exists the perturbation error term between the R and R, whose distribution benefits to choose the regularized parameter [17].Making full use of the two constraints, we formulate the recovery problem of the covariance matrix as the following optimization problem: , which determines the quality of fitting term.According to the covariance fitting criteria [17,18,26], the vectorization form of the perturbation error term ΔR satisfies the asymptotic normal (AsN) distribution as follows: After prewhitening processing, we obtain We can further derive the following formulation: where Asχ 2 M 2 represents the asymptotic chi-square distribution with M 2 degrees of freedom.According to (11), it is quite obvious how to choose the parameter η.We choose the parameter η large enough such that the inequality Unfortunately, this problem in ( 8) is nonconvex optimization problem and it is also a NP-hard problem.To sidestep the nonconvexity, we exploit convex relaxation technique, i.e., we relax the rank norm of Toep P u to the nuclear norm of Toep P u for approximating the low-rank matrix.Therefore, the convex relaxation form for (8) can be represented as min u,σ Combining ( 5) and ( 12), the optimization problem is rewritten as It is worth noting that there is noise variance variable that needs to be optimized in model (13).In fact, we can eliminate this variable by estimating the noise variance in advance, but the estimation of the nonuniform noise variance is not easy to attain directly, especially in the underdetermined DoA estimation circumstances while the sample covariance matrix no longer contains a low-rank component spanning the noise subspace.We can use two methods to handle this case, the first method is that we treat noise variance as a variable to be optimized in the optimal problem such as (13).The second method, we get rid of noise variables by linear transform firstly and then solve the optimization problem; in this paper, we focused on the latter.
It is interesting to observe that the noise term has M nonzero elements which contain in the diagonal element of the covariance matrix and correspond to the unknown noise variance σ 2 n .Based on this observation, we carry out the linear transform to remove the element of the noise variance terms.Mathematically, this operation can be formulated as International Journal of Antennas and Propagation where J is an M M − 1 × M 2 linear transform selecting matrix and can be represented as where 16 1 column vector with one at the ith position and zeros elsewhere.J m ∈ ℝ M 2 ×M .By applying the above operation, the unknown noise variance is eliminated effectively and easily.This elimination operation avoids the estimation of noise variance, therefore facilitating a noise-free low-rank covariance matrix reconstruction.Through the above derivation, ( 13) can further reformatted as From (17), we observe that there is only one vector u that needs to be optimized by the above processing.At the same time, the prewhitening matrix and regularization parameter should be changed.Combining the distribution of the perturbation error term with (14) yields Applying the prewhitening processing [19], we have where the prewhitening matrix C −1/2 is Hermitian square Like the choosing parameter η, the parameter η can be determined by η = χ 2 p M M − 1 .

Iteratively
Reweighted Nuclear Norm Minimization.We utilize the nuclear norm 12) instead of the rank norm rank Toep P u = σ Toep P u 0 = ∑ i g 0 σ i Toep P u in (8), where g 0 x = sign x , x ≥ 0, g 1 x = x, and x ≥ 0. The advantage is effective computable and the best convex approximation, but the drawback is the worst fitting for rank approximation.We can observe that there is a large gap between the rank and nuclear norm for robust and exact reconstruction of low-rank covariance matrix.Therefore, it is essential to find a surrogate function of rank function which not only have the better approximation but also needs a low computational complexity.Some methods utilize a class of nonconvex surrogate functions to approximate the rank function, such as the works in [27,28].
Among the nonconvex approximation functions, we find f x = 1 − e −x/ϵ closely matches g 0 x , so we use this function h ϵ x = f x = 1 − e −x/ϵ as the approximation function of g 0 x , i.e., rank norm.Combine the model (17), the following general nonconvex optimization model is proposed: where h ϵ σ Toep P u = ∑ i f σ i Toep P u denotes the rank approximating function.We expect that the proposed nonconvex optimization model can bridge the gap between the rank norm model ( 8) and nuclear norm (12).
Since the model in ( 20) is nonconvex optimization problem, how to guarantee for convergence to a global minimum is a key problem.Now, we propose an iteratively reweighted nuclear norm minimization to obtain a suboptimal solution, which is based on the framework of the majorization-maximization (MM) method.The MM method is an iterative approach, it iteratively solves a sequence of weighted convex surrogate optimization problems whose weights are attained by the previous optimal value.Thus, by calculating the sequence of convex problems given in (22), we can minimize the optimization function in (20).For the k-th iteration, the weights u are calculated by the previous optimal variable u k−1 .Since h ϵ Toep P u = ∑ i f σ i Toep P u is concave when f x = g 1 x , we have Therefore, the optimization problem for the k-th iteration becomes where W ϵ k ≜ ∇h ϵ Toep P u k .In order to calculate W ϵ k , we start by introducing a key useful proposition as follows.
Proposition 1. Suppose that h ϵ σ Z is expressed as

6
International Journal of Antennas and Propagation matrices with the eigenvalue decomposition (EVD) Z = U diag σ Z U H , Z ⪰ 0, and h ϵ and f are concave and differentiable, then the gradient of f Z at Z is calculated as where η Z = ∇h ϵ σ Z denotes the gradient of h ϵ at σ Z .
In (22), Z = Toep P u k , for the sake of simplicity, we denote σ Toep P u k as σ k and apply the result in [27]; we obtain where U k diag σ k U H k is the EVD of Toep P u k .We observe that ϵ is a key parameter, which controls the closeness between f σ and g 0 σ or g 1 σ .In particular, f σ is close to g 0 σ when ϵ is small, but it is subjected to many local minima.On the other hand, a large value of ϵ causes f σ toward the convex g 1 σ , whereas with poor approximation quality to g 0 σ .In order to avoid getting local minima, the algorithm starts with a large value of ϵ, and then the value ϵ decreases gradually during the iteration.
It is worth mentioning that this strategy is adopted in model (22).We define the weight matrix W ϵ k in each iteration, which uses the solution of the previous iteration for avoiding local minima.The convergence analysis of the proposed method for ( 20) is given as the following theorem, which is proven in the appendix.Theorem 1. Denote u k as the k-th iteration solution in (22), after that, the optimal variable u k satisfies the following properties: (i) h ϵ Toep P u k is monotonically decreasing and bounded as k tends to ∞ (ii) The optimal variable u k converges to a local minimum of (20)

3.3.
Steps of the Proposed Method.To sum up, the steps of the proposed method are presented as follows.Once the solution u of model ( 20) is obtained, we can determine the number of According to classical Vandermonde decomposition lemma [28], we can estimate the DoAs using positive semidefinite Toeplitz matrices Toep P u .So we can apply the grid-free DoA estimation method, such as Prony's method to efficiently estimate θ and σ 2 s .Firstly, by solving a linear system involving only û, we can obtain θ from zeros of a polynomial.Then, σ 2 s is solved by a polynomial rooting-based estimation method.
In the proposed method, firstly, it recovers the complete noiseless matrix Toep P u in advance under the nonuniform noise background, which is attained by solving an iteratively weighted nuclear norm minimization optimization.Compare with the ULAs which have the same number of physical elements, the proposed method forms a larger aperture ULA and gets more degree of freedoms, so it has the potential ability to estimate the DoA in the underdetermined case.Moreover, exploiting the grid-free DoA estimation method to estimate the DoAs, which do not need to discretize the angle grid, avoids the grid mismatch in the grid-based sparse method.

Simulation Results
In this section, we evaluate the DoA estimation performance of the proposed method and compare the proposed method with the others: the state-of-the-art method, spatial smoothing root MUSIC (SS-root-MUSIC) [7], covariance matrix reconstruction approach with root MUSIC (CMR-root-MUSIC) [22], He et al.'s method [23], OGSBL-NA [25], and OGVSBL-PVA [24].To provide a benchmark for evaluating the effectiveness of these methods, the Cramer-Rao lower bound (CRLB) for nonuniform noise given in [23] is also included though it is a lower bound for unbiased estimators.
In the following simulations, we consider a two-level nested array with optimal configuration [8], whose total element number is six.The first ULA subarray contains three elements with interelement spacing d 1 = d; the second ULA subarray contains three elements with spacing d 2 = 4d 1 , where d is equal to half wavelength.We assume that all the signals are far-field equal power uncorrelated signals; the additive noise is spatially nonuniform which is modeled as an independent complex Gaussian random vector with zero mean, and its covariance matrix is R n = diag σ 2 n , since the input signal-to-noise ratio (SNR) can be defined as , 25 where M = N 1 + N 2 is the number of NA elements, σ 2 s denotes the signal source power, and σ 2 n,m denotes the noise power in difference elements.For the performance metric of DoA estimation, the root-mean-square error (RMSE) is selected for evaluating the estimation accuracy of the DoAs, which is defined as , 26 where θ k,l denote the estimated θ r k of the k-th signal source from the l-th Monte Carlo trial.
In the first experiment, we investigate the resolution ability under the overdetermined case with the nested array geometries, respectively.We consider K = 5 spatial uncorrelated sources which are uniformly distributed within −45 °to 45 °.The SNR is set to 0 dB, and 500 snapshots are collected.The nonuniform noise is σ 2 n = 10,5 6,8 5,4 2,9 5,0 5 .In He et al.'s method, we set the grid size is 0.5 °. Figure 3 shows the normalized power of the five targets for the nested array, respectively.We use the dashed lines to represent the true DoAs.It is observed that the proposed method can exactly resolve all targets among all methods, and it has the most superior estimation performance since the DoA estimations are closest to the true target directions.The reason is that the proposed method makes the best of the features of covariance matrix, i.e., low rank, Toeplitz structure, and PSD, and approach the true rank by the iterative weighted idea.The OGVSBL-PVA and the CMR-root-MUSIC are a little worse than the proposed method, so we confirmed that the proposed method achieved higher DoA estimation accuracy.
In the second experiment, we investigate the detectable capability of all methods for the underdetermined case.We consider K = 7 narrowband uncorrelated sources whose number is greater than the number of physical elements, the signal source directions are uniformly distributed within the range from −π/4 to −π/7, and the nonuniform noise power is assumed that σ 2 n = 10,5 6,8 5,4 2,9 5,0 5 ; the SNR is set to 0 dB, and the number of snapshots is 500.Figure 4 displays the normalized power of all methods, respectively.It is observed that the OGVSBL-PVA method and the proposed method can detect all the signal sources exactly.He et al.'s method has a certain deviation; this method uses the denoising operation and covariance matrix fitting to estimate the DoAs, which is based on the on-grid model.From Figure 4, it is easy to observe that the DoA estimation of the proposed method is closest to the true DoAs among all the methods; the reason is that the NA array provides more DOFs than the number of physical elements, and the low-rank covariance matrix approximation is more exact than the covariance matrix reconstruction method by exploiting the iteratively reweighted nuclear norm minimization; moreover, the proposed method is not limited by the grid size like the traditional grid-based sparsity method.Thus, the proposed methods have the superior estimation performance in all methods.
In the third experiment, in order to further illustrate the superior performance of the proposed method, the RMSE of DoA estimation performance versus the angular separation between the source is examined for all methods.In the experiment, the DoA of the first source is fixed and equal to −5.3 °, whereas the DoA of the second source is varied and the varied step size is equal to 0.1 times beam width.When the angle separation is larger than 0.2 times beam width, the performance of the proposed method and CMR-root-MUSIC is achieving the CRLB; since the covariance matrix reconstruction method is based on the low-rank property, it circumvents the traditional grid search, which is the gridless method.The OGSBL-NA method and OGVSBL-PVA method also have the better DoA performance due to the linear approximation manifold vector.As shown in Figure 5, compared with other methods, the performance of the proposed method demonstrates the better performance even when the target is very close.
In the following experiments, unless otherwise specified, the element number in nested array is set to six; seven uncorrelated equal power signals are considered which are uniformly distributed between −π/4 and −π/7.In this case, when grid size is set to 0.5 °, for He et al.'s method, the DoAs do not fall into the predefined angular grid, i.e., grid mismatch.The nonuniform noise power is assumed that σ 2 n = 10,5 6,8 5,4 2,9 5,0 5 .All simulation results were obtained by averaging over 500 independent Monte Carlo trials.
In the fourth experiment, we examine the RMSE of DoA estimation performance with regard to the SNR for different methods.Figure 6 shows the RMSE of DoA estimation performance as a function of SNR when the number of snapshots is 500.Figure 7 shows that the RMSE versus the number of snapshots with the SNR is fixed to 0 dB.As we can see in Figures 6  and 7, the DoA estimation performance of the proposed method is improved with the SNR and increasing snapshot number, since the covariance matrix approximation is more accurate with larger SNR and more snapshots.The proposed method and the CMR-root-MUSIC outperform OGSBL-NA, 9 International Journal of Antennas and Propagation OGVSBL-PVA, and the other methods.CMR-root-MUSIC method and the proposed method provide almost similar estimation performance, while the proposed method is slightly better than CMR-root-MUSIC.Since the proposed method not only exploits the full noiseless covariance matrix but also utilizes the iteratively reweighted strategy to enhance the sparsity of matrix rank resulting in better covariance matrix approximation, we can see that in most of the cases    the part of the noiseless covariance matrix due to the denoising operation and when the DoAs do not fall into the predefined angular grid, the performance will degrade.OGSBL-NA method and OGVSBL-PVA method both improve the DoA performance by linear approximation manifold vector, but the precise depends on the order of the Taylor expansion.It needs to point that the proposed method considerably enhances the sparsity of the covariance matrix due to the iteratively weighted nuclear norm strategy; moreover, the search-free method is free of grid restriction.Thus, we can see that the proposed method has the best performance of DoA estimation accuracy in all tested methods.In the last experiment, we measure the DoA estimation performance in the nonuniform noise environment.We define the worst noise power ratio (WNPR) WNPR = σ 2 nmax /σ 2 nmin [29], where σ 2 nmax and σ 2 nmin denote the maximum and the minimum noise power, respectively.The RMSE of DoA estimation performance for all methods is shown in Figure 8 as a function of the WNPR.The SNR is fixed to 0 dB.We vary the WNPR from 20 to 120 by changing the noise power of the first element from 10 to 60.It is observed that the proposed method performs better than the other methods at all WNPR levels; this can be explained by the fact that using denoising operation mitigates the influence of the nonuniform noise and enhances the robustness of the proposed method; furthermore, by iteratively reweighed idea, the proposed method can exactly approximate the rank of covariance matrix, and search-free DoA method is not subject to the grid constraint.

Conclusion
In this paper, we addressed the direction finding problem in the presence of unknown nonuniform noise over nested array.We proposed a gridless DoA estimation method based on the low-rank covariance matrix approximation.In order to migitate the effect of nonuniform noise, we eliminated the nonuniform noise variables on the diagonal of the covariance matrix by linear transform.In order to recover the noise-free covariance matrix, we proposed an iterative reweighted nonconvex optimization method, whose regularization parameter is determined by the covariance fitting criteria.Finally, we exploited the search-free DoA estimation method to estimate the exact DoA parameters.Numerical simulations were carried out to verify the effectiveness of the proposed method, compared with the state-of-the-art methods, which achieves more accurate DoA estimation performance for the nonuniform noise environment in the presence of grid bias.

Appendix
Proof.Assume h ϵ Toep P u = ∑ i f σ i Toep P u .If function f is twice differentiable, strictly concave, for any bounded k and 0 ≤ x ≤ k, f ′′ ≤ −m k ′ < 0, and then h ϵ Toep P u is strictly concave.For any bounded u, v u ≠ v , there exists some λ > 0 such that With the aid of the result on the h ϵ Toep P u = ∑ i h ϵ σ i Toep P u ≥ 0, we can draw the first property.In order to derive the second property, we prove that the optimal variable u k is convergent.By applying classical Vandermonde decomposition lemma [28] on h ϵ Toep u , we have which implies that lim k⟶∞ Toep P u k+1 − u k = 0. Thus, u k is convergent.In order to prove that u k converges to a local minimum, it is noted that the models (20) and (22) have the same Karush-Kuhn-Tucker (KKT) conditions when the optimal variable u k converges, i.e., u k ⟶ u * as k ⟶ ∞.The reason is that the gradients of the objective functions and the constraints in the rank norm are the same as in nuclear norm models when the MM method 11 International Journal of Antennas and Propagation converges.We conclude that u * tend to be a local minimum of (20) since the objective function in (20) monotonically decreases.This completes the proof.

2 Toep P u k+1 − u k 2 F ≥ λ 2 Toep P u k+1 − u k 2 F A 6 For all k ≥ 0 summing the inequality above, we can attain λ 2 〠Toep P u k+1 − u k 2 F
h ϵ Toep P u k − h ϵ Toep P u k+1 ≥ ∇h ϵ Toep P u k Toep P u k − u k+1 * + λ +∞ k=0 ≤ h ϵ Toep P u 0 , A7 He et al.'s method showcases the worst performance with the largest error; the reason is that He et al.'s method only utilizes