Alternating projections gridless covariance-based estimation for DOA

Direction-of-arrival (DOA) estimation is localizing several sources arriving at an array of sensors. It is an important problem in a wide range of applications, including radar, sonar, etc. Compressive sensing (CS) based DOA estimation, which promotes sparse solutions, has advantages over traditional DOA estimation methods. [1, 2] DOAs exist on a continuous angular domain, and gridless CS can be employed. [3, 4, 5] We propose a gridless sparsity-promoting DOA estimation method and apply it to co-prime arrays, which can resolve more sources than the number of sensors.

CS-based DOA estimation exploits the framework of CS, which promotes sparse solutions, for DOA estimation and has a high-resolution capability, deals with single-snapshot data, and performs well with coherent arrivals. [2, 3] To estimate DOAs in a continuous angular domain, non-linear estimation of the DOAs is linearized by using a discretized angular-search grid of potential DOAs (“steering vectors”). Grid-based CS has a basis mismatch problem [6, 7] when true DOAs do not fall on the angular-search grid. To overcome the basis mismatch, gridless CS [8, 9] has been utilized for DOA estimation. [3, 4, 5, 10, 11]

Gridless SPICE (GLS) [5, 12], one of the off-grid sparse methods, is a gridless version of sparse iterative covariance-based estimation (SPICE) [13]. GLS re-parameterizes the data covariance, or sample covariance matrix (SCM), using a positive semi-definite (PSD) Toeplitz matrix, and finds the lowest rank Toeplitz matrix which fits the SCM. The Toeplitz-structured SCM is related to a Fourier-series, which is composed of harmonics. [14, 15] GLS-based DOA estimation retrieves DOA-dependent harmonics from the SCM parameter. [5] The GLS solver uses a semi-definite programming problem (SDP), which is infeasible in practice for high-dimensional problems.

Alternating projections (AP) algorithm [16, 17] has been introduced to solve matrix completion [18, 19, 20, 21] and structured low rank matrix recovery [22, 23] and consists of projecting a matrix onto the intersection of a linear subspace and a nonconvex manifold. Atomic norm minimization (ANM) [6, 8] solves gridless CS and is equivalent to a recovery of a Toeplitz-structured low rank matrix [24]. AP based on ANM has been applied to gridless CS for DOA estimation. [25, 26]

We propose AP-based GLS for gridless CS for DOA estimation. GLS reconstructs a DOA-dependent SCM matrix, which is a Toeplitz-structured low rank matrix and has a PSD matrix in its constraint. AP-GLS solves the reconstruction of the Toeplitz-structured low rank matrix by using a sequence of projections onto the following sets: Toeplitz set, rank-constraint set, and PSD set.

Co-prime arrays are introduced for DOA estimation and offer the capability of identifying more sources than the number of sensors. [27] Sparse Bayesian learning (SBL) deals with co-prime arrays without constructing a co-array based covariance matrix and shows accurate DOAs identifying more sources than the number of sensors. [28, 29] We apply AP-GLS to co-prime arrays and show that AP-GLS with co-prime arrays estimates more DOAs than the number of sensors.

We study the performance of AP-GLS with co-prime arrays for single- and multiple-snapshot data, incoherent and coherent sources, and when the number of sources exceeds the number of sensors.

Consider the ULA case and assume incoherent sources. (GLS is robust to source correlations. [5, 12, 13]) In the noiseless case, the SCM $\mathbf{R}^{\ast}\in\mathbb{C}^{M\times M}$ is given by

where $\mathbf{Y}^{\ast}$ is noise-free data and $p_{k}>0$, $k=1,\dotsc,K$ is the power of sources, i.e., $p_{k}=c_{k}^{2}$. The SCM $\mathbf{R}^{\ast}$ is a (Hermitian) Toeplitz matrix,

where $\mathbf{r}\in\mathbb{C}^{M}$. Moreover, ${\mathbf{R}^{\ast}}$ is PSD and has rank $K$. A PSD Toeplitz matrix of rank $K<M$ can be uniquely decomposed (Vandermonde decomposition) [5, 6, 30] as

where $\mathbf{A}=[\mathbf{a}(\theta_{1})\ \dotsc\ \mathbf{a}(\theta_{K})]\in\mathbb{C}^{M\times K}$.

GLS uses a SCM-related parameter $\mathbf{R}\in\mathbb{C}^{M\times M}$, which is a rank-$K$ PSD Toeplitz matrix, and fits the parameter $\mathbf{R}$ to SCM $\tilde{\mathbf{R}}=\mathbf{Y}\mathbf{Y}^{\mathsf{H}}/L\in\mathbb{C}^{M\times M}$. The covariance fitting is implemented, in the case of $L\geq M$ whenever $\tilde{\mathbf{R}}$ is non-singular, by minimizing the criterion, [5, 12, 13]

In the case of $L<M$, when $\tilde{\mathbf{R}}$ is singular, the following criterion is used instead, [5, 12]

GLS is achieved using the following optimization,

where $\mathbf{R}\succeq 0$ denotes $\mathbf{R}$ is a PSD matrix and $\mathbf{Z}\in\mathbb{C}^{M\times M}$ is a free variable. Consider the case of $\mathbf{R}=\mathbf{R}^{\ast}$, then $\mathrm{tr}(\mathbf{R})=M\sum_{k=1}^{K}p_{k}$. Defining $\mathrm{tr}(\mathbf{Z})=M\sum_{k=1}^{K}p_{k}$, the objective in (8), divided by $2M$, equals,

Note that, in ANM, [6, 8] minimizing $\sum_{k=1}^{K}p_{k}=\sum_{k=1}^{K}c_{k}^{2}$ is equivalent to minimizing the atomic norm,

The atomic norm is a convex relaxation of the atomic $l_{0}$ norm, [6]

Minimizing the atomic $l_{0}$ norm is equivalent to minimizing rank of $\mathbf{R}^{\ast}=\mathbf{Y}^{\ast}\mathbf{Y}^{\ast\mathsf{H}}/L$. [5, 6] Summarizing, the term $\mathrm{tr}(\mathbf{R})=\sum_{k=1}^{K}p_{k}$ is the nuclear norm, used as a convex relaxation of $\mathrm{rank}(\mathbf{R})$.

By using the rank minimization in (8), the resulting optimization is as follows,

For the coprime array, we use the row-selection matrix $\mathbf{\Gamma}_{\Omega}\in\{0,1\}^{M\times M_{\Omega}}$, i.e.,

where $\mathbf{Y}$ is data of full $M$-element ULA and the Moore-Penrose pseudo-inverse $\mathbf{\Gamma}_{\Omega}^{\dagger}$. The optimization for the coprime array is given as,

where $\tilde{\mathbf{R}}_{\Omega}=\mathbf{Y}_{\Omega}\mathbf{Y}_{\Omega}^{\mathsf{H}}/L\in\mathbb{C}^{M_{\Omega}\times M_{\Omega}}$ and $\mathbf{R}_{\Omega}=\mathbf{\Gamma}_{\Omega}\mathbf{R}\mathbf{\Gamma}_{\Omega}^{\mathsf{T}}\in\mathbb{C}^{M_{\Omega}\times M_{\Omega}}$. To minimize $\mathrm{rank}(\mathbf{R})$, $\mathbf{R}$ is calculated,

We suggest alternating projections to reconstruct Toeplitz-structured low rank matrix in (12) and (14). AP-GLS involves the following sets: Toeplitz set, positive semi-definite (PSD) set, and rank-constraint set.

We consider a ULA with $M=16$ elements, half-wavelength spacing, and a co-prime array with $M=8$ elements with $M_{1}=5$ and $M_{2}=2$ (elements 1, 3, 5, 6, 7, 9, 11, 16).

The signal-to-noise ratio (SNR) is defined, $\text{SNR}=10$ $\log_{10}[\mathbb{E}\{\lVert\mathbf{A}\mathbf{s}_{l}\rVert\}_{2}^{2}/\mathbb{E}\{\lVert\mathbf{e}_{l}\rVert\}_{2}^{2}]$, where $\mathbf{s}_{l}\in\mathbb{C}^{K}$ and $\mathbf{e}_{l}\in\mathbb{C}^{M}$, $l=1,\dotsc,L$, are the source amplitude and the measurement noise for the $l$th snapshot. The root mean squared error (RMSE) is, $\text{RMSE}=\sqrt{\mathbb{E}\left[\frac{1}{K}\sum_{k=1}^{K}\left(\hat{\theta}_{k}-{\theta}_{k}\right)^{2}\right]}$, where $\hat{\theta}_{k}$ and ${\theta}_{k}$ represent estimated and true DOA of the $k$th source.

We consider an $M=8$ co-prime array and $K=4$ stationary sources at DOAs $[-60,-3,3,50]^{\circ}$ with snapshot-varying magnitudes $[12,20]$ dB in four scenarios, see Fig. 1. Conventional beamforming (CBF), SBL [28, 29], and AP-GLS are compared. CBF cannot distinguish close two DOAs $[-3,3]^{\circ}$. AP-GLS solves the single snapshot case and resolve the close arrivals. We also consider the DOA performance with a coherent sources due to multipath arrivals. AP-GLS still shows accurate DOAs.

Co-prime arrays can estimate more sources than the number of sensors, see Fig. 2. We consider the same co-prime array and $K=10$ stationary sources uniformly distributed in $[-60,60]^{\circ}$, with $L=20$ and SNR 20 dB. The histogram shows the distribution of the DOA estimates of AP-GLS.

The DOA performance is evaluated with the RMSE versus SNR, see Fig. 3. RMSE larger than 10 times the median is outlier and eliminated. We consider $K=4$ sources, same as in Fig. 1 but with equal strengths. Cramér-Rao bound (CRB) [32], MUSIC and MUSIC with co-array interpolation (MUSIC-I) [33] are also compared. Compared to full ULA cases, AP-GLS has a bounded error even with high SNR, which is come from the fact that $\mathbf{R}$ is recovered from its submatrix $\mathbf{\Gamma}_{\Omega}\mathbf{R}\mathbf{\Gamma}_{\Omega}^{\mathsf{T}}$.

The DOA performance is evaluated with the RMSE versus number of sources $K$, see Fig. 4. We consider the same co-prime array and for each case, $K$ equal strength sources are generated randomly in $[-65,65]^{\circ}$, with $L=20$ and SNR 20 dB. AP-GLS has higher estimation accuracy than CBF and estimating more sources than the number of sensors.

We introduced alternating projections based gridless sparse iterative covariance-based estimation for direction-of-arrival estimation that is gridless and promotes sparse solutions. Numerical evaluations indicated that the proposed method shows a favorable performance even with single-snapshot data and coherent arrivals. For co-prime array data, the proposed algorithm resolved more sources than the number of sensors.