A Sparse-Based Off-Grid DOA Estimation Method for Coprime Arrays

Recently, many sparse-based direction-of-arrival (DOA) estimation methods for coprime arrays have become popular for their excellent detection performance. However, these methods often suffer from grid mismatch problem due to the discretization of the potential angle space, which will cause DOA estimation performance degradation when the target is off-grid. To this end, we proposed a sparse-based off-grid DOA estimation method for coprime arrays in this paper, which includes two parts: coarse estimation process and fine estimation process. In the coarse estimation process, the grid points closest to the true DOAs, named coarse DOAs, are derived by solving an optimization problem, which is constructed according to the statistical property of the vectorized covariance matrix estimation error. Meanwhile, we eliminate the unknown noise variance effectively through a linear transformation. Due to finite snapshots effect, some undesirable correlation terms between signal and noise vectors exist in the sample covariance matrix. In the fine estimation process, we therefore remove the undesirable correlation terms from the sample covariance matrix first, and then utilize a two-step iterative method to update the grid biases. Combining the coarse DOAs with the grid biases, the final DOAs can be obtained. In the end, simulation results verify the effectiveness of the proposed method.


Introduction
Direction of arrival (DOA) estimation of multiple far-field narrowband sources is a vital and interesting topic in a signal processing field, which has wide application in many areas such as radar, sonar, radio astronomy and wireless communications, etc. [1][2][3][4]. As we all know, the DOA estimation problem can be solved efficiently by classical subspace-based methods such as multiple signal classification (MUSIC) [5] and estimation of signal parameters via rotational invariance technique (ESPRIT) [6] in conditions where the receiving arrays are mostly uniform arrays, such as uniform linear array (ULA) [7], uniform rectangular array (URA) [8] and uniform circular array (UCA) [9]. The inter-element spacing of these receiving arrays is required to be less than or equal to the half-wavelength of received signal. Therefore, the detection performance of these DOA estimation methods is greatly confined by the number of physical sensors. For example, a maximum of M − 1 sources can be detected when these algorithms are applied to a ULA with M sensors. Accordingly, more sensors are needed in order to detect more sources and acquire the desirable estimation accuracy. However, a large number of sensors will increase the hardware cost and the difficulty of array calibration in practical applications. To overcome this challenge, some non-uniform array geometries, called sparse arrays, have been proposed, such as minimum redundancy arrays (MRAs) [10], nested arrays (NAs) [11,12] and coprime arrays (CPAs) [13][14][15]. Though MRAs can obtain more degrees of freedom (DOFs) through constructing an augmented covariance matrix, they have these methods are all applied on the traditional uniform linear array and they do not utilize the increased DOFs provided by the difference coarray of coprime arrays. Thus, they cannot detect more sources than sensors. As mentioned before, virtual array aperture extension can be described effectively by the vectorized covariance matrix model. However, it is not rational to directly apply most sparse-based methods to this model as those methods are usually designed to work in the raw data domain. Due to the superiority of coprime arrays, some off-grid DOA estimation methods have been proposed for coprime arrays. In reference [33], an off-grid DOA estimation method based on the framework of sparse Bayesian learning was proposed for coprime arrays. The predefined grid points are directly updated by iteratively decreasing a surrogate function majorizing the given objective function. Sun et al. [34] proposed an iterative method to obtain the final DOA estimation by gradually amending the offset vector with introducing a convex function majorizing the objective function. Both of the methods in references [33,34] need to use the gradient descent method, in which the gradient descent coefficient is difficult to choose. If the gradient descent coefficient is too large, the DOA estimation accuracy will be reduced. On the contrary, if the gradient descent coefficient is too small, the convergence speed of the algorithms will be limited.
To deal with the above difficulties, a sparse-based off-grid DOA estimation method for coprime arrays is proposed in this paper. The proposed method includes two processes, i.e., coarse estimation process and fine estimation process. In the coarse estimation process, we consider the estimation error caused by finite snapshots effect and remove the redundant elements in the covariance vector (i.e., vectorized covariance matrix). In addition, we circumvent estimating the noise variance by using a linear transformation to remove the noise-related part in the covariance vector. This is due to the fact that the noise variance is not so easy to evaluate especially under the underdetermined scenarios. Through a series of linear transformations, the covariance matrix estimation error can be normalized to an identity matrix. Then, according to the statistical property of the vectorized covariance matrix estimation error, DOA estimation can be cast as a problem of recovering a nonnegative sparse vector with an extended steering vector. The grid points corresponding to the non-zero elements in the recovered sparse vector are coarse DOA estimation results for the true signals, which are taken as the predefined grid points required for the fine estimation process. In practical applications, the number of snapshots is finite, which will cause the sample covariance matrix deviating from the actual covariance matrix largely [35]. Therefore, in the fine estimation process, the sample covariance matrix is modified firstly and then the grid biases are derived by applying a two-step iterative technique to the vectorized modified covariance matrix. The final estimated DOAs can be obtained by integrating the coarse DOA estimation results and the grid biases.
The rest of the paper is organized as follows. The far-field narrowband signal model for coprime arrays is briefly introduced in Section 2. Section 3 formulates the proposed sparse DOA estimation method. Section 4 provides several numerical simulation results to show the superior estimation performance for the proposed method and Section 5 concludes this paper.
Notations: In this paper, italic lower-case (upper-case) bold characters are used to denote vectors (matrices). In particular, I K denotes the K × K identity matrix. Sets are denoted by capital letters in blackboard boldface. Specifically, the symbols R, C and Z represent the sets of real numbers, complex numbers and integers, respectively. (·) * , (·) T , (·) H and (·) −1 imply conjugate, transpose, conjugate transpose and inverse, respectively. vec (·) denotes the vectorization operator. The symbols •, and ⊗ represent Hadamard product, Khatri-Rao product and Kronecker product, respectively. diag (a) stands for a diagonal matrix whose diagonal elements are taken from the given vector a, while diag (A) represents a column vector whose elements are the diagonal elements of the given matrix A. (·) means taking the real part of a complex value. The l 1 -norm and l 2 -norm are respectively denoted by · 1 and · 2 . tr (·) is the trace operator. E [·] denotes the statistical expectation operator. CN (a, B) denotes a complex Gaussian distribution with mean vector a and covariance matrix B.

Signal Model
As illustrated in Figure 1, consider the coprime arrays which are composed of two uniform linear subarrays. One includes N sensors with inter-element spacing Mλ/2, while the other owns 2M sensors with inter-element spacing Nλ/2, where λ is the wavelength of the received signal. As described in Section 1, M and N are coprime integers. Without loss of generality, we assume that M is less than N. The two subarrays share the zeroth sensor, which is also the reference sensor. Thus, the total number of sensors in the coprime arrays is 2M + N − 1, which are located at where L is sorted in ascending order with l 1 = 0 and l 2M+N−1 = (2M − 1) N. Suppose K far-field narrowband signals are incident on the coprime arrays from directions (θ 1 , θ 2 , · · · θ K ), where θ k ∈ − π 2 , π 2 , k = 1, 2, · · · , K is the direction of the k-th signal. Then, the received signal at time t can be expressed as In Equation (2), s (t) = [s 1 (t) , s 2 (t) , · · · , s K (t)] T ∈ C K×1 is the signal vector at time t (0 ≤ t ≤ T), which is independent from the signal vectors at other time instances. T is the number of snapshots. n (t) = [n 1 (t) , n 2 (t) , · · · , n 2M+N−1 (t)] T is the noise vector at time t, the elements of which are assumed to be temporally and spatially independent and identically distributed (i.i.d.) random variables which satisfy the complex Gaussian distribution with zero mean and variance σ 2 . When the number of sampled snapshots is sufficiently large, the signal vectors can be assumed to be independent from the noise vectors. A = [a (θ 1 ) , a (θ 2 ) , · · · , a (θ K )] ∈ C (2M+N−1)×K is the array manifold matrix, where a (θ k ) is the steering vector for source k and it can be written as a (θ k ) = e jπl 1 sin θ k , e jπl 2 sin θ k , · · · , e jπl 2M+N−1 sin θ k T . ( Since the number of sampled snapshots is finite in practical applications, the sample covariance matrix can be approximated asR

Proposed Sparse-Based Off-Grid Direction of Arrival Estimation Method
For the sake of solving the problem that the incident signals are not always located on the predefined grid points, an off-grid DOA estimation method is proposed in this section. The method includes two parts: coarse estimation process and fine estimation process. The coarse estimation process is to determine the nearest grid point to the direction of each incident signal, while the fine estimation process is to find the bias between the nearest grid point and the true DOA.

Coarse Estimation Process
One of the tasks for the coarse estimation process is to estimate the positions of the grid points closest to the incident signals. In addition, it also provides the parameters needed for the fine estimation process.
Under the situations where the number of sampled snapshots is sufficiently large, the received data covariance matrix can be denoted as Here, R s =E s (t) s H (t) = diag σ 2 1 , σ 2 2 , · · · , σ 2 K is the theoretical signal covariance matrix and σ 2 k (k = 1, 2, · · · , K) denotes the power of the k-th signal. By vectorizing the covariance matrix R xx , the covariance vector can be obtained as T . e i denotes a column vector whose elements are all zeros except that the i-th element is one. r can be seen as the received data from the virtual array with extended array aperture, where the virtual sensors are located at The sample covariance matrix is as described in Equation (4). Owing to the finite snapshots, the estimation error between the actual and sample covariance matrix needs to be considered. Thus, vectorizing the sample covariance matrix, we can get where ε is the vectorized covariance matrix estimation error and it satisfies an asymptotic Gaussian distribution [36,37] The number of extended virtual sensors in Equation (6) is (2M + N − 1) 2 , while the number of non-redundant virtual sensors is 3MN + M − N. These duplicate virtual sensors only provide unnecessary redundant information. Thus, for the sake of saving computation complexity, we remove the repeated rows ofÃ corresponding to the redundant virtual sensors and sort the remaining rows to construct a new virtual array manifold matrix B. Similarly, the same transformation is performed onr to obtain the corresponding received data vector z, i.e., where e is a column vector whose elements are all zeros except that the element at the center is one. It is obvious thatε is different from ε. The relationship betweenε and ε will be given later.
We defineD as the set, the elements of which are the locations of the virtual sensors corresponding to array manifold matrix B. It can be directly derived by subtracting the duplicates from the difference D and sorting the remaining entries in ascending order and can be denoted as is the virtual steering vector for the k-th signal with the i-th element e jπd i sin θ k ,d i ∈D.
According to [33], we can get the relationship between ε andε: For the convenience of notation expression, let is a selection matrix, the i-th row of which contains all zeros but a single one at the (l is also a selection matrix, the elements of which are all zeros in the i-th row except a single one at the d Let ϑ = (ϑ 1 , ϑ 2 , · · · , ϑ L ) be the grid points sampled from − π 2 , π 2 with appropriate step size, where L is the total number of grid points and L K in general. Then, the virtual received signal model in Equation (10) can be sparsely represented as ×L is the overcomplete dictionary matrix andp is the extension of p with zeros filling in the positions where no signals stay in. As mentioned before, σ 2 e has only one non-zero element at its own center. Therefore, to avoid calculating the noise variance, we remove the corresponding row in the received data vector z by a linear transformation, so that Equation (14) can be transformed to a noiseless received signal model, which can be expressed asz where J ∈ R (3MN+M−N−1)×(3MN+M−N) can be denoted as Combining Equations (9) with (12) and (15), we can get From Equations (15) and (17), we can derive where W = JF −1 R T xx ⊗ R xx JF −1 T T and W −1/2 is the Hermitian square root of W −1 . Then, Here, Asχ 2 (3MN + M − N − 1) denotes asymptotic chi-square distribution with (3MN + M − N − 1) DOFs. Thus, coarse DOA estimation results can be derived by solving the optimization problem: . According to experience, β is usually set as 10 −4 .
After attainingp, the grid points corresponding to the K maximum values ofp are the coarse DOA estimation results, called coarse DOAs, which are denoted as θ coarse = θ (1)

Fine Estimation Process
The following is the fine estimation process where the grid biases can be derived according to the coarse DOAs. First of all, we expand the sample covariance matrixR xx in Equation (4) aŝ Comparing (5) with (21), we can easily find that the expansion ofR xx in (21) contains four terms while R xx in (5) only contains two. The first two terms in (21) can be viewed as the estimates of R xx , which correspond to the signal and noise components, respectively. The last two terms in (21) can be seen as the estimates of the correlation parts between signal and noise vectors. Generally, the signal vector s (t) is assumed to be independent from the noise vector n (t). Therefore, when the number of snapshots is sufficiently large, the last two terms in (21) are equal to zeros in theory. However, in practical applications, the sampled snapshots available are finite so that the undesirable correlation terms between signal and noise vectors exist, which may have significant values. The existence of these correlation terms will cause the estimate of the received data covariance matrix deviating from the actual one and further result in DOA estimation performance degradation. Accordingly, in order to improve the DOA estimation performance, we should estimate and remove these correlation terms first to construct a new modified sample covariance matrix. From the signal model (2), we can derive the estimate of signal vector by minimizing the residual between observations and estimators, i.e., Here,Â = a θ (1)
Using the Least Square (LS) method to solve the optimization problem in (22), we can obtain the estimate of signal vector s (t) asŝ Then, the estimate of noise vector n (t) can be represented aŝ where PÂ =Â ÂÂ H −1Â H is an estimation for the projection matrix of signal component and P ⊥ A = I − PÂ is an estimation for the projection matrix of noise component. Combining Equations (23) with (24), the third term in (21) can be written as The fourth term in (21) is the conjugate transposition of the third one. Thus, the modified sample covariance matrix can be expressed asR where γ is the scale factor with a value between zero and one [35]. The value of γ depends on the estimation accuracy of the undesirable correlation terms. Ideally, γ should be equal to one when there are no errors in the estimates of the correlation terms. However, estimation errors are inevitable. Therefore, γ takes a value close to one when V in (25) has a small estimation error and a small value conversely if the estimate of V is inaccurate. The modified covariance vectorr can be generated by vectorizingR in (26). By extracting the entries in positions corresponding to the setD fromr, we can attain the covariance vectorz without redundant information. Here, the extended virtual array is consistent with that in (11). Generally, no matter how dense the grid points are, the true signals are impossible to always fall on the predefined grid points. Therefore, it is necessary to introduce off-grid problem for improving the DOA estimation performance. The coarse DOAs are the grid points closest to the true signals. Let δ k denote the grid bias for the k-th signal. Then, the grid biases for all K narrowband signals can be expressed as ∆ = diag (δ), where δ = [δ 1 , δ 2 , · · · , δ K ] T . The interval between two adjacent grid points in the coarse estimation process is r g , which determines the range of the grid bias is − r g 2 , r g 2 . By first-order Lagrangian approximation, the steering vector of the kth signal can be approximated as . Consequently, the off-grid model can be denoted as ×K is the derivative ofÃ with respect to θ coarse .p ∈ R K×1 is a nonnegative signal power vector.
Based on (27), a two-step iterative method is proposed for solving the grid biases. The first step is to fix the grid bias vector δ and solve the optimization problem below: The unconstrained optimization problem [38] of (28) can be expressed as where η 0 is a regularization parameter that can be adjustable to trade off the sparsity ofp and the least-squares error between the observations and estimators. The superscript (j + 1) represents the (j + 1)-th iteration. The optimization problem in (29) can be solved by CVX toolbox [39]. Then, turning to the second step, fix the obtained nonnegative vectorp (j+1) and update the grid bias vector δ by To solve the optimization problem about δ, we first expand (30) as Removing the terms that are not related to δ in (31) Let the derivative of Equation (33) with respect to δ be equal to zero. Then, we can get the updated grid bias vector δ (j+1) as As described before, the grid bias δ k (k = 1, 2, · · · , K) is constrained in − r g 2 , r g 2 , so we have After the iteration converges, the final estimate of the grid bias vector is obtained and can be denoted as δ f inal . Consequently, the final DOA estimation results can be expressed aŝ

Summary of the Steps about the Proposed Method
Synthesizing the foregoing analysis, the steps of the proposed method can be summarized as follows: Step 1: Compute the sample covariance matrixR xx according to Equation (4); Step 2: By vectorizingR xx , getr according to Equation (8). Then, remove the repeated entries ofr and sort the remaining entries to get z as Equation (10); Step 3: Represent z sparsely according to Equation (14) and get the noiseless received signal model according to Equations (15) and (16);

and (35);
Step 9: Terminate the iterative process if convergence criteria δ (j+1) − δ (j) δ (j) ≤ τ is satisfied or the number of iterations exceeds the maximum one. Otherwise, return to Step 7 and continue the iteration process. Output: The final grid bias vector δ f inal and DOA estimation resultsθ = θ coarse + δ f inal .

Numerical Simulations
In this section, we perform several simulation experiments to demonstrate the excellent performance of the proposed method for DOA estimation, and also compare it with SS-MUSIC [16], low rank matrix denoising (LRD) [17], LASSO [23] and off-grid sparse Bayesian inference (OGSBI) [28]. In particular, for estimation accuracy comparision, the Cramer-Rao bound (CRB) for coprime arrays [40] is included. The coprime arrays considered here are composed of ten sensors with two coprime integers M = 3 and N = 5. The positions of sensors in the first subarray are [0, 3,6,9,12]

Detection Performance
In this subsection, the detection performance of the proposed method is compared with that of the methods OGSBI, LASSO and LRD. All of these methods can detect more sources than sensors with using coprime arrays. In theory, LRD can only identify at most 17 sources under the coprime arrays set above. Therefore, for convenience of comparison, we consider 17 randomly distributed narrowband uncorrelated sources here, which are located at −50.88 From Figure 2, we can see that OGSBI, LASSO and the proposed method can detect all 17 sources successfully. However, there are still some small spurious peaks in Figure 2a for OGSBI and Figure 2c for LASSO. Furthermore, from Figure 2b, it is easy to find that LRD fail to detect all sources. In its spatial spectrum, several peaks are missed and some of the existing spectral peaks deviate from the true signals a lot. Therefore, it can be concluded that the detection performance of the proposed method is superior to the other three algorithms.

Resolution Ability
In this subsection, the resolution ability of the proposed method is still compared with that of the three algorithms mentioned in Section 4.1. We test the resolution abilities of these algorithms by detecting two closely spaced sources which are located at −31.25 • and −29.64 • , respectively. Here, the SNR is set as 0 dB and the number of sampled snapshots is set as 800. Additionally, the scale factor γ is identical with that in Section 4.1.
From Figure 3, it is obvious to find that there is only one peak in the spatial spectrums for LRD and LASSO, which means that both of them can not identify the two signals successfully. In addition, from Figure 3d, we can find that the proposed method identifies the two signals successfully and its estimated DOAs about the two signals are −31.05 • and −29.3 • , respectively. Though there are errors between the estimated DOAs and the true DOAs, the errors are very small and within an acceptable range. From Figure 3a, it can be seen that OGSBI can also identify the two signals, but the estimation error for OGSBI is larger than the proposed method. Therefore, it can be concluded that the resolution ability of the proposed method is the best among all four of the algorithms.

Estimation Accuracy
In this subsection, the estimation accuracy of each algorithm is reflected by its root-mean-square error (RMSE), which is defined as where Q is the total number of Monte Carlo trials,θ k,q represents the estimate of the k-th source in the q-th trial. For each scenario, we conduct 1000 Monte Carlo trials to get the corresponding RMSE. In the first experiment of this subsection, we fix the number of snapshots at 800 and vary SNR from −10 dB to 20 dB. When SNR is less than or equal to 2 dB, γ is set as 0.5. Otherwise, γ is set as 0.7. Figure 4 shows the RMSE curves of the five algorithms and CRB as a function of SNR, respectively. From Figure 4, we can see that the estimation accuracy is improved with the increase of SNR for all algorithms. It is obvious that the estimation performance of the proposed method outperforms the other four methods. The reason is that spatial smoothing technique causes the extended DOF loss for SS-MUSIC and LRD. Though LASSO can utilize all of the extended DOFs, the assumption that incident signals must fall on the predefined grid points causes its estimation performance degradation. The algorithm OGSBI needs to update the noise variance in each iteration, which makes it sensitive to the noise. In addition, the covariance matrix estimation error and the correlation terms between signal and noise vectors are not considered in OGSBI, which is another reason that the proposed method has better estimation performance than OGSBI.
For the second experiment, the SNR is fixed at 0dB and RMSE is a function of snapshots, the number of which varies from 10 to 1500. When the number of snapshots is less than or equal to 800, γ is set as 0.5. Otherwise, γ is set as 0.8. The RMSE values for all algorithms and CRB versus snapshots are shown in Figure 5.
It is obvious from Figure 5 that the estimation accuracy is improved with the increase of snapshots for all algorithms. For the same reason as described above, the estimation performance of SS-MUSIC, LRD, LASSO and OGSBI is inferior to the proposed method. It can be seen from Figure 5 that when the number of snapshots changes from 10 to 1500, the estimation accuracy of the proposed method is getting better and better than the other four algorithms and finally outperforms the other four algorithms notably by a large margin.
In addition, to have an intuitive understanding on the computation complexity, we compute the average execution time of each algorithm on an Intel Core i7-7700@3.6GHz, 16G RAM PC. Under the same conditions, it takes 0.006 s for SS-MUSIC, 0.670 s for LRD, 0.285 s for LASSO, 1.343 s for OGSBI and 1.599 s for the proposed method. It can be found that SS-MUSIC has the shortest running time. The algorithm OGSBI and the proposed method are a little time-consuming. However, it is worth noting that the simulation results show that the estimation performance of the proposed method is the best among all five of the algorithms. According to the characteristics of each algorithm, they can be applied in different scenarios, respectively. When applied in scenarios where the real-time demand is not very high but the estimation accuracy demand is high, the proposed method is the preferred choice. Conversely, the other several comparison algorithms can be used in the scenarios with high real-time requirements but low estimation accuracy requirements. As we all know, the trade-off between the estimation accuracy and the computation complexity has always been an intractable problem for nearly all algorithms. Therefore, reducing the computation complexity on the basis of the proposed method will be a research direction for us in the future.

Conclusions
In this paper, we have proposed a sparse-based off-grid DOA estimation method for coprime arrays. By using coprime arrays, we can obtain extended virtual array aperture and increased DOFs. Accordingly, the proposed method can detect more sources than sensors. In practical applications, the existence of the estimation error between the sample and actual covariance matrix is inevitable due to finite snapshots. According to the statistical property of the vectorized covariance matrix estimation error, we construct an optimization problem about the signal power vector in the coarse estimation process so as to obtain the positions of the grid points closest to the true DOAs. In the fine estimation process, we first remove the undesirable correlation terms between the signal and noise vectors to get the modified sample covariance matrix. Then, for solving the grid mismatch problem, an off-grid model with vectorized modified sample covariance matrix is constructed. By applying the two-step iterative method to the off-grid model, the grid biases can be estimated precisely. Combining the coarse estimation results with the estimated grid biases, the final DOA estimation results can be accurately acquired. By considering the covariance matrix estimation error and the correlation terms between signal and noise vectors, the proposed method can achieve better DOA estimation performance than those traditional algorithms that do not consider the finite snapshots effect. In addition, the introduction of an off-grid model solves the grid mismatch problem existing in the traditional on-grid model successfully. Numerical simulation experiments confirmed the superiority and effectiveness of the proposed method.