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

Authors: Weijie Tan; Xi’an Feng
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

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.


---

## Body

## 1. 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–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 oMN 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–11] and the Toeplitz completion method [12–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–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 Section2, 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 vectorx 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, σiX denote the i-th largest eigenvalue value, Ik denotes a K×K dimensional identity matrix, ·2 denotes the l2 norm, and Pr· denotes probability density function.


## 2. Array Model for Nested Array

We consider a two-level nested array consisting of two ULA subarrays, where the first ULA subarray hasN1 elements with the interelement spacing d1=d, while the second ULA subarray has N2 elements with the interelement spacing N1+1d1; the total number of element is M=N1+N2, as shown in Figure 1.Figure 1
The nested array structure.Let the first element of the array as the phase reference point, the position set of element in ascending is denoted as(1)P=p1,⋯,pM.Assume thatK uncorrelated narrow-band far-field sources from θ=θ1,⋯,θKT impinge on the nested array. The array output at time t can be represented as
(2)xt=Aθst+nt,where xt∈ℂM×1 is the output of nested array. st=s1t,s2t,⋯,sktT is the K×1 receive signal source vector, and nt∈ℂM×1 denotes the additive white Gaussian noise vector. The M×K array manifold matrix in array model (2) is defined by
(3)A=aθ1,aθ2,⋯,aθK,where aθk=ejπp1sinθk,ejπp2/sinθk,⋯,ejπpM/sinθk is the array steering vector and the element position is defined in (1). We note that each position pi 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 setP 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 aP~θk of the ULA can be given as aP~θk=ej2πP~1d/λsinθk,⋯,ej2πP~M~d/λsinθk.According to the index setP of the NA, we design a selection matrix Γ∈0,1M×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=ΓaP~θk. Then, the nested array model can be expressed as
(4)xt=Aθkst+nt=ΓAP~st+ΓnP~t.Under ideal condition, the array covariance matrix of the outputxt of nested array can be calculated as
(5)R=limN⟶∞1N∑n=1NxtnxHtl=ARsAH+diagσn2=ΓAP~RsAP~HΓT+Γdiagσn,P~2ΓT=ΓToepP~uΓT+Γdiagσn,P~2ΓT,where N is the number of collected snapshot, Rs=diagσs,12,⋯,σs,K2, and ToepP~u=AP~RsAP~H for some uP~∈ℂM~×1. Furthermore, the noise-free covariance matrix ToepP~u is a Hermitian Toeplitz matrix, which is determined by M complex numbers.
(6)ToepP~u=u1u2⋯uM~u2∗u1⋯uM~−1⋮⋮⋱⋮uM~∗uM~−1∗⋯u1.It is evident thatrankToepP~u=K≤M~−1, and ToepP~u is positive definite, i.e., ToepP~u⪰0. For nested array, Toepu holds that
(7)Toepum,n=∑k=1Kσs,k2amθkan∗θk=∑k=1Kσs,k2ej2πd/λm−nsinθk.LetD=m−n:m,n∈P,m≥n⊂P~. For the NA, D is a redundancy array which can be regarded as a coarray defines on the NA. The matrix Toepu contains all elements uP~ of explicitly by (6). It implies that the relation ToepP~u⟶Toepu is one-to-one.A coarray of nested array with six elements is shown in Figure2. When the positions of nested array are P=0,1,2,3,6,10, the coarray of nested array composes of M~+1=11 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].Figure 2
Weight coefficient of NA coarray with six elements.One of the aims of this paper is to recover the noise-free covariance matrixToepP~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 ToepP~u.


## 3. 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 matrixToepP~u. On one hand, we note that the covariance matrix ToepP~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=1NxtnxHtn/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:
(8)minu,σrankToepP~us.t.C^−1/2vecΔR22≤ηToepu⪰0,where C^=1/NRT⊗R represents the prewhitening weighted matrix, ΔR=R^−R is the perturbation error term, and η is an upper bound of C^−1/2vecΔR22, 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:
(9)vecΔR~AsN0M2,1,1NRT⊗R.After prewhitening processing, we obtain(10)C^−1/2vecΔR~AsN0,IM2.We can further derive the following formulation:(11)C^−1/2vecΔR22~Asχ2M2,where Asχ2M2 represents the asymptotic chi-square distribution with M2 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 C^−1/2vecΔR22≤η is satisfied with a high probability 1−p, e.g., 0.999, PrC^−1/2vecΔR22≤η=1−p, which is calculated by the matlab function chi2inv1−p,M2, so we have PrχM22≤η=p and η=χp2M2.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 ToepP~u to the nuclear norm of ToepP~u for approximating the low-rank matrix. Therefore, the convex relaxation form for (8) can be represented as
(12)minu,σToepP~u∗s.t.C^−1/2vecΔR22≤ηToepu⪰0.Combining (5) and (12), the optimization problem is rewritten as
(13)minu,σToepP~u∗s.t.C^−1/2vecR^−ΓToepP~u+diagσn2ΓT22≤ηToepP~u⪰0.Remark1.
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 hasM nonzero elements which contain in the diagonal element of the covariance matrix and correspond to the unknown noise variance σn2. 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
(14)JvecΔR=JvecR^−ΓToepP~uΓH+Γdiagσn2ΓH=JvecR^−JvecΓToepP~uΓH,where J is an MM−1×M2 linear transform selecting matrix and can be represented as
(15)J=J1,J2,⋯,JM−1T,where
(16)Jm=em−1M−1+2,em−1M−1+3,⋯,emM−1,m=1,2,…,M−1,with eii=m−1M−1+2,⋯,mM−1 being an M2×1 column vector with one at the ith position and zeros elsewhere. Jm∈ℝM2×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
(17)minuToepP~u∗s.t.C~−1/2JvecR^−JvecΓToepP~uΓH22≤η^ToepP~u⪰0.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
(18)JvecΔR~AsN0MM−1,1,1NJRT⊗RJT.Applying the prewhitening processing [19], we have
(19)C~−1/2JvecΔR~AsN0MM−1,1,IMM−1,where the prewhitening matrix C~−1/2 is Hermitian square root of C~−1, C~=JRT⊗RJT/N, and IMM−1 is an MM−1×MM−1 identity matrix. Like the choosing parameter η^, the parameter η^ can be determined by η^=χp2MM−1.


### 3.2. Iteratively Reweighted Nuclear Norm Minimization

We utilize the nuclear normToepP~u∗=∑iσiToepP~u=∑ig1σiToepP~u in (12) instead of the rank norm rankToepP~u=σToepP~u0=∑ig0σiToepP~u in (8), where g0x=signx, x≥0, g1x=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 findfx=1−e−x/ϵ closely matches g0x, so we use this function hϵx=fx=1−e−x/ϵ as the approximation function of g0x, i.e., rank norm. Combine the model (17), the following general nonconvex optimization model is proposed:
(20)minuhϵToepP~us.t.C~−1/2JvecR^−JvecΓToepP~uΓH22≤η^ToepP~u⪰0,where hϵσToepP~u=∑ifσiToepP~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 uk−1. Since hϵToepP~u=∑ifσiToepP~u is concave when fx=g1x, we have
(21)hϵToepP~u≤hϵToepP~uk+∇hϵToepP~uToepP~u−uk∗.Therefore, the optimization problem for thek-th iteration becomes
(22)minuWkϵToepP~u∗s.t.C~−1/2JvecR^−JvecΓToepP~uΓH22≤η^ToepP~u⪰0,where Wkϵ≜∇hϵToepP~uk. In order to calculate Wkϵ, we start by introducing a key useful proposition as follows.Proposition 1.
Suppose thathϵσZ is expressed as hϵσZ=∑ifσiZ, where Z is positive semidefinite matrices with the eigenvalue decomposition (EVD) Z=UdiagσZUH, Z⪰0, and hϵ and f are concave and differentiable, then the gradient of fZ at Z is calculated as
(23)∇hϵZ=UdiagηZUH,where ηZ=∇hϵσZ denotes the gradient of hϵ at σZ.In (22), Z=ToepP~uk, for the sake of simplicity, we denote σToepP~uk as σk and apply the result in [27]; we obtain
(24)Wkϵ=Ukdiag∇hϵσkUkH,where UkdiagσkUkH is the EVD of ToepP~uk.We observe thatϵ is a key parameter, which controls the closeness between fσ and g0σ or g1σ. In particular, fσ is close to g0σ when ϵ is small, but it is subjected to many local minima. On the other hand, a large value of ϵ causes fσ toward the convex g1σ, whereas with poor approximation quality to g0σ. 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 Wkϵ 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.
Denoteuk as the k-th iteration solution in (22), after that, the optimal variable uk satisfies the following properties:
(i)
hϵToepP~uk is monotonically decreasing and bounded as k tends to ∞
(ii)
The optimal variableuk 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 solutionu of model (20) is obtained, we can determine the number of sources K^ by the rank of ToepP~u. Given ToepP~u^, the remaining task is to estimate the directions of the signal θ^. According to classical Vandermonde decomposition lemma [28], we can estimate the DoAs using positive semidefinite Toeplitz matrices ToepP~u. So we can apply the grid-free DoA estimation method, such as Prony’s method to efficiently estimate θ^ and σ^s2. Firstly, by solving a linear system involving only u^, we can obtain θ^ from zeros of a polynomial. Then, σ^s2 is solved by a polynomial rooting-based estimation method.Algorithm 1: Gridless DoA estimation using iteratively reweighted covariance matrix reconstruction.
1:Input: X∈ℂM×K.
2:Initialize: k=0, u0=0M×1, ϵ, γ, ε=10−4, Itermax
3: Calculate the sample covariance matrixR^S of NA using
the observe snapshotsxΩt;
4:while k≤Itermax do
5:  Update the weight matrixWkϵ by the equation
6:  Updateuk+1 by the optimization problem (22), (24);
7:if the convergence condition uk+1−uk≤ε is satisfied
then
8:   break
9:else
10:ϵ=ϵ/γ, k=k+1;
11:end if
12:end while
13: Calculate the covariance matrixR^=Toepu^;
14: Using search-free DoA estimation method to estimate the
angle parameters.In the proposed method, firstly, it recovers the complete noiseless matrixToepP~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.


## 3.1. Nuclear Norm Minimization Optimization

As mentioned above, the key problem is recovery of the noise-free covariance matrixToepP~u. On one hand, we note that the covariance matrix ToepP~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=1NxtnxHtn/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:
(8)minu,σrankToepP~us.t.C^−1/2vecΔR22≤ηToepu⪰0,where C^=1/NRT⊗R represents the prewhitening weighted matrix, ΔR=R^−R is the perturbation error term, and η is an upper bound of C^−1/2vecΔR22, 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:
(9)vecΔR~AsN0M2,1,1NRT⊗R.After prewhitening processing, we obtain(10)C^−1/2vecΔR~AsN0,IM2.We can further derive the following formulation:(11)C^−1/2vecΔR22~Asχ2M2,where Asχ2M2 represents the asymptotic chi-square distribution with M2 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 C^−1/2vecΔR22≤η is satisfied with a high probability 1−p, e.g., 0.999, PrC^−1/2vecΔR22≤η=1−p, which is calculated by the matlab function chi2inv1−p,M2, so we have PrχM22≤η=p and η=χp2M2.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 ToepP~u to the nuclear norm of ToepP~u for approximating the low-rank matrix. Therefore, the convex relaxation form for (8) can be represented as
(12)minu,σToepP~u∗s.t.C^−1/2vecΔR22≤ηToepu⪰0.Combining (5) and (12), the optimization problem is rewritten as
(13)minu,σToepP~u∗s.t.C^−1/2vecR^−ΓToepP~u+diagσn2ΓT22≤ηToepP~u⪰0.Remark1.
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 hasM nonzero elements which contain in the diagonal element of the covariance matrix and correspond to the unknown noise variance σn2. 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
(14)JvecΔR=JvecR^−ΓToepP~uΓH+Γdiagσn2ΓH=JvecR^−JvecΓToepP~uΓH,where J is an MM−1×M2 linear transform selecting matrix and can be represented as
(15)J=J1,J2,⋯,JM−1T,where
(16)Jm=em−1M−1+2,em−1M−1+3,⋯,emM−1,m=1,2,…,M−1,with eii=m−1M−1+2,⋯,mM−1 being an M2×1 column vector with one at the ith position and zeros elsewhere. Jm∈ℝM2×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
(17)minuToepP~u∗s.t.C~−1/2JvecR^−JvecΓToepP~uΓH22≤η^ToepP~u⪰0.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
(18)JvecΔR~AsN0MM−1,1,1NJRT⊗RJT.Applying the prewhitening processing [19], we have
(19)C~−1/2JvecΔR~AsN0MM−1,1,IMM−1,where the prewhitening matrix C~−1/2 is Hermitian square root of C~−1, C~=JRT⊗RJT/N, and IMM−1 is an MM−1×MM−1 identity matrix. Like the choosing parameter η^, the parameter η^ can be determined by η^=χp2MM−1.


## 3.2. Iteratively Reweighted Nuclear Norm Minimization

We utilize the nuclear normToepP~u∗=∑iσiToepP~u=∑ig1σiToepP~u in (12) instead of the rank norm rankToepP~u=σToepP~u0=∑ig0σiToepP~u in (8), where g0x=signx, x≥0, g1x=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 findfx=1−e−x/ϵ closely matches g0x, so we use this function hϵx=fx=1−e−x/ϵ as the approximation function of g0x, i.e., rank norm. Combine the model (17), the following general nonconvex optimization model is proposed:
(20)minuhϵToepP~us.t.C~−1/2JvecR^−JvecΓToepP~uΓH22≤η^ToepP~u⪰0,where hϵσToepP~u=∑ifσiToepP~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 uk−1. Since hϵToepP~u=∑ifσiToepP~u is concave when fx=g1x, we have
(21)hϵToepP~u≤hϵToepP~uk+∇hϵToepP~uToepP~u−uk∗.Therefore, the optimization problem for thek-th iteration becomes
(22)minuWkϵToepP~u∗s.t.C~−1/2JvecR^−JvecΓToepP~uΓH22≤η^ToepP~u⪰0,where Wkϵ≜∇hϵToepP~uk. In order to calculate Wkϵ, we start by introducing a key useful proposition as follows.Proposition 1.
Suppose thathϵσZ is expressed as hϵσZ=∑ifσiZ, where Z is positive semidefinite matrices with the eigenvalue decomposition (EVD) Z=UdiagσZUH, Z⪰0, and hϵ and f are concave and differentiable, then the gradient of fZ at Z is calculated as
(23)∇hϵZ=UdiagηZUH,where ηZ=∇hϵσZ denotes the gradient of hϵ at σZ.In (22), Z=ToepP~uk, for the sake of simplicity, we denote σToepP~uk as σk and apply the result in [27]; we obtain
(24)Wkϵ=Ukdiag∇hϵσkUkH,where UkdiagσkUkH is the EVD of ToepP~uk.We observe thatϵ is a key parameter, which controls the closeness between fσ and g0σ or g1σ. In particular, fσ is close to g0σ when ϵ is small, but it is subjected to many local minima. On the other hand, a large value of ϵ causes fσ toward the convex g1σ, whereas with poor approximation quality to g0σ. 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 Wkϵ 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.
Denoteuk as the k-th iteration solution in (22), after that, the optimal variable uk satisfies the following properties:
(i)
hϵToepP~uk is monotonically decreasing and bounded as k tends to ∞
(ii)
The optimal variableuk 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 solutionu of model (20) is obtained, we can determine the number of sources K^ by the rank of ToepP~u. Given ToepP~u^, the remaining task is to estimate the directions of the signal θ^. According to classical Vandermonde decomposition lemma [28], we can estimate the DoAs using positive semidefinite Toeplitz matrices ToepP~u. So we can apply the grid-free DoA estimation method, such as Prony’s method to efficiently estimate θ^ and σ^s2. Firstly, by solving a linear system involving only u^, we can obtain θ^ from zeros of a polynomial. Then, σ^s2 is solved by a polynomial rooting-based estimation method.Algorithm 1: Gridless DoA estimation using iteratively reweighted covariance matrix reconstruction.
1:Input: X∈ℂM×K.
2:Initialize: k=0, u0=0M×1, ϵ, γ, ε=10−4, Itermax
3: Calculate the sample covariance matrixR^S of NA using
the observe snapshotsxΩt;
4:while k≤Itermax do
5:  Update the weight matrixWkϵ by the equation
6:  Updateuk+1 by the optimization problem (22), (24);
7:if the convergence condition uk+1−uk≤ε is satisfied
then
8:   break
9:else
10:ϵ=ϵ/γ, k=k+1;
11:end if
12:end while
13: Calculate the covariance matrixR^=Toepu^;
14: Using search-free DoA estimation method to estimate the
angle parameters.In the proposed method, firstly, it recovers the complete noiseless matrixToepP~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.


## 4. 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 d1=d; the second ULA subarray contains three elements with spacing d2=4d1, 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 Rn=diagσn2, since the input signal-to-noise ratio (SNR) can be defined as
(25)SNR=10log101M∑m=1Mσs2σn,m2,where M=N1+N2 is the number of NA elements, σs2 denotes the signal source power, and σn,m2 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)RMSE=1LK∑k=1K∑r=1Lθkr−θ^k,l21/2,where θ^k,l denote the estimated θkr 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 considerK=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 σn2=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.Figure 3
Normalized spatial spectra in nonuniform noise case forM=6, K=5, N=500, and SNR=0dB.In the second experiment, we investigate the detectable capability of all methods for the underdetermined case. We considerK=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 σn2=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.Figure 4
Normalized spatial power in nonuniform noise case forM=6, K=7, N=500, and SNR=0dB.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 Figure5, compared with other methods, the performance of the proposed method demonstrates the better performance even when the target is very close.Figure 5
RMSEs of DoA estimation versus angle separation forM=6, N=500, and SNR=0dB.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 σn2=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. Figure6 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, 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 He et al.’s method showcases the worst performance with the largest error; the reason is that He et al.’s method only utilizes 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.Figure 6
RMSEs of DoA estimation versus SNR forM=6, K=7, and N=500.Figure 7
RMSEs of DoA estimation versus snapshot number forM=6, K=7, and SNR = 0 dB.In the last experiment, we measure the DoA estimation performance in the nonuniform noise environment. We define the worst noise power ratio (WNPR)WNPR=σn2max/σn2min [29], where σn2max and σn2min 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.Figure 8
RMSEs of DoA estimation versus WNPR forM=6, N=500, and SNR=0dB.


## 5. 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.



---

*Source: 7657898-2019-03-18.xml*