Sparsity-Based Two-Dimensional DOA Estimation for Co-Prime Planar Array via Enhanced Matrix Completion

Abstract

In this paper, the two-dimensional (2-D) direction-of-arrival (DOA) estimation problem is explored for the sum-difference co-array (SDCA) generated by the virtual aperture expansion of co-prime planar arrays (CPPA). Since the SDCA has holes, this usually causes the maximum virtual aperture of CPPA to be unavailable. To address this issue, we propose a complex-valued, sparse matrix recovery-based 2-D DOA estimation algorithm for CPPA via enhanced matrix completion. First, we extract the difference co-arrays (DCA) from SDCA and construct the co-array interpolation model via nuclear norm minimization to initialize the virtual uniform rectangular array (URA) that does not contain the entire rows and columns of holes. Then, we utilize the shift-invariance structure of the virtual URA to construct the enhanced matrix with a two-fold Hankel structure to fill the remaining empty elements. More importantly, we apply the alternating direction method of the multipliers (ADMM) framework to solve the enhanced matrix completion model. To reduce the computational complexity of the traditional vector-form, sparse recovery algorithm caused by the Kronecker product operation between dictionary matrices, we derive a complex-valued sparse matrix-recovery model based on the fast iterative shrinkage-thresholding (FISTA) method. Finally, simulation results demonstrate the effectiveness of the proposed algorithm.

1. Introduction

Direction-of-arrival (DOA) estimation is a fundamental issue for array signal processing and has been widely applied in radar, sonar, and navigation [1,2,3]. Various methods have been proposed for DOA estimation, such as classical subspace-based MUSIC and ESPRIT methods [4,5,6,7,8,9]. However, these methods rely on uniform linear arrays (ULA) or uniform rectangular arrays (URA), leading to mutual coupling and redundancy problems between sensors. Recently, a novel coprime array (CPA) [10,11] has been proposed and received extensive attention. Compared with the ULA, the CPA has a virtual aperture expansion capability by calculating its difference co-array (DCA), improving angle estimation performance. In addition, the distance between the sensors of the CPA is no longer limited to the half-wavelength, which can alleviate the mutual coupling effect between the physical sensors.

In order to exploit the virtual aperture of the CPA, many DOA estimation algorithms have been proposed in [12,13,14,15,16,17,18]. A spatial smoothing-based MUSIC algorithm [12] and a sparse recovery-based method [13] have been proposed which utilize the consecutive parts of the DCA and discard the non-consecutive parts for DOA estimation. However, none of these methods utilize the non-consecutive part of the DCA, which usually causes the virtual aperture loss of the DCA, reducing the DOA estimation performance. In order to make full use of the virtual aperture of the CPA, co-array interpolation-based algorithms [14,15] were proposed to fill the holes in the DCA. In [16], the holes in the DCA are filled by minimizing the nuclear norm of the covariance matrix of the received signal, generating a consecutive virtual ULA to improve angle estimation performance. In [17,18], the covariance matrix with the Toeplitz structure corresponding to the interpolated virtual array was reconstructed through atomic norm minimization, improving the DOA estimation accuracy. In [19,20], co-array interpolation algorithms based on truncated nuclear norm regularization were proposed to interpolate the covariance matrix with the Toeplitz structure corresponding to a virtual ULA, thereby improving the virtual aperture of the CPA. A 2-D DOA estimation algorithm was proposed for coprime parallel arrays via sparse representation in [21], extending the advantages of CPA. Although this algorithm implies sparse representation to improve the performance of DOA estimation, it does not utilize the maximum virtual aperture of the coprime parallel array, so there is an inherent aperture loss. In [22], a 2-D DOA fast estimation of multiple signals with matrix completion theory in a coprime planar array was proposed. This algorithm uses the matrix completion technique to complete the part holes of the virtual URA, which limits the number of DOFs. In addition, the missing elements in the entire rows and columns of the virtual URA are not accurately filled, which reduces recovery performance. In [23], a partial spectral search algorithm was proposed for the coprime planar array (CPPA), which improves DOA estimation accuracy by exploiting the coprime property to eliminate ambiguity. To benefit from the virtual aperture expansion of CPPA, Shi et al. in [24] proposed a sparsity-based, two-dimensional (2-D) DOA estimation algorithm for the sum-difference co-array (SDCA) of CPPA. However, this algorithm only utilizes the consecutive parts of the SDCA and discards the non-consecutive parts, resulting in an inherent loss of virtual aperture. Furthermore, co-array interpolation-based methods [14,15,16,17] cannot be used to fill the holes in SDCA. The reason is that the holes of SDCA have empty elements of the entire rows and columns, so the nuclear norm minimization-based algorithm [25] fails to recover the missing elements.

This paper proposes a sparse matrix-recovery-based 2-D DOA estimation algorithm via enhanced matrix completion to fill the holes in the SDCA. Using the coprime property, we first establish the SDCA of CPPA by rearranging the cross-correlation matrix of two coprime planar sub-arrays. Since the SDCA has missing elements in the entire rows and columns, the nuclear norm minimization-based matrix completion algorithm cannot be utilized to fill the holes. To solve this problem, we interpolate the covariance matrices of the DCA with the Toeplitz structure, thereby initializing the SDCA that does not contain the entire rows or columns of empty elements. Then, we use the shift-invariance structure of the SDCA to construct the enhanced matrix with a two-fold Hankel structure to fill all the empty elements accurately. Subsequently, we construct the nuclear norm, minimization-based, enhanced matrix completion model. Then, we utilize the multi-block, alternating direction method of the multipliers (ADMM) [26] framework for solving this model to reduce the computational burden caused by the increase in enhanced matrix dimensions. More importantly, the proposed enhanced matrix completion model can obtain the maximum virtual aperture of the CPPA, improving the angle estimation accuracy. Finally, we also derive a complex-valued, sparse matrix-recovery model based on the fast iterative shrinkage-thresholding (FISTA) method [27] reducing the computational complexity of the traditional vector-form, sparse recovery algorithm.

The main contributions of the proposed algorithm are as follows.

• We extract the DCA from the SDCA and use the co-array interpolation algorithm to fill the holes in the DCA, thereby constructing a virtual URA model that does not contain the missing elements of the entire rows and columns.
• To accurately complete all the empty elements of the virtual URA, we construct the enhanced matrix-completion model. Furthermore, we derive an enhanced matrix-completion algorithm based on the ADMM framework, reducing the computational complexity of the matrix completion.
• We derive a complex-valued, sparse matrix-recovery algorithm based on the FISTA method, which avoids the Kronecker product operation between dictionary matrices, reducing the computational complexity of the traditional vector-form, sparse recovery algorithm.

The rest of this paper is organized as follows. In Section 2, we present the coprime planar array signal model, establish an enhanced matrix-completion model, and derive a matrix-completion algorithm based on the ADMM framework. In addition, we derive a complex-valued, sparse matrix-recovery algorithm based on the FISTA method. In Section 3, we demonstrate the simulation results. Discussion and conclusions are provided in Section 4 and Section 5, respectively.

Notations: We use lowercase bold letters for vectors, uppercase bold letters for matrices, and italics for scalars. ${\left(·\right)}^T$, ${\left(·\right)}^{*}$ and ${\left(·\right)}^H$ represent transpose, complex conjugate and conjugate transpose, respectively. ${\left(·\right)}^{-1}$ is matrix inversion. ⊗ and ∘ represent the Kronecker product and Hadamard product, respectively. $\mathrm{diag}\left(·\right)$ and $\mathrm{vec}\left(·\right)$ are the diagonalization and vectorization operations. ${\mathbf{I}}_N$ denotes the $N\times N$ identity matrix, and ${\mathbf{J}}_N$ represents the $N\times N$ exchange matrix, with ones on its anti-diagonal and zeros elsewhere. $\nabla \left(·\right)$ denotes the conjugate gradient of a complex-valued vector or matrix. $〈\mathbf{A},\mathbf{B}〉$ denotes the inner product of $\mathbf{A}$ and $\mathbf{B}$. $\mathrm{tr}\left(·\right)$ refers to the trace of the matrix. The subscript ${\left(·\right)}_{+}$ refers to a scalar $u_{+}>0$, then $u_{+}=u$ and is equal to zero otherwise. ${∥·∥}_1$, ${∥·∥}_2$, and ${∥·∥}_F$ are the $l_1$-norm, $l_2$-norm, and Frobenius norm, respectively.

2. Method

In this section, we first describe the signal model for CPPA. Then, we define the concept of SDCA by extracting the cross-correlation matrix of CPPA. Since there are entire rows and columns of missing elements in the SDCA, two co-array interpolation optimization problems are proposed to initialize the holes in the DCA to generate a virtual URA that does not contain the entire rows and columns of empty elements. Subsequently, we construct a nuclear norm, minimization-based, enhanced matrix-completion model to recover all missing elements in the virtual URA. Finally, a complex-valued, sparse matrix-recovery algorithm based on the FISTA method is derived to estimate the 2-D DOA of targets.

2.1. Signal Model

A CPPA consists of two planar sub-arrays with $2M\times 2M$ and $N\times N$ physical sensors, respectively, where M and N are coprime integers. As shown in Figure 1, the CPPA contains a total of $T=4M^2+N^2-1$ physical sensors, and the sensor position sets of the first and second planar sub-array are ${\mathbb{P}}_1$ and ${\mathbb{P}}_2$, respectively. Therefore, the sensor position set of the CPPA is $\mathbb{P}$ and can be expressed as

$$\begin{array}{cc}\hfill \mathbb{P}& ={\mathbb{P}}_1\cup {\mathbb{P}}_2\hfill \\ \hfill & =\left\{\left(m_x{d}_1,m_y{d}_1\right)\right\}\cup \left\{\left(n_x{d}_2,n_y{d}_2\right)\right\}\hfill \end{array}$$

(1)

where $0\le m_x,m_y\le 2M-1$, $0\le n_x,n_y\le N-1$, $d_1=Nd$ and $d_2=Md$ are the inter-element spacing of the first sub-array and the second sub-array along the x and y directions, respectively, and $d=\lambda /2$, $\lambda$ is the signal wavelength.

Assume there are K uncorrelated narrowband far-field signals impinging the CPPA from direction $\left(\mathbf{\theta },\mathbf{\phi }\right)=\left[\left({\theta }_1,{\phi }_1\right),\left({\theta }_2,{\phi }_2\right),\cdots ,\left({\theta }_K,{\phi }_K\right)\right]$, where ${\theta }_k$ and ${\phi }_k$ denote the elevation and azimuth of the kth signal. We present the main assumptions of the signal model in

Table 1. The main assumptions of the signal model.

The structure of array	$\mathit{M}>1$ and $\mathit{N}>1$ are co-prime integers
The statistical properties of the signal	follow the independent Gaussian distribution
The statistical properties of the noise	complex-valued additive white Gaussian noise
The inter-element spacing of the first sub-array	$d_1=Nd$
The inter-element spacing of the second sub-array	$d_2=Md$
The transform frequency domain variable u	$sin\left(\theta \right)cos\left(\phi \right)$
The transform frequency domain variable v	$sin\left(\theta \right)sin\left(\phi \right)$

.

After vectorization, the received signals of the first planar sub-array in the vector form of the qth snapshot can be expressed as

$${\mathbf{x}}_1\left(q\right)=\mathbf{As}\left(q\right)+{\mathbf{n}}_1\left(q\right)$$

(2)

where $\mathbf{s}\left(q\right)={\left[s_1\left(q\right),s_2\left(q\right),\cdots ,s_K\left(q\right)\right]}^T$ is the complex-valued signal vector, ${\mathbf{n}}_1\left(q\right)$ is the complex-valued white Gaussian noise with zero mean and variance ${\sigma }^2$, $\mathbf{A}={\mathbf{A}}_x\otimes {\mathbf{A}}_y$, ${\mathbf{A}}_x=\left[{\mathbf{a}}_x\left(u_1\right),{\mathbf{a}}_x\left(u_2\right),\cdots ,{\mathbf{a}}_x\left(u_K\right)\right]$, ${\mathbf{A}}_y=\left[{\mathbf{a}}_y\left(v_1\right),{\mathbf{a}}_y\left(v_2\right),\cdots ,{\mathbf{a}}_y\left(v_K\right)\right]$. ${\mathbf{a}}_x\left(u_k\right)$ and ${\mathbf{a}}_y\left(v_k\right)$ are the steering vectors along the x and y directions, respectively, and can be expressed as

$${\mathbf{a}}_x\left(u_k\right)={\left[1,e^{-j\pi N{u}_k},\cdots ,e^{-j\pi (2M-1)N{u}_k}\right]}^T$$

(3)

$${\mathbf{a}}_y\left(v_k\right)={\left[1,e^{-j\pi N{v}_k},\cdots ,e^{-j\pi (2M-1)N{v}_k}\right]}^T$$

(4)

where $u_k=sin\left({\theta }_k\right)cos\left({\phi }_k\right)$ and $v_k=sin\left({\theta }_k\right)sin\left({\phi }_k\right)$ represent the transform frequency domain variables. Similarly, the received signal model of the second planar sub-array can be given as

$${\mathbf{x}}_2\left(q\right)=\mathbf{Bs}\left(q\right)+{\mathbf{n}}_2\left(q\right)$$

(5)

where ${\mathbf{n}}_2\left(q\right)$ is the complex-valued, white Gaussian noise vector, $\mathbf{B}={\mathbf{B}}_x\otimes {\mathbf{B}}_y$, ${\mathbf{B}}_x=\left[{\mathbf{b}}_x\left(u_1\right),{\mathbf{b}}_x\left(u_2\right),\cdots ,{\mathbf{b}}_x\left(u_K\right)\right]$, ${\mathbf{B}}_y=\left[{\mathbf{b}}_y\left(v_1\right),{\mathbf{b}}_y\left(v_2\right),\cdots ,{\mathbf{b}}_y\left(v_K\right)\right]$. ${\mathbf{b}}_x\left(u_k\right)$ and ${\mathbf{b}}_y\left(v_k\right)$ can be represented as

$${\mathbf{b}}_x\left(u_k\right)={\left[1,e^{-j\pi M{u}_k},\cdots ,e^{-j\pi (N-1)M{u}_k}\right]}^T$$

(6)

$${\mathbf{b}}_y\left(v_k\right)={\left[1,e^{-j\pi M{v}_k},\cdots ,e^{-j\pi (N-1)M{v}_k}\right]}^T$$

(7)

2.2. The Virtual Ura with Sdca

Similar to [24], the difference co-array (DCA) between two planar sub-arrays of a CPPA can be defined as

$$\begin{array}{cc}\hfill {\mathbb{D}}^1& =\left\{p\mid p=p_1-p_2,p_1\in {\mathbb{P}}_1,p_2\in {\mathbb{P}}_2\right\}\hfill \\ \hfill & =\left\{N{m}_x+N{m}_y-\left(M{n}_x+M{n}_y\right)\right\}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\mathbb{D}}^2& =\left\{p\mid p=p_2-p_1,p_1\in {\mathbb{P}}_1,p_2\in {\mathbb{P}}_2\right\}\hfill \\ \hfill & =\left\{M{n}_x+M{n}_y-\left(N{m}_x+N{m}_y\right)\right\}\hfill \end{array}$$

(9)

where $0\le m_x,m_y\le 2M-1$, $0\le n_x,n_y\le N-1$. Then, we can define the sum-difference co-array (SDCA) as follows.

$$\mathbb{D}={\mathbb{D}}^1\cup {\mathbb{D}}^2$$

(10)

From Equations (8) and (9), the maximum virtual aperture position of ${\mathbb{D}}^1$ and ${\mathbb{D}}^2$ in x-direction or y-direction are $\left(2MN-N\right)$ and $-\left(2MN-N\right)$. Similar to [24], the continuous lag ranges of ${\mathbb{D}}^1$ and ${\mathbb{D}}^2$ are $\left[-N+1,MN+M-1\right]$ and $\left[-\left(MN+M-1\right),N-1\right]$. Therefore, we have

$$\begin{array}{c}2MN-N-\left(MN+M-1\right)\hfill \\ =MN-N-M+1\hfill \\ =\left(M-1\right)\left(N-1\right)>0\hfill \end{array}$$

(11)

where M and N are coprime integers greater than 1. From Equation (11), there are always entire rows and columns of missing elements for SDCA of CPPA with $M>1$ and $N>1$.

The cross-correlation matrix of the first and second planar sub-arrays can be expressed as

$$\begin{array}{cc}\hfill {\widehat{\mathbf{R}}}_{12}& =1/Q\sum _{q=1}^Q{\mathbf{x}}_1\left(q\right){\mathbf{x}}_2^H\left(q\right)\hfill \\ \hfill & =\mathbf{A}{\mathbf{R}}_s{\mathbf{B}}^H+{\mathbf{N}}_{12}\hfill \end{array}$$

(12)

where Q is the total number of snapshots, ${\mathbf{R}}_s=\mathrm{diag}\left({\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_K^2\right)$, ${\sigma }_k^2$ is the power of the kth signal and ${\mathbf{N}}_{12}$ is the remaining residual term and can be expressed as

$$\begin{array}{cc}\hfill & {\mathbf{N}}_{12}=\frac{\mathbf{A}}{Q}\sum _{q=1}^Q\mathbf{s}\left(q\right){\mathbf{n}}_2^H\left(q\right)+\hfill \\ \hfill & \sum _{q=1}^Q{\mathbf{n}}_1\left(q\right){\mathbf{s}}^H\left(q\right)\frac{{\mathbf{B}}^H}{Q}+\frac{1}{Q}\sum _{q=1}^Q{\mathbf{n}}_1\left(q\right){\mathbf{n}}_2^H\left(q\right)\hfill \end{array}$$

(13)

Then, we calculate the cross-correlation matrix of the second and the first planar sub-array, which can be given as

$$\begin{array}{cc}\hfill {\widehat{\mathbf{R}}}_{21}& =1/Q\sum _{q=1}^Q{\mathbf{x}}_2\left(q\right){\mathbf{x}}_1^H\left(q\right)\hfill \\ \hfill & =\mathbf{B}{\mathbf{R}}_s{\mathbf{A}}^H+{\mathbf{N}}_{21}.\hfill \end{array}$$

(14)

where ${\mathbf{N}}_{21}$ is the remaining residual term and can be written as

$$\begin{array}{cc}\hfill & {\mathbf{N}}_{21}=\frac{\mathbf{B}}{Q}\sum _{q=1}^Q\mathbf{s}\left(q\right){\mathbf{n}}_1^H\left(q\right)+\hfill \\ \hfill & \sum _{q=1}^Q{\mathbf{n}}_2\left(q\right){\mathbf{s}}^H\left(q\right)\frac{{\mathbf{A}}^H}{Q}+\frac{1}{Q}\sum _{q=1}^Q{\mathbf{n}}_2\left(q\right){\mathbf{n}}_1^H\left(q\right).\hfill \end{array}$$

(15)

By extracting and rearranging ${\widehat{\mathbf{R}}}_{12}$ and ${\widehat{\mathbf{R}}}_{21}$, we can obtain a matrix $\mathbf{Z}\in {\mathbb{C}}^{n\times n}$ corresponding to a virtual URA, where $n=2\left(2MN-N\right)+1$ and its elements are located at the lag positions in the $\mathbb{D}$, as shown in Figure 2. Furthermore, we initialize the holes in $\mathbf{Z}$ to zeros. From Figure 2, there are many missing elements in $\mathbf{Z}$. If we discard these discontinuous parts, it usually causes the loss of the virtual aperture of $\mathbf{Z}$. We introduce a matrix completion model to fill all holes to avoid virtual aperture loss.

2.3. Co-Array Interpolation Initializes the Virtual Ura

Since there are missing elements in the entire rows and columns in $\mathbf{Z}$, the popular matrix completion methods cannot recover these holes [25].

The main reason is that the empty elements of the entire rows and columns will make the singular vector of $\mathbf{Z}$ more sparse, thereby reducing the recovery performance. To solve this problem, we first extract a DCA from the y direction in $\mathbf{Z}$, thereby obtaining a virtual ULA ${\mathbf{y}}_1$ from $\left(0,-\left(2MN-N\right)\right)d$ to $\left(0,2MN-N\right)d$, where the holes at positions $\left(0,-\left(2MN-N-1\right)\right)d$ and $\left(0,2MN-N-1\right)d$ are initialized to zeros. Similarly, we can also extract a DCA from the x direction in $\mathbf{Z}$, getting a virtual ULA ${\mathbf{y}}_2$ from $\left(-\left(2MN-N\right),0\right)d$ to $\left(2MN-N,0\right)d$, where the holes at positions $\left(-\left(2MN-N-1\right),0\right)d$ and $\left(2MN-N-1,0\right)d$ are initialized to zeros.

Thus, ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$ can be represented as

$${\mathbf{y}}_1=\overline{\mathbf{A}}s+{\mathbf{e}}_1$$

(16)

$${\mathbf{y}}_2=\overline{\mathbf{B}}s+{\mathbf{e}}_2$$

(17)

where ${\mathbf{e}}_1$ and ${\mathbf{e}}_2$ are error terms, $\overline{\mathbf{A}}=\left[\overline{\mathbf{a}}\left(u_1\right),\overline{\mathbf{a}}\left(u_2\right),\dots ,\overline{\mathbf{a}}\left(u_K\right)\right]$, $\overline{\mathbf{B}}=\left[\overline{\mathbf{b}}\left(v_1\right),\overline{\mathbf{b}}\left(v_2\right),\dots ,\overline{\mathbf{b}}\left(v_K\right)\right]$. $\overline{\mathbf{a}}\left(u_k\right)$ and $\overline{\mathbf{b}}\left(v_k\right)$ can be represented as

$$\begin{array}{c}\overline{\mathbf{a}}\left(u_k\right)=\left[e^{-j\pi (L-1)u_k},0,e^{-j\pi (L-3)u_k},\dots ,e^{-j\pi u_k}\right.\\ {\left.1,e^{j\pi u_k},\dots ,e^{j\pi (L-3)u_k},0,e^{j\pi (L-1)u_k}\right]}^T\end{array}$$

(18)

$$\begin{array}{c}\overline{\mathbf{b}}\left(v_k\right)=\left[e^{-j\pi (L-1)v_k},0,e^{-j\pi (L-3)v_k},\dots ,e^{-j\pi v_k},\right.\\ {\left.1,e^{j\pi v_k},\dots ,e^{j\pi (L-3)v_k},0,e^{j\pi (L-1)v_k}\right]}^T\end{array}$$

(19)

where $L=(n+1)/2$.

To fill the holes of the virtual ULA, we usually divide ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$ into L sub-arrays to compute a smooth covariance matrix with Hermitian positive semidefinite (PSD) Toeplitz structure. According to [17], we can build two covariance matrices of virtual signals ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$ and they can be expressed as

$${\tilde{\mathbf{R}}}_1=\left[\begin{array}{cccc}{〈{\mathbf{y}}_1〉}_L& {〈{\mathbf{y}}_1〉}_{L+1}^{*}& \cdots & {〈{\mathbf{y}}_1〉}_{2L-1}^{*}\\ {〈{\mathbf{y}}_1〉}_{L+1}& {〈{\mathbf{y}}_1〉}_L& \cdots & {〈{\mathbf{y}}_1〉}_{2L-2}^{*}\\ ⋮& ⋮& \ddots & ⋮\\ {〈{\mathbf{y}}_1〉}_{2L-1}& {〈{\mathbf{y}}_1〉}_{2L-2}& \cdots & {〈{\mathbf{y}}_1〉}_L\end{array}\right]$$

(20)

$${\tilde{\mathbf{R}}}_2=\left[\begin{array}{cccc}{〈{\mathbf{y}}_2〉}_L& {〈{\mathbf{y}}_2〉}_{L+1}^{*}& \cdots & {〈{\mathbf{y}}_2〉}_{2L-1}^{*}\\ {〈{\mathbf{y}}_2〉}_{L+1}& {〈{\mathbf{y}}_2〉}_L& \cdots & {〈{\mathbf{y}}_2〉}_{2L-2}^{*}\\ ⋮& ⋮& \ddots & ⋮\\ {〈{\mathbf{y}}_2〉}_{2L-1}& {〈{\mathbf{y}}_2〉}_{2L-2}& \cdots & {〈{\mathbf{y}}_2〉}_L\end{array}\right]$$

(21)

${\tilde{\mathbf{R}}}_1$ and ${\tilde{\mathbf{R}}}_2$ have holes corresponding to the zeros of the virtual array signals ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$. Therefore, we build the binary matrix $\mathbf{D}\in {\mathbb{R}}^{L\times L}$ to distinguish virtual array signal positions and holes, where the elements of $\mathbf{D}$ corresponding to non-zero signal positions are ones, and the rest are zeros.

To fill the empty elements of ${\tilde{\mathbf{R}}}_1$ and ${\tilde{\mathbf{R}}}_2$, we formulate the nuclear norm minimization optimization problem as

$$min{∥{\mathbf{R}}_1∥}_{*},\mathrm{s}.\mathrm{t}.{∥\left({\tilde{\mathbf{R}}}_1-{\mathbf{R}}_1\right)\circ \mathbf{D}∥}_F\le {\delta }_1$$

(22)

$$min{∥{\mathbf{R}}_2∥}_{*},\mathrm{s}.\mathrm{t}.{∥\left({\tilde{\mathbf{R}}}_2-{\mathbf{R}}_2\right)\circ \mathbf{D}∥}_F\le {\delta }_2$$

(23)

where ${\delta }_1$ and ${\delta }_2$ are the thresholds of the fitting error, and ${∥·∥}_{*}$ represents the nuclear norm of a matrix, which is the relaxation of the rank norm of the matrix.

Since the optimization problems (22) and (23) are convex, we can solve them directly by employing the interior point method [28]. Then, we can obtain the estimated optimal covariance matrices ${\widehat{\mathbf{R}}}_1$ and ${\widehat{\mathbf{R}}}_2$ corresponding to ${\tilde{\mathbf{R}}}_1$ and ${\tilde{\mathbf{R}}}_2$. Subsequently, we extract the elements ${\widehat{\mathbf{R}}}_1\left(1,L-2\right)$ and ${\widehat{\mathbf{R}}}_1\left(L-2,1\right)$, filling the holes at positions $\left(0,-\left(L-2\right)\right)d$ and $\left(0,L-2\right)d$ in matrix $\mathbf{Z}$. Similarly, we use ${\widehat{\mathbf{R}}}_2\left(L-2,1\right)$ and ${\widehat{\mathbf{R}}}_2\left(1,L-2\right)$ to fill the holes of matrix $\mathbf{Z}$ at positions $\left(-\left(L-2\right),0\right)d$ and $\left(L-2,0\right)d$. After the above interpolation operation, we can initialize the matrix $\mathbf{Z}$ to obtain a matrix $\overline{\mathbf{Z}}$ that does not contain the missing elements of the entire rows and columns. In addition, we demonstrate the structure of the virtual signal position set $\overline{\mathbb{D}}$ corresponding to $\overline{\mathbf{Z}}$, as shown in Figure 3.

From Figure 3, $\overline{\mathbf{Z}}$ can be represented as a single-snapshot signal model received by a virtual URA and can be written as

$$\overline{\mathbf{Z}}=\left(\overline{\mathbf{A}}{\mathbf{R}}_s{\overline{\mathbf{B}}}^T+\mathbf{E}\right)\circ \mathbf{S}$$

(24)

where $\mathbf{E}\in {\mathbb{C}}^{n\times n}$ is the error term, and $\mathbf{S}\in {\mathbb{R}}^{n\times n}$ is a binary matrix whose elements correspond to holes in $\overline{\mathbf{Z}}$ as zeros, and the rest are ones.

2.4. Admm-Based Enhanced Matrix Completion

After the above co-array interpolation processing, there are no missing elements in the entire rows and columns in $\overline{\mathbf{Z}}$, which prompts us to use matrix completion to fill the remaining empty elements. Furthermore, we assume that $\overline{\mathbf{Z}}$ has a low-rank property because the number of targets is usually smaller than the dimension of $\overline{\mathbf{Z}}$. Based on these conditions, the matrix completion model can be expressed as

$$\underset{\mathbf{Y}}{min}{\parallel \mathbf{Y}\parallel }_{*},\mathrm{s}.\mathrm{t}.{∥U_{\mathrm{\Omega }}(\mathbf{Y}-\overline{\mathbf{Z}})∥}_F\le \alpha$$

(25)

where $\alpha$ is the threshold of the fitting error, $\mathrm{\Omega }$ is a set of positions of the known entries of $\overline{\mathbf{Z}}$, and the matrix operator $U_{\mathrm{\Omega }}$ can be defined as

$${\left(U_{\mathrm{\Omega }}\left(\mathbf{Y}\right)\right)}_{ij}=\left\{\begin{array}{cc}{\mathbf{Y}}_{ij},& \mathrm{if}(i,j)\in \mathrm{\Omega }\hfill \\ 0,& \mathrm{otherwise}\hfill \end{array}\right.$$

(26)

The optimization problem (25) does not use the TNNR-based matrix completion method proposed in [24]. The reason is that the TNNR method needs to know the number of targets, but, in practical problems, the number of targets is usually unknown. In addition, the optimization problem (25) is convex and can be solved by the interior point method [28]. However, the optimization model (25) does not take the inherent shift-invariant structure of $\overline{\mathbf{Z}}$, reducing the accuracy of matrix completion. What is worse, when the number of targets is larger than the dimension of $\overline{\mathbf{Z}}$, this model always fails to recover the completion matrix $\mathbf{Y}$.

In order to solve these problems, we propose an ADMM-based, enhanced matrix-completion algorithm. The proposed method makes full use of the shift-invariant structure of $\overline{\mathbf{Z}}$ to form the enhanced matrix ${\mathbf{M}}_e$ with a two-fold Hankel structure, thereby improving the performance of matrix completion. In addition, we also developed a fast ADMM algorithm to reduce the computational complexity of the NNM-based algorithm and named the algorithm EMaC-ADMM. According to [29], the enhanced matrix ${\mathbf{M}}_e$ corresponding to the matrix $\overline{\mathbf{Z}}$ can be defined as

$${\mathbf{M}}_e=\left[\begin{array}{cccc}{\mathbf{M}}_0& {\mathbf{M}}_1& \cdots & {\mathbf{M}}_{n_1-k1}\\ {\mathbf{M}}_1& {\mathbf{M}}_2& \cdots & {\mathbf{M}}_{n_1-k_1+1}\\ ⋮& ⋮& ⋮& ⋮\\ {\mathbf{M}}_{k_1-1}& {\mathbf{M}}_{k_1}& \cdots & {\mathbf{M}}_{n_1-1}\end{array}\right]$$

(27)

where $k_1=\frac{n+1}{2}$ is a pencil parameter. The ith block Hankel matrix ${\mathbf{M}}_i\in {\mathbb{C}}^{k_2\times \left(n-k_2+1\right)}$ contained in ${\mathbf{M}}_e$ can be expressed as

$${\mathbf{M}}_i=\left[\begin{array}{cccc}{\overline{\mathbf{Z}}}_{(i,0)}& {\overline{\mathbf{Z}}}_{(i,1)}& \cdots & {\overline{\mathbf{Z}}}_{\left(i,n-k_2\right)}\\ {\overline{\mathbf{Z}}}_{(i,1)}& {\overline{\mathbf{Z}}}_{(i,2)}& \cdots & {\overline{\mathbf{Z}}}_{\left(i,n-k_2+1\right)}\\ ⋮& ⋮& ⋮& ⋮\\ {\overline{\mathbf{Z}}}_{\left(i,k_2-1\right)}& {\overline{\mathbf{Z}}}_{\left(i,k_2\right)}& \cdots & {\overline{\mathbf{Z}}}_{(i,n-1)}\end{array}\right]$$

(28)

where $k_2=\frac{n+1}{2}$ is another pencil parameter.

Then, we can recover the missing elements of the noise-free enhanced matrix ${\mathbf{M}}_e$ by forming the matrix completion model as follows.

$$\begin{array}{c}min{∥{\mathbf{Y}}_e∥}_{*}\hfill \\ s.t.{∥U_{{\mathrm{\Omega }}_e}\left({\mathbf{Y}}_e\right)∥}_F={∥U_{{\mathrm{\Omega }}_e}\left({\mathbf{M}}_e\right)∥}_F\hfill \end{array}$$

(29)

where ${\mathrm{\Omega }}_e$ refers to the set of known element positions of the enhanced matrix ${\mathbf{M}}_e$, ${\mathbf{Y}}_e\in {\mathbb{C}}^{m\times m}$ is the enhanced form of $\mathbf{Y}\in {\mathbb{C}}^{n\times n}$, and $m=k_1{k}_2$. However, in practical applications, the existence of noise is inevitable. Therefore, the enhanced matrix completion model (29) for bounded white Gaussian noise can be expressed as

$$\underset{{\mathbf{Y}}_e}{min}{∥{\mathbf{Y}}_e∥}_{*},s.t.{∥U_{{\mathrm{\Omega }}_e}\left({\mathbf{Y}}_e-{\mathbf{M}}_e\right)∥}_F\le \delta$$

(30)

where $\delta >0$ is a parameter to constrain the noise. The optimization problem (30) can be efficiently solved by the interior point method [28]. However, the interior point method takes up large memory when completing large-scale matrices, reducing the running speed of the algorithm.

Given the superior convergence performance and fast convergence speed of the ADMM framework, we propose an EMaC-ADMM method to solve the optimization problem (30). Therefore, we rewrite the optimization problem (30) as

$$\begin{array}{c}\underset{{\mathbf{Y}}_e,{\mathbf{N}}_e,\mathbf{W}}{min}{∥{\mathbf{Y}}_e∥}_{*}\hfill \\ s.t.{\mathbf{Y}}_e+{\mathbf{N}}_e=\mathbf{W}\hfill \\ U_{{\mathrm{\Omega }}_e}\left(\mathbf{W}\right)=U_{{\mathrm{\Omega }}_e}\left({\mathbf{M}}_e\right)\hfill \\ {∥U_{{\mathrm{\Omega }}_e}\left({\mathbf{N}}_e\right)∥}_F\le \sigma \hfill \end{array}$$

(31)

where ${\mathbf{N}}_e\in {\mathbb{C}}^{m\times m}$ refers to the enhanced form of the noise matrix $\mathbf{N}\in {\mathbb{C}}^{n\times n}$, and $\mathbf{W}\in {\mathbb{C}}^{m\times m}$ is an additional variable.

The augmented Lagrangian function of (31) is defined as

$$\begin{array}{c}{L}_{\beta }\left({\mathbf{Y}}_e,{\mathbf{N}}_e,\mathrm{\Lambda }\right)={∥{\mathbf{Y}}_e∥}_{*}+<\mathrm{\Lambda },{\mathbf{Y}}_e+{\mathbf{N}}_e-\mathbf{W}>+\frac{\beta }{2}{∥{\mathbf{Y}}_e+{\mathbf{N}}_e-\mathbf{W}∥}_F^2\hfill \end{array}$$

(32)

where $\mathrm{\Lambda }\in {\mathbb{C}}^{m\times m}$ are the Lagrange multipliers (dual variable), and $\beta >0$ denotes the penalty parameter. We utilize the ADMM framework to iteratively update the primal optimal variables ${\mathbf{Y}}_e$, ${\mathbf{N}}_e$, $\mathbf{W}$, and the dual variable $\mathrm{\Lambda }$ by finding the saddle point of the augmented Lagrange. In the kth iteration, the specific execution process of ADMM can be expressed as

$${\mathbf{Y}}_e^{k+1}=\underset{{\mathbf{Y}}_e}{argmin}{L}_{\beta }\left({\mathbf{Y}}_e,{\mathbf{N}}_e^k,{\mathbf{W}}^k,{\mathrm{\Lambda }}^k\right)$$

(33)

$${\mathbf{N}}_e^{k+1}=\underset{{\mathbf{N}}_e}{argmin}{L}_{\beta }\left({\mathbf{Y}}_e^{k+1},{\mathbf{N}}_e,{\mathbf{W}}^k,{\mathrm{\Lambda }}^k\right)$$

(34)

$${\mathbf{W}}^{k+1}=\underset{U_{{\mathrm{\Omega }}_e}\left(\mathbf{W}\right)=U_{{\mathrm{\Omega }}_e}\left({\mathbf{M}}_e\right)}{argmin}{L}_{\beta }\left({\mathbf{Y}}_e^{k+1},{\mathbf{N}}_e^{k+1},\mathbf{W},{\mathrm{\Lambda }}^k\right)$$

(35)

$${\mathrm{\Lambda }}^{k+1}={\mathrm{\Lambda }}^k+\beta U_{{\mathrm{\Omega }}_e}\left({\mathbf{Y}}_e^{k+1}+{\mathbf{N}}_e^{k+1}-{\mathbf{W}}^{k+1}\right)$$

(36)

For problems (33) to (36), ADMM updates ${\mathbf{Y}}_e^{k+1}$, ${\mathbf{N}}_e^{k+1}$, ${\mathbf{W}}^{k+1}$ and ${\mathrm{\Lambda }}^{k+1}$ in an alternate iterative manner, where the dual variable ${\mathrm{\Lambda }}^{k+1}$ is updated by the gradient ascent method.

For sub-problem (33), by ignoring the constant term, ${\mathbf{Y}}_e^{k+1}$ can be updated as

$${\mathbf{Y}}_e^{k+1}=arg\underset{{\mathbf{Y}}_e}{min}{∥{\mathbf{Y}}_e∥}_{*}+\frac{\beta }{2}{∥{\mathbf{Y}}_e+{\mathbf{N}}_e^k-{\mathbf{W}}^k+\frac{{\mathrm{\Lambda }}^k}{\beta }∥}_F^2$$

(37)

The convex optimization problem (37) can be efficiently solved by the singular value shrinkage operator [30], and we have

$${\mathbf{Y}}_e^{k+1}=\mathbf{U}\mathrm{diag}\left(max\left\{{\sigma }_i-\frac{1}{\beta },0\right\}\right){\mathbf{V}}^H$$

(38)

where ${\sigma }_i$, $\mathbf{U}$ and $\mathbf{V}$ are the ith singular value, left and right singular vectors of ${\mathbf{M}}_e-{\mathbf{N}}_e^k-\frac{{\mathrm{\Lambda }}^k}{\beta }$, respectively.

Given ${\mathbf{Y}}_e^{k+1}$ and ignoring the irrelevant terms of ${\mathbf{N}}_e$, sub-problem (34) can be rewritten as

$$\begin{array}{c}{\mathbf{N}}_e^{k+1}=arg\underset{{\mathbf{N}}_e}{min}{∥{\mathbf{N}}_e-\mathbf{C}∥}_F^2\hfill \\ s.t.{∥U_{{\mathrm{\Omega }}_e}\left({\mathbf{N}}_e\right)∥}_F\le \delta \hfill \end{array}$$

(39)

where $\mathbf{C}=\left({\mathbf{W}}^k-{\mathbf{Y}}_e^{k+1}-{\mathrm{\Lambda }}^k/\beta \right)$. The solution to the optimization problem (39) can be divided into two cases. In the first case, when ${∥U_{{\mathrm{\Omega }}_e}\left(\mathbf{C}\right)∥}_F\le \delta$, the solution of (39) can be given directly

$${\mathbf{N}}_e^{k+1}=U_{{\mathrm{\Omega }}_e}\left(\mathbf{C}\right)$$

(40)

For the second case, if ${∥U_{{\mathrm{\Omega }}_e}\left(\mathbf{C}\right)∥}_F>\delta$, the constraints of the optimization problem (39) can be rewritten as equality constraints [31] and can be expressed as

$$\begin{array}{c}{\mathbf{N}}_e^{k+1}=arg\underset{{\mathbf{N}}_e}{min}{∥{\mathbf{N}}_e-\mathbf{C}∥}_F^2\hfill \\ s.t.{∥U_{{\mathrm{\Omega }}_e}\left({\mathbf{N}}_e\right)∥}_F=\delta \hfill \end{array}$$

(41)

therefore, the solution of optimization problem (41) can be given

$${\mathbf{N}}_e^{k+1}=\frac{\delta U_{{\mathrm{\Omega }}_e}\left(\mathbf{C}\right)}{{∥U_{{\mathrm{\Omega }}_e}\left(\mathbf{C}\right)∥}_F}$$

(42)

Fixing ${\mathbf{Y}}_e^{k+1}$ and ${\mathbf{N}}_e^{k+1}$ and ignoring the irrelevant terms of $\mathbf{W}$, (35) can be rewritten as

$${\mathbf{W}}^{k+1}=\underset{U_{{\mathrm{\Omega }}_e}\left(\mathbf{W}\right)=U_{{\mathrm{\Omega }}_e}\left({\mathbf{M}}_e\right)}{argmin}\frac{\beta }{2}{∥\mathbf{W}-({\mathbf{Y}}_e^{k+1}+{\mathbf{N}}_e^{k+1}+\frac{{\mathrm{\Lambda }}^k}{\beta })∥}_F^2$$

(43)

for (43), ${\mathbf{W}}^{k+1}$ has a closed-form solution and can be expressed as

$${\mathbf{W}}^{k+1}={\mathbf{Y}}_e^{k+1}+{\mathbf{N}}_e^{k+1}+\frac{{\mathrm{\Lambda }}^k}{\beta }$$

(44)

Then, we add constraints on $\mathbf{W}$ and fix the known elements of ${\mathbf{M}}_e$, and finally update ${\mathbf{W}}^{k+1}$ as

$${\mathbf{W}}^{k+1}=U_{{\mathrm{\Omega }}_e^c}\left({\mathbf{W}}^{k+1}\right)+U_{{\mathrm{\Omega }}_e}\left({\mathbf{M}}_e\right)$$

(45)

where ${\mathrm{\Omega }}_e^c$ refers to the set of positions of the unobserved elements of the enhanced matrix ${\mathbf{M}}_e$.

The specific execution process of EMaC-ADMM is summarized in Algorithm 1. EMaC-ADMM utilizes the potential shift-invariance structure of $\overline{\mathbf{Z}}$ to improve the matrix completion performance. In addition, it also employs the ADMM framework to reduce the computational complexity of large-scale matrix completion problems caused by matrix enhancement.

Algorithm 1 EMaC-ADMM

1: Input:${\mathbf{M}}_e$, ${\mathrm{\Omega }}_e$, ${\epsilon }_1$, $\delta$, $\beta >0$; 2: Initialization:$k=1$, ${\mathbf{Y}}_e^1={\mathbf{M}}_e$, ${\mathbf{N}}_e^1=\mathbf{0}$, ${\mathbf{W}}_e^1={\mathbf{M}}_e$, ${\mathrm{\Lambda }}^{k+1}={\mathbf{M}}_e$. 3: repeat 4:    Update ${\mathbf{Y}}_e^{k+1}$ via (37). 5:    Update ${\mathbf{N}}_e^{k+1}$ via (39). 6:    Update ${\mathbf{W}}^{k+1}$ via (43). 7:    Update the Lagrange multiplier ${\mathrm{\Lambda }}^{k+1}$ via (36). 8:    $k=k+1$. 9: until ${∥{\mathbf{Y}}_e^{k+1}-{\mathbf{Y}}_e^k∥}_F<{\epsilon }_1$. 10: Output ${\widehat{\mathbf{Y}}}_e={\mathbf{Y}}_e^k$.

2.5. Sparse Matrix-Recovery-Based 2-D Doa Estimation

After the enhanced matrix completion in the previous sub-section, we can recover the completion matrix $\widehat{\mathbf{Y}}$ from the enhanced form ${\widehat{\mathbf{Y}}}_e$, and it can be represented as

$$\widehat{\mathbf{Y}}=\overline{\mathbf{A}}{\mathbf{R}}_s{\overline{\mathbf{B}}}^T+\mathbf{E}$$

(46)

where $\mathbf{E}$ is the error term.

To estimate the 2-D DOA of targets, we establish the following $l_1$-norm-based sparse matrix recovery model.

$$\underset{\mathbf{X}}{min}\left\{\frac{1}{2}{∥\widehat{\mathbf{Y}}-\tilde{\mathbf{A}}\mathbf{X}{\tilde{\mathbf{B}}}^T∥}_F^2+\rho {∥\mathbf{X}∥}_1\right\}$$

(47)

where $\rho$ is the regularization parameter, $\mathbf{X}\in {\mathbb{C}}^{I\times J}$ is an unknown sparse matrix, $\tilde{\mathbf{A}}=\left[\overline{\mathbf{a}}\left(u_1\right),\overline{\mathbf{a}}\left(u_2\right),\dots ,\overline{\mathbf{a}}\left(u_I\right)\right]$ and $\tilde{\mathbf{B}}=\left[\overline{\mathbf{b}}\left(v_1\right),\overline{\mathbf{b}}\left(v_2\right),\dots ,\overline{\mathbf{b}}\left(v_J\right)\right]$ are overcomplete dictionary matrices in the transform domain $u\in \left[0,1\right]$ and $v\in \left[-1,1\right]$, respectively, and I and J represent the number of sampling grids.

The optimization problem (47) can be transformed into a sparse recovery problem in vector form and solved by the interior point method [28]. However, the Kronecker product operation between dictionary matrices will significantly increase the dimension of the dictionary matrix, which increases the computational burden of the vector-form, sparse recovery algorithm.

To reduce the computational burden of the traditional vector-form, sparse recovery algorithm, we utilize the FISTA method [32] to perform sparse matrix recovery. However, the FISTA method only works on real-valued data and cannot be directly applied to complex-valued data. Therefore, we derive a complex-valued FISTA algorithm named (CV-FISTA). The specific derivation process is as follows.

Similar to [32], we rewrite the optimization problem (47) as

$$\underset{\mathbf{X}}{min}H\left(\mathbf{X}\right)+\rho F\left(\mathbf{X}\right)$$

(48)

where $H\left(\mathbf{X}\right)=\frac{1}{2}{∥\widehat{\mathbf{Y}}-\tilde{\mathbf{A}}\mathbf{X}{\tilde{\mathbf{B}}}^T∥}_F^2$ is a smooth convex function, $F\left(\mathbf{X}\right)={∥\mathbf{X}∥}_1$ is a non-smooth convex function. To solve the optimization problem (48), we apply a quadratic approximation of $H\left(\mathbf{X}\right)$ at a given point $\mathbf{P}$.

$$\begin{array}{c}{Q}_L\left(\mathbf{X},\mathbf{P}\right)=H\left(\mathbf{P}\right)+\mathrm{tr}\left(\nabla H{\left(\mathbf{P}\right)}^T\left(\mathbf{X}-\mathbf{P}\right)\right)+\frac{L_f}{2}{∥\mathbf{X}-\mathbf{P}∥}_F^2+\rho F\left(\mathbf{X}\right)\hfill \end{array}$$

(49)

where $L_f={∥\tilde{\mathbf{B}}\otimes \tilde{\mathbf{A}}∥}_F^2$ is the Lipschitz constant. After some simplification, we can rewrite (49) as

$$Q_L\left(\mathbf{X},\mathbf{P}\right)=\frac{L_f}{2}{∥\mathbf{X}-U\left(\mathbf{P}\right)∥}_F^2+\rho {∥\mathbf{X}∥}_1+C\left(\mathbf{P}\right)$$

(50)

where $C\left(\mathbf{P}\right)$ is a constant term concerning $\mathbf{P}$, and $U\left(\mathbf{P}\right)$ can be expressed as

$$U\left(\mathbf{P}\right)=\mathbf{P}-\frac{1}{L_f}{\tilde{\mathbf{A}}}^H\left(\tilde{\mathbf{A}}\mathbf{X}{\tilde{\mathbf{B}}}^T-\widehat{\mathbf{Y}}\right){\tilde{\mathbf{B}}}^{*}$$

(51)

Proof: please see Appendix A.

Ignoring the constant term, we rewrite the optimization problem (50) as

$$arg\underset{\mathbf{X}}{min}{Q}_L\left(\mathbf{X},\mathbf{P}\right)=arg\underset{\mathbf{X}}{min}\left\{\frac{L_f}{2}{∥\mathbf{X}-U\left(\mathbf{P}\right)∥}_F^2+\rho {∥\mathbf{X}∥}_1\right\}$$

(52)

The optimization problem (52) can be solved by the soft threshold method [33] and can be expressed as

$$arg\underset{\mathbf{X}}{min}{Q}_L\left(\mathbf{X},\mathbf{P}\right)=\mathrm{soft}\left(U\left(\mathbf{P}\right),\frac{\rho }{L_f}\right)$$

(53)

where $\mathrm{soft}\left(·\right)$ can be defined as

$$\mathrm{soft}\left(\mathbf{X},c\right)={\mathbf{X}}_{ij}{\left(1-\frac{c}{\left|{\mathbf{X}}_{ij}\right|}\right)}_{+}$$

(54)

where $c>0$ is a constant and ${\mathbf{X}}_{ij}$ is the $\left(i,j\right)$th element of matrix $\mathbf{X}$.

We give the specific execution process of the CV-FISTA algorithm in Algorithm 2. After several iterations, we can get the sparse recovery matrix $\widehat{\mathbf{X}}$, and the azimuth and elevation angles can be obtained by calculating

$${\widehat{\phi }}_k={tan}^{-1}\left({\widehat{v}}_k/{\widehat{u}}_k\right),{\widehat{\theta }}_k={sin}^{-1}\sqrt{{\widehat{u}}_k^2+{\widehat{v}}_k^2}$$

(55)

where ${\widehat{u}}_k$ and ${\widehat{v}}_k$ are the estimated parameters to obtain the kth target from $\widehat{\mathbf{X}}$.

Algorithm 2 CV-FISTA

1: Input: the completion matrix $\widehat{\mathbf{Y}}$; overcomplete dictionary matrices $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$; tolerance ${\epsilon }_2$; the regularization parameter $\rho >0$. 2: Initialization: $k=1$, ${\mathbf{X}}_1=\mathbf{0}$, ${\mathbf{X}}_0=\mathbf{0}$, $t_1=1$, $t_0=1$. 3: repeat 4:    ${\mathbf{P}}_k={\mathbf{X}}_k+\frac{t_{k-1}-1}{t_k}\left({\mathbf{X}}_k-{\mathbf{X}}_{k-1}\right)$. 5:    $U\left({\mathbf{P}}_k\right)={\mathbf{P}}_k-\frac{1}{L_f}{\tilde{\mathbf{A}}}^H\left(\tilde{\mathbf{A}}{\mathbf{P}}_k{\tilde{\mathbf{B}}}^T-\widehat{\mathbf{Y}}\right){\tilde{\mathbf{B}}}^{*}$. 6:    ${\mathbf{X}}_{k+1}=soft\left(U\left({\mathbf{P}}_k\right),\frac{\rho }{L_f}\right)$. 7:    $t_{k+1}=(1+\sqrt{4t_k^2+1})/2$. 8:    $k=k+1$. 9: until ${∥{\mathbf{X}}_{k+1}-{\mathbf{X}}_k∥}_F<{\epsilon }_2$. 10: Output $\widehat{\mathbf{X}}={\mathbf{X}}_k$.

2.6. Complexity Analysis

We analyze the computational complexity of the proposed algorithm and compare it with state-of-the-art algorithms. The computational complexity of the proposed algorithm mainly includes two parts: the first part is matrix completion, and the second part is sparse matrix recovery. For the first part, the computational complexity of co-array interpolation initialization and enhanced matrix completion is about $\mathcal{O}\left(L^3{K}_1+m^3{K}_2\right)$, where m is the dimension of the enhanced matrix ${\mathbf{Y}}_e$, $K_1$ is the number of iterations of co-array interpolation, and $K_2$ is the number of iterations of the EMaC-ADMM algorithm. Then, we analyze the computational complexity of sparse matrix recovery, mainly caused by matrix multiplication, which is approximate $\mathcal{O}\left(I^3{K}_3\right)$, where I is the number of the searching grid and $K_3$ is the number of iterations of the CV-FISTA algorithm.

Therefore, the total computational complexity of the proposed algorithm is about $\mathcal{O}\left(L^3{K}_1+m^3{K}_2+I^3{K}_3\right)$. Subsequently, we analyze the computational complexity of the TNNR algorithm [20], which is mainly caused by singular value decomposition, about $\mathcal{O}\left(r{n}^2{K}_4\right)$, where $K_4$ is the number of iterations, and r is the rank of the matrix,. Then, the CV-FISTA method is used to estimate the angle of targets. Therefore, the final computational complexity of the TNNR algorithm is about $\mathcal{O}\left(r{n}^2{K}_4+I^3{K}_3\right)$.

Next, we analyze the computational complexity of the partial spectral search (PSS) algorithm [23] and the iterative sparse recovery (ISR) algorithm [24]. The computational complexity of the PSS algorithm is mainly caused by spectral search and eigenvalue decomposition, which is about $\mathcal{O}\left(I\left(\frac{{\left(2M\right)}^4}{N^2}+\frac{N^4}{{\left(2M\right)}^2}\right)+{\left(2M\right)}^3+N^3\right)$. Since I is much larger than M and N, the total computational complexity of the PSS algorithm is about $\mathcal{O}\left(I\left(\frac{{\left(2M\right)}^4}{N^2}+\frac{N^4}{{\left(2M\right)}^2}\right)\right)$. The computational complexity of the ISR algorithm is approximately $\mathcal{O}\left(r^3{I}^3{K}_5+I^3\right)$, where $K_5$ is the number of iterations. We clearly describe the computational complexity of the above algorithms in

Table 2. Computation complexity comparison of different algorithms.

Proposed algorithm	$\mathcal{O}\left({\mathit{L}}^3{\mathit{K}}_1+{\mathit{m}}^3{\mathit{K}}_2+{\mathit{I}}^3{\mathit{K}}_3\right)$
TNNR in [20]	$\mathcal{O}\left(r{n}^2{K}_4+I^3{K}_3\right)$
PSS in [23]	$\mathcal{O}\left(I\left(\frac{{\left(2M\right)}^4}{N^2}+\frac{N^4}{{\left(2M\right)}^2}\right)\right)$
ISR in [24]	$\mathcal{O}\left(r^3{I}^3{K}_5+I^3\right)$

.

Finally, we compare the computational complexity of the above algorithms versus the number of sensors in Figure 4, where M is fixed at 2, $N=\left[3,5,7,9\right]$. As shown in Figure 4, the computational complexity of the proposed algorithm is higher than the TNNR algorithm. The cause is that, in the process of matrix completion, the proposed algorithm needs the co-array interpolation initialization process, so it has higher computational complexity than the TNNR algorithm. The PSS algorithm has the lowest computational complexity because it does not require iterative operations. When $I\gg M,N,L$, and m, the computational complexity of the proposed algorithm is less than the ISR algorithm. However, as the total number of sensors increases, the proposed algorithm’s computational complexity increases faster than the other tested algorithms. Therefore, when the number of sensors is large, the computational complexity of the proposed algorithm is higher than that of the ISR algorithm. The reason is that the proposed algorithm’s enhanced matrix further expands the dimension of the virtual URA, causing the complexity to increase.

3. Results

In our simulations, we consider a CPPA with $M=2$ and $N=3$. The parameters $\delta =1$, ${\delta }_1=1$, ${\delta }_2=1$, $\beta =0.5$, and $\rho =1e^2$. After matrix completion processing, we can finally obtain a virtual URA with virtual sensor positions from $\left(-9,-9\right)d$ to $\left(9,9\right)d$. We compare the 2-D DOA estimation performance of the proposed algorithm with several recently proposed algorithms, including the TNNR algorithm [20], the ISR algorithm [24] and the PSS algorithm [23]. In addition, the Cramer–Rao bound (CRB) [24] is used as a reference. Finally, the root mean square error (RMSE) is utilized to compare the angle estimation performance of the above test algorithm and can be defined as

$$\mathrm{RMSE}=\sqrt{\frac{1}{TK}\sum _{t=1}^T\sum _{k=1}^K\left[{\left({\phi }_k-{\widehat{\phi }}_{k,t}\right)}^2+{\left({\theta }_k-{\widehat{\theta }}_{k,t}\right)}^2\right]}$$

(56)

where $\left({\widehat{\phi }}_{k,t},{\widehat{\theta }}_{k,t}\right)$ denotes the estimate of kth target $\left({\phi }_k,{\theta }_k\right)$ in the tth Monte Carlo trial, $T=100$ is the number of Monte Carlo trials.

In the first experiment, we compare the number of degrees of freedom (DOF) of the proposed algorithm with the other tested algorithms. After the enhanced matrix completion, the proposed algorithm has 361 DOFs corresponding to the matrix $\widehat{\mathbf{Y}}$. Since the ISR algorithm utilizes the consecutive part of the SDCA, it has 175 DOFs. A standard URA consisting of 24 physical sensors has 24 DOFs. In addition, the PSS algorithm has only 9 DOFs. Therefore, the proposed algorithm has the maximum DOFs compared to other algorithms. Figure 5a depicts the contour map of the normalized spectrum of the proposed algorithm, where the signal-to-noise ratio (SNR) is 0 dB, the number of snapshots is 500, and the searching grid is 0.002 for $u\in \left[0,1\right]$ and $v\in \left[-1,1\right]$. As shown in Figure 5a, the proposed algorithm can resolve all 27 targets using only 24 physical sensors. Figure 5b shows the contour map of the normalized spectrum of the TNNR algorithm. From Figure 5b, the TNNR algorithm can resolve up to 16 targets. This is because the rank is limited by the dimension of the matrix $\mathbf{Z}$, thereby reducing the maximum number of resolved targets of the TNNR algorithm. Figure 5c depicts the standard URA 2-D DOA estimation, where the 2-D unitary ESPRIT algorithm [7] is utilized. From Figure 5d, the standard URA can resolve up to eight targets under low SNR conditions. Figure 5d shows the 2-D DOA estimation of the ISR algorithm. Since the ISR algorithm only utilizes 15 DOFs for the initialization process of u or v, the maximum number of targets detected is limited. Therefore, the virtual aperture expansion ability of the proposed algorithm is better than that of the TNNR algorithm, ISR algorithm, PSS algorithm, and standard URA.

In the second experiment, we evaluate the estimation performance of all tested algorithms for uncorrelated targets. Figure 6 shows the RMSE of the angle estimate versus SNR, where two uncorrelated targets are located at $\left({\theta }_1,{\phi }_1\right)=\left({10}^{\circ },{10}^{\circ }\right)$ and $\left({\theta }_2,{\phi }_2\right)=\left({30}^{\circ },{30}^{\circ }\right)$, the number of snapshots is 10, and the SNR varies from 0 dB to 10 dB. As shown in Figure 6, when the SNR is greater than or equal to 4 dB, the proposed algorithm has the best estimation performance compared to the other test algorithms. The reason is that the proposed algorithm utilizes virtual aperture extension and sparse matrix recovery. However, when the SNR is lower than 4 dB, the performance of the ISR algorithm is better than that of the proposed algorithm because the matrix completion accuracy of the proposed algorithm decreases, leading to decrease in the estimation accuracy of the proposed algorithm. The performance of the TNNR algorithm degrades when the SNR is high. The cause is that the TNNR algorithm cannot complete the missing elements of the entire rows and columns in $\mathbf{Z}$, which significantly affects the sparse matrix recovery accuracy and reduces the DOA estimation performance. The PSS algorithm has the worst estimation performance due to not utilizing virtual aperture expansion. Furthermore, as the SNR increases, the estimated performance of all tested algorithms becomes relatively flat, which is limited by the pre-defined sampling grid.

In the third experiment, we compare the estimated performance of all tested algorithms versus the number of snapshots, where we assume two uncorrelated targets located at $\left({\theta }_1,{\phi }_1\right)=\left({10}^{\circ },{10}^{\circ }\right)$ and $\left({\theta }_2,{\phi }_2\right)=\left({30}^{\circ },{30}^{\circ }\right)$, the SNR is 0 dB, and the number of snapshots varies from 10 to 50. As demonstrated in Figure 7, when the number of snapshots exceeds 20, the proposed algorithm can obtain the optimal estimation performance due to the virtual aperture expansion and sparse matrix recovery. When the number of snapshots is less than or equal to 20, the ISR algorithm is better than the other test algorithms because it uses an iterative sparse recovery algorithm to improve the angle estimation performance. Under the condition of a low SNR and low snapshots, the TNNR algorithm outperforms the proposed algorithm. The explanation is that the TNNR algorithm has better completion accuracy than the proposed algorithm for empty elements that are not entire rows and columns. Therefore, the proposed algorithm has better estimation performance than the other test algorithms under mild conditions.

In the fourth experiment, we evaluate the estimation performance of all tested algorithms for closely spaced targets. We assume two closely spaced targets are located at $\left({\theta }_1,{\phi }_1\right)=\left({10}^{\circ },{10}^{\circ }\right)$ and $\left({\theta }_2,{\phi }_2\right)=\left({10}^{\circ }+∆,{10}^{\circ }+∆\right)$, where $∆$ varies from $2^{\circ }$ to ${10}^{\circ }$. Figure 8 shows the RMSE versus the angular separation, where the number of snapshots is 30, and the SNR is fixed at 0 dB. As shown in Figure 8, all the test algorithms fail when the angular separation is less than or equal to $4^{\circ }$. The rationale is that the two closely spaced targets will cause the rank loss of the cross-correlation matrix between the sub-arrays, thereby reducing the angle estimation performance. When the angular separation is more significant than $4^{\circ }$, the proposed algorithm shows optimal estimation performance compared to the ISR, TNNR, and PSS algorithms. The reason is that the proposed algorithm obtains the largest virtual aperture by the enhanced matrix completion technique, thereby improving the estimation performance. However, the ISR algorithm only utilizes part of the consecutive lags, so the performance of this algorithm is lower than that of the proposed algorithm. Since the TNNR algorithm cannot complete the empty elements of the entire rows and columns, its performance is lower than that of the proposed algorithm. However, it has more virtual apertures than the ISR and PSS algorithms, so the estimation performance is better than these two algorithms. Since the PSS algorithm exploits the co-prime property to remove ambiguity and does not exploit virtual aperture expansion, it has the worst estimation performance compared to the proposed algorithm and the ISR algorithm.

In the fifth experiment, we evaluate the RMSE of the proposed algorithm and the TNNR algorithm versus the percentage of observations. We set that the received data of a $19\times 19$ URA is ${\mathbf{X}}_{out}$, where the inter-element spacing is half a wavelength, and the SNR is set to 10 dB and $Q=1$. Suppose two targets are located at $\left({\theta }_1,{\phi }_1\right)=\left({10}^{\circ },{10}^{\circ }\right)$ and $\left({\theta }_2,{\phi }_2\right)=\left({30}^{\circ },{30}^{\circ }\right)$, and the searching grid is set to 0.002 for $u\in \left[0,1\right]$ and $v\in \left[0,1\right]$. We randomly select $\left[50\%,60\%,70\%,80\%,90\%\right]$ elements of ${\mathbf{X}}_{out}$ as available observation elements. Figure 9 shows RMSE versus the percentage of observations. It can be seen from Figure 9 that, when there are no missing elements in the entire row or column, the TNNR algorithm is slightly better than the proposed algorithm because the TNNR algorithm knows accurately the number of targets. However, the number of targets in practical problems is usually unknown, so the proposed algorithm has more practical application value.

In the sixth experiment, we compare the RMSE and CPU time of all tested algorithms versus the number of $N=\left[3,5,7,9\right]$, where M is fixed to 2. We assume two targets are located at $\left({\theta }_1,{\phi }_1\right)=\left({10}^{\circ },{10}^{\circ }\right)$ and $\left({\theta }_2,{\phi }_2\right)=\left({30}^{\circ },{30}^{\circ }\right)$, where the SNR is 0 dB, the number of snapshots is 50, and the searching grid is 0.0025 for $u\in \left[0,1\right]$ and $v\in \left[0,1\right]$. Figure 10 shows RMSE versus the total number of sensors. From Figure 10, as the number of sensors increases, the performance of all the tested algorithms increases. However, when the number of sensors is greater than 64, the performance of all the tested algorithms tends to be stable because the sampling grid limits further improvement in performance. Compared with the other tested algorithms, the proposed algorithm has the best performance because the proposed algorithm utilizes co-array interpolation initialization and enhanced matrix completion techniques. Figure 11 shows CPU time versus the total number of sensors, where the CPU time of all tested algorithms is based on the AMD R7-5800H processor.

From Figure 11, the CPU time of all the tested algorithms increases as the number of sensors increases. When the total number of sensors is less than 64, the CPU time of the proposed algorithm is higher than that of the PSS and TNNR algorithms and lower than that of the ISR algorithm. When the total number of sensors exceeds 64, the proposed algorithm has the highest CPU time due to the proposed algorithm’s co-array interpolation and enhanced matrix completion. Furthermore, the proposed algorithm has better estimation accuracy than other tested algorithms under certain conditions. Therefore, the proposed algorithm has a better trade-off between computational complexity and estimation accuracy.

4. Discussion

To exploit the virtual aperture of CPPA, we propose a sparsity-based, 2-D DOA estimation algorithm via enhanced matrix completion and co-array interpolation techniques. The first experiment in Section 3 shows that the proposed algorithm has the most significant DOF compared to the ISR, TNNR, and PSS algorithms. Therefore, the proposed algorithm can distinguish more targets than the ISR, TNNR, and PSS algorithms, consistent with the conclusions of [24].

In this study, we consider the estimation performance of the proposed method versus SNR and the number of snapshots. According to Figure 6 and Figure 7, we find that the proposed algorithm has a better estimation performance than the ISR, TNNR, and PSS algorithms when the SNR is greater than 4 dB and the number of snapshots is greater than 20. The reason is that the proposed algorithm utilizes the virtual aperture expansion of CPPA and the sparse matrix recovery. In addition, when performing matrix completion, we utilize the shift-invariant structure of the virtual URA, thereby increasing the matrix completion accuracy, consistent with the ideas of [29]. The ISR algorithm performs better when the SNR is less than 4 dB and the number of snapshots is less than 20. The cause is reduced accuracy of the proposed enhanced matrix completion under low SNR and low snapshot conditions, resulting in a basis mismatch phenomenon when performing sparse matrix recovery, which reduces the estimation performance. Therefore, the DOA estimation performance of the proposed method is not only controlled by the virtual aperture, but also related to the SNR and the number of snapshots. From Figure 10 and Figure 11, the proposed algorithm’s estimated performance and CPU time increase as the number of sensors increases. However, when the number of sensors exceeds 64, the proposed algorithm has a higher CPU time than the other tested algorithms. Therefore, the proposed algorithm is suitable for scenarios with a small number of sensors. Furthermore, the performances of the proposed algorithm, the ISR algorithm, the TNNR algorithm, and the PSS algorithm become relatively flat as the SNR and the number of snapshots and sensors increase. Due to the pre-defined sampling interval limiting further improvement in estimation performance, this is consistent with the conclusions of [17].

However, the proposed method also has certain limitations. When performing virtual aperture expansion, we do not consider the position and phase errors of the sensors, which are unavoidable in practical applications. If there are sensor position and phase errors, this will significantly reduce the matrix completion accuracy of the proposed algorithm, thereby reducing the estimation performance. Therefore, in future research, we need to consider the position and phase errors of the sensors and construct a matrix completion model with errors to further improve the robustness of matrix completion.

5. Conclusions

This paper proposes a complex-valued, sparse matrix-recovery-based 2-D DOA estimation algorithm for CPPA through enhanced matrix completion. First, we use the co-array interpolation technique to reconstruct two covariance matrices with the Toeplitz structure to initialize the holes of the SDCA, thereby constructing a virtual URA that does not contain missing elements of the entire rows and columns. Subsequently, an enhanced matrix completion model based on nuclear norm minimization is developed to complete all empty elements in the virtual URA, avoiding aperture loss. To reduce the computational burden of matrix completion, we develop a fast iterative algorithm based on the ADMM framework. Finally, we derive a complex-valued, sparse matrix-recovery algorithm based on the FISTA method to estimate the 2-D DOA of targets to reduce the computational complexity of the traditional vector-form, sparse recovery algorithm. The simulation results demonstrate the superior performance of the proposed algorithm.